A two-channel approach to the average retarding force of metals for slow singly ionized projectiles

Based on the fundamental momentum-transfer theorem [Phys. Rev. Lett. 15, 11 (1965)] a novel contribution to the retarding force of metallic systems for slow intruders is derived. This contribution is associated with sudden charge-changing cycles during the path of projectiles. The sum of the novel and the well-known conventional contributions, both expressed in terms of scattering phase shifts, are used to discuss experimental data obtained for different targets. It is found that our two-channel modeling, with two nonlinear channels, improves the agreement between several data and theory and thus, as predictive modeling, can contribute to the desired convergence between experimental and theoretical attempts on the retarding force.


I. INTRODUCTION AND MOTIVATION
According to a basic book on quantum mechanics by Landau [1] one of the most important quantities in the interaction of heavy charged projectiles with fixed atoms is the average energy loss. This time-independent quantity, a kind of deposited energy, is an observable and due to conservation laws its measurement is feasible in experiments. Thus in this subfield of nature (physics, human therapy) the real challenge resides in the convergence of measurements and theories. Their interplay, a continuos one over a century, fertilizes the developments on both side of approaches which can result in a transferable knowledge [2].
The present contribution is dedicated to the case where singly ionized projectiles interact with constituents of metallic targets. The main challenge addressed here is to find a reasonable combination of the quantum statistical and atomistic aspects of the energy loss process in real targets. As motivation, on which our new attempt is partially based, we start with an established result. The well-known, conventional form [3][4][5] for the stopping power (written in Hartree atomic units, where e 2 = m e = = 1) of a homogeneous degenerate electron gas (characterized by Fermi velocity v F ) for heavy intruders is given by where v and v e ∈ [0, v F ] are the projectile and system-electron velocities. The relative velocity v r is expressed by these kinematical variables as v 2 r = (v 2 + v 2 e − 2v e v x). Clearly, in an interpretation based on independent-electron scattering off a heavy projectile moving with constant velocity v, the remaining task resides in a common two-body interaction V (r), in order to perform the statistical averaging over a Fermi-Dirac distribution with In this scattering interpretation the analysis is based on the concept of asymptotic states in the infinite past and future, i.e., involving large time differences. Sudden processes in time, like a local charge-change in metals, requires a refined approach on associated transition amplitudes in time-dependent perturbation theory. In kick-like sudden [1] processes one may use predetermined states as a complete set to treat the matrix elements in strong (but Such an implementation rests on averaging of quantum mechanical time-dependent energy differences over certain time scales in order to define a force-like quantity as stationary observable in stopping [6,7]. In these large-scale simulations specification of the initial conditions is required to real-time propagation. For instance, in [7] two options were considered for helium in aluminum target. In the first one, the screened atom was included in the determination of the static initial state. In the second one, the initial condition was set up by adding an α-particle and thus producing a sudden change in the external potential.
In both cases the authors control only the initial state and not the subsequent dynamics which is given by the time-dependent single-particle equations within TDDFT. Therefore, a smoothed evolving picture, without fast local charge-changing processes, is employed.
First, as concretization of motivation, we integrate Eq.(1) by using models for the momentum transfer cross section in order to get useful information to phenomenology made in Section II after Eqs. (6). Namely, we take the form of σ tr (v r ) = 4πA α /v α r , in which α = 2 and α = 4. By straightforward quadrature we obtain [5] from Eq.(1) for these models The above model results are, of course, in agreement with expected limits (dE/dx) = at v → 0 and v → ∞, respectively. Earlier, careful theoretical analysis [8] performed within an adiabatic framework on velocity-dependence states that the next term beyond the v-proportional one is at least second order in velocity.
Our closed expressions for v ≤ v F are in harmony with this important statement. Furthermore, a certain weighted combination of our two expressions at v ≤ v F would result in an almost perfect v-proportionality. That, at this point simple mathematical, observation will become a more transparent and physical one in Section II, where we extend the theory on the average retarding force beyond the common fixed-potential approximation by considering physically reasonable force matrix elements as independent channel contributions.
The rest of this paper is organized as follows. Section II is devoted to the theory and the discussion of the results obtained. The last Section contains a short summary and few dedicated comments. As above, we use atomic units throughout this work.

II. RESULTS AND DISCUSSION
We begin this section by outlining few important elements of stationary scattering theory.
According to basic rules of quantum mechanics on expectation values of operators, one should consider the force matrix element [9] between orthonormal components of a scattering state to characterize the associated momentum transfer. Applying this quantum mechanical theorem, where σ tr (v r ) ∝ ∞ l=0 (l + 1)[I 1 (l, v r )] 2 , one has [10][11][12] for the matrix elements We stress that this remarkable identity rests on those states which refer to the scattering Schrödinger wave equation with v 2 r /2 energy and V (r) external field. However, with partial waves based on V (r), but with a net Coulomb field ∆V c (r) = −q/r in space of V (r) we get and with unperturbed (u) partial wave components the corresponding result becomes Here we used the spherical Bessel functions of the first kind, i.e., the components of an unperturbed plane-wave (momentum) state, instead of self-consistent radial functions. These forms in Eqs. (4)(5) are based on the fact that in cases with abrupt perturbations the original stationary system has no time [1] to relax to the stationary state of a new Hamiltonian.
We will consider these amplitudes as the proper ones when there is a sudden change in the self-consistent V (r), as in the case of charge-changing (q = 1) processes generated by the binary interaction with fixed lattice ions. This charge-change results in an excess bare field ∆V c (r) = −1/r. The square of [I 2 (l, v r ) − I (u) 2 (l, v r )] 2 can characterize, in a quantum mechanical interpretation, an extra (kick-like) momentum transfer due to the sudden change in the external field. That square is, in fact, a regularized transition probability. Such a regularization is needed since both I 2 and I (u) 2 would give divergent results after l-summation. This regularized channel gives (at q = 0) a novel form for the associated cross section to which a simple trigonometrical identity (1 − cos α) 2 = 4 sin 4 (α/2) is employed.
Before our quantitative analysis, we continue with phenomenology. There are important differences between Eq.(6) and Eq. (2), i.e., between σ (2) tr (v r ) and the conventional one given by Eq.(2) and denoted from here by σ (1) tr (v r ). The kinematical prefactors show that the new term vanishes faster at large scattering wave number v r . However, at v r ≃ 1 values, which represents at small intruder velocity the range of the Fermi velocity of metals, and at δ ≃ π for a dominating phase shift, the new term can become the dominating one.
with the integrated characteristics found with separated model cross sections in the Introduction, signals that a velocity-proportionality in the stopping power holds, practically upto v ≃ v F from below. Now, we turn to the quantitative part of this Section. We will determine numerically the two quantities, denoted by Q (1) (v F ) and Q (2) (v F ), by which the low-velocity stopping power of metals (a system of an electron gas and lattice ions) takes a friction-like form where the two coefficients (when q = 0) are given by the following expressions Our summation of two channel cross sections in Eq.(7) resembles, mathematically, to the well-known [13] rule in potential scattering with a spin-orbit interaction term where we sum the direct (non-spin-flip) and spin-flip partial differential cross sections for electron scattering for any spin orientation before scattering. There, the integrated cross sections, needed to observables, are obtained by integrating over all scattering angles. Remarkably, the spin-flip part depends on an amplitude difference, similarly to our regularized difference.
. Phase shifts, based on the orbital version of DFT [19][20][21], are employed. See the text for further details. We stress at this point that we employ to summation in Eq. (7) an a priori equal-weight assumption. In reality, i.e., at channeling-like conditions in metals, the impact parameterdependent closest approach of intruders and lattice ions [7,[14][15][16] may influence that assumption. In more simple terms, our present weighting would refer to random-collision situations. Nonequal weighting might be based on certain probabilistic inputs [17] to sum two nonlinear channel. But, such inputs need, in our modeling, an additional justification, since one can not apply stationary linear-response ideas to a sudden effect.  [19][20][21]. Both Q (1) and Q (2) are oscillating functions, but the oscillation in Q (2) (r s , Z 1 ) is a moderate one. This term in the sum [Q (1) + Q (2) ] can make an important modulation in the sharp Z 1 -oscillation of Q (1) , especially around the minima of this latter. For Z 1 = 1, we take our values for Q (1) obtained within the explicit version [5] of DFT. There a single Euler equation is solved in an iterative self-consistent way. That calculation does not consider a doubly-populated weakly bound (extended) state around an embedded proton in an electron gas, in harmony with experimental facts, obtained by positive muons, on the nonexistence of muonium in metals.
Despite this, there is a perfect agreement [5] with Q (1) results obtained from the implicit, orbital-based, DFT. This agreement signals that it is the short-range part of proton-screening which needs a nonlinear treatment. In simple terms, that range is the most important one to determine the first few phase shifts. Our Q (2) = 0 values for Z 1 = 1 are in accord with a screened-proton picture without bound state, where there is no charge-changing cycle, thus q = 0 during the motion of the projectile. Since the experimental data, obtained with low velocity proton projectiles for Al and Ni targets, are in very reasonable harmony [15,16,18] with nonlinear theory [19][20][21] based solely on Q (1) (v F ), we have a transferable knowledge in this case. A desired convergence between the two sides of understanding is achieved.
For all other Z 1 ≥ 2 we consider, for velocities v ≤ v F , the q = 1 value as the most plausible one. This conservative value seems to be a realistic one with singly ionized intruders.
Higher q values might have relevance when there is a large electronic overlap between clouds of colliding atoms. We believe that such a partial channel with q > 1 would need more energetic, head-on-like collisions. Theoretically, it would be interesting to model the transition from our discrete-q modeling of charge-changing cycles to the pioneering [22] quasiclassical work where the retarding force (the observable) is related to an electron-density-flux constructed from the statistical Thomas-Fermi theory of atoms. With a transition-study one might arrive at a deeper understanding of an integrated (classical trajectory Monte Carlo) approach [23] for energy loss and capture processes.
For the channeling simulation our Q (1) also gives a very reasonable agreement with [6,7].
We stress that [24] uses α-particle (He ++ ) as projectile and an atom-centered optimized Gaussian basis set to model the energy transfer. The observed agreements (see, the discussion at Fig. 2) are especially remarkable in the light of careful experiment [25] performed on electron emission from aluminum. There a perfect linearity in the velocity of helium ions, with v ≤ 0.6, was obtained, and thus a quantitative agreement with [6] was concluded.
Notice, in the spirit of discussion made already in [7] for proton and helium intruders, tr (v F )]. We will return to this experiment, at the discussion of Fig. 3 which is devoted to an important comparison for Z 1 > 2.
Related to our prediction with [Q (1) + Q (2) ] > Q (1) values for the average retarding force in metals, we turn to a brief discussion of data [14] obtained for an other free-electron-like material, Mg. For this metal v F ≃ 0.7 and the experiment with He + and proton projectiles was performed for v ≥ v F . It was found that the ratio (R) of stopping powers with these intruders becomes about two, in contrast to a ratio of about unity which is based on Q (1) values of self-consistent DFT. Our novel approach would result in a ratio (R > 2) which is not in contradiction with experimental suggestion. As support, we note that in [26], i.e., in surface experiment, the helium per proton stopping ratio was found to be always higher than unity, even for r s (z) > 3. Clearly, the desired convergence of theories and experiments requires further studies for Mg (r s ≃ 2.7) and, say, for Ca (r s ≃ 3) as well, within large-scale TDDFT simulations with proton and helium intruders at v ≤ v F . In the field of high-energy-density plasmas, our q-mediated enhancement in stopping power may contribute to the proper determination of the ignition threshold [28] in a deuteriumtritium-alpha energy deposition process. There, via a plausible postulation, the theoretical underprediction of stopping data was associated [28] with ion-ion nuclear scattering. Now, we illustrate our novel results by three Figures. In Fig.1, for r s = 1.5 of the Wigner-Seitz parameter, we plot the dimensionless ratios of , ratios of nonlinear quantities. One might consider [18,29] these ratios as a kind of effective charge. This Figure   reflects, in a highly phenomenological manner, that the so-called Z 1 oscillations may get important modulations especially around the minima of the conventional R 1 ratio. Notice that the such-defined ratios are square roots of physical magnitudes. This mathematical operation has a smoothing character (c.f., Fig. 3) with renormalized oscillating functions. corrections [30], is used. In such an approximate theory, the stopping is proportional to Z 2 1 . The corresponding form, employed in a cornerstone work [20] as well, is given by where a = πv F is an abbreviation. Its value is about a ≃ 2.83 for Al. This expression is a particular realization of the modeling made in the Introduction with σ tr (v r ) ∝ (1/v 4 r ). Our present results (at Z 1 = 2) correspond to the green solid curve [Q (1) + Q (2) ], and the green dashed curve [Q (1) ]. The red curve with dots, and the black curve with dots, are taken from Fig. 5b of [7]. They refer to off-channeling and channeling conditions, respectively.
Notice that an earlier TDDFT result of [6] (not shown here) agrees precisely with the black curve. There is a fortuitous similarity between the RPA result for the homogeneous electron gas and results plotted via curves by solid green and dash-dotted red with dots. Neither the screening-treatment nor the scattering-description of RPA is correct to a nonlinear situation.
Experimental [15] data, used already in TDDFT to comparison [7], are plotted here by a dashed magenta curve with filled triangles. This curve signals a two-slope behavior with linearities in projectile velocity. Remarkably, a quite similar, i.e., two-slope, behavior was found in [16] for Ni (r s ≃ 1.8) with singly-ionized He + intruder. There, a comparison with TDDFT results [31] is made, by using 1.15 as multiplying factor for the simulation-based results. As we already discussed above, we can image such a two-slope behavior within the present theoretical framework with certain, presumably closest-approach-dependent [32], finer tuning of our two nonlinear channels. A complete convergence is still not achieved.
The two black squares, for α-projectile, are taken from Fig. 5b of [24], for our velocity range. They are based on TDDFT with an optimized, localized atomistic, Gaussian basis set. We speculate that, for extended systems with slow ions, the screening action of the metallic electron gas needs a further consideration. Moreover, singly ionized He + intruder, instead of He ++ , might be more close to the experimental situation. Their systematic and statistical errors are analyzed in the experimental work [15]. Finally, the green dotted curve refers to Eq.(10) with Z 1 = 2. See the text for further details. . In Fig. 3 we plot the observable quantities, (1/v)(dE/dx), as a function of Z 1 . The experimental data (black dots and triangles) are used earlier [20] to a comparison with Q (1) , which is denoted here by a dashed curve. It was stated in this pioneering work that there is a substantial disagreement with data in magnitude, particularly around the minimum. Our new result, [Q (1) + Q (2) ], is denoted by a solid curve. Notice that data symbols, without error bars, refer to v = 0.411 (dots) and v = 0.826 (triangles). The target is the frequently used prototype of free-electron metal, aluminum. By inspection, one can observe an essential improvement in agreement between data and the novel approach. Here we return to experiment in [26], i.e., to the above-mentioned surface experiment. There, although with somewhat smaller deviations from the conventional Q (1) [v F (z)]-type scaling, also a systematic upward enhancement in stopping power was established. In the present two-channel modeling, such an enhancement can be associated with a Q (2) -proportional contribution. in a pioneering work leaded by Echenique [20]. They are denoted by solid circles and triangles. .

III. SUMMARY AND COMMENTS
In this theoretical paper we have investigated the problem of average retarding force of metallic targets for slow projectiles. Beyond the well-known contribution, established for electron-slow-intruder scattering in a degenerate electron gas, a novel contribution is derived which is associated with charge-change cycles due to lattice ions. Our main closed result is given by Eq. (6), which is, in the terminology of this sub-field of physics, a nonlinear form, similarly to the more conventional one in Eq. (2). These forms are implemented here by standard phase shifts obtained by applying the orbital-based DFT to screening in an electron gas. For helium projectiles, we made comparisons with selected results of experiments [15] and large-scale simulations [6,7] in TDDFT for a free-electron-like metal, Al. With our novel contribution to the retarding force, the agreement with these results is improved.
A challenging problem in ratio-data-interpretation [14] for Mg is discussed as well. Our two-channel-based result is in reasonable agreement with data at around v = v F ≃ 0.7.
The conventional theoretical estimation is not in agreement with experiment. A detailed comparison with data for Z 1 ∈ [1,18] in Al is made, and improved agreement is found.
Based on these agreements, we suggest further efforts within TDDFT along these lines.
Notice that a recent adiabatic modeling [29], motivated by experiment in [18], also results in remarkable deviations from a simple modeling with Q (1) (v F ). There the density inhomogeneity was considered, via lattice-atom-volume averaging of Q (1) [v F (r)], as a modulating effect. Such an averaging was applied successfully [33] for stopping of swift Z 1 = ±1 in order to discuss the charge-sign effect in Si. In the theoretical modeling [29] a remarkable similarity with experimental effective charges (defined at Fig. 1)  i.e., Q (1) -dependent, theoretical estimation. Both seem to be, a priori, relevant in reality.
Their proper weights and interplay need future investigations. Cases with self-irradiated condition [34,35] could be especially important in this (q = 0) respect. For instance, Ni ions in Ni target [35]. In such a symmetric case we can image (for a metal) even q = 2 to our Q (2) channel. For r s ≃ 2, Ni ion with its Z 1 = 28 represents the second minimum in Z 1 -oscillation [19,21]. In our modeling we get [Q (1) + Q (2) ] ≃ [0.28 + (q 2 × 0.72)].
At q = 2 and v = 1, one arrives at [(Q (1) + Q (2) )/Q (1) ] ≃ 11, thus the corresponding stopping power would change steeply to about (dE/dx) ≃ 160 eV /A. For transition metals, that show a high electronic stopping power [35], the spin-flip process needs a thorough investigation. The electron spin direction is no longer conserved during electron-atom collision.
One has to consider the total angular momentum j = (l + s) operator in order to construct a complete set of spin-angle functions which are needed to expansions. We left this exciting sub-problem in stopping theory with a new (spin) degree of freedom to future studies.
Notice that at high ion velocity, our new term would scale as (q/v) 2 with respect to the conventional, i.e., Bethe-like [1], leading one [7,24]. There, a term with [q(v)/v] 2 can give a slowly vanishing enhancement. Thus the high-velocity limit, first of all under self-irradiated condition [35], also requires further investigations. The precise relevance of permutationbased, similarity-aided level crossing [22,36] behind higher q(v)-values seems to be an other interesting sub-problem in stopping theory. The Bethe-limit, especially for metals with their dense electron gas, is not a simple cumulative sum of isolated atomic contributions [37].
Based on the established capability of our new modeling for metals, we believe that the two-channel approach developed here can find application in other important fields as well. For instance, in the friction problem of diatomic molecules during their dissociative adsorption on metallic surfaces. There, based on an empirically motivated local-densityfriction approximation, a local Q (1) [v F (r)] is employed [38]. We argue here that transient electronic processes, due to dissociation in an electron gas, could be related to Q (2) in Eq. (6). For instance, the case of N, with its Z 1 = 11, might be a good candidate, as our Table I suggests. We stress, however, that at high target-temperatures, the coupling to phonon modes, i.e., to quanta of lattice vibrations, can open a new [39] channel to inelastic processes, beyond the friction-like channel discussed in our present study for cold metals. Still, as Figure 1 of [39] signals, the proper magnitude of this latter channel could be important. Indeed, the so-far neglected [39] charge-transfer-type [related to Q (2) (v F )] processes, especially with highly reactive molecules, may have impact on conclusions.
We close with few general comments. The wave functions of the conventional orbitalbased DFT for embedded Z 1 are used [19][20][21] here to calculate the induced electron density.
That is the basic variable of the underlying variational theory. The phase shifts are, therefore, auxiliary quantities [19]. Their sums over angular momentum quantum numbers always satisfy the associated neutrality condition of a self-consistent orbital-based approximation, i.e., the Friedel sum rule and the Levinson theorem [13] for local interactions. Since these are satisfied by construction for any form of a local many-body term in the Schrödinger-like equations, the physical quality of DFT results needs further, i.e., energetic, justifications.
But, in accord with closely related statements [19,21], the highly improved quantitative agreements with experimental facts justify, a posteriori, our phase-shift-based two-channel modeling with a novel term for the retarding force. Generally, and in accord with Landau basic attempt [36] for fermi liquids, a modeling is good if it contains few adjustable elements, agrees with several observations, and makes controllable predictions. We stress, finally, that the truly exciting theoretical problem of interparticle-interaction, i.e., correlated motion of electrons, is considered in stopping calculation only at the mean-field level. However, at least for a prototype two-particle correlated model system, recent exact result [40] for the energy shift in time-dependent (passing) perturbations indicates that proper independent modes, rather than effective single-particle states, could pave the path to future developments.