An arbitrary order time-stepping algorithm for tracking particles in inhomogeneous magnetic fields

TL;DR

This paper introduces an improved high-order Boris-SDC algorithm for particle tracking in magnetic fields, enhancing accuracy and stability while reducing computational costs compared to the standard Boris method.

Contribution

The authors develop a new Boris-SDC algorithm using GMRES-based convergence acceleration, achieving arbitrary order accuracy with better energy stability and efficiency.

Findings

Better accuracy at lower computational cost

Enhanced long-term energy stability

Effective in magnetic mirror and equilibrium scenarios

Abstract

The Lorentz equations describe the motion of electrically charged particles in electric and magnetic fields and are used widely in plasma physics. The most popular numerical algorithm for solving them is the Boris method, a variant of the St\"ormer-Verlet algorithm. Boris' method is phase space volume conserving and simulated particles typically remain near the correct trajectory. However, it is only second order accurate. Therefore, in scenarios where it is not enough to know that a particle stays on the right trajectory but one needs to know where on the trajectory the particle is at a given time, Boris method requires very small time steps to deliver accurate phase information, making it computationally expensive. We derive an improved version of the high-order Boris spectral deferred correction algorithm (Boris-SDC) by adopting a convergence acceleration strategy for second order problems based on the Generalised Minimum Residual (GMRES) method. Our new algorithm is easy to implement as it still relies on the standard Boris method. Like Boris-SDC it can deliver arbitrary order of accuracy through simple changes of runtime parameter but possesses better long-term energy stability. We demonstrate for two examples, a magnetic mirror trap and the Solev'ev equilibrium, that the new method can deliver better accuracy at lower computational cost compared to the standard Boris method. While our examples are motivated by tracking ions in the magnetic field of a nuclear fusion reactor, the introduced algorithm can potentially deliver similar improvements in efficiency for other applications.

Figures

Tables

Table 1. Table 1 : Energy error for Solev’ev equilibrium after 100000 100000 100000 time steps of length 1   ns times 1 nanosecond 1\text{\,}\mathrm{ns} .

Iterations $k$	BGSDC w. Picard	BGSDC w. Sweeps	Newton-BGSDC
1	$7.0\times {10}^{-6}$	$1.0\times {10}^{-6}$	$7.4\times {10}^{-1}$
2	$3.7\times {10}^{-7}$	$3.8\times {10}^{-7}$	$5.1\times {10}^{-6}$
3	$3.8\times {10}^{-7}$	$3.8\times {10}^{-7}$	$3.8\times {10}^{-7}$

Equations

$${\mathbf{v}}^{n+1/2}={\mathbf{v}}^{n-1/2}+\alpha \Delta t\frac{{\mathbf{v}}^{n+1/2}+{\mathbf{v}}^{n-1/2}}{2}\times {\mathbf{B}}^n$$

$${\Vert {\mathbf{v}}^{n+1/2}\Vert }^2={\Vert {\mathbf{v}}^{n-1/2}\Vert }^2$$

$$\frac{1}{2}\left(\mathbf{f}({\mathbf{x}}_n,{\mathbf{v}}_n)+\mathbf{f}({\mathbf{x}}_{n+1},{\mathbf{v}}_{n+1})\right)\mapsto \mathbf{f}({\mathbf{x}}_{n+1/2},{\mathbf{v}}_{n+1/2})$$

$${\mathbf{v}}_{m+1}^{k+1}={\mathbf{v}}_m^{k+1}+\Delta {\tau }_m\mathbf{f}({\mathbf{x}}_{m+1/2}^{k+1},{\mathbf{v}}_{m+1/2}^{k+1})-\Delta {\tau }_m\mathbf{f}({\mathbf{x}}_{m+1/2}^k,{\mathbf{v}}_{m+1/2}^k)+{\mathbf{IF}}_m^{m+1}$$

$${\mathbf{v}}_{n+1/2}={\mathbf{v}}_n+\frac{\Delta t}{2}\mathbf{f}({\mathbf{x}}_n,{\mathbf{v}}_n),$$

Full text

Answers to reviewers

JCOMP-D-18-02028 “An arbitrary order time-stepping algorithm for tracking particles in inhomogeneous magnetic fields”

We would like to thank the three reviewers for their detailed comments on the manuscript. In response to the issues raised, we have made major revisions to the paper, highlighted in colour:

• •

Modifications made in reply to issues raised by Reviewer #1 are marked in red.

• •

Modifications made in reply to issues raised by Reviewer #2 are marked in green.

• •

Modifications made in reply to issues raised by Reviewer #3 are marked in blue.

Please note that we have also removed or shortened substantial parts of the paper in the revised form in response to issue (i) raised by Reviewer #3.

1 Reviewer #1.

1. Page 3: It would be helpful to readers not familiar with SDC methods, to add a few sentences here describing SDC methods. Perhaps before the paragraph starting ”This paper presents…”, or by expanding the paragraph below which starts ”Boris-GMRES-SDC (below,…”.

Answer**.**

Two sentences have been added to briefly characterise SDC, including proper attribution to Dutt et al who introduced the method (p. 3, lines 60–63). In the original version this reference was only given in the methodology section, but after the reviewer’s comments we think it deserves to be mentioned in the introduction.

2. Page 4: ”the electric field generated through particle interactions is negligible”. This is not true in all cases, though it can be a useful approximation. Close to the edge of tokamaks, for example, radial electric fields play an important role. I suggest qualifying this statement, to say that it is a reasonable approximation in many cases, fast particle transport in the core region usually being one of them.

Answer**.**

We agree that this statement was way too general. It has been rephrased to point out clearly that we focus on tracking fast particles in the core region (see page 4, lines 85–89). Please see also the answer to the first issue raised by Reviewer #2, who also pointed out this issue. While we still neglect particle interactions, we now include the effect of $\boldsymbol{\mathbf{E}}\times\boldsymbol{\mathbf{B}}$ drift in the second example by including a weak, external electric field.

3. Page 4: ”a real reactor where the field is given by measured data”. For experimental scenarios the magnetic field would be given by a solution to the Grad-Shafranov equation, fitted to measurements of the magnetic field made around the boundary of the plasma: The field is not directly measured in the plasma. I suggest modifying to ”a real reactor where the field is given by a numerical solution to the Grad-Shafranov equation.”

Answer**.**

This was indeed formulated a bit sloppy. Changed as suggested to “in contrast to a real reactor where the field is given by a numerical solution to the Grad-Shafranov equation fitted to experimental data” (see page 4, lines 94–95).

4. Page 5: Please make it clearer that Algorithm 1 corresponds to unstaggered Boris (equation 2 + ”Boris trick”)

Answer**.**

The caption for Algorithm 1 has been amended to make this clearer (see page 5).

5. Page 7: Surely the reason why 6d x 6d equations can be solved rather than 6Md x 6Md is that there is no coupling between the particles: the system trivially separates into equations for each particle separately. This is not a feature of Boris-SDC as implied by ”Instead, Boris-SDC solves…”

Answer**.**

This was not made very clear in the original paper: we were referring to the fact that, for a collocation method, the equations for the $M$ stages are all coupled. SDC decouples those independent of whether particles are interacting or not. The independence of particles leads to another form of decoupling though: each of the $M$ nonlinear systems consisting of $6d$ equations decouples into blocks of size $6$, one for each particle. This is, however, not required by GMRES-SDC. Note that this sentence no longer appears in the revised manuscript.

6. Page 17: The location of the table, before the definition of the quantities it contains, is a bit confusing. Perhaps the definition of the symbols could be in the table, the table caption, or the table moved to be closer to the definitions.

Answer**.**

Table 1 and Table 2 in the original manuscript have been merged together and are now Table 1 on page 15 in the revised manuscript.

7. Page 17: ”different than the ones for the numerical examples” – could this be made a bit clearer, since Figure 1 is a kind of numerical example. I suggest ”different than the parameters used to quantitatively evaluate the numerical method, which are discussed in section 3.1.1”.

Answer**.**

This has been rephrased and is hopefully clearer now, see page 14, lines 199–201.

8. Page 18: In equation 47, I suggest a modulus around the gradient of B e.g. $|\nabla B|$

Answer**.**

Fixed (see what is now Eq. (30) on page 15).

9. Page 20 (Fig 2),21. The saturation in the error at around $10^{-3}$ is only seen in the GMRES-SDC results. Is this really because the original Boris method has not reached that error level, or due to a problem/bug in the SDC implementation? Is it possible to either a) Run the non-stg Boris method with higher number of RHS evaluations, to show that it also saturates at the same level; or b) show that $10^{-3}$ is the error you would expect given the finite value of $\varepsilon$?

Answer**.**

We extended the range of time steps for staggered Boris in the Figure to demonstrate that it saturates at the same level as BGSDC. Note that, in order to maintain reasonable runtimes for very small $\Delta t$, we had to reduce $t_{\textrm{end}}$, see Table 1 (left) on p. 15, but this did not have a visible impact on the results.

10. Page 21: $\varepsilon\gg 0$ is a bit ambiguous because it does not give an indication of the scale; is 1e-2 ”large”? Equation 47 (p18) specifies $\varepsilon\ll 1$ as the adiabatic limit, so 1e-2 may still be said to be in the adiabatic regime.

Answer**.**

We agree. This has been rephrased, see p. 16, line 209. Furthermore, the headings have been renamed from ”Near adiabatic regime” and ”Non-adiabatic regime” to more neutral ”Scenario 1”, ”Scenario 2” with the order of magnitude of $\varepsilon$ given explicitly (see p. 15, line 205 and p. 17, line 233).

11. Page 25: Figure 5 caption please add ”for Boris-SDC (top) and GMRES-SDC (bottom)” as the reader (e.g. me) may miss the different captions initially.

Answer**.**

Added as suggested by the reviewer, see Fig. 5 on page 21 in revised manuscript.

12. Page 27, table 3: Please add units to the quantities given e.g Wb, m, T

Answer**.**

Fixed, see Table 2 on p. 23 in revised manuscript.

13. Page 31 ”compare for the classification” is unclear. ”compare to”?

Answer**.**

Replaced with ”following the classification by Gear”.

14. Page 32 ”our result suggest” $\rightarrow$ ”our results suggest”

Answer**.**

Fixed.

Reviewer #2.

The following assertion is made on page 4, toward the end of the manuscript’s introduction: “While GMRES-SDC can be applied to problems where an electric field is present, we focus here on scenarios motivated by nuclear fusion reactors. There, particle motion is governed by a strong magnetic field and the electric field generated through particle interactions is negligible.” It is true that the strong magnetic field in magnetic fusion reactors plays a dominant role in the single particle motion, causing it to gyrate rapidly about field lines. However, no modern kinetic simulation of these reactors (see e.g. XGC, GENE, GEM, ORB5, GYRO, and others) would dream of neglecting the self-consistent electric field arising from particle interactions. The reasons for this are myriad, but a particularly noteworthy one is the E-cross-B drift induced by electric fields. While this is indeed small compared to typical particle gyration velocities in tokamaks and stellarators, it gives rise to drift motion perpendicular to the magnetic field that has important physical consequences on time-scales long compared to the gyroperiod. I am not an expert on the fast ion motion that is mentioned as an application for this work – if it is the case that for fast ions the electric field may be reasonably neglected, it should be emphasized that this is a very special case and a more detailed justification for neglecting the electric field in this case should be given.

Answer**.**

We have extended this paragraph to make clearer that we focus on fast ion tracking and that the electric field cannot be neglected for many other processes in fusion reactors, see p. 4, lines 85–89 (please note these changes are coloured red, as the issue was also raised by Reviewer #1). For fast ion tracking, ignoring effects of the E-field is common: we have added two references ([3], [20] in the revised manuscript), which state “fast ion redistribution in real and velocity space due to the presence of radial electric fields […] is assumed to be negligible”, “A toroidal electric field can be imposed upon the model fast ion population, […]; however, this typically results in only $\approx$ 3% of the ohmic heating power being diverted into the fast ion population” and “We also neglect the effect of the electric field on the fast ion.” (for the sake of brevity, only the references but not the quotes are included in the revised paper).

Moreover, testing in scenarios featuring identically zero electric field considerably simplifies the energy conservation properties of an integration scheme (more on this in the notes below) – indeed, in this case a particle’s total energy is simply its kinetic energy, independent of particle location. Because of these concerns, I feel strongly that a test case with non-zero electric field should be added to the manuscript before it can be accepted. The Penning trap case in [17] would suffice, although a case with spatially varying magnetic field as well would be preferable for reasons the authors outline in the manuscript.

Answer**.**

We agree with reviewer’s comment. We’ve added non-constant radial electric field to Solov’ev example and updated all the results in section 3.2. The given fields create $E\times B$ drift and slightly changes the trajectories. Updated energy plots represent total energy ($E_{kin}+E_{pot}$) error for BGSDC as well as for stg Boris method. This has not affected the basic outcome, though: BGSDC is competitive over staggered Boris if precision of the order of millimetres is desired.

The authors point out that the “unstaggered” velocity-Verlet scheme of equation (2) and the “Boris algorithm” of equation (3) have quite similar properties, which is true. However, there is an important distinction that is particularly relevant when no electric field is present, as in the current test cases in the manuscript. In the Boris method (3), the property that the magnetic field does no work in the continuum is preserved exactly in the discrete. In the case when the electric field vanishes, this leads equation (3) to feature exact energy conservation. This is easy to show: the Boris velocity update with no electric field is just

[TABLE]

Dotting this expression with the average of the two velocity values immediately leads to

[TABLE]

Thus, when evaluated at the half-step (as is typically done in PIC codes), the magnitude of a particle’s velocity – and thus its kinetic energy – is constant in time, thereby conserving energy identically in the absence of an electric field. This distinction should at least be noted in the text,

Answer**.**

This is now mentioned explicitly in the discussion after Eq. (4) on page 5, lines 102–103.

but it also raises interesting questions: Could a GMRES-SDC scheme like the one in the manuscript built around the staggered Boris scheme share this property? If not, why? Additionally, the unstaggered velocity Verlet scheme on which the present scheme is built can be modified to share this exact conservation property by simply modifying

[TABLE]

in equation (2b) where quantities at the half-step are defined as the average of the corresponding quantities at times n and n+1. The Boris trick may still be used here to eliminate the apparent implicitness. It is not immediately obvious to me what other consequences this may have on the scheme – it is clearly still 2nd order, but does it feature energy drift? - but it would be an interesting avenue to explore.

Answer**.**

This is indeed an intriguing idea that we had not thought of. While a thorough investigation is beyond the scope of this reply, we implemented a basic version of a Boris-SDC with an velocity update based on the implicit mid-point rule instead of trapezoidal rule, as suggested by the reviewer. Interestingly, this does improve long-term energy error compared to standard Boris-SDC, see Fig. 2. However, the mid-point based update causes problems when trying to apply the GMRES acceleration. The resulting velocity sweep reads

[TABLE]

compared to Eq. (23) in the original manuscript (the revised manuscript no longer contains this equation). This no longer allows to express a sweep as the preconditioned iteration given by Eq. (16) in the revised manuscript, because the update for $\boldsymbol{\mathbf{x}}$ now uses values $\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{x}}_{m},\boldsymbol{\mathbf{v}}_{m})$ at full steps whereas the update of $\boldsymbol{\mathbf{v}}$ uses values $\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{x}}_{m+1/2},\boldsymbol{\mathbf{v}}_{m+1/2})$ at half-steps. While there are probably ways to apply the GMRES accelerator to the mid-point based formulation, this is, in our opinion, beyond the scope of this reply. Given that BGSDC (trapezoidal rule based Boris-SDC with GMRES acceleration) outperforms Boris-SDC-TP, we argue that for the moment, BGSDC still is the most competitive variant.

However, we agree with the reviewer that this is a very interesting direction for future research and have added a footnote about the mid-point based variant to the revised manuscript, see pages 8 and 9.

I do not understand equation (4). It appears to me that the particle’s velocity is advanced in time by 1.5 time-steps using this procedure (from n-1/2 to n + 1) but the force f on it is multiplied by a total factor of 1 delta t. Am I missing something? Should the factor of 1/2 on the f in (4a) be removed? In this case, this is just (3) but with an extra step on the end to sync position and velocity? This sync step is not typically used in applications, to the best of my knowledge.

Answer**.**

Yes, the reviewer is correct that (4) is simply a reformulation of (3) which also provides velocity information at full steps which is sometimes useful. While the factor $\frac{1}{2}$ is correct, there was indeed a mistake: (4a) should read

[TABLE]

with $\boldsymbol{\mathbf{v}}_{n}$ instead of $\boldsymbol{\mathbf{v}}_{n-1/2}$ as written incorrectly in the manuscript. While the sync step may not always be required, the kick-drift-kick formulation is a common way of deriving the Boris method by way of Strang splitting, see, for example, the following paper by Patacchini and Hutchinson in Journal of Computational Physics in 2009: https://doi.org/10.1016/j.jcp.2008.12.021.

The magnetic mirror field in equation (45) is problematic. As written, it is not divergence free. Looking at the authors’ code on GitHub (the authors are to be applauded for making this available), it appears that there is a typo in the z-component. The minus sign should be changed to a plus. Even after fixing this, though, this field represents an approximation to the true mirror field due, in the simplest case, to two circular current loops, and this approximation is only valid near the center of the mirror. For instance, the standard quantity of the “mirror ratio” that determines which particles are confined is undefined for the field in equation (45), since the magnitude of the field in unbounded. This field should either be replaced with the true field from two current rings (a sum of expressions found in Simpson et al., Simple analytic expressions for the magnetic field of a circular current loop (2001), or the near-axis approximation thereof) or it should be clarified that the given field is a mathematical model featuring similar characteristics to a magnetic mirror rather than a mirror itself.

Answer**.**

We’ve fixed the minus sign in the paper. The reviewer is correct that, for this example, we use a simple mathematical model that approximates the magnetic field between the two coils of a mirror trap. However, we assume that the maximum magnetic field is on the coil $B_{max}=B(z_{0})$ and the critical pitch angle $\varphi_{min}$ can be evaluated using $B_{max}$ and $B_{0}$ (magnetic field at the start position). We only study particle’s who angle is $\varphi_{0}$ is greater than $\varphi_{min}$ and which will stay inside the trap. Explanations including a reference the Simpson et al. (2001) have been added, see p. 13 above and below Eq. (13).

In the Solov’ev equilibrium example, it would be useful to report the gyrofrequency in units of inverse nanoseconds. This quantity varies in space, of course, but an order of magnitude estimate would aid the reader in interpreting the results.

Answer**.**

The gyrofrequency is now stated explicitly in the caption of Table 3, page 23.

The name “GMRES-SDC” has already been used by Huang et al in reference [18]. It is confusing to re-use it here. I would suggest either sticking to “Boris-GMRES-SDC” for the entire manuscript or coming up with a different shorthand – e.g. BGS – for the new scheme.

Answer**.**

We followed the suggestion by the reviewer but, to remain in line with acronyms for other variants of SDC, kept the SDC part. The method is now abbreviated as BGSDC (for Boris-GMRES-SDC) throughout the revised paper.

Near the start of section 3, there appears “Although most of today’s efforts focus on Tokamak design…”. It would be wise to include stellarators in a description of today’s efforts.

Answer**.**

Following comments by Reviewer #3, to streamline the paper in order to reduce its lengths we have removed this part. Given that the paper’s focus is not on reactor design, we felt it was unnecessary.

Bullet point with further typographical notes.

Answer**.**

Those have all been fixed or the highlighted equations/paragraphs no longer appear in the revised manuscript.

Reviewer # 3.

(i) With more than 30 pages, the paper is too long for the novel contributions. The first 12 pages appear to be mostly a recollection of the predecessor [17] in different but quite standard formulation. Pages 11 and 12 rephrase the SDC differently, essentially in the way of [17] (to which the reader is related for details), and do not really help the understanding. The novel methodology is presented in pages 13-16 only. The remaining pages are devoted to interesting and instructive numerical tests, with a rather verbose presentation. I think the text could be shortened significantly to, say, around 24 pages, and I believe this would be beneficial for the paper.

Answer**.**

We have removed the whole subsection 2.2, which was indeed not critical for understanding the method. After introducing the collocation formulation, the revised manuscript now proceeds directly to the matrix formulation of Boris-SDC and BGSDC. Since removing 2.2 also removed the component-wise formulation of the Boris-SDC sweep, we have added the very similar equations used to apply the preconditioner within BGSDC instead, which seems more relevant to the paper since this is what is actually implemented in the provided code.

We have cut some superfluous comments and text at other places. When using final option instead of preprint in the LaTeX template, the paper is now 23 pages long, 2 of which are references. While this is probably still on the longer side, we think that trimming down the description of the method further would make the paper difficult to read while removing numerical examples would weaken the evidence for the method’s usefulness.

(ii) The authors advocate a particular three-step handling of mildly nonlinear problems, by (a) computing a provisional approximation with a Boris integrator, (b) solving a linearized collocation problem with a GMRES-accelerated Boris-SDC iteration, and (c) performing a couple of Picard iterations for the nonlinear collocation problem. Unfortunately, they do not give a satisfactory rationale for their design. Immediate questions that arise are:

Why Picard iterations and not Boris sweeps? (On p 24 they state that Picard iterations are cheaper and can be parallelized, but then they probably give less accurate results. The trade-off is unclear.)

Answer**.**

We tested a variant where GMRES-SDC is followed by fully nonlinear Boris-SDC sweeps instead of Picard iterations. Table 1 shows the resulting errors. Except for the $k=1$ case, there doesn’t seem to be much benefit. This is not enough of an improvement to justify the higher computational cost. However, there may well be situations (e.g., if the nonlinearity is stronger), where full sweeps outperform Picard iterations. We have revised the corresponding paragraph on page 11, lines 164–171 to acknowledge this.

Why not perform the GMRES part (realizing a truncated Newton step) at the end? The benefit would be that close to the solution Newton convergence is much faster than SDC or Picard convergence, such that the linearization error is probably negligible. I’d guess that (a) provisional solution (b) Boris-SDC sweeps (c) GMRES would be more efficient. I may be wrong, but this topic is not discussed at all.

Answer**.**

Following up on the reviewer’s remarks, we implemented a fully nonlinear Newton version of BGSDC, using GMRES with a fixed number of iterations as an inner solver. For the problems studied in the paper, we found it to be not competitive with BGSDC + Picard. Table 1 shows the final energy error for BGSDC with Picard iterations and the Newton-version for 1000 steps for the Solev’ev equilibrium with $\Delta t=5\times 10^{-9}$. The first column indicates the number of iterations $k$: BGSDC uses $k$ iterations in GMRES, followed by $k$ Picard iterations whereas Newton-BGSDC uses $k$ out iterations and $k$ GMRES iterations in every linear solve. Note that this makes Newton-BGSDC much more computationally costly, since it effectively uses $k^{2}$ GMRES iterations whereas BGSDC+Picard uses $k$ GMRES + $k$ Picard iterations. Nevertheless, for $k=1$ and $k=2$, Newton-BGSDC delivers less precision than BGSDC with Picard. For $k=3$ and, as reference, $k=12$, both methods have converged and deliver the error of the underlying collocation method. Similar results were found for other parameters. While Newton-BGSDC may be competitive for other problems, e.g. with strong nonlinearities, using Picard iterations provides better results for the problems studied in the paper.

We have modified the paragraph on page 10, lines 152–153 to reflect this.

Similarly, why Picard at the end and not just a second round of GMRES on a new linearization point? This could exploit the fast Newton convergence close to the solution, and GMRES-SDC converges probably faster than Picard.

Answer**.**

This relates to the answer to the previous point: the issue with Picard iterations are not their convergence speed but their very limited domain of convergence. However, for the accurate starting value provided by GMRES, this is not a problem and we therefore use Picard iterations as the cheapest way to adjust for the nonlinearity.

Algorithm 1: It’s hard to understand the essence of Boris’ trick from the algorithm. I’d rather have a few lines showing why $v=a+v\times B$ implies $v=(a+a\times B+(a^{T}B)B)/(1+B^{T}B)$.

Answer**.**

We agree that the geometric motivation for the trick is difficult to understand from the algorithm. However, our intent was mainly to provide the necessary information for a quick implementation so that a reader can reproduce the results, not necessarily to derive Boris’ trick in full. To help a reader who is interested in the motivation behind the algorithm, we have added a reference to Section 4-4 in the book by Birdsall and Langdon, where a detailed geometric derivation of Boris’ trick can be found.

p 24, last line: The argument that Picard is cheaper than Boris sweeps is valid, but does not match the effort metric of rhs evaluation used throughout the paper. This should probably be put into perspective more clearly.

Answer**.**

This has been rephrased to be clearer, see page 20, lines 311–314.

Minor issues in eq (4a), eq (7), eq (13), eq (14).

Answer**.**

These have either been fixed or the corresponding equations no longer appear in the revised manuscript.