Viscous Current-Induced Forces

Viscous current-induced forces Vladimir U. Nazarov,1Tchavdar N. Todorov,2and E. K. U. Gross1 1Fritz Haber Research Center of Molecular Dynamics, the Hebrew University of Jerusalem, Institute of Chemistry, Israel∗ 2Queen’s University Belfast, School of Mathematics and Physics, Belfast, UK We study the motion (translational, vibrational, and rotational) of a diatomic impurity immersed in an electron liquid and exposed to electronic current. An approach based on the linear response time-dependent density functional theory combined with the Ehrenfest dynamics leads to a system of linear algebraic equations, which account for the competing and counteracting effects of the currentinduced force (electron wind) and the electronic friction. These forces, by means of the dynamic exchange-correlation kernel fxc(r,r′, ω), include the electronic viscosity contribution. Starting from the ground state at the equilibrium inter-nuclear distance and applying a current pulse, we observe three phases of the motion: (I) acceleration due to the prevalence of the current-induced force, (II) stabilization upon balancing of the two forces, and (III) deceleration due to the friction after the end of the pulse. The viscous contribution to the force largely increases the acceleration (deceleration) at the first (third) phase of the process. For the aluminium HEG electron density, we find this correction to amount to up to 70% of the total electron wind and friction effects. The current densities in atomic wires can exceed those in macroscopic conductors by many orders of magnitude [1]. Current flow in a lightbulb causes heating, light emission, and electromigration [2], leading to the eventual failure of the current-carrying element. So what should we expect in the microscopic world of atomic wires? This question has prompted intense research into current-driven dynamics in nanoscale conductors for over 25 years, resulting in simulation techniques of great sophistication [3, 4]. Of central importance in these studies is the mean force exerted by the current on individual atomic nuclei. Techniques typically rely on static density functional theory (DFT) or self-consistent tight binding (TB) method for the calculation of this all-important quantity under non-equilibrium open-boundary conditions (see Ref. 5 and references therein). But it is known that static DFT/TB miss key dynamical effects, which have been shown to be of importance for related phenomena, such as electronic stopping [6], bulk impurities resistivity [7], and nanoscale conductance [8–11]. These dynamical corrections, akin to electron viscosity, have never been investigated for current-induced forces, to the best of our knowledge. The aim of this Letter is to bridge this gap. We develop a formalism based on the frequencydependent kernel of time-dependent DFT (TDDFT) and apply it to the electron-wind force on impurities in jellium. The upshot is to correct the current-induced forces and the friction effect by up to 70%. This viscous correction to the bare wind force is far from being of academic interest alone. When combined with thermal activation, this increase in the wind force can result in very significant changes to the impurity electromigration rates. A related phenomenon where these corrections become of central interest are non-conservative forces under current ∗vladimir.nazaro[email protected] and the waterwheel effect [12], an application we propose to study in the near future. We consider two classical nuclei, of atomic numbers Z1 and Z2, at positions R1and R2, respectively, immersed in an otherwise homogeneous electron gas (HEG) of density parameter rs, where ¯n−1=4 3πr3 s, ¯nbeing the HEG density (atomic units are used throughout). Equilibrium configuration.– In the ground state (GS) at the equilibrium inter-nuclear distance, the forces on each nucleus vanish Fα=−Z∇Rα Zα |r−Rα|n0(r;R1,R2)dr −∇Rα ZαZβ |Rα−Rβ|=0, α = 1,2, β = 2,1, (1) where n0(r;R1,R2) is the nuclear-position-dependent GS electronic density. Time-dependent perturbation.– We apply to the equilibrium state a weak external electric field δEext(t) = δEext(ω)e−iωt.(2) To first order in the perturbation (2), the current density induced in the system is δjind,i(r, ω) = c iω Zˆχij(r,r′, ω) [δEext,j(ω) + 2 X γ=1 (δRγ(ω)·∇′)∇′ j Zγ |r′−Rγ|#dr′, (3) where iand jare Cartesian indices, and summation over the doubly repeated index jis implied. In Eq. (3), ˆχij(r,r′, ω) is the tensorial current density response function [13], in which a parametric dependence on R1,2is implied. δRγ(ω) is the displacement, to first order in Eext(ω), of nucleus γfrom its equilibrium position. The two terms in the square brackets in Eq. (3) stand for the bare external field and that of the Coulomb charges of the displaced nuclei, respectively. arXiv:2403.07081v1 [cond-mat.mtrl-sci] 11 Mar 2024 2 Based on Eq. (3), in Appendix A we show that the electric field due to the dynamical redistribution of the electrons is δEe(r, ω)= ω2 ω2−ω2 p δEext(ω)+ 4πc ω2∇−2∇∇iZˆχij(r,r′, ω) ×"δEext,j (ω) + 2 X γ=1 (δRγ(ω)·∇′)∇′ j Zγ |r′−Rγ|#dr′, (4) where ωp=√4π¯nis the plasma frequency of the HEG. As the next step, using a relation between the tensorial ˆχof the TDCDFT and the scalar χof the TDDFT, we rewrite Eq. (4) as (see Appendix B) δEe(r, ω) = ω2 ω2−ω2 p δEext(ω) + 1 ω2−ω2 p∇Z1 |r−r′′|[χ(r′′,r′, ω)−χ(r′′,r′,0)](δEext(ω)·∇′)V0(r′)dr′dr′′ +∇Z1 |r−r′′|χ(r′′,r′, ω)"2 X γ=1 (δRγ(ω)·∇′)Zγ |r′−Rγ|#dr′dr′′, (5) where V0(r) = − 2 X γ=1 Zγ |r−Rγ|(6) is the bare static potential by the nuclei. The advantage of Eq. (5) over Eq. (4) is that the former relies on the more tractable scalar density response function χ. The last step is to evaluate the forces on each nucleus. To the first order in the perturbation, they read δFα(ω) = ZαδEe(Rα, ω) −Zα[(δRα(ω)·∇Rα)+(δRβ(ω)·∇Rβ)]∇Rα Zβ |Rα−Rβ| +Zα(δRα(ω)·∇Rα)∇RαZn0(r) |Rα−r|dr, β =α. (7) Indeed, the first term on the RHS of Eq. (7) is the field of Eq. (5) at the unshifted position of the α-th nucleus times its charge. The second line stands for the change of the bare force from another nucleus due to the displacements of both of them. The third line accounts for the change of the force in the potential of the groundstate electronic distribution due to the displacement of the nucleus α. Equations (5) - (7) constitute the main result of our theory. On the equal footing, they include current-induced forces and the electronic friction [5, 14], and the electron-electron interaction is fully accounted for through, the formally exact, density response function χ(r,r′, ω) [15]. Within the Ehrenfest dynamics, Eqs. (5) - (7) are complemented with the Newton’s equations −Mαω2δRα(ω) = δFα(ω),(8) thus closing the system of equations to be solved to find the displacements δRα(ω). The most challenging element of the described scheme is the determination of the density response function χ(r,r′, ω;R1,R2) of the system of the HEG plus two nuclei in their respective positions (here we indicate the parametric dependence on the latter explicitly). In order to avoid significant computational difficulties while still keeping the essential physics, in the implementation of our theory we resort to the Weak electrons-impurity interaction approximation.– This approximation amounts to replacing the density response function in Eq. (5) by its HEG counterpart χh(|r−r′|, ω) while, in the equilibrium condition (1), taking n0(r;R1,R2) = Zχh(|r−r′|,0)V0(r′)dr′.(9) Then the force and the potential energy at the given separation dbetween the nuclei are evaluated to (see Appendix C) F=Z1Z2 d28Zdq χh(q, 0) q3(sin qd−qd cos qd)+1,(10) U=Z1Z2 d8Zdq χh(q, 0) q3sin qd + 1,(11) where we have Fourier transformed χh(r,0) to the wavevector variable q. In Fig. 1, the force and the potential energy of Eqs. (10) and (11) are plotted versus the internuclear separation d, for HEG of density parameter rs= 2. We note that, in the weak interaction approximation, the equilibrium separation is independent on the nuclear charges, being a function, via χh, of the HEG density only. Instead of the external field δEext(ω), we introduce the current density δj(ω) in the HEG as it would be in the absence of the impurity, the two quantities being related by δj(ω) = iωω2 p 4π(ω2−ω2 p)δEext(ω),(12) 3 −1 −0.5 0 0.5 1 4 6 8 10 12 14 16 18 20 rs= 2 a.u. F,U/(Z1Z2)×103(a.u.) d(a.u.) U F FIG. 1. The potential energy U(11) and the force F(10) versus the distance dbetween two point charges Z1and Z2in the HEG of the density parameter rs= 2 a.u. The equilibrium distance between the charges (marked with an arrow) is d= 5.17 a.u. χh(q, 0) used was obtained with Eq. (17) and the MCP07 approximation for fxc(q, 0).[16] which is a consequence of the Drude formula. We, furthermore, introduce the coordinate of the center of mass (c.m.) δRc(ω) = M1δR1(ω) + M2δR1(ω) Mc , Mc=M1+M2, the relative coordinate δRr(ω) = δR2(ω)−δR1(ω),(13) and the corresponding velocities δVc(ω) and δVr(ω), which, finally, leads to a system of coupled equations of motion in terms of the latter velocities −ωδVc(ω)=−4π(Z1+Z2)ω Mcω2 p δj(ω) + 2 πMcωZdqq q4χh(q, ω)−χh(q, 0)Z2 1+Z2 2+ 2Z1Z2e−iq·dδj(ω) ¯n−δVc(ω)·q −2 πM2 cωZdqq q4χh(q, ω)−χh(q, 0)M1Z2 2−M2Z2 1+Z1Z2(M1−M2)e−iq·d(δVr(ω)·q), (14) −ωδVr(ω) = −4πω ω2 pZ2 M2−Z1 M1δj(ω) −2 πMcωZdqq q4χh(q, ω)−χh(q, 0)(δVr(ω)·q)M1Z2 2 M2 +M2Z2 1 M1−2Z1Z2eiq·d +2 πω Zdqq q4χh(q, ω)−χh(q, 0)δj(ω) ¯n−δVc(ω)·qZ2 2 M2−Z2 1 M1 +1 M2−1 M1Z1Z2eiq·d +2 πω 1 M1 +1 M2Zdqq q4Z1Z2χh(q, 0)eiq·d[δVr(ω)·q]−1 ω1 M1 +1 M2[δVr(ω)·∇d]∇d Z1Z2 d. (15) While the c.m. motion and the relative one are coupled in Eqs. (14)-(15), we show that the motion in the direction parallel to dand that in the perpendicular direction are independent (see Appendix D) This allows us to study the two geometries separately. Individual terms in Eqs. (14)-(15) can be attributed a meaning as follows. The first terms on the RHS of both equations are due to the direct field, i.e., the external field (2) screened in HEG. The second terms stand for the combined action of the electron wind and the friction, in Eq. (14), and the friction only in Eq. (15). In the limit ω→0, a factor ∂χ(q, ω)/∂ω|ω=0 appears under the integrals. The latter is familiar from theories of the stopping power for ions and the impurity resistivity [6, 7]. At finite ω, these terms still keep their meaning, although they are quantitatively changed according to the magnitude of ω. The third terms are the cross ones, accounting for the c.m. motion dependence on the relative one and vice verse. The last two terms in Eq. (15) stand for the elastic restoring force in the relative motion, as becomes clear noting that iδVr(ω)/ω =δRr(ω). This force gives rise to vibrations in the system shown in Fig 4. When the motion 4 is decomposed into the parallel and perpendicular to d components, this force survives in the parallel part only [see Eqs. (D8′) and (D10′) of Appendix D]. Equations (14)-(15) solve our problem in the case of a monochromatic field (2). Furthermore, utilizing the linearity of the present approach, by means of the forward and inverse Fourier transforms, we can construct the solution for an arbitrary current pulse ¯ δj(t). We choose a pulse of rectangular shape ¯ δj(t)=[H(T−t)−H(−t)] j0,(16) where H(t) is the Heaviside step function and Tis the duration of the pulse. We have conducted calculations using Eqs. (14)-(16) for the impurity comprised of a proton and a deuteron in a HEG of rs= 2. In Figs. 2 and 3, we present the time evolution of the velocity of the c.m. of the impurity, in the direction parallel to the axis dand perpendicular to it, respectively. In Figs. 4 - 5 the same is shown for the relative velocities. The c.m. motion (Figs. 2 and 3) 0 5 10 15 20 25 30 35 0 2 4 6 8 10 12 14 16 18 1/¯n rs= 2 a.u. T Vck/j0(a.u.) t×10−4(a.u.) MCP07 RPA FIG. 2. Velocity of the center-of-mass of the impurity comprised of a proton and deuteron in HEG of rs= 2 a.u. A rectangular current pulse is applied in the direction parallel to the impurity’s axis. At t < 0, the system is in its ground state, at the equilibrium separation between the nuclei. The solid (MCP07) and dashed (RPA) lines show results of calculations by Eqs. (14) – (15) with the use of fh xc(q, ω) of Ref. 16 and the random phase approximation [fh xc(q, ω) = 0], respectively. The pulse is schematically shown with a rectangle in short-dashed lines. are qualitatively similar in the two geometries. Upon the application of the pulse at t= 0, acceleration of the impurity as a whole under the electron wind takes place. The acceleration gradually slows down as the counteracting friction force grows with the increase of the velocity, until the c.m. velocity stabilizes at a value of j0/¯n, i.e., at the electronic drift velocity. Upon the pulse switching off at t=T, the electron wind force disappears, leading to a fast slowing down of the impurity under the friction force alone. 0 5 10 15 20 25 30 35 0 2 4 6 8 10 12 14 16 18 1/¯n rs= 2 a.u. T Vc⊥/j0(a.u.) t×10−4(a.u.) MCP07 RPA FIG. 3. Same as Fig. 2, but the pulse is applied in the direction perpendicular to the impurity axis. In Fig. 4, the relative velocity evolution in the direction parallel to the impurity axis is shown. A dominating feature here is the vibrational motion around the equilibrium separation between the two nuclei. The vibrations are, however, of attenuating amplitudes, after the kicks of the pulse switching on and off. This attenuation is due to the coupling between the relative and c.m. motions in Eqs. (14)-(15): when the velocity of the c.m. approaches saturation around the middle of the pulse (see Fig. 2), the electron wind ceases to support the vibrations, while the friction persists, leading to the decay of the oscillations. The frequency and the damping of the oscillations are worked out in Appendix E. Figure 5 also shows the evolution of the relative velocity, but in the geometry with the current perpendicular to the impurity axis. This regime corresponds to the rotation of the impurity around its c.m. Similar to vibrations, the velocity of the rotation quickly increases after the pulse switching on and off, but it falls off in the region of the stabilization (see Appendix E). In our theory, a factor of the primary role is the coupling between the motion of the c.m. of the impurity and that of the nuclei relative to each other. This is due to the mediation by the medium and, obviously, it is absent in the motion in vacuum. One of the consequences of the coupling is that, in the middle of the pulse duration, when the stabilization of the c.m. velocity is reached, the DC current does not support vibrations and rotation, because the impurity is moving as a whole with the saturation velocity j0/¯n. However, the pulse switching on or off constitutes kicks on the system, as a result all frequencies get involved, leading to the commencement/resumption of vibrations and rotation, which die out afterwards due to the friction (see Figs. 4 and 5). In the determination of χh(q, ω) in Eqs. (14) – (15), 5 −1.5 −1 −0.5 0 0.5 1 1.5 0 2 4 6 8 10 12 14 16 18 rs= 2 a.u. Vrk/j0(a.u.) t×10−4(a.u.) MCP07 RPA FIG. 4. Same as Fig. 2, but for the velocity of the nuclei relative to each other. −10 −5 0 5 10 0 2 4 6 8 10 12 14 16 18 rs= 2 a.u. T Vr⊥/j0(a.u.) t×10−4(a.u.) MCP07 RPA FIG. 5. Same as Fig. 4, but the pulse is applied in the direction perpendicular to the impurity’s axis. we rely on the linear response TDDFT equality [15] 1 χh(q, ω)=1 χh s(q, ω)−4π q2−fh xc(q, ω),(17) where χh s(q, ω) is the Kohn-Sham single particle density response function and fh xc(q, ω) is the dynamic xc kernel. While the former is known exactly and analytically by the Lindhard formula [17], the knowledge of the latter is limited to approximations. In our calculations we use the constraint-based fh xc(q, ω) of Ref. 16 (MCP07). Results are compared with those in the random phase approximation (RPA) [i.e., fh xc(q, ω) = 0]. From Figs. 2 – 5 we conclude that the electronic viscosity, i.e., the frequency-dependence of fxc, plays an important role in the motion. In the c.m. motion (Figs. 2 and 3), it leads to an increase in the acceleration due the electron wind at the first phase of the motion, and to a faster deceleration upon the end of the pulse, the latter being due to friction. The maximal effect of the viscosity can be observed in the case of the rotational motion (Fig. 5), where, in the stabilization region (phase II), it leads to a decrease of the velocity by up to 70%. In conclusion, we have studied the motion of a diatomic impurity in an electron liquid under the action of an electronic current. A consistent linear response TDDFT approach combined with Ehrenfest dynamics has been utilized. In contrast to applications of this method to the related problems of the electronic stopping power and the impurity resistivity in metals, we find the role of the viscosity of the electron liquid to be of major importance. In particular, in the case of the rotational motion, the viscosity contribution comprises up to 70% of the electron wind and the electronic friction effects. Fundamental implications for the electromigration theory are envisaged. ACKNOWLEDGMENTS This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and information programme (Grant Agreement No. ERC-2017-AdG-788890). [1] T. H. Todorov, J. Hoekstra, and A. P. Sutton, Currentinduced forces in atomic-scale conductors, Philosophical Magazine B 80, 421 (2000). [2] R. S. Sorbello, Theory of electromigration, in Solid State Physics, Solid State Physics, Vol. 51, edited by H. Ehrenreich and F. Spaepen (Academic Press, 1998) pp. 159– 231. [3] J.-T. L¨u, S. Leitherer, N. R. Papior, and M. Brandbyge, Ab initio current-induced molecular dynamics, Phys. Rev. B 101, 201406 (2020). [4] M. Stamenova, P. Stamenov, and T. Todorov, Phonon and magnon jets above the critical current in nanowires with planar domain walls, Phys. Rev. Lett. 131, 206302 (2023). [5] W. Dou and J. E. Subotnik, Perspective: How to understand electronic friction, The Journal of Chemical Physics 148, 230901 (2018). [6] V. U. Nazarov, J. M. Pitarke, Y. Takada, G. Vignale, and Y.-C. Chang, Including nonlocality in the exchangecorrelation kernel from time-dependent current density functional theory: Application to the stopping power of electron liquids, Phys. Rev. B 76, 205103 (2007). [7] V. U. Nazarov, G. Vignale, and Y.-C. Chang, Dynamical many-body corrections to the residual resistivity of metals, Phys. Rev. B 89, 241108(R) (2014). 6 [8] M. Koentopp, K. Burke, and F. Evers, Zero-bias molecular electronics: Exchange-correlation corrections to landauer’s formula, Phys. Rev. B 73, 121403(R) (2006). [9] N. Sai, M. Zwolak, G. Vignale, and M. Di Ventra, Dynamical corrections to the dft-lda electron conductance in nanoscale systems, Phys. Rev. Lett. 94, 186810 (2005). [10] J. Jung, P. Bokes, and R. W. Godby, Phys. Rev. Lett. 98, 259701 (2007). [11] N. Sai, M. Zwolak, G. Vignale, and M. Di Ventra, Sai et al. reply:, Phys. Rev. Lett. 98, 259702 (2007). [12] D. Dundas, E. J. McEniry, and T. N. Todorov, Currentdriven atomic waterwheels, Nature Nanotechnology 4, 99 (2009). [13] G. Vignale and W. Kohn, Current-density functional theory of linear response to time-dependent electromagnetic fields, in Electronic Density Functional Theory: Recent Progress and New Directions, edited by J. F. Dobson, G. Vignale, and M. P. Das (Plenum, New-York, 1998) pp. 199–216. [14] W. Dou, G. Miao, and J. E. Subotnik, Born-oppenheimer dynamics, electronic friction, and the inclusion of electron-electron interactions, Phys. Rev. Lett. 119, 046001 (2017). [15] E. K. U. Gross and W. Kohn, Local density-functional theory of frequency-dependent linear response, Phys. Rev. Lett. 55, 2850 (1985). [16] A. Ruzsinszky, N. K. Nepal, J. M. Pitarke, and J. P. Perdew, Constraint-based wave vector and frequency dependent exchange-correlation kernel of the uniform electron gas, Phys. Rev. B 101, 245135 (2020). [17] G. F. Giuliani and G. Vignale, Quantum Theory of the Electron Liquid (Cambridge University Press, Cambridge, 2005). Appendix A: Derivation of Eq. (4) By Maxwell’s equations ∇×δE(r, ω) = iω cδH(r, ω),(A1) ∇×δH(r, ω) = −iω cδE(r, ω)−4π cδj(r, ω),(A2) we can write δE(r, ω) = −4π iω δj(r, ω) + c2 ω2∇×[∇×δE(r, ω)].(A3) Since, for a uniform field, Eext(ω) = −4π iω jext(ω),(A4) we can rewrite Eq. (A3) as δE(r, ω) = δEext(ω)−4π iω δjind(r, ω) + c2 ω2∇×[∇×δE(r, ω)].(A5) where jind(r, ω) = j(r, ω)−jext(ω).(A6) Using a vector analysis identity, Eq. A5 can be written as δE(r, ω) = δEext(ω)−4π iω δjind(r, ω) + c2 ω2∇[∇·δE(r, ω)] −c2 ω2∇2δE(r, ω).(A7) or, with account of Eq. (A2), δE(r, ω) = δEext(ω)−4π iω δjind(r, ω)−4πc2 iω3∇[∇·δjind(r, ω)] −c2 ω2∇2δE(r, ω).(A8) Adding and subtracting −4π iω δjind(ω), we have δE(r, ω) + c2 ω2∇2δE(r, ω) = δEext(ω)−4π iω δjind(ω)−4π iω [δjind(r, ω)−δjind(ω)] −4πc2 iω3∇[∇·δjind(r, ω)].(A9) and then δE(r, ω) = δEext(ω)−4π iω δjind(ω)−(1 + c2 ω2∇2)−14π iω [δjind(r, ω)−δjind(ω)] + 4πc2 iω3∇[∇·δjind(r, ω)],(A10) 7 where we have used the fact that, on a uniform field, (1 + c2 ω2∇2)−1acts as the identity operator. Then, in the limit c→ ∞, δE(r, ω) = δEext(ω)−4π iω δjind(ω)−4π iω ∇−2∇[∇·δjind(r, ω)].(A11) Since quantities averaged over the whole space are those of the HEG, we have δE(r, ω) = δEext(ω) + δEind(ω)−4π iω ∇−2∇[∇·δjind(r, ω)] =δE(ω)−4π iω ∇−2∇[∇·δjind(r, ω)] = 1 1−ω2 p ω2 δEext(ω)−4π iω ∇−2∇[∇·δjind(r, ω)].(A12) Equation (A12) together with Eq. (3) lead us to Eq. (4). Appendix B: Derivation of Eq. (5) Using the relation [7] ω2−ω2 pZˆχij(r,r′, ω)dr′=Zˆχik(r,r′, ω)∇′ k∇′ jV0(r′)dr′+ω2 cn0(r)δij,(B1) V0(r) given by Eq. (6), we have from Eq. (4) δE(r, ω) = 1 1−ω2 p ω2 δEext(ω) +4πc ω2(ω2−ω2 p)∇−2∇Z∇iˆχik(r,r′, ω)∇′ k∇′ jV0(r′)dr′+ω2 c∇in0(r)δijδEext,j (ω) +4πc ω2∇−2∇Z∇iˆχij(r,r′, ω)∇′ j"2 X γ=1 (δRγ(ω)·∇′)Zγ |r′−Rγ|#dr′. (B2) With the use of the relation between the tensorial ˆχof the TDCDFT and the scalar one χof the TDDFT [7] χ(r,r′, ω) = −c ω2∇iˆχij(r,r′, ω)∇′ j,(B3) Eq. (B2) can be rewritten as δE(r, ω) = 1 1−ω2 p ω2 δEext(ω)−4π ω2−ω2 p∇∇−2Zχ(r,r′, ω)(δEext(ω)·∇′)V0(r′)dr′−(δEext(ω)·∇)n0(r) −4π∇∇−2Zχ(r,r′, ω)"2 X γ=1 (δRγ(ω)·∇′)Zγ |r′−Rγ|#dr′, (B4) or, with the use of the static sum-rule, Zχ(r,r′,0)∇′V0(r′)dr′=∇n0(r),(B5) Eq. (B4) can be written as δE(r, ω) = 1 1−ω2 p ω2 δEext(ω)−4π ω2−ω2 p∇∇−2Z[χ(r,r′, ω)−χ(r,r′,0)](δEext(ω)·∇′)V0(r′)dr′ −4π∇∇−2Zχ(r,r′, ω)"2 X γ=1 (δRγ(ω)·∇′)Zγ |r′−Rγ|#dr′, (B6) which proves Eq. (5). 8 Appendix C: Derivation of Eqs. (10) and (11) Rewriting the force of Eq. (1) with the substitution of n0(r) in the weak interaction regime (9), we obtain Fα=1 (2π)3Z∇r Zα |r−Rα|χh(q, 0)ei(r−r′)·qZα |r′−Rα|+Zβ |r′−Rβ|drdr′dq−∇Rα ZαZβ |Rα−Rβ|, α = 1,2, β = 2,1, (C1) or, after the integration by parts, Fα=−1 (2π)3ZiqZα |r−Rα|χh(q, 0)ei(r−r′)·qZα |r′−Rα|+Zβ |r′−Rβ|drdr′dq−∇Rα ZαZβ |Rα−Rβ|, α = 1,2, β = 2,1. (C2) or Fα=−1 (2π)3ZiqZα |r|χh(q, 0)ei(r−r′)·qZα |r′|drdr′dq −1 (2π)3ZiqZα |r−Rα|χh(q, 0)ei(r−r′)·qZβ |r′−Rβ|drdr′dq−∇Rα ZαZβ |Rα−Rβ|, α = 1,2, β = 2,1, (C3) where the integration variables substitutions has been performed in the first term. Obviously, the latter evaluates to zero, and we have Fα=−1 (2π)3ZiqZα |r−Rα|χh(q, 0)ei(r−r′)·qZβ |r′−Rβ|drdr′dq−∇Rα ZαZβ |Rα−Rβ|, α = 1,2, β = 2,1,(C4) In the last equation, integrations over rand r′are carried out explicitly, producing Fα=−(4π)2ZαZβ (2π)3Ziq q4χh(q, 0)e−id·qdq−∇d ZαZβ d=0, α = 1,2, β = 2,1.(C5) Equation (10) is retrieved by the integration in Eq. (C5) over angular variables of the vector q. Finally, Eq. (11) is obtained by the formula U(d) = ∞ Z d F(d)dd. (C6) Appendix D: Proof of the independence of the motions parallel and perpendicular to the impurity’s axis within the weak interaction approximation For an arbitrary vector b=b∥+b⊥, where b∥and b⊥are vectors parallel and perpendicular to d, respectively, equalities hold Zq(b·q)dΩq=4π 3q2b,(D1) Zq(b∥·q)eid·qdΩq=4π d3qb∥[2qd cos qd + (q2d2−2) sin qd],(D2) Zq(b⊥·q)eid·qdΩq=4π d3qb⊥(sin qd −qd cos qd),(D3) where the integrations are taken over the full solid angle of the vector q. In Eqs. (14) and (15), let us decompose δj(ω) = δj∥(ω) + δj⊥(ω),(D4) δVc(ω) = δVc∥(ω) + δVc⊥(ω),(D5) δVr(ω) = δVr∥(ω) + δVc⊥(ω).(D6) 9 Then, with the use of Eqs. (D1)-(D3), −ωMcδVc∥(ω)=−4π(Z1+Z2)ω ω2 p δj∥(ω) +8 ωZdq χh(q, ω)−χh(q, 0)1 3(Z2 1+Z2 2) + 2Z1Z2 q3d3[2qd cos qd + (q2d2−2) sin qd]"δj∥(ω) ¯n−δVc∥(ω)# −8 McωZdq χh(q, ω)−χh(q, 0)1 3(M1Z2 2−M2Z2 1)+(M1−M2)Z1Z2 q3d3[2qd cos qd + (q2d2−2) sin qd]δVr∥(ω). (D7) −ωδVr∥(ω) = −4πω ω2 pZ2 M2−Z1 M1δj∥(ω) +8 ωZdq χh(q, ω)−χh(q, 0)1 3Z2 2 M2−Z2 1 M1+1 M2−1 M1Z1Z2 q3d3[2qd cos qd + (q2d2−2) sin qd]"δj∥(ω) ¯n−δVc∥(ω)# −8 McωZdq χh(q, ω)−χh(q, 0)1 3M1Z2 2 M2 +M2Z2 1 M1−2Z1Z2 q3d3[2qd cos qd + (q2d2−2) sin qd]δVr∥(ω) +8 ω1 M1 +1 M2Zdqχh(q, 0)Z1Z2 q3d3[2qd cos qd + (q2d2−2) sin qd]δVr∥(ω)−2Z1Z2 ωd31 M1 +1 M2δVr∥(ω). (D8) −ωMcδVc⊥(ω)=−4π(Z1+Z2)ω ω2 p δj⊥(ω)+ 8 ωZdqχh(q, ω)−χh(q, 0)1 3(Z2 1+Z2 2)+ 2Z1Z2 q3d3(sin qd−qd cos qd)δj⊥(ω) ¯n−δVc⊥(ω) −8 McωZdq χh(q, ω)−χh(q, 0)1 3(M1Z2 2−M2Z2 1)+(M1−M2)Z1Z2 q3d3(sin qd −qd cos qd)δVr⊥(ω), (D9) −ωδVr⊥(ω) = −4πω ω2 pZ2 M2−Z1 M1δj⊥(ω) +8 ωZdq χh(q, ω)−χh(q, 0)1 3Z2 2 M2−Z2 1 M1+1 M2−1 M1Z1Z2 q3d3(sin qd −qd cos qd)δj⊥(ω) ¯n−δVc⊥(ω) −8 McωZdq χh(q, ω)−χh(q, 0)1 3M1Z2 2 M2 +M2Z2 1 M1−2Z1Z2 q3d3(sin qd −qd cos qd)δVr⊥(ω) +8 ω1 M1 +1 M2ZdqZ1Z2χh(q, 0) 1 q3d3(sin qd −qd cos qd)δVr⊥(ω) + Z1Z2 ωd31 M1 +1 M2δVr⊥(ω). (D10) Finally, we note that, because dis the equilibrium inter-nuclear distance nullifying the force of Eq. (10), Eqs. (D8) and (D10) can be rewritten as −ωδVr∥(ω) = −4πω ω2 pZ2 M2−Z1 M1δj∥(ω) +8 ωZdq χh(q, ω)−χh(q, 0)1 3Z2 2 M2−Z2 1 M1+1 M2−1 M1Z1Z2 q3d3[2qd cos qd + (q2d2−2) sin qd]"δj∥(ω) ¯n−δVc∥(ω)# −8 McωZdq χh(q, ω)−χh(q, 0)1 3M1Z2 2 M2 +M2Z2 1 M1−2Z1Z2 q3d3[2qd cos qd + (q2d2−2) sin qd]δVr∥(ω) +8 ω1 M1 +1 M2Zdqχh(q, 0) sin qdZ1Z2 qd δVr∥(ω). (D8′)