Tailoring the Structural and Optical Properties of Germanium Telluride Phase-Change Materials by Indium Incorporation

Chalcogenide phase-change materials (PCMs) based random access memory (PCRAM) enter the global memory market as storage-class memory (SCM), holding great promise for future neuro-inspired computing and non-volatile photonic applications. The thermal stability of the amorphous phase of PCMs is a demanding property requiring further improvement. In this work, we focus on indium, an alloying ingredient extensively exploited in PCMs. Starting from the prototype GeTe alloy, we incorporated indium to form three typical compositions along the InTe-GeTe tie line: InGe3Te4, InGeTe2 and In3GeTe4. The evolution of structural details, and the optical properties of the three In-Ge-Te alloys in amorphous and crystalline form, was thoroughly analyzed via ab initio calculations. This study proposes a chemical composition possessing both improved thermal stability and sizable optical contrast for PCM-based non-volatile photonic applications.

PCMs can be switched rapidly and reversibly between their amorphous and crystalline phases via Joule heating induced by electrical or optical pulses [1,35]. The notable contrast in either electrical resistivity or optical reflectivity between each phase is utilized to encode digital information [1]. Several demanding requirements, such as high programming speed, good thermal stability, low power consumption, stable property contrast window and long cycling endurance, have to be well satisfied for high-performance PCRAM. Germanium chalcogenides, in particular, GeTe and GeTe-Sb 2 Te 3 pseudo-binary compounds (GST), especially Ge 2 Sb 2 Te 5 [36], are one of the most successful material families that could meet these challenging requirements simultaneously. Doping and alloying are frequently used to tailor the material properties for faster speed and/or better retention temperature, targeting different application scenarios [37][38][39][40][41][42][43][44].
For decades, indium has been an important alloying element used in rewritable optical data storage products [45]. The flagship PCM is AgInSbTe [45][46][47][48]. Recently, indium-alloyed GeTe [49][50][51][52] and GST [53,54] were reported, and their enhanced amorphous stability makes them suitable candidates for high-temperature PCRAM applications. In addition, indium forms a unique PCM In 3 SbTe 2 [55][56][57][58][59] that exhibits metallic behavior in its crystalline phase, but semiconducting behaviors in its amorphous phase, in contrast to conventional PCMs, which remain semiconducting during memory programming. It has been suggested that even InTe could also be a potential PCM for non-volatile electronics [53,60]. In this work, we focus on the InTe-GeTe (IGT) tie line, in particular, the three stoichiometric compositions, namely InGe 3 Te 4 , InGeTe 2 and In 3 GeTe 4 . By performing thorough ab initio calculations and chemical bonding analyses, we elucidate the role of indium in altering the structural and optical properties of GeTe.

Computational Details
We performed ab initio molecular dynamics (AIMD) simulations based on density functional theory (DFT) to generate melt-quenched amorphous structures [61]. The secondgeneration Car-Parrinello method [62] as implemented in CP2K package [63] was employed along with Perdew-Burke-Ernzerhof (PBE) functional [64] and the Goedecker pseudopotentials [65]. The canonical NVT ensemble was used and the time step was set at 2 fs. Vienna Ab-initio Simulation Package (VASP) [66] was employed to relax the amorphous structures and crystalline counterparts, prior to the calculations of electronic structure and optical response. For VASP calculations, we applied the PBE functional and projector augmented-wave (PAW) pseudopotentials [67]. The energy cutoff for plane waves was set at 500 eV. Chemical bonding analyses were conducted with the LOBSTER code [68][69][70]. Crystal orbital Hamilton populations (COHP) were applied to separate the covalent interactions into bonding (positive −COHP) and antibonding (negative −COHP) contributions. Bader charges were calculated to evaluate the atomic charge transfer in the structures [71]. Frequency-dependent dielectric matrix was calculated within the independent-particle approximation without considering local field effects and many body effects, which proved to be adequate to account for the optical contrast between crystalline and amorphous PCMs [72][73][74]. The absorption α(ω) and reflectivity R(ω) can be calculated from the dielectric functions [75]: where ε 1 and ε 2 are the real and imaginary parts of the dielectric function. n and k are the refractive index and extinction coefficient, which can be calculated from the dielectric functions: Van der Waals correction based on the Grimme's D3 method was considered in all AIMD and DFT calculations [76,77]. All the electronic structures, chemical bonding and optical properties were calculated using relaxed structures at zero K with VASP. For standard calculations, only gamma point was used to sample the Brillouin zone of the supercell models, while a 3 × 3 × 3 k-point mesh was used to converge the optical response calculations. For statistics, we built three crystalline and three amorphous models for each composition.

Results and Discussion
As reported in Ref. [78], a single-phase rock-salt structure was obtained over a wide compositional range of 8-75 mole% InTe in IGT at ambient conditions, in which Te atoms occupied one sublattice, while Ge and In atoms shared the other one. The three IGT compositions considered in this work, namely InGe 3 Te 4 , InGeTe 2 and In 3 GeTe 4 , fall in this compositional range, and were expected to take the rock-salt structure. To account for the compositional disorder of Ge and In atoms on the cation-like sublattice, we built 3 × 3 × 3 supercells (216 atoms in total) and distributed Ge and In atoms using a quasirandom number generator. Three independent models were considered for each composition. Each supercell model was fully relaxed with respect to both atomic coordinates and cell volume by DFT calculations. The relaxed cell-edge of crystalline (c-) InGe 3 Te 4 , InGeTe 2 and In 3 GeTe 4 is 18.21, 18.28 and 18.55 Å, respectively. The corresponding unit cell lattice parameters, 6.07, 6.09 and 6.18 Å, are in good agreement with experimental values (5.97, 6.00 and 6.06 Å) [78]. The relaxed structure of c-InGeTe 2 is shown in Figure 1a and the structures of the other two compositions are shown in Figure S1. Nanomaterials 2021, 11, x FOR PEER REVIEW 3 of optical properties were calculated using relaxed structures at zero K with VASP. F standard calculations, only gamma point was used to sample the Brillouin zone of t supercell models, while a 3 × 3 × 3 k-point mesh was used to converge the optical respon calculations. For statistics, we built three crystalline and three amorphous models for ea composition.

Results and Discussion
As reported in Ref. [78], a single-phase rock-salt structure was obtained over a wi compositional range of 8-75 mole% InTe in IGT at ambient conditions, in which Te atom occupied one sublattice, while Ge and In atoms shared the other one. The three IGT com positions considered in this work, namely InGe3Te4, InGeTe2 and In3GeTe4, fall in this com positional range, and were expected to take the rock-salt structure. To account for t compositional disorder of Ge and In atoms on the cation-like sublattice, we built 3 × 3 × supercells (216 atoms in total) and distributed Ge and In atoms using a quasi-rando number generator. Three independent models were considered for each compositio Each supercell model was fully relaxed with respect to both atomic coordinates and c volume by DFT calculations. The relaxed cell-edge of crystalline (c-) InGe3Te4, InGeT and In3GeTe4 is 18.21, 18.28 and 18.55 Å, respectively. The corresponding unit cell latti parameters, 6.07, 6.09 and 6.18 Å, are in good agreement with experimental values (5.9 6.00 and 6.06 Å) [78]. The relaxed structure of c-InGeTe2 is shown in Figure 1a and t structures of the other two compositions are shown in Figure S1. The relaxed crystalline supercells were then used to generate amorphous (a-) mode following a melt-quench protocol [61]. The supercell models were quickly heated to very high temperature to remove the crystalline order. After randomization at 3000 K f 15 ps, the models were quenched down to and equilibrated at 1200 K, above the meltin point of IGT alloys (∼550-750 °C) [49] for 30 ps. Amorphous models were then generat by quenching the liquids down to 300 K with a cooling rate of 12.5 K/ps. During th quenching process, we stopped the simulation after every 100 K, and the simulation b size was increased to reduce the internal stress. Within each temperature window, o NVT calculation was performed using a fixed box size. The model was equilibrated at 3 K for 30 ps. This density value was then used to generate two additional melt-quench amorphous models for each composition. All three amorphous models show The relaxed crystalline supercells were then used to generate amorphous (a-) models following a melt-quench protocol [61]. The supercell models were quickly heated to a very high temperature to remove the crystalline order. After randomization at 3000 K for 15 ps, the models were quenched down to and equilibrated at 1200 K, above the melting point of IGT alloys (∼550-750 • C) [49] for 30 ps. Amorphous models were then generated by quenching the liquids down to 300 K with a cooling rate of 12.5 K/ps. During this quenching process, we stopped the simulation after every 100 K, and the simulation box size was increased to reduce the internal stress. Within each temperature window, one NVT calculation was performed using a fixed box size. The model was equilibrated at 300 K for 30 ps. This density value was then used to generate two additional melt-quenched amorphous models for each composition. All three amorphous models showed consistently low pressure values below 3 kbar. The obtained cell edges of amorphous InGe 3 Te 4 , InGeTe 2 and In 3 GeTe 4 are 18.90, 19.08 and 19.32 Å, respectively. This increase in the cell edge of the amorphous phase is consistent with the trend observed in their crystalline counterparts.
Further optimization of the internal stress or the use of the NPT ensemble for melt-quench simulations could potentially lead to some numerical differences in the mass density, but is not expected to alter the amorphous structures much.
The amorphous structure of InGeTe 2 is shown in Figure 1a and the snapshots of the other two compositions are in Figure S1. The partial radial distribution functions (RDFs) of each atomic pair in the three amorphous compounds are shown in Figure 1b. The peak positions of the heteropolar bonds In-Te (2.87 Å) and Ge-Te (2.78 Å) are not varied with composition, whereas the homopolar or "wrong" bonds show small shifts. As developed in our previous work [79], the "bond-weighted distribution function (BWDF)" provides direct information on the length of chemical bonds in amorphous IGT alloys ( Figure S2). Despite the change in chemical composition, the bond length shows very similar values in amorphous IGT alloys, i.e., Ge-Te 3.20 Å, In-Te 3.40 Å, Ge-Ge 3.20 Å, In-In 3.25 Å and Ge-In 3.40 Å. In all three amorphous IGT alloys, Te-Te shows mostly antibonding interactions. These bond length values are used as cutoffs for the interatomic distance for the following structural analysis.
The angle distribution function (ADF) of the three amorphous structures (Figure 2a) shows that In and Ge atoms mainly form local motifs with the central bond angles ranging from 90 • to 109.5 • , which corresponds to the bond angles in octahedral and tetrahedral motifs, respectively. As the concentration of indium increases, the ADF peak for indium atoms clearly shifts toward 109.5 • , implying an increase in indium-centered tetrahedral motifs. We used bond order parameter q [80] to quantify the fraction of tetrahedral motifs in amorphous IGT alloys. Such parameters are frequently used for the structural analysis of amorphous PCMs [81][82][83]. As shown in Figure 2b, as indium concentration increases, the total fraction of In-and Ge-centered tetrahedral motifs (short as tetra-In and tetra-Ge) increases from 31.6% (a-InGe 3 Te 4 ), 33.5% (a-InGeTe 2 ) to 42.4% (a-In 3 GeTe 4 ) in relation to the total number of In and Ge atoms. Specifically, the fraction of tetra-In increases from 10.1%, 18.6% to 34.6%, while the fraction of tetra-Ge decreases from 21.5%, 14.9% to 7.8%. The fraction of tetrahedrons for a-InGeTe 2 is smaller than that reported in previous work (38 % in total, with 31.5% for tetra-In and 6.5% for tetra-Ge) [50], due to the different choice of cutoff values for the interatomic distance and the deviation in calculated density values (vdW interactions were included in the current work). The ratios of tetrahedral motifs in amorphous IGT alloys are all higher than that of their parent phase-GeTe, where 25-30% tetrahedral atoms are typically found in the rapidly quenched amorphous phase [84][85][86].
In addition to the increase in the total number of tetrahedral motifs from a-InGe 3 Te 4 to a-In 3 GeTe 4 , the local bonding configuration also shows a major difference. Despite the change in chemical compositions, nearly all the tetra-Ge atoms are bonded with at least one Ge or In atom, while the majority of tetra-In atoms are heteropolar-bonded in the three amorphous IGT alloys. To quantify the role of "wrong" bonds, we carried out a projected COHP (pCOHP) analysis. As shown in Figures 3 and S3, the pCOHP of heteropolar-bonded tetra-Ge atoms demonstrates sizable antibonding interactions right below Fermi energy (E F ), and the presence of wrong bonds (including Ge-Ge and Ge-In) largely reduces such antibonding contributions, stabilizing the tetrahedral motifs locally. In contrast, the pCOHP of tetra-In atoms with and without wrong bonds mostly shows bonding interactions below E F . These results are consistent with the bonding configuration in the two parent phases, a-GeTe [79] and a-InTe [83], though In-Ge bonds are present in all three amorphous IGT alloys. Since indium atoms do not require homopolar bonds to stabilize tetrahedral motifs, the ratio of tetrahedral units is increased with the indium concentration. As compared to their crystalline counterparts, where all In and Ge atoms are octahedrally bonded, the enlarged structural deviation will enhance the thermal stability of amorphous IGT alloys. This observation is consistent with experimental findings, as the crystallization temperature T x of doped InGeTe 2 [49] and undoped InTe thin films [60] is increased to ∼276 and ∼300 • C, respectively, as compared to that of GeTe T x ∼190 • C [87].
is increased to ∼276 and ∼300 °C, respectively, as compared to that of GeTe Tx ∼190 [87].  The calculated density of states (DOS) of the three IGT alloys in both crystalline a amorphous forms are shown in Figure 4a,b. Regarding the crystalline models, the over DOS profiles are quite similar, and all three alloys exhibit metallic features. By contra all three amorphous models are narrow-gap semiconductors. Statistical sampling yie consistent results ( Figure S4). The large difference in DOS between the crystalline a amorphous IGT results in a wide resistance contrast window for PCRAM applicatio    The calculated density of states (DOS) of the three IGT alloys in both crystalline an amorphous forms are shown in Figure 4a,b. Regarding the crystalline models, the over DOS profiles are quite similar, and all three alloys exhibit metallic features. By contra all three amorphous models are narrow-gap semiconductors. Statistical sampling yiel consistent results ( Figure S4). The large difference in DOS between the crystalline an amorphous IGT results in a wide resistance contrast window for PCRAM applicatio The calculated density of states (DOS) of the three IGT alloys in both crystalline and amorphous forms are shown in Figure 4a,b. Regarding the crystalline models, the overall DOS profiles are quite similar, and all three alloys exhibit metallic features. By contrast, all three amorphous models are narrow-gap semiconductors. Statistical sampling yields consistent results ( Figure S4). The large difference in DOS between the crystalline and amorphous IGT results in a wide resistance contrast window for PCRAM applications [49]. The Bader charge analysis (Figure 4c) details larger net charges for In atoms than for Ge atoms due to the difference in electronegativity. The bimodal feature of the charges of indium atoms is consistent with previous work [50], stemming from different local environments of indium atoms. The enlarged charge transfer in amorphous structures increases the probability of long-distance electromigration under the transient electrical field induced by programming pulses [88], which is detrimental to the cycling endurance of devices [89,90]. For RESET operations, the higher the melting temperature T m , the greater the power consumption. The melting temperature for IGT alloys has a "W" shape profile, according to the InTe-GeTe phase diagram, which shows that InGeTe 2 has a higher melting temperature T m (740 • C) than GeTe and InTe whose T m are 715 and 688 • C, respectively [49]. Interestingly, the two other compositions, InGe 3 Te 4 and In 3 GeTe 4 , demonstrated reduced T m (645 and 565 • C). Taking into account all these factors for practical applications, we would suggest keeping the IGT composition within the range of InGeTe 2 to In 3 GeTe 4 for balanced device performance.  [49]. The Bader charge analysis (Figure 4c) details larger net charges for In atoms than Ge atoms due to the difference in electronegativity. The bimodal feature of the charges indium atoms is consistent with previous work [50], stemming from different local en ronments of indium atoms. The enlarged charge transfer in amorphous structures creases the probability of long-distance electromigration under the transient electri field induced by programming pulses [88], which is detrimental to the cycling enduran of devices [89,90]. For RESET operations, the higher the melting temperature Tm, greater the power consumption. The melting temperature for IGT alloys has a "W" sha profile, according to the InTe-GeTe phase diagram, which shows that InGeTe2 has a high melting temperature Tm (740 °C) than GeTe and InTe whose Tm are 715 and 688 °C, resp tively [49]. Interestingly, the two other compositions, InGe3Te4 and In3GeTe4, demo strated reduced Tm (645 and 565 °C). Taking into account all these factors for practi applications, we would suggest keeping the IGT composition within the range of InGe to In3GeTe4 for balanced device performance. The enhanced amorphous stability is also useful for non-volatile photonic appli tions [91], yet the incorporation of indium makes IGT alloys metallic, which could aff the optical contrast between the amorphous and crystalline phase. The significant contr of ∼30% in the optical reflectivity of PCMs stems from a fundamental change in bondi nature from covalent to metavalent bonding (MVB) upon crystallization [92][93][94][95][96][97][98][99]. Ho ever, in comparison with GeTe and GST, which have three p electrons per site (a key f ture of MVB), InTe has a deficient number of p electrons, turning the rock-salt phase fro semiconducting to metallic. As a result, MVB in IGT alloys is expected to be weakened For verification, we carried out optical response calculations using the relaxed cr talline and amorphous IGT structures. We focused on the spectrum range from 400 to 16 nm, covering both the visible light region (∼400 nm to 800 nm) for optical displays [1 The enhanced amorphous stability is also useful for non-volatile photonic applications [91], yet the incorporation of indium makes IGT alloys metallic, which could affect the optical contrast between the amorphous and crystalline phase. The significant contrast of ∼30% in the optical reflectivity of PCMs stems from a fundamental change in bonding nature from covalent to metavalent bonding (MVB) upon crystallization [92][93][94][95][96][97][98][99]. However, in comparison with GeTe and GST, which have three p electrons per site (a key feature of MVB), InTe has a deficient number of p electrons, turning the rock-salt phase from semiconducting to metallic. As a result, MVB in IGT alloys is expected to be weakened.
For verification, we carried out optical response calculations using the relaxed crystalline and amorphous IGT structures. We focused on the spectrum range from 400 to 1600 nm, covering both the visible light region (∼400 nm to 800 nm) for optical displays [19][20][21] and the telecom wavelength bands (∼1500 to 1600 nm) for silicon-waveguideintegrated photonic applications [23][24][25][26]. As shown in Figure 5a, the optical absorption and reflectivity profiles vary slightly with the chemical compositions in the amorphous phase, while strong changes are found in the crystalline phase. For InGe 3 Te 4 and InGeTe 2 , sizable contrast in reflectivity between the crystalline and amorphous phases is observed over the whole spectrum, with an average value ∼20%. However, for In 3 GeTe 4 , the contrast in reflectivity nearly vanishes at around ∼900 nm, and it becomes greater below 700 nm, or above 1000 nm.  21] and the telecom wavelength bands (∼1500 to 1600 nm) for silicon-waveguide-integrated photonic applications [23][24][25][26]. As shown in Figure 5a, the optical absorption and reflectivity profiles vary slightly with the chemical compositions in the amorphous phase, while strong changes are found in the crystalline phase. For InGe3Te4 and InGeTe2, sizable contrast in reflectivity between the crystalline and amorphous phases is observed over the whole spectrum, with an average value ∼20%. However, for In3GeTe4, the contrast in reflectivity nearly vanishes at around ∼900 nm, and it becomes greater below 700 nm, or above 1000 nm. The three IGT crystals have a crossover at ∼1100 nm. For long wavelength or small photon energy (below ∼1 eV) regions, the optical excitation is mainly determined by states near EF. With increased DOS near EF by heavier indium alloying (Figure 4a), the absorption in the long wavelength region is enhanced. However, for the short wavelength or large photon energy (above ∼2 eV) region, the electronic states of a wider energy range would participate in the optical excitation. According to the projected DOS (Figure 5b), the valence states below EF are mainly contributed by Te atoms for all three crystals. However, indium alloying shifts the DOS peak in the conduction band to a higher energy range. Since fewer excited states of low energy could contribute to short wavelength excitation, c-In3GeTe4 shows the smallest absorption and reflectivity in the short wavelength region (Figure 5a).
We note that our optical calculations were performed using DFT-PBE functional with an independent-particle approximation, excluding local field effects and many body The three IGT crystals have a crossover at ∼1100 nm. For long wavelength or small photon energy (below ∼1 eV) regions, the optical excitation is mainly determined by states near E F . With increased DOS near E F by heavier indium alloying (Figure 4a), the absorption in the long wavelength region is enhanced. However, for the short wavelength or large photon energy (above ∼2 eV) region, the electronic states of a wider energy range would participate in the optical excitation. According to the projected DOS (Figure 5b), the valence states below E F are mainly contributed by Te atoms for all three crystals. However, indium alloying shifts the DOS peak in the conduction band to a higher energy range. Since fewer excited states of low energy could contribute to short wavelength excitation, c-In 3 GeTe 4 shows the smallest absorption and reflectivity in the short wavelength region (Figure 5a).
We note that our optical calculations were performed using DFT-PBE functional with an independent-particle approximation, excluding local field effects and many body effects. Therefore, the absolute values of optical profiles could vary if more advanced methods are employed. Nevertheless, the observation of weakened optical contrast due to heavier indium alloying should remain valid. Taking into account the enhanced crystallization temperature and the reduced melting temperature, we predict an optimal IGT composition within the ranges of InGeTe 2 and In 3 GeTe 4 for high-performance non-volatile photonics. To the best of our knowledge, thorough optical measurements of IGT alloys are still lacking. Therefore, we anticipate future experiments exploring the suitability of IGT alloys for optical and photonic PCM applications.

Conclusions
In summary, we have carried out systematic ab initio calculations for three typical compositions of indium incorporated GeTe compounds, InGe 3 Te 4 , InGeTe 2 and In 3 GeTe 4 , to elucidate the evolution of structural and optical properties along the InTe-GeTe tie line. Upon indium alloying, the crystalline phase of all the three alloys turns metallic, while their amorphous counterparts all show semiconducting features with narrow band gaps. This stark contrast in the electronic structure guarantees a large resistance window between amorphous and crystalline In-Ge-Te alloys for electrical PCRAM. Yet, too much indium should be avoided, because the stronger charge transfer could be harmful to cycling endurance due to electromigration. Regarding optical properties, both InGe 3 Te 4 and InGeTe 2 show sizable optical contrast between the crystalline and amorphous phases in the spectrum range from 400 nm to 1600 nm, covering both visible-light and telecom bands. Meanwhile, In 3 GeTe 4 shows a less robust contrast window, due to weakened MVB. Moreover, the increased indium concentration enlarges the ratio of tetrahedral motifs in the amorphous phase and consequently increases the structural barrier for crystallization. The InTe-GeTe phase diagram establishes that the melting temperature reaches minimum around In 3 GeTe 4 , indicating lowest power consumption for melt-quench amorphization. Taking all these factors into account, we suggest that the optimal chemical composition for In-Ge-Te alloys should be located in the range between InGeTe 2 and In 3 GeTe 4 , which could result in the most balanced device performance for PCM-based non-volatile electronic and photonic applications. Our work should serve as a stimulus for further investigations into indium-incorporated PCMs.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/nano11113029/s1. Figure S1: Atomic structures of crystalline and amorphous InGe 3 Te 4 and In 3 GeTe 4 (denoted as c/a-134 and 314); Figure S2: Bond population (left panel) and bond weighted distribution functions BWDFs (right panel) for each interatomic pair in the three amorphous compounds. The crossover values from positive to negative in the BWDFs represent cutoff values for bonding interactions, which are used for structural analysis. Further technical details about BWDF can be found in [79] of the main text; Figure S3: Projected COHP (pCOHP) for tetra-In and tetra-Ge motifs in a-InGe 3 Te 4 and a-In 3 GeTe 4 . Tetra-hedral motifs are grouped as the ones with only heteropolar bonds, denoted as In[In 0 Ge 0 Te 4 ] and Ge[Ge 0 In 0 Te 4 ], and the others with at least one wrong bond indicated as In[In x Ge y Te 4−x−y ] and Ge[Ge x In y Te 4−x−y ] (x or y ≥ 1, x + y ≤ 4). In the a-314 structure, Ge[Ge 0 In 0 Te 4 ] motif is absent; Figure S4: Density of states (DOS) for crystalline and amorphous structures of the three IGT compositions. Three models for each composition were built, which show consistent results.
Author Contributions: Investigation, visualization, X.W., X.S. and S.S.; writing-original draft preparation, X.W.; conceptualization, funding acquisition, writing-review and editing, W.Z. All authors have read and agreed to the published version of the manuscript.
Funding: This work was funded by National Natural Science Foundation of China (61774123) and 111 Project 2.0 (BP2018008).

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.