PhD Laurent L'homme

PhD David Dombrowicz
email: david.dombrowicz@pasteur-lille.fr

in vitro models of human mast cells: How to get more and better with induced pluripotent stem cells?

Keywords: Mast cell, Induced pluripotent stem cell (iPSC), mast cell

Mast cells (MCs) are tissue-resident innate immune cells mostly present in connective and mucosal tissues, as in those constituting the skin, lung and gastrointestinal tract. MCs play key roles in host defense against pathogens and venoms but are also directly involved in the development of several pathological conditions such as allergic reactions, mastocytosis, atherosclerosis and cancer [START_REF] Dudeck | Mast cells as protectors of health[END_REF] . MCs express several receptors including the high affinity receptor for IgE (FcRI), CD117/c-Kit, Mas-related G-protein coupled receptor member X2 (MRGPRX2), several Toll-like receptors (TLRs) and complement receptors, allowing them to sense and respond to several innate stimuli. According to the MC receptor involved, activation triggers the release of preformed mediators stored in granules and/or of neosynthetized mediators such as short lived prostaglandins and leukotrienes and, at a later time point, cytokines. The presence of cytoplasmic granules and the ability to degranulate within minutes allows MCs to respond faster than most immune cells upon physiological or pathological activation. MC granules contains a wide variety of molecules including bioactive amines (histamine), proteases (chymases, tryptases), anti-microbial peptides, proteoglycans and cytokines 2 .

Several MC-deficient mouse models, mice with constitutive granule compound-deficiency and mice with MC-specific deficiencies for target genes have been developed to study the role of MCs in vivo. These tools have substantially increased our knowledge on the role of MCs in physiological and pathophysiological conditions. Moreover, experimental studies in vitro have provided a better understanding of the molecular mechanisms underlying MC functions such as the signaling events involved in FcRI-dependent degranulation. Despite the strong improvement in cell separation techniques during the past decades, purification of human tissue-resident MCs remains unpractical due to the low number of MCs in tissues and the issues related to the obtention of large biopsies. As MCs differentiate in peripheral tissues from blood progenitors arising from bone marrow, blood has been used as primary source to generate human MCs in vitro, giving rise to cord/peripheral bloodderived MCs (CB/PB-MCs). MC cell lines have also been generated and several human cell lines (HMC-1, LAD, LUVA, ROSA, MCPV-1) are currently available. These models display specific advantages and limits (Table 1), and no model fully recapitulating human tissue-resident MC phenotype is currently available.

In this issue of the Journal of Allergy and Clinical Immunology, Luo et al. described a new protocol for the generation of MCs from human induced pluripotent stem cells (iPSCs), namely iPSC-derived MCs (iPSC-MCs) [START_REF] Luo | A novel approach for studying mast cell-driven disorders: Mast cells derived from induced pluripotent stem cells[END_REF] . Authors developed a four steps protocol where they first differentiate human iPSCs into hematopoietic stem cells (HSCs). The following steps are similar to classical protocols using HSCs from blood, in which HSCs are expanded and differentiated into MCs. Despite the additional step of iPSCs culture and its associated extra cost, the advantages of initiating the process from iPSCs are that these cells provide a self-renewing, stable and standardized material compared to blood progenitors. iPSCs are produced from somatic cells that are reprogrammed to become pluripotent and may be theoretically generated from any source, such as cells from patients harboring a specific mutation. In the past few years, iPSCs banks were generated and made accessible to the scientific community in order to share iPSC lines and improve the consistency between laboratories. Similar to MC cell lines, iPSC culture can be easily scaled up to provide a high number of MCs. iPSC-MCs may therefore be of high value for experiments requiring a large amount of cells with a high level of phenotypic homogeneity over time, such as for drug screening or mast cell activation tests (MAT). However, the greater reproducibility conferred by iPSC-MCs is achieved at the expense of genetic diversity. Primary human tissue MCs and in vivo experiments in animal models remain the gold standard methods to assess the biological relevance.

iPSC-MCs displayed all the characteristics of mature MCs such as the presence of abundant cytoplasmic granules containing tryptases, chymase and histamine or a high expression of the FcRI and CD117/c-Kit. Furthermore, iPSC-MCs properly degranulated in response to IgE and substance P, although at a slightly lower level than primary skin MCs. RNA sequencing analysis of iPSC-MCs and skin MCs indicated a high similarity between these two MC models but also revealed a non-negligible divergence in the transcriptomic signature. Among the differences, Luo et al. observed a less "mature" signature in iPSC-MCs and a lower expression of FcRI and MRGPRX2, probably underlying the reduced degranulation in response to IgE and substance P. Of note, MC cell lines and other MC models derived from progenitors, such as CB-MCs and PB-MCs, are also usually less mature than tissue MCs (Table 1), probably because the differentiation steps occurring physiologically in the tissue environment are extremely difficult to fully recapitulate in vitro. Moreover, MCs not only display some tissue-specificy in vivo but also different subpopulation coexist within the same tissue, with differences affecting granule composition 2 , receptor expression [START_REF] Fujisawa | Expression of Mas-related gene X2 on mast cells is upregulated in the skin of patients with severe chronic urticaria[END_REF] and response to stimuli 5 (Table 1). Various modifications in culture conditions were proposed to mimic these tissue-and subpopulation-specific phenotypes and favor the differentiation of MC with either a higher tryptase or chymase activity. Among them, in addition to stem cell factor (SCF) and interleukin (IL)-3 commonly used for MC differentiation, IL-6 or IL-4 were also employed to improve MC differentiation towards one of this phenotype. The culture conditions and the cytokines employed are critical for the yield and maturation of MCs and probably underlie the large variability observed between groups for the same MC model. Luo et al. used SCF, IL-3, IL-6 and, in a rather uncommon way, low-density lipoproteins (LDL) as a differentiation cocktail. In fact, the same authors had previously proposed LDL as a source of lipids to support proliferation and the high requirement of phospholipids for granule formation [START_REF] Schmetzer | A novel method to generate and culture human mast cells: Peripheral CD34+ stem cell-derived mast cells (PSCMCs)[END_REF] . In serum-free culture conditions, to counteract the effect of factors preventing MC proliferation and differentiation present in serum, Luo et al demonstrated that the addition of LDL strongly improved the yield of mature MC [START_REF] Schmetzer | A novel method to generate and culture human mast cells: Peripheral CD34+ stem cell-derived mast cells (PSCMCs)[END_REF] . While not being totally innovative compared to published protocols using blood progenitors, Luo et al. provide a straightforward protocol using only commercial and standardized reagents to generate phenotypically homogenous MCs from IPSC lines. iPSC-MCs generated with such a protocol are mature, functional, close to tissue-resident MCs and appropriated for studying both IgE-dependent and IgE-independent activation. Whether iPSC-MCs are also suitable for studying other MC functions such as phagocytosis, IgG-dependent activation or production of extracellular traps remains to be determined.

The limited amount of mature MCs in tissues has led to the development of a large variety of in vitro models. MC cell lines are easy-to-grow and provide cells with an homogeneous phenotype but show several abnormalities, such as immature phenotypes and the lack of key MC characteristics [START_REF] Passante | Mast Cell and Basophil Cell Lines: A Compendium[END_REF][START_REF] Guhl | Mast cell lines HMC-1 and LAD2 in comparison with mature human skin mast cells--drastically reduced levels of tryptase and chymase in mast cell lines[END_REF] (Table 1). CB/PB-MCs are a widely used model for human studies that show a highly similar phenotype with primary MCs 3,9 (Table 1). However, CB-MCs express less FcRI and show a reduced IgE-dependent activation than PB-MCs [START_REF] Iida | Selective downregulation of high-affinity IgE receptor (FcepsilonRI) alpha-chain messenger RNA among transcriptome in cord blood-derived versus adult peripheral blood-derived cultured human mast cells[END_REF] . Unsurprisingly, iPSC-MCs, that are differentiated in a similar manner as CB/PB-MCs, also exhibit an advanced maturity and functionality close to that from tissue-resident MCs (Table 1). Both CB/PB-MCs and iPSC-MCs are interesting models that phenocopy several features of tissue-resident cells. The main limitations of CB/PB-MCs are the donor's heterogeneity and the difficulty of upscaling MCs production. Both issues can be solved by the use of IPSCs. Thus, iPSC-MCs represent an interesting alternative when a large amount of phenotypically homogenous MCs is required or for long term studies, as they provide a stable material that should help to reduce inter-and intra-laboratory variability.

Table 1 :

 1 Main characteristics of the most common in vitro models of MCs. There is currently no study comparing all these models simultaneously in the same settings. It is important to note that several protocols exist for each model giving different level of MC maturation. This table provides a summary of commonly obtained results for each model.

		Tissue-resident MCs	Cell lines	CB/PB-MCs	iPSC-MCs by Luo et al.
	Availability of material	Limitation of biopsy (access/ethics, size, ...)	Easy to expand and maintain	Limitation of blood (access/ethics, volume, ...)	iPSCs easy to expand and maintain
	Donor heterogeneity	Yes	No	Yes	No
	Preparation time	4 to 24 hrs	Ready to use	21 to 63 days	52 days
	Typical range yield	0.1 to 1 x 10^6 per g of tissue.	> 10^6 (doubling time from 48 hrs to 2 weeks	0.5 to 5 x 10^6 per 50 mL of blood	0.5 to 0.8 x 10^6 per 40 iPSC colonies
	Cell surface
	markers
	CD117/c-Kit	++++	+++(+) (WT or mut)	++++	++++
	FcRI	++++	-(HMC-1) or ++(+)(+)	++(+)(+)	+++(+)
	MRGPRX2/MrgprB2	+(+)(+)(+)	+(+)(+)	+++	+++
	Granule
	composition
	Histamine	+++	+(+)	++++	++++
	Tryptases	++++	+(+)	++++	++++
	Chymases	-or ++++	-(HMC-1) or +(+)	++++	++++
	Degranulation in
	response to
	IgE-mediated receptor cross-linking	++++	-(HMC-1) or +(+)(+)	++(+)(+)	+++
	MRGPRX2/MrgprB2 receptor activation	+(+)(+)(+)	+(+)(+)	+++	+++

DD work is supported in part by grants from the ANR and the European Union: EGID ANR-10-LABX-46 and by Fondation pour la Recherche Médicale (Nouveaux défis en allergologie).


Jérôme Carrayrou
email: jerome.carrayrou@unistra.fr

Caroline Bertagnolli

Marwan Fahs

Algorithms for activity correction models for geochemical speciation and reactive transport modelling

Keywords: Geochemical modelling, instantaneous equilibrium chemistry, activity correction, outer fixed-point algorithm, Davis activity, Debye-Hückel activity, B-dot model

Reactive transport softwares are today one of the cornerstones of environmental research. They contain multiphysics with very complex algorithms, including flow, transport, chemical and sometimes heat transport, mechanical and/or biological algorithms. Because of this complexity, some parts of these algorithms still have not been sufficiently studied. In this work, we focus on algorithms for activity correction, a specific subset of equilibrium chemistry algorithms. We show that the most used algorithm (the inner fixed-point algorithm) or the most rigorous algorithm (the full Newton) might not be the most efficient, and we propose a new one, the outer fixed-point algorithm, which is more robust and faster than other algorithms.

INTRODUCTION

The problem of groundwater management has received increasing attention, and many tools have been developed to address this issue. One of these tools, reactive transport modeling, was first limited to laboratory experiments [START_REF] Lefèvre | Migration of strontium in clayey and calcareous sandy soil: Precipitation and ion exchange[END_REF] and then extended to the comprehension of problems in various fields [START_REF] Wanner | A chromate-contaminated site in southern Switzerland -Part 2: Reactive transport modeling to optimize remediation options[END_REF]. Reactive transport modeling is actually a mature research field that has produced important results in many environmental domains, such as water management, sea water intrusion [START_REF] Werner | Seawater intrusion processes, investigation and management: Recent advances and future challenges[END_REF], long-term nuclear waste storage [START_REF] De Windt | Intercomparison of reactive transport models applied to UO2 oxidative dissolution and uranium migration[END_REF], CO 2 sequestration [START_REF] Navarre-Sitchler | Elucidating geochemical response of shallow heterogeneous aquifers to CO2 leakage using high-performance computing: Implications for monitoring of CO2 sequestration[END_REF][START_REF] Pruess | Code intercomparison builds confidence in numerical simulation models for geologic disposal of CO2[END_REF] and heavy metal contamination [START_REF] Wanner | Benchmarking the simulation of Cr isotope fractionation[END_REF].

Accepted Article

This article is protected by copyright. All rights reserved.

Numerous reactive transport codes are available, and some review articles [START_REF] Pruess | Code intercomparison builds confidence in numerical simulation models for geologic disposal of CO2[END_REF][START_REF] Steefel | Reactive transport codes for subsurface environmental simulation[END_REF][START_REF] Yeh | A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components[END_REF][START_REF] Yaron | Contaminant geochemistry-a new perspective[END_REF][START_REF] Steefel | New directions in hydrogeochemical transport modeling: Incorporating multiple kinetic and equilibrium reaction pathways[END_REF][START_REF] Van Der Lee | Present state and future directions of modeling of geochemistry in hydrogeological systems[END_REF][START_REF] Steefel | Reactive transport modeling: An essential tool and a new research approach for the Earth sciences[END_REF] propose an overview of them.

General algorithms and/or numerical methods for geochemical speciation modelling are widely recognized as very challenging due to the non-linearity of the mathematic system, the very large range of magnitude of the chemical species concentrations and the equilibrium constants. Since first works [START_REF] Morel | Numerical method for computing equilibriums in aqueous chemical systems[END_REF][START_REF] Westall | MINEQL: a computer program for the calculation of chemical equilibrium composition of aqueous system[END_REF][16], many authors worked on improving these algorithms and/or numerical methods [START_REF] Zhang | Reactive Transport Modeling in Variably Saturated Media with MIN3P: Basic Model Formulation and Model Enhancements[END_REF][START_REF] Steefel | CRUNCHFLOW: Software for Modeling Multicomponent Reactive Flow and Transport[END_REF][START_REF] Steefel | Reactive transport codes for subsurface environmental simulation[END_REF][START_REF] Parkhurst | PhreeqcRM: A reaction module for transport simulators based on the geochemical model PHREEQC[END_REF][START_REF] Hoffmann | A general reduction scheme for reactive transport in porous media[END_REF][START_REF] Poulet | escript-RT: Reactive transport simulation in Python using escript[END_REF][START_REF] Xu | TOUGHREACT Version 2.0: A simulator for subsurface reactive transport under non-isothermal multiphase flow conditions[END_REF][START_REF] Lichtner | FLOTRAN users manual: two-phase nonsothermal coupled thermalhydrologic-chemical (THC) reactive flow and transport code version 2[END_REF][START_REF] Van Der Lee | Module-oriented modeling of reactive transport with HYTEC[END_REF][START_REF] Meeussen | ORCHESTRA: An Object-Oriented Framework for Implementing Chemical Equilibrium Models[END_REF][START_REF] Paz-García | Computing multi-species chemical equilibrium with an algorithm based on the reaction extents[END_REF][START_REF] Carrayrou | New efficient algorithm for solving thermodynamic chemistry[END_REF][START_REF] Nazer | Convergence acceleration of iterative sequences for equilibrium chemistry computations[END_REF]. An examination of these articles shows that all these codes include one or more activity correction models. Even though the different models of activity correction are usually well-detailed, the algorithmic method used to compute the activity coefficients and to incorporate these calculations into the entire chemistry algorithm is usually not given. We found mention of it only in the user's guide for CHESS [START_REF] Van Der Lee | Thermodynamic and mathematical concept of CHESS[END_REF] and PHREEQC [31] codes.

In this work, we will present the two algorithmic ways found in the literature to incorporate activity correction models into an entire geochemical speciation algorithm; propose a new one and compare these three algorithms. All these algorithms can be easily implemented into geochemical codes but they are not all equivalent and some leads to an increase of numerical instabilities. To discuss on activity correction algorithm, we need to briefly present an entire geochemical speciation algorithm.

Although some approaches have abandoned formulation based on mass action law [START_REF] Leal | Efficient chemical equilibrium calculations for geochemical speciation and reactive transport modelling[END_REF][START_REF] Allan | An overview of computational methods for chemical equilibrium and kinetic calculations for geochemical and reactive transport modeling[END_REF] or resolution based on a "Newton" type method [START_REF] Nazer | Convergence acceleration of iterative sequences for equilibrium chemistry computations[END_REF][START_REF] Amir | Jacobian Free Methods for Coupling Transport with Chemistry in Heterogenous Porous Media[END_REF][START_REF] Nazer | Sur la modélisation du transport réactif dans les réserves d'eaux potables. Université du Littoral Côte d'Opale[END_REF], we will apply our activity correction algorithms to this historical formulation, so as not to make it too cumbersome. Because we will focus on a specific part of the geochemical algorithm, we will not discuss, neither on the specific methods used to overpass the high non-linearity of the mathematical problem, nor on the ways to include specific chemical phenomena such as precipitation-dissolution, liquid-gas exchange or surface complexation modelling. This simplification made for the presentation of the model does not imply any loss of generality because the specificity neglected here are usually treated as an "external" loop in geochemical modelling [START_REF] Steefel | CRUNCHFLOW: Software for Modeling Multicomponent Reactive Flow and Transport[END_REF][START_REF] Carrayrou | New efficient algorithm for solving thermodynamic chemistry[END_REF][START_REF] Brassard | A feasible set for chemical speciation problems[END_REF][START_REF] Van Der Lee | Module-oriented modeling of reactive transport with HYTEC[END_REF] We first present the general concepts including fundamental equations (mass action law, activity correction model) and a simple formulation of the Newton method applied to geochemical speciation. The second part is devoted to the presentation of the two algorithms for activity correction found in literature and to the presentation of the new one. The specificities of these three algorithms will be explained on a tricky test case. The third part of this article is devoted to benchmarking the studded algorithms against 3 different and challenging test cases. We can then conclude and propose some recommendations to improve the robustness and efficiency of geochemical speciation codes.

GENERAL CONCEPTS

Fundamental equations

A general formulation of a chemical reaction leading to the formation of one of the Nc species (C i ) from the number Nx of chosen components X j depending on stoichiometric coefficients a i,j is written as:

, 1 Nx i i j j j C a X (1)
Components and species [START_REF] Morel | Numerical method for computing equilibriums in aqueous chemical systems[END_REF] may be named Primary and Secondary species [START_REF] Lichtner | Principles and practice of reactive transport modeling[END_REF] in some other works. We prefer the term "component" rather than "primary species". Saying "Primary species" implies that they are a subsets of the chemical species whereas the term "component" does not. By this way a component can be different from a chemical species. Some authors [START_REF] Lützenkirchen | The constant capacitance model and variable ionic strength: An evaluation of possible applications and applicability[END_REF][START_REF] Lützenkirchen | Parameter estimation for the triple layer model. Analysis of conventional methods and suggestion of alternative possibilities[END_REF] describe surface complexation phenomena by including the electrical potential at the surface as a component.

In this work, we use the electron e -as component to describe redox reaction in the iron-chromium test case.

Instantaneous equilibrium chemistry is usually described using two fundamental concepts: mass conservation equations and mass action laws. According to the classical formulation stated by Morel

and Morgan [START_REF] Morel | Numerical method for computing equilibriums in aqueous chemical systems[END_REF], mass conservation equations describe the conservation of the total concentrations of the components (T j ), and mass action laws describe the formation of each chemical species as a combination of the Nx chosen components.

Mass conservation equations are written using the species concentrations [C i ] to conserve the total concentration T j of each component:

, 1 Nc j i j i i T a C (2)
On the other hand, mass action laws are written using the species {C i }, the equilibrium constant K i and components {X j } of the activities:

Accepted Article

This article is protected by copyright. All rights reserved. (3)

To ensure the closure of the system, an activity model is used. The activity coefficient ( i ) allows determining the species activity from its concentration.

i i i CC ( 4
)
Several activity models have been developed, most of which use the ionic strength of the solution (I)

and the electric charge of each species z i :

2 1 1 2 Nc ii i I z C (5)
One of the most used activity models for reactive transport is the Davies model, where the activity coefficients are given by:

2 log 1 ii I Az BI I (6)
The parameters A and B are defined by: Because we found that the algorithmic behavior of the methods is quite independent on the activity model used, we only present here the Davies model and all the tests will be done using this one. Details on extended Debye-Hückel and B-dot model are given on appendix A-1.

Newton algorithm

The methodology for computing chemical equilibrium is well established [START_REF] Morel | Numerical method for computing equilibriums in aqueous chemical systems[END_REF][START_REF] Westall | MINEQL: a computer program for the calculation of chemical equilibrium composition of aqueous system[END_REF][START_REF] Carrayrou | New efficient algorithm for solving thermodynamic chemistry[END_REF][START_REF] Van Der Lee | Thermodynamic and mathematical concept of CHESS[END_REF]31,[START_REF] Marinoni | Thermodynamic equilibrium solutions through a modified Newton Raphson method[END_REF][START_REF] Steefel | Crunch -user's guide[END_REF] and is usually built on a Newton procedure. By incorporating the Nc mass action laws into the Nx conservation equation, one can obtain an Nx by Nx nonlinear system that must be iteratively solved. Here, we present the approach where mass action laws are written in a logarithmic form.

We define the logarithm of the components activity j  as: ln ln j j j j XX [START_REF] Steefel | Reactive transport codes for subsurface environmental simulation[END_REF] This can be rewritten as a matrix formulation:

ln ln X diag X (9)

Accepted Article

This article is protected by copyright. All rights reserved.

Where ln X diag is a diagonal matrix, size These can be rewritten by using the vector of species concentrations C :

ln ln ln C diag C K A (11)
Where

Combining both equations, we obtain:

exp ln ln T C diag T A K A (13)
This nonlinear system is iteratively solved using a Newton procedure. At the nth iteration, the error is given by:

exp ln ln n T n n C diag Y T A K A (14)
By computing the derivative of the error versus the logarithm of the components of the activity, we obtain the Jacobian matrix Z .

Historically, activity correction has been neglected when computing the Jacobian matrix. This approximation is justified if:

(i)
The ionic strength is small enough so that the activity coefficients can be assumed to be equal to one (ideal solution). (ii)

The changes in ionic strength are sufficiently small during one Newton iteration so that the activity coefficients can be assumed to be constant.

The complete computation of the Jacobian matrix leads to equation ( 15):

, , , 11 , , , 11 , , , 1

exp ln ln exp ln ln exp ln ln

Nc Nx j j k i j i i i h h ih kk Nc Nx j k i j i i i h h ih k Nx j k i j i i i h h ih Y Z a K a Z a K a Z a K a , 11 , , , 1 , , , 1 ln ln ln
ln Nc Nx i i i h h h k Nc j k i j i i i k k i k Nc i j k i j i i k i k Ka Z a C a Z a C a (15)

Accepted Article

This article is protected by copyright. All rights reserved.

This is because:

, , , , 11 exp exp Nx Nx ii i h h i k i h h i k i hh k i i KK a a a a c (16)
Neglecting the derivative of the activity coefficient (hypothesis (i) or (ii)), one obtains (17):

, , , 1
Nc j k i j i k i i Z a a c ( 17
)
Alternatively, with diag C a diagonal matrix, size

, Nc Nc and , i ii diag c C : T diag Z A C A (18)
Solving the linear system [START_REF] Steefel | Reactive transport codes for subsurface environmental simulation[END_REF] gives the step of the Newton method n :

n n n ZY (19)
Then, a new value of the activity component is given by:

1 n n n ( 20
)
The procedure is repeated until the error is sufficiently small. In this work, we check the convergence using relation [START_REF] Hoffmann | A general reduction scheme for reactive transport in porous media[END_REF] with the convergence criteria 12 [START_REF] Yaron | Contaminant geochemistry-a new perspective[END_REF] . Relation (21) defines a relative error Err that ensures an equivalent accuracy for both major and trace element.

, 1 max j Nc j i j i i Y Err T a C ( 21
)
If convergence is not reached until a prescribed number of iterations, 1,000 in this work, the procedure is stopped and the method fails to converge. The robustness of a method reports its ability to reach the convergence within the maximum prescribed iterations.

Accepted Article

This article is protected by copyright. All rights reserved.

POSSIBLE ALGORITHMS FOR ACTIVITY CORRECTION

We found two algorithms reported in the literature to handle activity correction: the inner fixedpoint and the full Newton algorithm. We suppose this point is poorly reported for several reasons:

(i) Activity coefficients do not change a lot and are usually close to one. Because they depends on the ionic strength of the solution, they are supposed to be less modified during the optimization loops. (ii)

Other numerical difficulties are much more complicated to handle. Among them: the large range of magnitude of equilibrium constants, species and components concentrations; important nonlinearities of the system due to stoichiometric coefficients in the mass action laws; very high condition number of the Jacobian matrix. (iii)

Activity correction can be included in the global geochemical algorithm by a very simple way (see inner fixed-point algorithm).

Nevertheless, by working on numerical instabilities on computing geochemical equilibrium, we found that points (i) and (ii) are not verified for some specific chemical system. We then develop a third algorithm to fix these problems: the outer fixed-point algorithm.

Inner fixed-point activity algorithm

This algorithm is the most commonly used algorithm, according to the few information obtained. It is used by CHESS codes [START_REF] Van Der Lee | Thermodynamic and mathematical concept of CHESS[END_REF], SPECY [START_REF] Carrayrou | Looking for some reference solutions for the reactive transport benchmark of MoMaS with SPECY[END_REF] and we found it be examining the FORTRAN code of IMPACT [START_REF] Jauzein | A flexible computer code for modelling transport in porous media: impact[END_REF].At each Newton loop, the activity coefficients are updated according to the new ionic strength (Figure 1). From a strictly mathematical point of view, this method is then a quasi-Newton method because the Jacobian matrix Z is an approximation of the derivative of the error using equation [START_REF] Zhang | Reactive Transport Modeling in Variably Saturated Media with MIN3P: Basic Model Formulation and Model Enhancements[END_REF].

Full Newton activity algorithm

One can find an approach in the PHREEQC code [31] for including the derivative of the activity coefficients in the Jacobian matrix. This lead to the Full Newton algorithm (Figure 1) where the Jacobian matrix is computed using equation [START_REF] Westall | MINEQL: a computer program for the calculation of chemical equilibrium composition of aqueous system[END_REF].

The derivative of the activity coefficients must be calculated. Because the activity coefficients are strongly dependent on the ionic strength, we write: ln ln ii kk I I [START_REF] Poulet | escript-RT: Reactive transport simulation in Python using escript[END_REF] and then compute two parts:

(i)
The derivative of the activity coefficients versus the ionic strength

ln i I
. This component is model-dependent and simple to compute regardless of which model is used. Details for Davies model are given on equation [START_REF] Nazer | Sur la modélisation du transport réactif dans les réserves d'eaux potables. Université du Littoral Côte d'Opale[END_REF] , details for extended Debye-Hückel and B-dot models are given on Appendix A-1.

Accepted Article

This article is protected by copyright. All rights reserved.

(ii)

The derivative of the ionic strength versus the logarithm of the activity components, which is independent on the activity correction model used.

We propose here two computation ways. A recursive formulation using equation ( 29) and an explicit formulation using equation [START_REF] Amir | Jacobian Free Methods for Coupling Transport with Chemistry in Heterogenous Porous Media[END_REF]. Both formulations lead to the same Jacobian matrix whatever the test case and only differs on the computing time. Except for CPU time comparison, we will use the explicit formulation.

Computation method 1: a recursive formulation

Computing the derivative of the ionic strength yields:

C I z z C a ( 23
)
We obtain a recursive formulation for the derivative of the activity coefficients: [START_REF] Lichtner | FLOTRAN users manual: two-phase nonsothermal coupled thermalhydrologic-chemical (THC) reactive flow and transport code version 2[END_REF] This recursive formulation can be rearranged as: [START_REF] Van Der Lee | Module-oriented modeling of reactive transport with HYTEC[END_REF] The vector k N contains the explicit portion of the equation. This vector differs depending on the component used to derive the activity coefficient through the stoichiometric coefficient

p ,k a . 2 , 1 ln 1 2 Nc i i p p k p p N z a C I (26)
The matrix M contains the coefficients of the implicit portion of the equation. This matrix is the same regardless of the component used in the derivation. A complete derivative of the activity coefficients is given by the following linear system:

ln ln k kk NM or ln k k Id M N (28)

Accepted Article

This article is protected by copyright. All rights reserved.

Computation method 2: an explicit formulation

Computing the derivative of the ionic strength yields: [START_REF] Leal | Efficient chemical equilibrium calculations for geochemical speciation and reactive transport modelling[END_REF] We obtain an explicit formulation for the derivative of the ionic strength:

This allows for an explicit formulation of the derivative of the activity coefficients:

Accepted Article

This article is protected by copyright. All rights reserved.

Derivative of the activity coefficients versus the ionic strength for the extended Debye-Hückel and the B-dot models can be found on Appendix A-1.

Outer fixed-point activity algorithm

Here, we propose this new algorithm. For a given set of values for the activity coefficients, the Newton procedure is iteratively run until convergence is achieved (Loop 1 on Figure 1). Once convergence is reached, the activity coefficients are updated according to the ionic strength computed at the time of convergence (Loop 2 on Figure 1). This procedure is repeated until no changes are detected neither on the component activities (Test convergence = YES on Figure 1) nor on the activity coefficients (Activity changed = YES on Figure 1). We then obtain a true Newton algorithm where the Jacobian matrix calculated by equation ( 17) is truly the derivative of the error.

Comparison of the algorithms on a trivial test case

To compare the activity correction algorithms without any interferences of the other numerical difficulties reported to the Newton procedure, we propose a test case without any chemical reactions (Table 1) called "test case with only activity correction" (details can be found appendix A-2). Any nonlinearity is only due to activity correction. The chemical system is composed of chloride Cl  ions, calcium 2 Ca  ions, aluminum 3 Al  ions and tin 4 Sn  ions. We neglect water dissociation and all chemical reactions. The details and equilibrium solutions are given in appendix A-2.

One should note that this test case is not chemically realistic. Moreover, its numerical value comes only from activity correction if the unknowns of the nonlinear system [START_REF] Morel | Numerical method for computing equilibriums in aqueous chemical systems[END_REF] are the logarithms of the activity components j  [START_REF] Steefel | Reactive transport codes for subsurface environmental simulation[END_REF]. Otherwise, if the unknowns are component concentrations j X , the problem becomes trivial and linear, and its solution is the total concentration 1,..., jj X T j Nx .

We present two scenarios for the test case with only activity correction: one with a low ionic strength and one with a high ionic strength. The objective is to determine the influence of the activity correction on the Newton procedure depending on the algorithm used. For the situation with a low ionic strength, this influence is expected to be negligible, whereas we expect a greater impact in the situation with a high ionic strength.

For the low ionic strength situation, the initial component activities are 5.0 10 -7 mol l -1 for all components. The ionic strength is 7.80 10 -6 mol l -1, and we obtain the species concentrations and activity values, which are given in Table 2. Also in Table 2, we show the first Newton steps   proposed by the fixed-point algorithms (inner and outer) and by the full Newton algorithm.

Accepted Article

This article is protected by copyright. All rights reserved.

The Jacobian matrices given by both fixed-point algorithms are the same: 5.00 10 0.00 0.00 0.00 0.00 5.00 10 0.00 0.00 0.00 0.00 5.00 10 0.00 0.00 0.00 0.00 5.00 10

Z                 (36)
The full Newton algorithm leads to the following derivative matrix of the activity coefficients: The following Jacobian matrix is also obtained:

Z ( 38
)
For the scenario with a low ionic strength, the Jacobian matrices are quite equivalent for the full Newton and fixed-point algorithms, leading to similar Newton steps (   ).

The high ionic strength scenario starts with an initial component activities equal to 5.0 10 -4 mol.l -1 for all components. The ionic strength equals 1.197 mol.l -1 . The species concentrations, activity values and first Newton steps   proposed by the fixed-point algorithm (inner and outer) and full Newton algorithm are given in Table 3.

The Jacobian matrices given by both fixed point algorithms are the same: 5.00 10 0.00 0.00 0.00 0.00 5.00 10 0.00 0.00 0.00 0.00 5.00 10 0.00 0.00 0.00 0.00 5.00 10

Z ( 39
)
The full Newton algorithm leads to the following derivative matrix of the activity coefficients:

Accepted Article

This article is protected by copyright. All rights reserved. 1. [START_REF] Carrayrou | Looking for some reference solutions for the reactive transport benchmark of MoMaS with SPECY[END_REF] The following Jacobian matrix is also obtained: 7.12 10 1. [START_REF] Steefel | CRUNCHFLOW: Software for Modeling Multicomponent Reactive Flow and Transport[END_REF] For the scenario with a high ionic strength, the Jacobian matrices Z are very different for the full Newton and fixed-point algorithms, leading to very different Newton steps (   ). Moreover, we find an increase in the condition number of matrix Z for the full Newton algorithm. The condition number of the matrix ( 41) is 21.6, whereas it equals 1 for matrix [START_REF] Lützenkirchen | The constant capacitance model and variable ionic strength: An evaluation of possible applications and applicability[END_REF].

The symmetry of the Jacobian matrix in equations ( 38) and ( 41) is specific to this test case. As shown in equation ( 18), the Z matrix for the fixed-point algorithms is symmetric, whereas equation [START_REF] Westall | MINEQL: a computer program for the calculation of chemical equilibrium composition of aqueous system[END_REF] proves that it is usually not symmetric for the full Newton algorithm.

We solved this test case with only activity correction using the three presented algorithms with the low and high ionic strength scenario. We plot on a graph (Figure 2) the evolution of the norm of the error (eq 16) versus the Newton iterations.

Comparing the three algorithms on the test case with only activity correction, one can see in Figure 2 that:

-The outer fixed-point algorithm requires the fewest Newton iterations to reach convergence, whereas the inner fixed-point algorithm requires the most iterations. The full Newton algorithm requires an intermediate number of Newton iterations. -For the three algorithms, obtaining the solution at a low ionic strength requires less Newton iteration than at a high ionic strength. This point is obvious: for this test case, activity corrections are the only nonlinearity of the problem, and they are less important at a low ionic strength than at a high ionic strength. -The outer fixed-point algorithm runs 3 minimization loops (Loop 2 on Figure 1) for situations with both low and high ionic strength. The first loops converge at 10 (low) and 23 (high) iterations; the second loop converges at 15 (low) and 29 (high) iterations. The third loop is the confirmation loop used to check that no changes in the ionic strength computation occur and then to confirm the convergence of the algorithm (Test convergence is YES and activity changed is YES on Figure 1)

Accepted Article

This article is protected by copyright. All rights reserved.

It seems that the outer fixed point algorithms converges at lower value than others but it's an illusion. Convergence criteria  is set to 10 -12 for the three algorithms (eq 19). This illusion is because the convergence rate of the outer fixed point algorithms is higher than other algorithms near the solution, leading to a sharper slope of it's curve on Figure 2).

Accepted Article

This article is protected by copyright. All rights reserved.

BENCHMARKING OF THE THREE ALGORITHMS

Testing procedure

It is well-known that the initial guess of the values of the components plays a critical role in the convergence of the Newton methods [START_REF] Steefel | Reactive transport modeling: An essential tool and a new research approach for the Earth sciences[END_REF][START_REF] Hoffmann | A general reduction scheme for reactive transport in porous media[END_REF]. To test several initial guesses, we generate a large number (30 000) of activity component values according to the following procedure: We then obtain 30 000 realizations of the optimization procedure using the same chemical test case but different initial guesses. To analyze this large amount of data, we construct a frequency graph of the number of Newton iterations needed to reach convergence. We plot graphs of the cumulative ratio of the resolutions that converge within a given number of Newton iterations. It gives the ratio of simulations that converge within a given number of Newton iteration. According to the graph, the algorithm that reaches a cumulative frequency of 1 is said to be robust. The algorithm that reaches a high cumulative frequency for a low number of Newton iterations is said to be fast. The difference between 1 and the highest cumulative frequency gives the failure ratio of the algorithm.

min max log 1 log log X rand X rand X (42)

Chemical test cases

We propose 5 chemical test cases including the only activity correction test case presented previously. This selection includes a quite easy test (phosphoric acid), a quite difficult one (gallic acid), one very hard (iron-chromium) and one of very large size (Morel-Morgan).

Phosphoric acid test case

This test presents reactions between phosphoric acid and salt water. It includes 4 components and 8 chemical species. We handle only acid-base reactions: water dissociation and the 3 phosphoric acid reactions. This test case is numerically a simple chemical test, which is interesting in this context because of the PO 4

3-species which leads to important activity correction due to its -3 electrical charge. A table including the stoichiometric coefficients, equilibrium constants, total concentrations and equilibrium solutions is given in appendix A-3.

Gallic acid test case

This test case was proposed by Brassard and Bodurtha [START_REF] Brassard | A feasible set for chemical speciation problems[END_REF]. It includes 3 components and 17 chemical species. It is a classical test case, and many difficulties in convergence have been reported while solving it by using Newton or Newton-like algorithms [START_REF] Carrayrou | New efficient algorithm for solving thermodynamic chemistry[END_REF][START_REF] Nazer | Convergence acceleration of iterative sequences for equilibrium chemistry computations[END_REF][START_REF] Marinoni | Thermodynamic equilibrium solutions through a modified Newton Raphson method[END_REF][START_REF] Machat | Comparison of linear solvers for equilibrium geochemistry computations[END_REF]. It should be selective according to the activity correction algorithm because of the presence of many highly charged

Accepted Article

This article is protected by copyright. All rights reserved.

species such as AlL 3 6-. A table including the stoichiometric coefficients, equilibrium constants, total concentrations and equilibrium solutions is given in appendix A-4.

Iron-chromium test case

This test case concerns the rehabilitation of chromium-contaminated industrial soil using an ironchromium reduction [START_REF] Wanner | A chromate-contaminated site in southern Switzerland -Part 2: Reactive transport modeling to optimize remediation options[END_REF]46]. Chromium (VI), which is the most toxic and mobile form of chromium, is reduced by iron (II) to yield chromium (III), which has a much lower solubility and is less toxic [47]. This test is reported to be a very difficult one [START_REF] Marinoni | Thermodynamic equilibrium solutions through a modified Newton Raphson method[END_REF][START_REF] Machat | Comparison of linear solvers for equilibrium geochemistry computations[END_REF], so here we use some favorable testing conditions to increase the convergence of the Newton algorithm. Its numerical difficulties are attributed to the redox reactions described that lead to a large range of stoichiometric coefficients (from -4 to +8) and equilibrium constants (from 10 -83.17 to 10

Morel-Morgan test case

This test has been presented in the article that founded Newton algorithm as a key tool of equilibrium computation [START_REF] Morel | Numerical method for computing equilibriums in aqueous chemical systems[END_REF]. It is one of the larger test case reported in the literature, using 52 components to generate 781 species. It includes highly charged species such as P 3 O 10 5-, a large range of stoichiometric coefficients (from -4 to +6) and equilibrium constants (from 10 -31.7 to 10 66 ). A table including the stoichiometric coefficients, equilibrium constants, total concentrations and equilibrium solutions is given in appendix A-6.

Frequency graphs Test case with only activity correction

The test case with only activity correction is a simple chemical test case. It makes sense only for studying the activity correction algorithms. It is solved by all the algorithms (see Table 5) within 150 Newton iterations (Figure 3). The fastest algorithm is the outer fixed-point algorithm, regardless of the ionic strength. Moreover, this algorithm shows a very low sensitivity to the ionic strength by resolving the low ionic strength case within 24 or 25 iterations and the high ionic strength case within 21 iterations regardless of the initial guess. The inner fixed-point and the full Newton algorithms are much more sensitive to the ionic strength, with significant increases in the number of iterations required to converge in the case with a high ionic strength. For this case, we find that the best algorithm is the outer fixed-point algorithm, and the inner fixed-point algorithm is the worst according to the number of Newton iterations. Taking the computing time of one Newton iteration into account (Table 4), we see that the full Newton algorithm is the slowest and the outer fixed-point algorithm is the fastest.

Accepted Article

This article is protected by copyright. All rights reserved.

point algorithm is the fastest in terms of CPU times for one Newton iteration, usually faster than or equivalent to the other algorithms in terms of the number of required Newton iterations and the most robust.

We have thus shown that the external fixed point algorithm provides greater robustness and speed than the other two. The influence of the choice of the activity correction algorithm is amplified in this work, because the general geochemical equilibrium resolution algorithm used here is extremely simple. The use of specific numerical methods to reduce or overcome the strong non-linearities of geochemical systems, as implemented in modern codes, will obviously make this influence less visible. Nevertheless, and according to the results presented here, we recommend the outer fixedpoint algorithm. This algorithm is the least time consuming for one Newton iteration, it usually requires the fewest number of iterations, and it is the most robust and least sensitive to the initial guess. Moreover, its implementation with existing geochemical algorithms is very simple and requires very few modifications.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

  . For example, precipitation-dissolution phenomena are usually solved as follow: (a) solving a first geochemical equilibrium with a given set of minerals. (b) Checking saturation of the existing and potential minerals. (c) Solving a new geochemical equilibrium including the existence of oversaturated (as calculated at step c) minerals and removing the undersaturated minerals. (d) Repeat (a), (b) and (c) to obtain a definitive set of minerals.

,

  Davies model, the activity coefficients are given in logarithmic form by:

  the description of the chemical test case to handle a representative range of concentrations.

Figure 1 :

 1 Figure 1: Inner fixed-point, full Newton and outer fixed-point activity algorithms

  80.90 ). Highly charged species such as Cr 3+ or Fe 2 (OH) 24+ will be challenging for activity coefficients algorithms. It includes 7 components and 39 species. A table including the stoichiometric coefficients, equilibrium constants, total concentrations and equilibrium solutions is given in appendix A-5.

Table 1 :

 1 46. Wanner C, Eggenberger U, Kurz D, et al (2012) A chromate-contaminated site in southern Switzerland -Part 1: Site characterization and the use of Cr isotopes to delineate fate and transport. Appl Geochem 27:644-654. http://dx.doi.org/10.1016/j.apgeochem.2011.11.009 47. Richard FC, Bourg ACM (1991) Aqueous geochemistry of chromium: A review. Water Res 25:807-816. https://doi.org/10.1016/0043-1354(91)90160-R Chemical table for the test case with only activity correction.

							Cl -	Ca 2+	Al 3+	Sn 4+	K
						Cl -	1	0	0	0	1
						Ca 2+	0	1	0	0	1
						Al 3+	0	0	1	0	1
						Sn 4+	0	0	0	1	1
	TOTAL (mol.l -1 )	9.10 -5	10 -5	10 -5	10 -5
	Initial low I (mol.l -1 )	5.10 -7	5.10 -7	5.10 -7	5.10 -7
	Initial high I (mol.l -1 )	5.10 -4	5.10 -4	5.10 -4	5.10 -4
	C	i		1 Nx j	, i j j a X (mol.l -1 )	10 -9	10 -9	10 -9	10 -9
	j T	Nc	, i j i a C (mol.l -1 )	5.10 -1	5.10 -1	5.10 -1	5.10 -1
			i	1

Table 2 :

 2 Initial values for the situation with a low ionic strength.

			Nx
				a
				, i j
		C K X			Nc	Nc
		i	i	j		1		2	1	2
						I		z C	I	z C
			j	1				i i	i i
					C C	2	i	1	2	i	1
					i i i
		(mol.l -1 )		Fixed-point	Full Newton
	Cl -	5.02 10 -7	0.997	8.41	8.41
	Ca 2+	5.07 10 -7	0.987	19.00	18.58
	Al 3+	5.15 10 -7	0.970	19.00	18.06
	Sn 4+	5.27 10 -7	0.948	19.00	17.34

Table 3 :

 3 Initial values for the situation with a high ionic strength.

			Nx
				a
				, i j
		C K X			Nc	Nc
		i	i	j		1		2	1	2
						I		z C	I	z C
			j	1				i i	i i
					C C	2	i	1	2	i	1
					i i i
		(mol.l -1 )		Fixed-point	Full Newton
	Cl -	7.12 10 -4	0.702	-0.82	-1.07
	Ca 2+	2.06 10 -3	0.243	-0.98	-1.78
	Al 3+	1.20 10 -2	4.16 10 -2	-0.98	-2.75
	Sn 4+	1.42 10 -1	3.51 10 -3	-0.98	-4.12

Table 4 :

 4 CPU times for 1 loop (ms) Outer fixed-point Inner fixed-point Full explicit Full recursive

	Accepted Article

This article is protected by copyright. All rights reserved.

Table 5 :

 5 Failure ratio after 1,000 iterations

		Only-I	Phosphoric acid Gallic acid	Iron-Chromium Morel-Morgan
	Inner fixed-point	0.00 %	0.00 %	0.00 %	1.62 %	17.91 %
	Outer fixed-point	0.00 %	0.00 %	0.00 %	0.00 %	0.00 %
	Full Newton	0.00 %	0.00 %	0.05 %	91.35 %	3.87 %

This article is protected by copyright. All rights reserved.

Acknowledgements

The authors acknowledge the French programme NEEDS for its financial support to the project NewSolChem. We greatly thanks the anonymous reviewers and Michel Kern for their helpful comments.

Phosphoric acid test case

. As said, this test case is quite easy to solve. We find that all the algorithms reach 100 % resolution (Table 5) in fewer than 70 iterations. Nevertheless, the presence of highly charged species such as PO [START_REF] De Windt | Intercomparison of reactive transport models applied to UO2 oxidative dissolution and uranium migration[END_REF] 3-allows to differentiate the behavior of the algorithms. The frequency graph shows that the full Newton algorithm usually converges within the fewest iterations and that the inner fixed-point algorithm requires the most iterations to converge. Nevertheless, the outer fixed-point algorithm provides very interesting results for this test case: it converges within 36 and 39 iterations regardless of the initial guess.

Gallic acid test case

The results for the gallic acid test case (Figure 5) confirm those of the phosphoric acid test case. The inner fixed-point algorithm requires the most iterations to converge (approximately 150). The full Newton algorithm usually requires fewer iterations (approximately 130) but sometimes requires many more iterations (200); sometimes it fails to converge at all (0.05 % failure rate, see Table 5).

The outer fixed-point algorithm converges with the fewest iterations (77 iterations maximum). Moreover, this algorithm gives the sharpest frequency graph, indicating that it is not sensitive to the initial guess.

Iron-chromium test case

The iron-chromium test case is the strongest test used here. It is reported [START_REF] Marinoni | Thermodynamic equilibrium solutions through a modified Newton Raphson method[END_REF][START_REF] Machat | Comparison of linear solvers for equilibrium geochemistry computations[END_REF] to generate very high condition numbers. The range of initial guesses is chosen with the goal of keeping this test reasonable. Nevertheless, one can see (Figure 6) that it is very hard for the full Newton algorithm to converge, with a 91.35 % failure rate (Table 5). Notably, this full Newton algorithm sometimes results in faster resolutions (according to the number of Newton iterations) because it is the only algorithm that sometimes converges with fewer than 200 Newton iterations. Both the inner and outer fixedpoint algorithms converge most frequently within 230-240 iterations. The inner fixed-point algorithm reaches 80 % of its realizations after 240 iterations but needs up to 1000 iterations to complete the set and fails to converge at a rate of 1.62 % (Table 5). The outer fixed-point algorithm requires between 226 and 233 iterations to converge regardless of the initial point. Moreover, it always succeeds in solving this test case.

Morel-Morgan test case

The Morel-Mogran test case is the largest test used here. It can be seen on Figure 7 that the inner fixed-point algorithm solves quite fast 40 % of the simulations using less than 200 iterations each. But this algorithm needs more and more iterations to solve 40 % of the other simulation and fails to converge for 17.91 % of the cases (Table 5). The Full Newton algorithm produces more regular results succeeding to solve easiest simulation within less than 100 iteration and using all the 1000 iterations to solve some of the hardest simulations. It fails to converge at the rate of 3.87 % (Table 5). The Outer fixed point algorithm is very efficient for this test case too, solving all the cases within less than 300 Newton iteration each.

Accepted Article

This article is protected by copyright. All rights reserved.

Comparison of CPU times and failure ratios

As expected, the CPU times for one iteration loop increase as the number of species increases in the test case and as the complexity of the algorithm increases (see Table 4). We show that the outer fixed-point algorithm is the fastest and that the full Newton algorithm is the slowest. Nevertheless, CPU times for the outer and the inner fixed-point algorithms are very close and their differences can be considered as negligible. Moreover, the recursive formulation of the full Newton algorithm is the slowest. For this reason, we strongly recommend not using the recursive formulation, and we prefer the explicit formulation.

We find that the only algorithm that solves all the test cases with a 100 % success rate is the outer fixed-point algorithm.

CONCLUSION

Algorithms to handle activity correction into geochemical models have been poorly mentioned in the literature. In this work, we extensively described the 2 algorithms we found mentioned and presented a new one. We compared 3 algorithms based on their ability to handle activity correction in equilibrium chemistry solvers. They were compared on 5 chemical test cases selected according to their increasing numerical difficulty and/or size

The full Newton algorithm is the most integrated algorithm from a mathematical point of view.

Nevertheless, we found it to be the slowest and weakest algorithm. We suppose this algorithm increases the nonlinearity of the chemical system by injecting activity corrections into the mass action equations and conservation laws. It increases the condition number of the Jacobian matrix, as shown by comparing [START_REF] Lützenkirchen | The constant capacitance model and variable ionic strength: An evaluation of possible applications and applicability[END_REF] and [START_REF] Marinoni | Thermodynamic equilibrium solutions through a modified Newton Raphson method[END_REF]. It has been shown [START_REF] Marinoni | Thermodynamic equilibrium solutions through a modified Newton Raphson method[END_REF][START_REF] Machat | Comparison of linear solvers for equilibrium geochemistry computations[END_REF] that a condition number that is too high leads to inaccurate steps   [START_REF] Steefel | Reactive transport codes for subsurface environmental simulation[END_REF] in the Newton methods, leading to numerical difficulties or nonconvergence. Because chemical equilibrium computation is still a highly nonlinear problem, increasing its nonlinearity by injecting activity correction seems to be an inefficient choice.

The inner fixed-point algorithm includes an intermediate integration of activity correction into the Newton loop. Both loops, Newton for the mass action equations and conservation laws and fixedpoint for activity correction, run together. In this way, changes induced by activity correction disturb the Newton minimization. This point explains the convergence difficulties of the inner fixed-point algorithm when activity correction becomes important.

The outer fixed-point algorithm proposes a complete separation between the Newton and activity correction loops. By this way, nonlinearity induced by activity correction cannot disturb the Newton convergence, and the condition number of the Jacobian matrix is lower than that obtained by the full Newton algorithm. This leads to a more stable and robust algorithm. We found that the outer fixed-

Accepted Article

This article is protected by copyright. All rights reserved.


N Cavassilas

M Lannoo

M Luisier

Y Lee

M Bescond
email: marc.bescond@iis.u-tokyo.ac.jp

D Logoteta

Anharmonic phonon-phonon scattering for modeling of three-dimensional atomistic transport: an original quantum treatment

Anharmonic phonon-phonon scattering modeling of three-dimensional atomistic transport: An efficient quantum treatment

I. INTRODUCTION

The semiconductor fabrication technology exhibits nowadays the capability of manipulating materials down to the atomic level, therefore offering the opportunity to realize well-defined nanoscale structures [1][2][3][4][5] . Nanoscale systems are very attractive, since their use has been recognized as a powerful option to promote technology advances in several different application fields, ranging from nanoelectronics (e.g. transistors [START_REF] Suk | IEEE International Electron Devices Meeting[END_REF][START_REF] Mertens | IEEE Symposium on VLSI Technology[END_REF] and tunnel diodes [START_REF] Schmid | [END_REF] ) and nanophotonics (e.g. photovoltaics 9,10 ) to energy harvesting devices (e.g. thermoelectric devices 11,12 ).

The development of nanoscale devices requires the modeling, on a quantum mechanical basis, of the physics underlying the transport processes. Large efforts have been traditionally devoted to investigate and accurately describe electron transport mechanisms [START_REF] Mahan | Many-Particle Physics[END_REF][START_REF] Datta | Electronic Transport in Mesoscopic Systems[END_REF][START_REF] Haug | Quantum Kinetics in Transport and Optics of Semiconductors[END_REF][START_REF] Ferry | Transport in Nanostructures[END_REF][START_REF] Zhao | [END_REF][18][19] . More recently, however, an increasing number of researches [20][21][22][23][24][25] have emphasized the fundamental role played by thermal effects at the nanometric scale. Self-heating, for example, has been shown to considerably reduce performance 26 and lifetime [START_REF] Srinivasan | International Conference on Dependable Systems and Networks[END_REF] of nanoscale electronic and optoelectronic devices. Furthermore, nanostructures have been recognized as promising building blocks for thermoelectric applications, due to their small thermal conductivity [START_REF] Hicks | [END_REF] . Despite this increasing interest, the understanding of nanoscale quantum phonon transport is still lacking compared to the electron one, and further studies are needed to gain a detailed physical picture, particularly for what concerns the role of anharmonic phonon-phonon interactions.

Nanoscale thermal transport has been investigated by adopting several approaches (e.g. the linearized Boltzmann transport equation based on the Fermi golden rule 22,[29][30][31][32] , (non)-equilibrium molecular dynamics [33][34][35][36][37] and first-principles perturbation theory [38][39][40] ), able to capture quantum mechanical effects within different degrees of approximation.

Among them, the nonequilibrium Green's function (NEGF) framework has attracted wide interest, as it offers the opportunity to perform full quantum simulations, including many-body effects (e.g. anharmonic phonon interactions 20,21,25,41,42 ) via suitable self-energies.

From the numerical point of view, the main bottleneck in applying the NEGF formalism, usually comes from the iterative solution within the self-consistent Born approximation (SCBA) of the non-linear Dyson equation. This task, particularly heavy in the simulation of threedimensional realistic nanostructures, requires performing a large number of iterations in order to obtain accurate results for the physical quantities of interest. This is partly due to the intrinsic convergence properties of the SCBA iterative scheme, which provides a self-energy fulfilling the conservation laws (in particular, a Φ-derivable self-energy [43][44][45][46] ) only asymptotically.

In the last decade, several works, focusing on the modeling of the electron-phonon interaction, have investigated the possibility of avoiding such an expensive iterative scheme [47][48][49][50][51][52][53][54] . In this paper, we tackle this problem by considering a different kind of self-energy, which describes the dominant phonon-phonon anharmonic interaction 41 and is not a linear functional with respect to the phonon Green's function. We will show that the so-called Lowest Order Approximation (LOA) 49,50,52,54 , combined with a rescaling technique and the Padé approximation, can also be effectively used in this case to compute the chargeconserving expectation values of a generic linear operator at any order. Importantly, the specific case considered turns out to be just an example of a very wide class of self-energies to which the method is applicable.

In order to assess the numerical efficiency of the method against the SCBA, we simulated the anharmonic phonon transport in Si and Ge nanowires (NWs). They are a relevant example of 3D nanostructures with potentialities in a variety of application fields where the thermal transport plays a relevant role and needs to be accurately modeled. Several configurations have been considered, differing in the NW crystallographic orientation, the value of the temperature and the regularity of the NW surface. The atomistic simulations have been carried out within the NEGF framework by using the OMEN 25 code, which implements a modified-valenceforce field (VFF) model to describe the interatomic potential.

As key findings, we report that for all the considered structures at room temperature the third-order thermal currents obtained from LOA diagrams, with an analytical continuation of Padé approximants, can reconstruct the SCBA values and lighten the SCBA computational burden. For anharmonic effects at high temperature (500 K -900 K), it turns out that higher-order LOA thermal currents are required. In such a case, we suggest as a solution the application of the first-order Richardson extrapolation [55][START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory[END_REF][START_REF] Brezinski | [END_REF] to the sequence of the Padé approximants N -1/N . The rest of the paper is organized as follows. In Sec. II, we first describe the NEGF model for the anharmonic phonon transport. We then introduce the LOA rescaling technique for the case of the considered phononphonon interaction self-energy, and the application of the Richardson extrapolation to the Padé approximants N -1/N . In Section III, we present the simulation results obtained within the LOA and systematically compare them with the corresponding obtained within the SCBA. Conclusions and perspective for future developments are followed in Sec. IV.

II. THEORY

A. Phonon transport model

The phonon Hamiltonian in three dimensions reads as

H ph = l1 µ P 2 µ 2M µ + 1 2 l1,l2 µ,ν ûl1 µ ûl2 ν • ∂ 2 ∂R l1 µ ∂R l2 ν V harm + 1 3! l1,l2,l3 µ,ν,σ ûl1 µ ûl2 ν ûl3 σ • ∂ 3 ∂R l1 µ ∂R l2 ν ∂R l3 σ V anh + O(û 4 µ ), (1)
where the first, second, and third term refers to the kinetic, harmonic, and anharmonic contributions, respectively. The indices µ, ν, and σ denote the atomic postions, while the indices l 1 , l 2 , and l 3 run over the x, y, and z directions. R µ , M µ , and Pµ are the position, the mass, and the momentum operator, respectively, of the atom at the position µ. The operator ûl1 µ is the quantized displacement of the atom at the position µ from the equilibrium position along the direction l 1 . In this work, V harm and V anh are the VFF harmonic 58,59 and anharmonic 25 potential energies, respectively. Free boundary conditions are enforced on the external surface of the nanowire.

The second derivative of the VFF harmonic potential energy with respect to the atomic position µ in the l 1 direction and to the atomic position ν in the l 2 direction corresponds to the entries of the dynamical matrix Φ

Φ l1l2 µν = ∂ 2 ∂R l1 µ ∂R l2 ν V harm . (2)
The retarded (D R (ω)) and greater/lesser (D ≷ (ω)) phonon Green's functions satisfy the equations:

l1 µ M ν ω 2 δ l2l1 δ νµ -Φ l2l1 νµ -Π RB, l2l1 νµ (ω) -Π RS, l2l1 νµ (ω) D R, l1l3 µσ (ω) = δ l2l3 δ νσ , D ≷, l2l3 νσ (ω) = l 1 ,l 1 µ,µ D R, l2l1 νµ (ω) Π ≷B, l 1 l 1 µµ (ω) + Π ≷S, l 1 l 1 µµ (ω) D A, l 1 l3 µ σ (ω) (3)
where

D A (ω) = D R (ω)
† is the advanced phonon Green's function, ω is the phonon frequency, Π RB (ω) (Π ≷B (ω)) is the retarded (greater/lesser) open boundary self-energy 60 , and Π RS (ω) (Π ≷S (ω)) is the retarded (greater/lesser) anharmonic phonon-phonon scattering self-energy 25 . The kinetic and harmonic terms of the phonon Hamiltonian are therefore directly included in the phonon Green's function, while the anharmonic term is perturbatively treated by defining a scattering selfenergy.

In this work, we model the anharmonic interactions in terms of the lowest order three-phonon process depicted in Fig. 1. It consists in the decay of a high-energy phonon (ω + ω ) into two lower energy phonons (ω and ω ) or vice versa 25,61 . The corresponding greater/lesser self-energy turns out to be Φ-derivable [START_REF] Stefanucci | Nonequilibrium Many-Body Theory of Quantum Systems[END_REF] and reads as 25 :

Π ≷S, l 2 l 3 νσ (ω) = 2i l 1 ,l 1 ,l 1 ,l 1 µ ,µ ,µ ,µ ∞ -∞ dω 2π dV (3)l 2 l 1 l 1 νµµ × D ≷, l 1 l 1 µµ (ω + ω ) D ≶, l 1 l 1 µ µ (ω ) × dV (3)l 1 l 1 l 3 µ µ σ , (4)
with

dV (3)l 2 l 1 l 1 νµµ = ∂ 3 ∂R l 2 ν ∂R l 1 µ ∂R l 1 µ V anh . (5)
The retarded scattering self-energy can be derived from the lesser and greater self-energies as

Π RS, l 2 l 3 νσ (ω) = 1 2 Π >S, l 2 l 3 νσ (ω) -Π <S, l 2 l 3 νσ (ω) + iP dω 2π Π >S, l 2 l 3 νσ (ω ) -Π <S, l 2 l 3 νσ (ω ) ω -ω , (6)
where P is the Cauchy principal integral value. In order to exploit the efficient recursive algorithm of Ref. 63 to solve Eqs. (3), we resort to a diagonal approximation 25,54 of the scattering self-energy with respect to the atomic position indexes ν and σ. Although this approximation may lead to an underestimation of the interaction strength, in our model this is partly compensated through the calibration of the model parameters for the anharmonic force against available experimental data 25 . As a further approximation, we neglect the real part of the retarded self-energy, as it mainly contributes to an energy renormalization [START_REF] Svizhenko | [END_REF][65][START_REF] Frey | Proceedings of the 38th European Solid-State Device Research Conference (ESSDERC)[END_REF] . Finally, the phonon frequency shift due to the temperature-dependent lattice expansion is not taken into account.

The self-consistent solution of Eqs. ( 3), ( 4) and ( 6) provides the interacting Green's functions within the SCBA. The local thermal current I p,p+1 ph between two unit cells p and p + 1 can then be calculated as

I p,p+1 ph = 2π ν∈p σ∈p+1 l2,l3 × ∞ 0 dω ω Φ l2l3 νσ D <,l3l2 σν -D <,l2l3 νσ Φ l3l2 σν , (7)
where the atoms at the positions ν and σ belong to the unit cells p and p + 1, respectively. A unit cell is composed of N consecutive atomic layers along the transport direction, with N = 4 for the 100 direction and N = 2 for the 110 one.

B. Lowest Order Approximation

The approach routinely used in NEGF numerical simulations to compute an interacting Green's function D consists in iteratively solving, within the SCBA, the Dyson equation 13

D = d 0 + d 0 Π[DD]D, (8)
where d 0 denotes the non-interacting Green's function, Π is the self-energy that accounts for the interactions, and the simplified matrix notation [i = (r i , ω i ), d 0 ΠD = d2 d2 d 0 (1; 2) Π(2; 2 )D(2 ; 1 )] has been used. The product DD in the argument of Π highlights that the considered self-energy (4) depends on the square of the Green's function. The iteration scheme to solve Eq. ( 8) can be expressed as

D N = [d -1 0 -Π[D N -1 D N -1 ]] -1 , (9)
where D N denotes the Green's function approximation at the N th iteration and D 0 = d 0 . By performing a Taylor series expansion in powers of Π[D N -1 D N -1 ], we can also rewrite Eq. ( 9) as

D N = d 0 + d 0 Π[D N -1 D N -1 ]d 0 + d 0 Π[D N -1 D N -1 ]d 0 Π[D N -1 D N -1 ]d 0 + • • • . (10)
Moreover, by recursively expanding

D N -1 in terms of Π[D N -2 D N -2 ], Π[D N -3 D N -3 ], • • • , Π[d 0 d 0 ]
, we can rearrange the right hand side of Eq. ( 10) in the form of products and compositions of the basic term d 0 Π[d 0 d 0 ]:

D N = d 0 + d 0 Π[d 0 d 0 ]d 0 + d 0 Π[d 0 d 0 ]d 0 Π[d 0 d 0 ]d 0 + • • • . ( 11
)
The second term on the right hand side of Eq. ( 11) corresponds to the diagram of Fig. 1, while the other terms correspond to higher order diagrams. Here and in the following, the order of a term/diagram is intended as its order in Π, namely as the number of occurrences of Π in that term/diagram. For each N , the expansion (11) can be split in two parts: a perturbation series truncated at the order N , which includes all the possible (non-crossing) diagrams of order smaller or equal to N , and the remaining infinite sum, which misses some diagrams at any order. It can be shown 50 that, if the self-energy is Φ-derivable, the perturbation series part always provides a conserving approximation of the Green's function. On the contrary, due to the presence of the renormalizing higher order terms, the SCBA provides a conserving approximation only asymptotically, when the self-consistency is achieved. Differently from the SCBA, the LOA provides a sequence of expectation values of a generic linear operator by directly relying on the Green's function perturbation series. The N th element of the sequence, which is referred to as LOA of order N , is computed by approximating the Green's function with its perturbation series truncated at the order N . According to the considerations above, each element of the sequence fulfills the conservation laws. However, as the perturbation series has a finite convergence radius 53 , the sequence of expectation values can also diverge. In order to turn this sequence into a convergent one, analytical continuation techniques can be resorted.

By using these techniques, higher order effects missed in the truncated perturbation series can be included, therefore obtaining a renormalized approximation of the LOA expectation values. These latter can still fulfill the conservation laws, if the analytical continuation preserves the conserving character of the LOA. Essentially, while both the SCBA and the analytically continued LOA implement a convergent iterative scheme, only the latter provides conserving results at each iteration. In the context of electron-phonon interaction, this difference has been shown to result in a significant acceleration of the convergence of the LOA with respect to the SCBA 52,54 . In the following, we will show how, by using a rescaling technique, the LOA expectation values can be obtained from the non-conserving estimates provided by the SCBA during the iteration process.

The idea of the rescaling technique is to precondition the expression of Π[D N -1 D N -1 ] in order to obtain, as the result of the SCBA iteration, a "rescaled" version of the expansion (11), in which terms corresponding to diagrams of increasing order are divided by increasing powers of the same scaling factor λ. A rescaled expression for D N , in which terms of order i are divided by λ i , can be easily obtained by using in the SCBA algorithm a version of Π[D N -1 D N -1 ] in which any occurrence of the functional Π is divided by λ. If λ is suitably chosen, the nonconserving terms of order higher than N will be strongly suppressed relative to the conserving ones, thus making it possible to isolate them and build the corresponding LOA. We proceed now to show how the rescaling technique can be implemented in practice, by describing its application to the first two SCBA iterations.

By rescaling Π[d 0 d 0 ] by λ 1 , the first SCBA iteration gives

D 1 = d 0 + d 0 Π[d 0 d 0 ]d 0 λ 1 + d 0 Π[d 0 d 0 ]d 0 Π[d 0 d 0 ]d 0 λ 2 1 + • • • . (12)
For λ 1 large enough, the non-conserving terms with order higher than the first are strongly suppressed, and the LOA expectation value of a linear operator O can be extracted as

O 1 = O 0 + ∆O 1 = O 0 (d 0 ) + O(d 0 Π[d 0 d 0 ]d 0 ) = O 0 (d 0 ) + λ 1 [O(D λ1 1 ) -O(d 0 )],
where D λ1 1 denotes the first-order phonon LOA Green's function rescaled by λ 1 :

D λ1 1 = d 0 + d 0 Π λ1 1 d 0 , Π λ1 1 = Π [d 0 d 0 ] λ 1 , (13)
To proceed with the computation of the second order LOA, we build the rescaled version Π λ2 2 of the phonon self-energy Π[D 1 D 1 ]:

Π λ2 2 = 1 λ 2 Π D λ2 1 D λ2 1 , = 1 λ 2 Π (d 0 + 1 λ 2 ∆d 1 )(d 0 + 1 λ 2 ∆d 1 ) , = 1 λ 2 Π [d 0 d 0 ] + 1 λ 2 2 Π [2d 0 ∆d 1 ] + 1 λ 3 2 Π [∆d 1 ∆d 1 ] , (14)
where ∆d

1 = d 0 Π 1 d 0 = d 0 Π [d 0 d 0 ]
d 0 and the linearity of Π has been exploited. It can be checked that terms of order i are divided by λ i 2 . This remains true for any polynomial in Π λ2 2 , and thus for the rescaled Green's function D λ2

2 provided by the second SCBA iteration. Again, if the scaling factor λ 2 is appropriately chosen 54 , the (nonconserving) terms of order higher than the second are strongly suppressed, and the second-order LOA expectation value can be obtained as

O 2 = O 0 + ∆O 1 + ∆O 2 , with ∆O 2 = λ 2 2 [O(D λ2 2 ) -O(d 0 )] -λ 2 ∆O 1 .
This strategy can be generalized to the computation of the phonon Green's function LOA at any order. In particular, from a Green's function exact to the (N -1) th order, rescaled by a factor λ N

D λ N N -1 = d 0 + 1 λ N ∆d 1 + • • • + 1 λ N -1 N ∆d N -1 , (15)
where ∆ i denotes the term of order i of the Green's function perturbation series, we can obtain a rescaled self energy Π λ N N exact to the N th order as

Π λ N N = Π D λ N N -1 D λ N N -1 = Π       d 0 d 0 λ N + • • • + N -1 i=0 ∆d i ∆d N -1-i λ N N + N -2 i=0 ∆d i+1 ∆d N -1-i λ N +1 N + • • •       . ( 16
)
In this paper, a conserving analytical continuation technique relying on the Padé approximation is used 50,52 . We compute the LOA expectation values of the thermal current up to the third order, as our results indicate that the Padé approximants 0/1, 1/1 and 1/2 built from them can already provide good estimates of the SCBA converged current.

The rescaling technique has been presented for the case of the self-energy (6), which depends on the square of the phonon Green's function. However, it is easy to argue that it can be used exactly in the same form in the case of the composition Σ[G k ] of a linear self-energy Σ with any positive integer power k of the Green's function G. By going a small step further, we observe that, due to the linearity of Σ, the method can be applied to selfenergies of the form Σ[

M i=0 f i G k(i) ] = M i=0 f i Σ[G k(i) ],
where the function k associates i a positive integer. If the functional Σ is bounded, the previous equality holds also for M = ∞, and, by taking k(i) = i, we can conclude that the rescaling technique is applicable to all the conserving self-energies of the form Σ[f (G)], where f is an analytic function.

These considerations show that the case of the phononphonon scattering self-energy on which the paper focuses, while important, represents just a very simple example among the wide range of physical situations the method can deal with. The achievable numerical efficiency will obviously depend on the particular case at hand.

C. Richardson extrapolation

When the anharmonic scattering is strong, as in the case of temperatures much higher than room temperature (see below, subsection III B), the convergence of the series of the Padé approximants can be slow. In such cases, we propose to apply the Richardson extrapolation to the LOA expectation values analytically continued by the Padé approximants N -1/N .

The Richardson extrapolation is a method to accelerate the convergence of monotone series. If a sequence S N approaches monotonically a certain value S, we can write S N in an asymptotic form as 56

S N ≈ S + c 1 N + c 2 N 2 + c 3 N 3 + • • • , (17)
where c 1 , c 2 , • • • are unknown parameters. By considering only the first correction term for S N and S N +1

S N = S + c 1 N , S N +1 = S + c 1 N + 1 , (18)
the first-order Richardson extrapolation can be expressed as

S [1] N = (N + 1)S N +1 -N S N . (19)
We used this technique to accelerate the convergence of the sequence O N -1/N of the LOA expectation values of an operator O, analytically continued by the Padé approximants N -1/N . If, as N → ∞, O N -1/N converges monotonically, we can set S N = O N -1/N and apply Eq. ( 19). In the paper, the Richardson extrapolation is applied to the sequence elements O 0/1 and O 1/2 as

S [1] 1 = 2S 2 -S 1 = 2O 1/2 -O 0/1 , (20)
to obtain an improved estimate S

[1]

1 of the corresponding SCBA result.

III. RESULTS AND DISCUSSIONS

In this section, we present and discuss the simulation results for the thermal current in Si and Ge NWs, obtained by using either the LOA analytically continued with Padé approximants and possibly further improved by using the Richardson extrapolation or the SCBA. The SCBA results have been obtained within a tolerance of 1% on the current conservation. All the NWs are undoped. The resulting low electron density allows us to neglect the electron contribution to the thermal current and the electron-phonon scattering. In all the considered cases, the thermal current is sustained by a temperature gradient of 0.1 K, enforced at the NW extremities.

In subsection III A, we consider the case of NWs with a uniform cross-section at room temperature. In subsection III B, we discuss the results obtained by increasing the temperature up to 900 K. In these conditions, the increase in the phonon population entails a significant enhancement of the phonon-phonon scattering rate. Finally, as an example of geometry-driven mode-mixing enhancement, in subsection III C we address the case of NWs with discontinuous cross-section.

A. Uniform cross-section nanowires at room temperature

The atomic structure of the considered Si and Ge NWs is depicted in Fig. 2. The nanowires are 60 nm long and have a circular cross-section with a diameter of 3 nm. The Ge NW is oriented in the 110 direction, while for the Si NW both the 110 and 100 directions are considered. In the following, we refer as "ballistic" the results for the thermal current obtained with uncoupled phonon modes, namely in the harmonic approximation.

Fig. 3 (a) illustrates the values of thermal current along the 110 -oriented Si NW, in the ballistic regime and in the presence of anharmonic phonon-phonon scattering. For this latter case, the SCBA converged value and the LOA approximations up to the third order are reported. The results confirm that the rescaling technique is actually able to provide current-conserving expectation values at any order. However, the LOA power series diverges, due to the already significant anharmonic scattering strength at room temperature 25 . Analogously to the case of the previously investigated electron-phonon selfenergy 52,54 , an analytical continuation of the LOA series based on the Padé approximation has been resorted to obtain a convergent behavior. Fig. 3 (b) shows the Padé approximants built from the LOA values. The sequence of Padé approximants 0/1, 1/1 and 1/2 oscillates around the SCBA value, the best approximation of which is obtained by the 1/2 approximant with a relative error ε (= 100 × |I SCBA -I Pade |/I SCBA ) lower than 9.2%. This result is achieved with around half of the computational burden needed for the SCBA (6 against 10 SCBA iterations). The same behavior is observed for the 100oriented Si NW (not shown). In this case, the divergence of the LOA series is faster, which suggests stronger anharmonic interactions. The Padé approximant 1/2 is also in this case in good agreement with the SCBA value (ε = 10%).

Figure 4 illustrates the corresponding results obtained for the 110 -oriented Ge NW. The divergence of the LOA series is slower than for the Si NWs, which indicates a higher ballisticity (the lower thermal conductivity of Ge NWs 40,67-69 with respect to Si NWs is a consequence of the lower group velocity of the acoustic phonon branches 25 , which mainly sustain the thermal current). Figure 4 (b) shows that, once again, the 1/2 Padé approximant provides a good estimate of the SCBA current (ε = 10%), while resulting in an improvement of the numerical efficiency similar to that found for the Si NWs.

B. Uniform cross-section nanowires at higher temperatures

A different scenario is observed when temperatures much higher than 300 K are considered. Table I reports the LOA thermal current analytically continued by the 0/1 and 1/2 Padé approximants for the 110 -oriented Si NW considered in subsection III A, and the corresponding values obtained through a first-order Richardson extrapolation. The SCBA and ballistic current, and the number of SCBA iterations needed to obtain each result are also shown. The data are provided for temperatures ranging from 300 K to 900 K. It can be noticed that the current degradation with respect to the ballistic case δI = 100 × [1 -I SCBA /I bal ], where I SCBA and I bal are the SCBA and ballistic value of the current, respectively, significantly increases with the temperature: δI = 46% at 300 K, 64% at 500 K, 74% at 700 K and 81% at 900 K. This is a consequence of the enhancement of the anharmonic interactions at high temperature 25 . The accuracy of the Padé approximant 1/2 in retrieving the SCBA result progressively degrades with increasing temperature (ε = 17.4% at 500 K, 20.8% at 700 K, and 29.2% at 900 K), suggesting that higher-order LOA terms are needed. However, it can be seen that in all the cases the Richardson extrapolation obtained from the Padé approximants 0/1 and 1/2 significantly improves the estimate of the SCBA result (we find ε = 4.0% at 300 K, 0.0% at 500 K, 1.3% at 700 K, and 11.3% at 900 K). We remark that the use of the Richardson extrapolation does not require further iterations and therefore does not increase the computational burden. In this way, it is possible to considerably boost the numerical efficiency of the LOA with respect to the SCBA scheme. For example, Table I shows that approximating the anharmonic thermal current at T = 900 K within a relative error close to 10% is more than 7 times faster than computing the corresponding SCBA result.

We tested the first-order Richardson extrapolation also in the case of the self-energies describing electron-phonon interactions considered in Ref. 54. In all the cases in which the sequence of the N -1/N Padé approximants converged monotonically, we again obtained a significant acceleration.

C. Nanowires with discontinuous cross-section

The atomic structure of the NWs with discontinuous cross-section is sketched in Fig. 5. We consider a 100oriented Si NW with L NW = 60 nm and d NW = 4 nm at room temperature. The choice of the 100 orientation allows us to discuss a different convergence behavior, which we did not observe in the 110 orientation. The NW in Fig. 5 (a), which we will refer to as NWR1, has one recessed region of length L Rec = 10 nm. The recess is obtained by removing the outermost layer of atoms. The length of the left/right unrecessed parts is L Left = L Right = 25 nm. The NW in Fig. 5 (b), which we will refer to as NWR2, has two recessed regions of length L Rec1 = L Rec2 = 5 nm. The lengths of the left, center, and right unrecessed regions are L Left = 12.5 nm, L Center = 25 nm, and L Right = 12.5 nm, respectively. These structures can be considered as elemental cells from which one-dimensional nanoscale phononic crystals with tailored phonon band structures 70 could be assembled. We remark that, in our model, the discontinuities of the NW cross-section affect the self-energy in a qualitatively different way with respect to an increase of the temperature. This latter induces an increase of the anhar- monic scattering rate essentially by modifying the noninteracting phonon Green's function. Differently, in the case of geometrical irregularities, this increase is mainly due to a change in the interaction strength, through the derivatives of the anharmonic components of the potential energy (see Eq. ( 5)). Figure 6 illustrates the results obtained for the structure NWR1. Contrary to the cases previously discussed in subsections III A and III B, it can be noticed that the SCBA and the first-order LOA thermal currents are higher than the ballistic one. This behavior is similar to the one observed in the n-type 111 -oriented Si NWs in the case of electron-phonon scattering (see Table I in Ref. 54). Physically, the increase of the current in the presence of anharmonic scattering originates from the occurrence, in ballistic conditions, of localized phonon modes, induced by the cross-section nonuniformity. In the presence of mode-mixing, these modes can contribute to transport, therefore entailing an increase of the current. Figure 6 (b) shows that the best approximation of the SCBA current is provided by the 1/1 Padé approximant, able to almost perfectly recover the SCBA value (ε = 0.1%). This could suggest that in similar cases the sequence of N/N Padé approximants converges faster than the N -1/N one. The LOA approximation turns out to be one order of magnitude faster, as only 3 SCBA iterations are needed to compute the 1/1 Padé approximant, against the 32 iterations required to obtain the converged SCBA result. This behavior is of course not always observed for NWs with surface recesses. Figure 7 shows the results for the slightly more complex NWR2 structure. In this case, the anharmonic scattering induces a decrease of the thermal current, and the behavior is the same as that observed for NWs with uniform cross-section. Accordingly, the best approximation of the SCBA thermal current is provided by the 1/2 Padé approximant. The relative error with respect to the SCBA result is ε = 1.4%, and the numerical efficiency is improved by more than 3 times.

It is interesting to compare the results for the structure NWR1 with those obtained for the NW with uniform cross-section at T = 700 K. We observe that, while the total number of SCBA iterations slightly increase from 28 to 32, the relative error of the best Padé approximant with respect to the SCBA significantly improves from 20.8% to 1%. The different behavior can be attributed to the local character of the NW cross-section discontinuities. A localized scattering center is less efficient in triggering high order scattering processes with respect to a homogeneous increase of the temperature. This explains the higher accuracy of the analytically continued LOA in the case of the NWR1 structure. On the other hand, within the SCBA scheme, localization effects typically result in significant spatially localized fluctuations of the current, which need numerous iterations to be canceled out. The SCBA scheme cannot benefit from the same improvement as the LOA, because its numerical efficiency is limited by the non-automatic current conservation.

IV. CONCLUSION

We investigated the opportunity to adopt the LOA approach as a more efficient option with respect to the SCBA in the description of the anharmonic phononphonon scattering in thermal transport NEGF simulations. We found that the LOA, in the rescalingtechnique-based implementation, is able to provide the correct thermal current-conserving expectation values at any order. This method can be generalized to any conserving self-energy that analytically depends on the Green's function. We benchmarked the numerical efficiency of the LOA approach against the SCBA through atomistic phonon transport simulations in Si and Ge NWs. For constant cross-section NWs at room temperature, our results indicate that the Padé approximants built with the LOA currents up to the third order can provide high accuracy and a significant computational advantage with respect to the SCBA (typically, a reduction of the numerical burden by a factor 2). The largest improvement (up to a factor 10) of the computational efficiency relative to the SCBA has been observed in NWs with discontinuous cross-section. A similar improvement can be obtained also at temperatures substantially higher than the room temperature, if the first-order Richardson extrapolation is applied to the Padé N -1/N sequence in order to accelerate the convergence.

V. ACKNOWLEDGEMENT

This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) and funded by the NOODLES contract (ANR-13-NANO-0009) of the ANR-French National Research Agency. The author Y. Lee would like to thank the AMIDEX for financial support.

FIG. 1 .

 1 FIG. 1. Lowest order anharmonic phonon-phonon self-energy diagram, describing the anharmonic decay of a high-energy phonon (ω + ω ) into two lower energy phonons (ω and ω ).

FIG. 2 .

 2 FIG. 2. Top panel: atomic structure of the 100 -oriented Si or Ge NWs investigated in subsection III A. The length and the diameter of the nanowires are LNW = 60 nm and dNW = 3 nm, respectively. Bottom panel: cross-section view of (left) the 100 -oriented and (right) the 110 -oriented nanowires.

FIG. 3 .

 3 FIG. 3. (a) Room temperature ballistic, SCBA and LOA thermal currents up to the third order in the 110 -oriented Si NW. Inset: enlargement around the value of the SCBA current. (b) Padé approximants 0/1, 1/1 and 1/2 built from the LOA results. The ballistic and SCBA currents are also represented.

FIG. 4 .

 4 FIG. 4. (a) Room temperature ballistic, SCBA and LOA thermal currents up to the third order in the 110 -oriented Ge NW. (b) Padé approximants 0/1, 1/1 and 1/2 built from the LOA results. The ballistic and SCBA currents are also represented.

FIG. 5 .

 5 FIG. 5. (a) Illustration of the atomic structure (not to scale) of the 100 -oriented Si NW in the presence of a single recessed region (NWR1). (b) Illustration of the atomic structure (not to scale) of the 100 -oriented Si NW in the presence of two recessed regions (NWR2). For the sake of clarity, the Si atoms in the recessed regions are shown in yellow.

FIG. 6 .

 6 FIG. 6. (a) Room temperature ballistic, SCBA and LOA thermal currents up to the third order in the NWR1 structure. Inset: enlargement around the value of the SCBA current. (b) Padé approximants 0/1, 1/1 and 1/2 built from the LOA results. The ballistic and SCBA currents are also represented. Inset: enlargement around the value of the SCBA current.

FIG. 7 .

 7 FIG. 7. (a) Room temperature ballistic, SCBA and LOA thermal currents up to the third order in the NWR2 structure. Inset: enlargement around the value of the SCBA current. (b) Padé approximants 0/1, 1/1 and 1/2 built from the LOA results. The ballistic and SCBA currents are also represented.

TABLE I .

 I LOA thermal currents analytically continued by the Padé approximants 0/1 and 1/2 and the corresponding Richardson extrapolations for the 110 -oriented Si NW, at different temperatures. The number of SCBA iterations required to obtain each approximation, the relative difference with respect to the SCBA results, and the corresponding ballistic and SCBA currents are also reported.

		Ballistic	1 st LOA + Padé 0/1	3 rd LOA + Padé 1/2	1 st Richardson	SCBA	Temperature
	I ph [W]	5.16e-10	2.17e-10	2.53e-10	2.89e-10	2.78e-10	300 K
	ε [%]	85.6	21.9	9.0	4.0	0.0
	Iterations	0	1	6	6	10
	I ph [W]	5.71e-10	1.34e-10	1.71e-10	2.07e-10	2.07e-10	500 K
	ε [%]	175.8	35.3	17.4	0.0	0.0
	Iterations	0	1	6	6	19
	I ph [W]	5.89e-10	8.87e-11	1.22e-10	1.56e-10	1.54e-10	700 K
	ε [%]	282.5	42.4	20.8	1.3	0.0
	Iterations	0	1	6	6	28
	I ph [W]	5.97e-10	6.05e-11	8.14e-11	1.02e-10	1.15e-10	900 K
	ε [%]	419.1	47.4	29.2	11.3	0.0
	Iterations	0	1	6	6	44