Fast low-rank approximations of multidimensional integrals in ion-atomic collisions modelling
aa r X i v : . [ phy s i c s . c o m p - ph ] A p r NUMERICAL LINEAR ALGEBRA WITH APPLICATIONS
Numer. Linear Algebra Appl. Fast low-rank approximations of multidimensional integrals inion-atomic collisions modelling
M. S. Litsarev ∗ , I. V. Oseledets , Skolkovo Institute of Science and Technology, Novaya St. 100, Skolkovo, 143025 Moscow Region, Russia Institute of Numerical Mathematics, Gubkina St. 8, 119333 Moscow, Russia
SUMMARYAn efficient technique based on low-rank separated approximations is proposed for computation of three-dimensional integrals arising in the energy deposition model that describes ion-atomic collisions. Directtensor-product quadrature requires grids of size which is unacceptable. Moreover, several of suchintegrals have to be computed simultaneously for different values of parameters. To reduce the complexity,we use the structure of the integrand and apply numerical linear algebra techniques for the constructionof low-rank approximation. The resulting algorithm is faster than spectral quadratures in sphericalcoordinates used in the original DEPOSIT code. The approach can be generalized to other multidimensionalproblems in physics. Copyright c (cid:13) Received . . .
KEY WORDS: Low-rank approximation; 2D cross; Separated representation; Exponential sums; 3DIntegration; Slater wave function
1. INTRODUCTIONComputation of multidimensional integrals is one of the most time-consuming tasks in physics.Standard approaches either have high complexity or require sophisticated quadrature schemes.Already in a three-dimensional case the integrand may depend on many parameters and shouldbe computed many times, so the computational cost of the simplest tensor-product quadratures is animportant issue.One of the promising tools to reduce the dimensionality of the problem (and hence the numberof mesh points where the integrand should be calculated) is the usage of separation of variables in the integrand (see, for example, [1, 2]). The idea was known for a long time (for examples werefer the reader to the review [3]), but it has become practically useful after the fast algorithms ofdecompositions of functions in a separated form had been obtained in two- [4, 5], three- [6, 7] andmultidimensional cases [8, 9].Let F ( x, y ) be a function of two variables x, y where point ( x, y ) is in a certain rectangle [ a x , b x ] ⊗ [ a y , b y ] . The function is said to be in a separated form if it can be represented as a sum ofproducts of univariate functions: F ( x, y ) = q X τ =1 σ τ u τ ( x ) g τ ( y ) . (1)The minimal number q such that (1) exists is called separation rank . Direct generalizationof (1) to multivariate functions is referred to as a canonical polyadic (CP, also known as ∗ Correspondence to: [email protected] c (cid:13)
Prepared using nlaauth.cls [Version: 2010/05/13 v2.00]
CANDECOMP/PARAFAC) decomposition [3, 10]. The reader can find examples of applicationsof separated representations in multidimensional cases in [1, 11, 12, 13, 14, 15, 16, 17, 18, 19].For a function given in the separated form the integration is simplified a lot. Indeed,
Z Z F ( x, y ) dxdy = q X τ =1 σ τ Z b x a x u τ ( x ) dx Z b y a y g τ ( y ) dy, (2)and the problem is reduced to the computation of one-dimensional integrals, which can be computedusing fewer quadrature points than the original integral.In case of a discrete approximation of (1) F ( x i , y j ) ≈ q X τ =1 σ τ u τ ( x i ) g τ ( y j ) (3)with the error ε in the Frobenious norm the number q is called ε -rank. We shall assume the notionof ε -rank when the term rank will be mentioned in the text bellow. Expression (3) can be rewrittenin the matrix form A ≈ U Σ G ⊤ , (4)where A is an n × m matrix with elements A ij = F ( x i , y j ) , U is an n × q matrix with elements U iτ = u τ ( x i ) , G is an m × q matrix with elements G jτ = g τ ( y j ) and Σ is a q × q diagonal matrixwith elements σ τ on the diagonal. This is a standard low-rank approximation problem for agiven matrix. Provided that a good low-rank approximation exists, there are very efficient crossapproximation algorithms [4, 5] that need only O (( n + m ) q ) elements of a matrix to be computed.In this paper we describe how to apply this technique to speedup computations of three-dimensional integrals in the energy deposition model intended to describe ion-atomic collisions.This model was introduced by N. Bohr [20] and developed further by A. Russek and J. Meli [21],C.L. Cocke [22], and V.P. Shevelko at al. [23]. Theoretical development of the model is presentedin [23, 24, 25, 26]. Examples of calculations are reported in [27, 28, 29, 30]. Detailed description ofthe computer code DEPOSIT and user guide are given in [31] and its updated version based on theseparated representations is avialable in [32].The code DEPOSIT allows to calculate total and multiple electron loss cross sections andionization probabilities needed for estimation of losses and lifetimes of fast ion beams, backgroundpressures and pumping requirements in accelerators and storage rings. All of them are, in fact,functionals of the deposited energy T ( b ) ( b is the impact parameter of the projectile ion), which inturn is a three-dimensional integral over the coordinate space. To calculate any of these parametersone has to compute T ( b ) in all points of the b -mesh.In the original work [31] an advanced quadrature technique was used, and the computational timehas appeared to be much faster in a comparison with direct usage of uniform meshes. Calculation of the deposited energy T ( b ) for a given atomic shell in one point b took about several seconds.For complex ions it was about few hours on one processor core in total, that is not enough fast.To overcome this issue a fully scalable parallel variant of the algorithm was proposed, but thecomputational time was still large.We present an entirely different approach for the computation of T ( b ) in many points of the b -mesh based on the idea of separation of variables (1). An approximation of functions to be integratedby a sum of products of univariate functions allows to effectively decrease the dimensionality of theproblem. This requires active usage of numerical and analytical tools.The problem setting is as follows. The deposited energy T ( b ) is a three-dimensional integral T ( b ) = Z Z Z ∆ E ( x, z − b ) ρ ( x, y, z ) dxdydz. (5) It involves cylindrically symmetric function of two variables (the energy gain ∆ E during an ion-atomic collision) which is a smooth finite function and spherically symmetric function of threevariables (electron density in a Slater-type approximation) which decays exponentially. For details Copyright c (cid:13)
Numer. Linear Algebra Appl. (0000)
Prepared using nlaauth.cls
DOI: 10.1002/nla and definitions we refer the reader to Appendix A. Previous approach [31] used tensor-productquadratures in spherical coordinates. We use very fine uniform meshes and low-rank approximation.In Section 2.1 the Slater-type function of three variables is decomposed by the exponential sumsapproach [33, 34]. The integral is immediately reduced to a two-dimensional one of a simplerstructure. In Section 2.2 for a function of two variables we use the pseudo-skeleton decomposition ofmatrices [35, 36, 37] computed via a variant of the incomplete cross approximation algorithm [4].Combining these two representations we obtain then an efficient algorithm with O ( n ) scaling, incomparison with O ( n ) complexity for direct integration over a three-dimensional mesh. We shownumerically that the function in question can be well-approximated by a separable function. Thus,the approximation can be computed in O ( n ) time, where n is the number of grid points in onedimension. The computation of T ( b ) on the whole b -mesh takes less then one minute instead ofseveral hours and total speedup of the program is about ∼ times. Illustrative examples are givenin Section 2.3. All the equations related to the physical model are written in atomic units.2. NUMERICAL PROCEDURE In this section we present analytical expansion of the spherically symmetric electron density ρ γ ( r ) inCartesian coordinates as a sum of separable functions. We use this decomposition later in Section 2.2for analytical integration in one dimension.A three-dimensional electron density ρ γ ( r ) is taken in a Slater-type approximation ρ γ ( r ) = C γ r α γ e − β γ r , r = p x + y + z , (6)where integer parameter α γ and real C γ and β γ correspond to one atomic shell labeled by γ (seealso Appendix B for description of parameters). For the density ρ γ ( r ) the following normalizationcondition occurs Z ∞ ρ γ ( r ) dr = N γ , (7)where N γ is the number of electrons in a γ -shell. For the sake of simplicity index γ will be skippedand only one shell will be considered in equations bellow.For a function ρ ( x, y, z ) defined in (6) the separation of variables (1) in Cartesian coordinates canbe done analytically [33, 34, 38, 39]. The main idea is to approximate the Slater density function bya sum of Gaussians ρ ( r ) ≈ K X k =0 λ k e − η k r . (8) Once the approximation (8) is computed, the separation of variables in Cartesian coordinates isimmediately done ρ ( x, y, z ) ≈ K X k =0 λ k e − η k x e − η k y e − η k z . (9)The technique for computation of nodes λ k and weights η k is based on the computation of inverseLaplace transform. Let us consider a function f αβ ( t ) such that its Laplace transform is a function F αβ ( s ) : F αβ ( s ) = ρ ( √ s ) C = (cid:0) √ s (cid:1) α e − β √ s = Z ∞ e − st f αβ ( t ) dt, (10) where α and β are parameters of Slater density (6). Inverse Laplace transform f αβ ( x ) can becomputed analytically for known F αβ ( s ) (10). In Appendix B we present explicit expressions forfunctions f αβ ( t ) corresponding to function (10) for integer α and real positive β . Copyright c (cid:13)
Numer. Linear Algebra Appl. (0000)
Prepared using nlaauth.cls
DOI: 10.1002/nla
Once function f αβ ( t ) in expression (10) is known, the integral (10) can be computed numericallyand approximated by a quadrature formula ρ ( r ) ≈ C K X k =0 w k e t k f αβ ( e t k ) e − r e tk . (11)Here w k and e t k are quadrature weights and nodes, respectively. The procedure to compute weightsand nodes was taken from the paper [34]. For the reader’s convenience we give the formula and itsderivation in Appendix C.Comparision with equation (8) gives λ k = C w k e t k f αβ ( e t k ) , η k = e t k . (12)It appears that not many quadrature points (at fixed r ) are required to achieve the accuracy of theexpansion of order − in the Chebyshev norm due to the hyper-exponential decay. Practically,such an accuracy is quite enough for the physical meaning of the model. However, r is going tobe very small, as soon as the finer grids (i.e. for a large number of nodes) are required for higher precision. Estimation of upper bound K in sum (11) which has the logarithmic dependence is givenin Appendix D. T ( b ) . In this section we describe a numerical scheme based on the cross decomposition of two dimensionalintegrand. Dimensionality reduction (from three to two dimensions) is achieved by means of theseparated representation of Slater density obtained in the previous Section.Discretization of one-dimensional integrals in (2) by some quadrature formula with nodes x i ∈ [ a x , b x ] , i = 1 , . . . , n , y j ∈ [ a y , b y ] , j = 1 , . . . , m and weights w ( x ) i , w ( y ) j , leads to the approximation Z Z F ( x, y ) dxdy ≈ q X τ =1 σ τ n X i =1 w ( x ) i u τ ( x i ) m X j =1 w ( y ) j g τ ( y j ) . (13)To get the decomposition (3) we apply algorithm [4, 40] implemented in [32].A three-dimensional integral T ( b ) defined in (5) can be reduced to a two-dimensional integral bymeans of the decomposition (9) T ( b ) = K X k =0 λ k Z Z Z ∆ E ( x, z − b ) e − η k x e − η k y e − η k z dxdydz (14)and analytical evaluation of one-dimensional Gaussian integral Z ∞−∞ e − η y dy = r πη , (15) T ( b ) = √ π K X k =0 λ k √ η k Z Z ∆ E ( x, z − b ) e − η k x e − η k z dxdz. (16)Suppose that ∆ E ( x, z − b ) has been decomposed as follows ∆ E ( x, z − b ) ≈ q X τ =1 σ τ u τ ( x ) g τ ( z − b ) . (17) Then the integration (16) can be reduced to a sequence of one-dimensional integrations. T ( b ) = √ π K X k =0 λ k √ η k q X τ =1 σ τ I τk J τk ( b ) , (18) Copyright c (cid:13)
Numer. Linear Algebra Appl. (0000)
Prepared using nlaauth.cls
DOI: 10.1002/nla I τk = Z b x a x u τ ( x ) e − η k x dx, (19) J τk ( b ) = Z b y a y g τ ( z − b ) e − η k z dz. (20)For the numerical approximation of integrals (19) and (20) we use quadrature formula with uniformquadrature nodes (although any suitable quadrature can be used) I τk = X i ∈ Ω xǫ ( k ) w ( x ) i u τ ( x i ) e − η k x i , Ω xǫ ( k ) = { i | e − η k x i > ǫ } (21) x i = − a x + i h x , ≤ i ≤ N x , h x = a x /N x , (22) J τk ( b ) = X j w ( z ) j g τ ( z j − b ) e − η k z j , (23) z j = − a z + j h z , ≤ j ≤ N z , h z = a z /N z . (24)We sample the impact parameter b (which can take only positive values) with the same step h z b l = l h z , ≤ l ≤ N z . (25)This allows us to introduce a new variable ˜ z = z − b discretized as follows ˜ z k = − a z + k h z , ≤ k ≤ N z , (26)and such that for the boundary conditions (24), (25), (26) z j − b l = ˜ z j − l + N z . (27)Thus, the approximation problem (17) reduces to the low-rank approximation of the extended (2 N x + 1) × (3 N z + 1) matrix ∆ E ( x i , ˜ z j ) = r X τ =1 σ τ u τ ( x i ) g τ (˜ z j ) . (28)This should be done only once (using the cross approximation algorithm), and the finalapproximation of integral (23) reads J τk ( b l ) = X j ∈ Ω ˜ zǫ ( k ) w (˜ z ) j g τ (˜ z j − l + N z ) e − η k ˜ z j . (29)The calculation of T ( b ) can be summarized in the following algorithm. for every γ -shell of the projectile ion do compute the decomposition (8) for ρ ( r ) compute the cross approximation for the matrix ∆ E ( x i , ˜ z j ) defined in (28) for k = 0 . . . K do for τ = 1 . . . q do compute the integral I τk defined in (21) for every b l required do for k = 0 . . . K do for τ = 1 . . . q do compute the integral J τk ( b l ) defined in (29) compute T γ ( b l ) , equation (18) Copyright c (cid:13)
Numer. Linear Algebra Appl. (0000)
Prepared using nlaauth.cls
DOI: 10.1002/nla
Finally, the question is how to compute T ( b ) for many different values of b . The direct summationrequires O ( N ) operations. But a closer look reveals that it is in fact a discrete convolution, whichcan be computed in linear cost. Nevertheless, due to the exponential decay of factors e η k x i and e η k ˜ z j in sums (21) and (29) correspondingly, there is a very small number of terms, such that a directsummation after the truncation is much faster. This can be easily seen from values of parameter θ x defined as θ x = P Kk =0 N (Ω xǫ ( k ))(2 N x + 1)( K + 1) , (30)where N ( X ) is the cardinality of set X . For details, please, see Table II and discussion in thefollowing Section 2.3. Once the analytical expansion (9) is obtained, the full calculation of T ( b ) consists of two steps (seealgorithm in Section 2.2): calculation of the cross decomposition of integrand (16) on the extended ˜ z -mesh (26) and computation of all integrals I τk and J τk ( b l ) for all b l from (25). There is no need inall the time recalculation of the cross approximation (28) for the full computation of these integrals.After it is computed (that is σ τ , u τ ( x i ) and g τ (˜ z j ) are known), the calculation of I τk and J τk ( b l ) starts.In Table I we present times T cross and ranks q calculated for the energy gains ∆ E γ ( x, ˜ z ) corresponding to two ion-atomic collision examples ( Au + O and U + Xe ). Values of T cross are reported for different sizes of ( x, ˜ z ) -mesh and the relative accuracy of the cross decomposition.One can find that T cross is about − ∼ second in order of magnitude. Given that thefull computation of integrals I τk and J τk ( b l ) for all b l takes approximately · · − = 25 seconds (the average value of column T S in Table II for b values), we can conclude thatthe cross approximation becomes then a pre-computing stage with a tiny contribution to the totalcomputational time.An important parameter in sum (28) is the rank value q . It determines the complexity of thealgorithm (the smaller q , the better). As it follows from the numerical experiments, ranks of ∆ E γ ( x, ˜ z ) decomposition are small (Table I) against the mesh size. It means that for real physicalsystems the cross decomposition is a very prominent tool. It allows to decrease the complexity ofthe problem in practice from O ( n ) elements to O ( q · n ) elements where q ≪ n .In Table II we present the program speedup for every atomic shell. In quadrature sum (21)all terms less then ǫ = 10 − were thrown out for every x i due to the exponential decay offactors e − η k x i for every fixed η k . Parameter θ x is defined as a relative number of terms in totalsummation of I τk over all k and i (equation (21) and (30)), above the threshold ǫ . As it can be seenfrom the table, there are only a few percent of terms to be summed, which considerably reduces thesum and speeds up the total calculation. Parameter θ x decreases when going from top to bottom ofthe third column (for one system). That is why the time T S also decreases while both ranks K and q increase. The same situation occurs for sum (29) over ˜ z j .Another important question is the full accuracy of the computation. As it was mentioned above, itconsists of two contributions acquired from the cross approximation and the quadrature summations.In equations (21) and (29) values σ τ , u τ ( x i ) and g τ (˜ z j ) are approximated with the cross accuracy ε c ,while the quadrature sum is calculated with an error ε i . In Table III and Table IV we provide anumerical example demonstrating the actual error for the integral T ( b ) for small enough mesh sizes.We consider the worst case b = 0 , when the integrand has the most singular behavior (because theSlater density has no higher derivatives at r = 0 ). As it follows from our results the scheme holdsthe third order for the Simpson rule, even in the worst singular case. In other cases b > it holds thefourth order. These results show that the claimed accuracy is achieved and the numerical schemeholds the order of integration up to the accuracy of the approximated function.Finally, we can conclude that usage of the technique based on separated representations (18) allows to decrease the total time of computing of T ( b ) by a factor of ∼ in comparison to theprevious version. In practice, computational time is reduced from several hours to one minute orless on the same hardware. Copyright c (cid:13)
Numer. Linear Algebra Appl. (0000)
Prepared using nlaauth.cls
DOI: 10.1002/nla
Table I. Ranks q of the decomposition (28) calculated by incomplete cross approximation algorithm [4]for the energy gain ∆ E ( x, ˜ z ) . Two cases are considered: collision of Au ions with an Oxigen atom ata collision energy E = 6 . MeV/u and collision of U ions with a Xenon atom at a collision energy E = 2 . MeV/u. The number of x -mesh points are taken equal to N + 1 , the number of ˜ z -mesh pointsis taken equal to N + 1 in correspondence to equations (22) and (26), a x = a z = 8 . The relative error ofthe approximation in the Frobenius norm ε and N are placed in the bottom line as a couple ( ε, N ) . Theycorrespond to different numerical tests. Calculations were carried out on . GHz Intel Core i5 processor.Column T cross corresponds to the time the cross algorithm takes. The numerical results confirm almostlinear scaling of the approach in N . System γ -Shell q T cross (sec) q T cross (sec) q T cross (sec) Au + O df
13 0 .
21 21 0 .
24 2.41 sp
13 0 .
21 21 0 .
24 2.40 d
14 0 .
19 22 0 .
26 2.58 sp
16 0 .
25 24 0 .
29 2.64 sp
17 0 .
28 25 0 .
30 2.71 sp
17 0 .
27 25 0 .
30 2.70 U + Xe sp
14 0 .
20 22 0 .
26 2.17 df
15 0 .
23 24 0 .
27 2.58 sp
17 0 .
28 25 0 .
30 2.69 d
17 0 .
27 25 0 .
30 2.77 sp
17 0 .
27 25 0 .
30 2.75 sp
17 0 .
26 25 0 .
30 2.71 sp
17 0 .
26 25 0 .
30 2.69( − , ) ( − , ) ( − , )These results were obtained by using of our implementation of the cross approximation algorithmand the low-rank format of the Slater density (11). The latest version of the cross decompositioncode implemented in C ++ can be downloaded from the link [40].3. CONCLUSIONS AND FUTURE WORKA new and efficient technique for computation of three-dimensional integrals based on low-rankand separated representations is proposed for the energy deposition model. Due to a general formof the integrand, which is a product of two functions with cylindrical and spherical symmetries,this methodology can be applied to many types of integrals having similar structure. Such anapproach significantly reduces computational time and allows to achieve a given accuracy (because it uses the cross approximation and quadrature summations). The general concept can be appliedto more complicated models (like ion-molecular collisions with electron loss and charge-changingprocesses) that lead to multidimensional integrals. For the multidimensional case we plan to use fastapproximation techniques based on the tensor train (TT) format [9]. ACKNOWLEDGEMENT
This work was partially supported by Russian Science Foundation grant 14-11-00659.
A. PHYSICAL MODEL
In this section we introduce functions and parameters involved into the definition of deposited energyintegral (5) T γ ( b ) = Z Z Z ∆ E γ ( p ) ρ γ ( x, y, z ) dxdydz, (31) Copyright c (cid:13)
Numer. Linear Algebra Appl. (0000)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nla
Table II. Timings to compute T ( b ) at fixed b are presented for two cases: the DEPOSIT code (old) T D and the code based on separated representations (18) T S . Collision systems are the same as in Table I.In the third column K labeles maximum value of index in expansion (9). The mesh (54) for radialdensity is taken with a t = − , b t = 45 , h t = ( b t − a t ) / . The meaning of parameter θ x is explainedin Section 2.3. Calculations were carried out for the relative accuracy of the cross decomposition ε = 10 − and [ − , ⊗ [ − , mesh with × points. The last column shows speedup of the program. System γ -Shell K θ x × − T S ( × − sec) T D (sec) T D /T S Au + O df
74 9 . .
94 3 .
89 4904 sp
69 6 . .
92 3 .
83 7783 d
73 4 . .
59 3 .
88 10803 sp
72 3 . .
81 3 .
82 10032 sp
107 1 . .
42 3 .
86 15921 sp
209 0 . .
24 3 .
88 3120 U + Xe sp
62 13 . . .
94 3904 df
70 6 . .
05 3 .
90 6444 sp
67 5 . .
00 3 .
94 7883 d
71 3 . .
88 3 .
92 10113 sp
70 3 . .
52 3 .
90 11062 sp
105 1 . .
99 3 .
87 19451 sp
207 0 . .
04 3 .
88 3723
Table III. Convergence of integrals T ( b ) for df , sp and d shells of Au + O at . MeV/u on two-dimensional mesh [ − , ⊗ [ − , with (2 N + 1) × (3 N + 1) points for different N (first column),see equations (22) and (26). Simpson weights w xi , w (˜ z ) j are used in (21) and (29). For radial density themesh (54) is taken with parameters a t = − , b t = 45 , h t = ( b t − a t ) / . Calculations were carried out forthe relative accuracy of the cross decomposition ε c = 10 − and fixed value of parameter b = 0 . . Lastcolumn shows the rellative error ε i . Extrapolated value of integral is calculated by Romberg method (withRichardson extrapolation) [41], p. 161. The order of scheme p e is defined by Aitken rule [41], p. 344. γ -Shell N q T ( b ) p e ε i df . . · − . . · − . .
07 7 . · − . .
01 4 . · − . .
00 3 . · − . .
02 2 . · − extrapolated . sp . . · − . . · − . .
06 8 . · − . .
01 5 . · − . .
00 3 . · − . .
01 2 . · − extrapolated . d . . · − . . · − . .
28 5 . · − . .
17 6 . · − . .
11 8 . · − . .
06 9 . · − . .
04 1 . · − . .
03 1 . · − extrapolated . Copyright c (cid:13)
Numer. Linear Algebra Appl. (0000)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nla
Table IV. Convergence of integrals T ( b ) for sp , sp and s shells of Au + O at . MeV/u on two-dimensional mesh [ − , ⊗ [ − , with (2 N + 1) × (3 N + 1) points for different N (first column),see equations (22) and (26). Simpson weights w xi , w (˜ z ) j are used in (21) and (29). For radial density themesh (54) is taken with parameters a t = − , b t = 45 , h t = ( b t − a t ) / . Calculations were carried out forthe relative accuracy of the cross decomposition ε c = 10 − and fixed value of parameter b = 0 . . Lastcolumn shows the rellative error ε i . Extrapolated value of integral is calculated by Romberg method (withRichardson extrapolation) [41], p. 161. The order of scheme p e is defined by Aitken rule [41], p. 344. γ -Shell N q T ( b ) p e ε i sp . . · − . . · − . .
69 6 . · − . .
91 7 . · − . .
98 9 . · − . .
99 1 . · − . .
00 1 . · − . .
00 1 . · − . .
01 2 . · − . .
01 4 . · − extrapolated . sp . . · − . . · − . .
64 4 . · − . .
91 5 . · − . .
98 7 . · − . .
99 9 . · − . .
00 1 . · − . .
00 1 . · − . .
00 1 . · − . .
00 2 . · − . .
02 3 . · − extrapolated . sp . . · − . . · − . .
31 7 . · − . .
40 8 . · − . .
29 9 . · − . .
17 1 . · − . .
09 1 . · − . .
05 1 . · − . .
02 2 . · − . .
01 2 . · − . .
01 3 . · − extrapolated . p = p ( z − b ) + x , (32)where the integration is done over the whole coordinate space. Integral (31) is written for an atomic shellwith principal quantum number n and orbital quantum number l labeled by γ = nl .The energy gain ∆ E γ ( p ) is a kinetic energy deposited to the projectile’s γ -shell by the field U ( R ) of thetarget atom. ∆ E γ ( p ) = ∆ E <γ ( p ) n f ( v γ − ϑ ) + ∆ E >γ ( p ) n f ( ϑ − v γ ) , (33) n f ( x ) = 1 e − kx + 1 , (34) Copyright c (cid:13)
Numer. Linear Algebra Appl. (0000)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nla where k = 3 by default and is an input parameter of the model. Expression (33) consists of two termscorresponding to low ∆ E <γ ( p ) = ξ γ p + ζ γ X i =1 φ i F ( χ i p ) (35)and high ∆ E >γ ( p ) = ν γ p + µ γ X i =1 φ i F ( χ i p ) ! , (36)energies. Function F ( x ) = xK ( x ) = Z ∞ e − √ x + y dy (37)is defined in terms of the modified Bessel function of the second kind K ( x ) (see [42], p. 375, Eq. 9.6.2)with asympthotical behavior ([42], p. 378, Eq. 9.8.3 and p. 379, Eq. 9.8.7) F ( x → + ) = 1 + (cid:16)
12 ln x . (cid:17) x + O ( x ) . (38)Fitting parameters φ i , χ i of the atomic field U ( R ) are obtained from the Dirac-Hartree-Fock-Slatercalculations [43]. The distance between the center of the field and the projectile electron of γ -shell is labeledby R . Atomic field is given in the Yukawa potential form U ( R ) = − ZR X i =1 φ i e − χ i R , (39)where Z is the nuclear charge of the target atom.Parameter ϑ is a relative velocity of the projectile. It is related with the collision energy of the ion-atomicsystem. The rest parameters v γ , ξ γ , ζ γ , ν γ and µ γ concern to fundamental properties of the projectile ionand the target atom (such as binding energy, average velocity, mean radius of the shell, atomic radius andcharge). They should be considered as positive constants in the integral (31) for a given ion-atomic systemfor all b . Detailed description of how to calculate them can be found in [31]. Examples of input files with theoriginal code can be downloaded from link [44]. B. INVERSE LAPLACE TRANSFORM SOURCES
For integer α and real positive β inverse Laplace transform f αβ ( t ) of F αβ ( s ) from equation (10) may becalculated analytically and expressed via the Kummer’s confluent hypergeometric function M ( a, b ; z ) ([42],chapter 13) as follows f αβ ( t ) = M (cid:16) α , , − β t (cid:17) t α Γ (cid:0) − α (cid:1) − β M (cid:16) α , , − β t (cid:17) t α Γ (cid:0) − α (cid:1) , (40)where M ( a, b ; z ) = 1 + ab z
1! + a ( a + 1) b ( b + 1) z
2! + . . . (41)and Γ( x ) is the Gamma function .Below we present f αβ ( t ) explicitly for most practically usefull cases ( α = 0 , , . . . , ). Due to thedifference of normalization conditions in spherical and Cartesian coordinates for the Slater density (6) ρ ( r ) = N γ (2 β ) µ +1 Γ(2 µ + 1) r µ e − βr , (42)parameter α in (6) is related to parameter µ in (42) as follows α = 2 µ − . (43)The number of electrons in the shell γ is labeled as N γ . Parameter µ is greater or equal to unity. It is aninteger or half-integer depending on the principal quantum number n and the orbital quantum number l Copyright c (cid:13)
Numer. Linear Algebra Appl. (0000)
Prepared using nlaauth.cls
DOI: 10.1002/nla of the atomic shell. Details can be found in [45, 46]. For example, µ s = 1 , α = 0 ; µ sp = 2 , α = 2 ; µ d = 3 . , α = 5 . Finally, f β ( t ) = g (cid:0) t/β (cid:1) √ π β , g ( t ) = e − t t / (44) f β ( t ) = g (cid:0) t/β (cid:1) √ π β , g ( t ) = − e − t t / (cid:16) − t (cid:17) (45) f β ( t ) = 3 g (cid:0) t/β (cid:1) √ π β , g ( t ) = − e − t t / (cid:16) − t (cid:17) (46) f β ( t ) = 3 g (cid:0) t/β (cid:1) √ π β , g ( t ) = e − t t / (cid:16) − t + 43 t (cid:17) (47) f β ( t ) = 15 g (cid:0) t/β (cid:1) √ π β , g ( t ) = e − t t / (cid:16) − t + 415 t (cid:17) (48) f β ( t ) = 15 g (cid:0) t/β (cid:1) √ π β , g ( t ) = − e − t t / (cid:16) − t + 4 t − t (cid:17) (49) f β ( t ) = 105 g (cid:0) t/β (cid:1) √ π β , g ( t ) = − e − t t / (cid:16) − t + 45 t − t (cid:17) (50)For practical reasons higher values of α are not necessarily due to the limitation of shell filling with electrons. C. QUADRATURE FORMULA FOR THE LAPLACE INTEGRAL
To obtain the decomposition (8) for given α and β we make a substitution s → s into the equation (10) F αβ ( s ) = s α e − βs = Z ∞ e − s x f αβ ( x ) dx, (51)then introduce another variable x = e t in the integral F αβ ( s ) = s α e − βs = Z ∞−∞ e − s e t + t f αβ ( e t ) dt. (52)The function under the integral (52) has exponential decay both in the spatial and frequency domains.Therefore the truncated trapezoidal rule gives the optimal convergence rate. It turns out to the finalapproximation of the form F αβ ( s ) ≈ K X k =0 w k e t k f αβ ( e t k ) e − s e tk (53)with trapezoidal weights w k and the integrand values in the nodes e − s e tk + t k f αβ ( e t k ) . The Gaussian-likepart is split out in a separate factor in correspondence with decomposition (8). Parameters of the formula t k = a t + kh t , h t = ( b t − a t ) /K (54)have to be selected in such a way that the resulting quadrature formula approximates the integral for a widerange of parameter s . Typically, the choice a t & − , b t . , and K ∼ gives good accuracy ( ≤ − ).As an example, in Table II the required number of terms in sum (53) is presented. Accurate error analysiscan be found in [34]. Copyright c (cid:13)
Numer. Linear Algebra Appl. (0000)
Prepared using nlaauth.cls
DOI: 10.1002/nla K ( N ) FOR THE SLATER DENSITY SERIES
In this section we estimate an upper bound K in sum (11) which has the logarithmic dependence on themesh size N . We follow the proof given in paper [47] for Lemma 4.3. In this Lemma a Slater-type function ρ ( √ s ) from eq. (10) is considered for α = 0 . Below, this function is considered for integer α .The integrand in (10) after the change of variables t = e x P αβ ( x ) = e − se x + x f αβ ( e x ) (55)has the following decay on the real axis (skipping the numerical factor before the exponent) P αβ ( x ) ≈ e − se x − c x , as x → ∞ , c = α + 1 − ( α mod 2)2 , (56) P αβ ( x ) ≈ e − β e | x | + c | x | , as x → −∞ , c = α + 12 . (57)This immediately implies expression (5.3) from [47] for b = min( β , s ) , b is taken in the notation of [47].Following then the statement I of Lemma 4.3 we may conclude, that K = O ( | log ε | ( | log ε | + log 1 /b )) withthe error ε of the approximation and b ∼ /N . REFERENCES1. Beylkin G, Mohlenkamp MJ. Numerical operator calculus in higher dimensions.
Proc. Nat. Acad. Sci. USA (16):10 246–10 251, doi:10.1073/pnas.112329799.2. Khoromskaia V, Khoromskij BN, Schneider R. Tensor-structured factorized calculation of two-electron integrals ina general basis. SIAM J. Sci. Comput. (2):A987–A1010, doi:10.1137/120884067.3. Kolda TG, Bader BW. Tensor decompositions and applications. SIAM Review (3):455–500, doi:10.1137/07070111X.4. Tyrtyshnikov EE. Incomplete cross approximation in the mosaic–skeleton method. Computing (4):367–380, doi:10.1007/s006070070031.5. Bebendorf M. Approximation of boundary element matrices. Numer. Mathem. (4):565–589, doi:10.1007/pl00005410.6. Oseledets IV, Savostianov DV, Tyrtyshnikov EE. Tucker dimensionality reduction of three-dimensional arrays inlinear time. SIAM J. Matrix Anal. Appl. (3):939–956, doi:10.1137/060655894.7. Rakhuba MV, Oseledets IV. Fast multidimensional convolution in low-rank formats via cross approximation. SIAMJ. Sci. Comput.
SIAM J. Sci. Comput. (5):3744–3759, doi:10.1137/090748330.9. Oseledets IV. Tensor-train decomposition. SIAM J. Sci. Comput. (5):2295–2317, doi:10.1137/090752286.10. Harshman RA. Foundation of the parafac procedure: Model and conditions for an explanatory multimode factoranalysis. UCLA Working papers in phonetics :1–84.11. Oseledets IV. Constructive representation of functions in low-rank tensor formats. Constr. Appr. (1):1–18,doi:10.1007/s00365-012-9175-x. URL http://pub.inm.ras.ru/pub/inmras2010-04.pdf .12. Beylkin G, Mohlenkamp MJ. Algorithms for numerical analysis in high dimensions. SIAM J. Sci. Comput. (6):2133–2159.13. Hackbusch W, Khoromskij BN. Low-rank Kronecker-product approximation to multi-dimensional nonlocaloperators. I. Separable approximation of multi-variate functions. Computing (3-4):177–202, doi:10.1007/s00607-005-0144-0.14. Hackbusch W, Khoromskij BN. Low-rank Kronecker-product approximation to multi-dimensional nonlocaloperators. II. HKT representation of certain operators. Computing (3-4):203–225, doi:10.1007/s00607-005-0145-z.15. Khoromskij BN. Fast and accurate tensor approximation of multivariate convolution with linear scaling indimension. J. Comp. Appl. Math. (11):3122–3139, doi:10.1016/j.cam.2010.02.004.16. Khoromskij BN. O ( d log n ) –Quantics approximation of N – d tensors in high-dimensional numerical modeling. Constr. Appr. (2):257–280, doi:10.1007/s00365-011-9131-1.17. Khoromskij BN. Tensor-structured preconditioners and approximate inverse of elliptic operators in R d . Constr.Approx
SIAM J. Sci. Comput. (4):3002–3026, doi:10.1137/080730408.19. Khoromskij BN, Khoromskaia V, Chinnamsetty SR, Flad HJ. Tensor decomposition in electronic structurecalculations on 3D Cartesian grids. J. Comput. Phys. (16):5749–5762, doi:10.1016/j.jcp.2009.04.043.20. Bohr N. On the decrease of velocity of swiftly moving electrified particles in passing through matter.
PhilosophicalMagazine Series 6 (178):581–612, doi:10.1080/14786441008635432.21. Russek A, Meli J. Ionization phenomena in high-energy atomic collisions. Physica (2):222–243.22. Cocke CL. Production of highly charged low-velocity recoil ions by heavy-ion bombardment of rare-gas targets. Phys. Rev. A
Sep 1979; :749–758, doi:10.1103/PhysRevA.20.749.Copyright c (cid:13) Numer. Linear Algebra Appl. (0000)
Prepared using nlaauth.cls
DOI: 10.1002/nla
23. Shevelko VP, Litsarev MS, Tawara H. Multiple ionization of fast heavy ions by neutral atoms in the energydeposition model.
Journal of Physics B: Atomic, Molecular and Optical Physics (11):115 204.24. Shevelko VP, Kato D, Litsarev MS, Tawara H. The energy-deposition model: electron loss of heavy ions incollisions with neutral atoms at low and intermediate energies. Journal of Physics B: Atomic, Molecular and OpticalPhysics (21):215 202.25. Song MY, Litsarev MS, Shevelko VP, Tawara H, Yoon JS. Single- and multiple-electron loss cross-sections forfast heavy ions colliding with neutrals: Semi-classical calculations. Nuclear Instruments and Methods in PhysicsResearch Section B: Beam Interactions with Materials and Atoms (14):2369 – 2375.26. Shevelko VP, Litsarev MS, Song MY, Tawara H, Yoon JS. Electron loss of fast many-electron ions colliding withneutral atoms: possible scaling rules for the total cross sections.
Journal of Physics B: Atomic, Molecular andOptical Physics (6):065 202.27. Shevelko V, Litsarev M, Stohlker T, Tawara H, Tolstikhina I, Weber G. Electron loss and capture processes incollisions of heavy many-electron ions with neutral atoms. Atomic Processes in Basic and Applied Physics , SpringerSeries on Atomic, Optical, and Plasma Physics , vol. 68, Shevelko V, Tawara H (eds.). 2012; 125 – 152.28. Litsarev MS, Shevelko VP. Multiple-electron losses of highly charged ions colliding with neutral atoms.
PhysicaScripta (T156):014 037.29. Tolstikhina IY, Shevelko VP. Collision processes involving heavy many-electron ions interacting with neutral atoms.
Physics-Uspekhi (3):213.30. Tolstikhina IY, Litsarev MS, Kato D, Song MY, J-S Y, Shevelko V. Collisions of be, fe, mo and w atoms and ionswith hydrogen isotopes: electron capture and electron loss cross sections. Journal of Physics B: Atomic, Molecularand Optical Physics (3):035 206.31. Litsarev MS. The deposit computer code: Calculations of electron-loss cross-sections for complex ions collidingwith neutral atoms. Computer Physics Communications (2):432–439.32. Litsarev MS, Oseledets IV. The deposit computer code based on the low-rank approximations.
Computer PhysicsCommunications (10):28012802.33. Beylkin G, Monz´on L. On approximation of functions by exponential sums.
Appl. Comput. Harm. Anal. (1):17–48, doi:10.1016/j.acha.2005.01.003.34. Beylkin G, Monz´on L. Approximation by exponential sums revisited. Appl. Comput. Harm. Anal. (2):131–149, doi:10.1016/j.acha.2009.08.011.35. Tyrtyshnikov EE. Mosaic-skeleton approximations. Calcolo (1):47–57, doi:10.1007/BF02575706.36. Goreinov SA, Tyrtyshnikov EE, Zamarashkin NL. A theory of pseudo–skeleton approximations. Linear AlgebraAppl. :1–21, doi:10.1016/S0024-3795(96)00301-1.37. Goreinov SA, Zamarashkin NL, Tyrtyshnikov EE. Pseudo–skeleton approximations by matrices of maximumvolume.
Mathematical Notes (4):515–519, doi:10.1007/BF02358985.38. Hackbusch W, Braess D. Approximation of x by exponential sums in [1 , ∞ ] . IMA J. Numer. Anal. (4):685–697.39. Gavrilyuk IP, Hackbusch W, Khoromskij BN. Tensor-product approximation to the inverse and related operators inhigh-dimensional elliptic problems. Computing https://bitbucket.org/appl_m729/schur_cross2d .41. Stoer J, Bulirsch R.
Introduction to Numerical Analysis . 3rd ed. edn., Springer, 2002.42. Abramowitz M, Stegun IA.
Handbook of Mathematical Functions with Formulas, Graphs, and MathematicalTables . tenth dover printing, tenth gpo printing edn., Dover: New York, 1972.43. Salvat F, Martnez JD, Mayol R, Parellada J. Analytical dirac-hartree-fock-slater screening function for atoms (z=1-92).
Phys. Rev. A :467–474, doi:10.1103/PhysRevA.36.467.44. C++ code: Deposit 2014 ; URL https://bitbucket.org/appl_m729/code-deposit .45. Slater J. Quantum theory of atomic structure . International series in pure and applied physics, McGraw-Hill, NewYork, 1960.46. Shevelko VP, Vainshtein LA.
Atomic physics for hot plasmas . Institute of Physics Pub., 1993.47. Khoromskij B. Structured rank-(r1, . . . ,rd) decomposition of function-related tensors in rd.
Comput. Methods Appl.Math. :194–220.Copyright c (cid:13) Numer. Linear Algebra Appl. (0000)
Prepared using nlaauth.clsnlaauth.cls