High frequency limit of spectroscopy

We consider an arbitrary quantum mechanical system, initially in its ground-state, exposed to a time-dependent electromagnetic pulse with a carrier frequency $\omega_0$ and a slowly varying envelope of finite duration. By working out a solution to the time-dependent Schr\"odinger equation in the high-$\omega_0$ limit, we find that, to the leading order in $\omega_0^{-1}$, a perfect self-cancellation of the system's linear response occurs as the pulse switches off. Surprisingly, the system's observables are, nonetheless, describable in terms of a combination of its linear density response function and nonlinear functions of the electric field. An analysis of jellium slab and jellium sphere models reveals a very high surface sensitivity of the considered setup, producing a richer excitation spectrum than accessible within the conventional linear response regime. On this basis, we propose a new spectroscopic technique, which we provisionally name the Nonlinear High-Frequency Pulsed Spectroscopy (NLHFPS). Combining the advantages of the extraordinary surface sensitivity, the absence of constraints by the traditional dipole selection rules, and the clarity of theoretical interpretation utilizing the linear response time-dependent density functional theory, NLHFPS has a potential to evolve into a powerful characterization method for nanoscience and nanotechnology.


I. INTRODUCTION
In optical spectroscopy, systems of interest are exposed to light, and their response allows us to explore their structure and composition. A significant part of spectroscopy involves linear effects, such as when light is absorbed or scattered off material targets, allowing their imaging and characterization, teaching us almost solely about dipole-allowed transitions. Nonlinear spectroscopy methods go beyond this limitation, studying otherwise hidden or dark changes applicable to a large variety of systems and processes. 1 Examples of nonlinear spectroscopy include the second-order harmonic generation (SHG) approach, used to study interfaces and adsorbed molecules and serves as high-resolution optical microscopy in biological systems, 2 multiphoton excitation fluorescence (MPEF), as well as various Raman scattering methods. 3 The use of nonlinear spectroscopies, especially in surface and nano-sciences, is growing due to their high interfacial sensitivity. However, results in nonlinear spectroscopies are often challenging to interpret since their description involves much more sophisticated theoretical techniques as compared to their linear counterparts. 4 This article studies the high-frequency limit of the electronic response, singling out a pathway which leads to a major simplification in the description of nonlinear spectroscopies, as long as the observables are analysed after the field acting on a system dies out. We find that the nonlinear behaviour of the system observables is expressible in terms of the linear electron density response function, the latter occurring on the time-scale of the pulse's a) Electronic mail: nazarov.vu@mipt.ru b) Electronic mail: roi.baer@huji.ac.il enveloping shape. By this, we present an approach to the problem of the nonlinear electronic response in the case of high-frequency pulses, which turns out no more theoretically and computationally demanding than the solution of the conventional linear response problem. Specifically, the well-developed methods of the linear response timedependent density functional theory 5 (TDDFT) can be readily invoked, expanding the reach of the latter to the realm of the nonlinear physics.
We validate our theory numerically using the exactly solvable hydrogen atom system propagating under a time-dependent field. Then we consider applications to nano-films and nano-dots, which demonstrate the power of the proposed method by revealing the modes in the excitation spectra of these systems, latent when probed within the linear regime. Finally, we present an example of molecular spectroscopy showing dipole-forbidden transitions.

II. FORMALISM
We consider a many-electron system subject to the time-dependent (TD) modulated periodic potential. We are concerned with solving the Schrödinger equation (in the following, atomic units are used unless indicated otherwise) where the unperturbed Hamiltonian iŝ N and v ext (r) being the number of electrons and the external (electron-nuclear Coulomb) potential, respectively, arXiv:2101.09467v3 [cond-mat.mes-hall] 4 Aug 2022 and the harmonic perturbation is enveloped with the po-tentialŴ For simplicity, we assume that the time-dependence in the pulse potential W (r, t) factorizes, i.e., where C(t) is the pulse envelope and W (r) determines the coordinate dependence of the potential, although, extensions to more general forms of the potential are straightforward.
Our principal result, the proof of which is postponed until Appendix A and the Supplemental Material, is an expression for the probability amplitude to find, after the end of the pulse, the system in its excited state Ψ α =0 In Eq. (5), E α are the eigenenergies of the system, C 2 (ω) is the Fourier transform of the square of the envelope function is the electron density operator, and n = 4 and 2, in the case of the uniform applied electric field and all other cases, respectively. Corresponding F n (r) are Finally, Eq. (5) holds to the leading non-vanishing order ω −n 0 in each of the cases. Further developments (see Appendix A) show, that the time-dependent oscillations in the electron density after the end of the pulse are given by δn(r, t > T ) = 1 iω n 0 e −iωt C 2 (ω)Imχ(r, r , ω)F n (r )dr dω, (10) where χ(r, r , ω) is the linear density response function of the interacting electron system. 6 Furthermore, to the leading order in ω −1 0 , we find for the total energy absorbed by the system during the pulse action ∆E = − π 4ω 2n 0 ω| C 2 (ω)| 2 F n (r)Imχ(r, r , ω)F n (r )dωdrdr .
Clearly, the case of the uniform electric field (n = 4) is relevant to the problem of the illumination by light. Although, strictly speaking, the latter should be described with the transverse vector potential A z (t − x/c), the usual practice is, neglecting the retardation, to reduce the problem to that with the homogeneous A z (t) and then, by the gauge transformation, to the equivalent problem with the scalar potential (7). 7 Apart from the lower bound on the frequency, inherent to our high-frequency asymptotic theory, ω 0 ω low , the neglect of the retardation imposes a standard upper bound ω 0 ω high = c/d, where c is the velocity of light, and d is the size of the system. Another case, n = 2, is relevant to processes with the excitation by longitudinal fields, such, e.g., as with moving charges. 8 This is promising for the construction of TDDFT of the stopping power of matter for fast ions beyond the adiabatic approximation for the exchangecorrelation potential, which theory now exists in the lowvelocity limit only. 9,10 Importantly, in Eqs. (10) and (11) we witness a hybridization of linear and quadratic response quantities: the linear density-density response function is multiplied by the quadratic frequency envelop C 2 (ω). In the illustrative calculations below, we will see that such hybridization leads to interesting effects.
In the field of the light-matter interactions, the application of the acceleration-frame method of Kramers and Henneberger (KH) 11,12 has led to a great many advancements in the theory. [13][14][15][16][17][18] Instructively, our formulas above can be re-derived in an alternative way using the KH method, as it is shown in Appendix B. However, this is possible to do in the case of the uniform field only (n = 4), since this case is inherent within the KH formalism.

A. Hydrogen atom
We now investigate how the high-frequency limit is approached as the frequency increases by calculating a precisely solvable system, namely, the hydrogen atom. First, assuming an atom, initially in its ground state, is subjected to the doubly modulated Gaussian pulse with a spherically symmetric quadrupole potential we numerically time-propagate the Scrödinger equation (1). In the pulse (12), the carrier frequency ω 0 serves to set the scene for the high-frequency regime, while the second frequency ω couples the pulse to the excitations in the system. Upon the end of the pulse, we look at the populations of the excited states, plot them in Fig. 1 versus the enveloping function frequency ω (the second frequency), and compare with the asymptotic limit. The latter, according to Eqs. (5), (9), and (12) is given by where φ n,s (r) are the hydrogenic s-orbitals and n are the corresponding eigenenergies, and we have restricted the comparison to the transitions to the s-states only. We note that the spherically symmetric quadrupole po-  Excitation probablity Excitation probability [the modulus squared of Eq. (5)] upon the end of the pulse of Eq. (12), from the ground-state of the hydrogen atom to a number of its excited s-states. The solid black line is the asymptotic limit of Eq. (13). Spectra at finite frequencies are obtained by the numerical propagation of the TD Schrödinger equation (1). The parameters of the pulse used were σ = 15 a.u. and W0 = 0.125 a.u. tential (12) is purely model one, which we use to demonstrate the convergence of the numerical solution of the Schrödinger equation to the asymptotic solution (5) for a non-uniform field (n = 2).
Similarly, in the case of the uniform field (n = 4), we propagate the system under the potential 19 For the hydrogen atom where Y lm (θ, φ) are spherical harmonics. Evidently, only transitions from the ground state to s-and d-states are possible, which have the following amplitudes where, in Eq. Figures 1 and 2 demonstrate the convergence, with the growth of ω 0 , of the excitation processes' outcome to their ω 0 → ∞ limits of Eqs. (5), for the cases of the quadrupole and dipole exciting potentials, respectively. Remarkably, in the quadrupole (Fig. 1) and the dipole (Fig. 2) cases, the asymptotic regime is approached in very different ways: in the former case, the peaks' positions and shape change dramatically with the frequency growth, while in the latter, the amplitude of the peaks varies monotonously only. In the dipole case, the convergence, with respect to peaks' amplitudes, is very slow, and it is not reached at practically achievable values of ω 0 . 21 At the same time, the excitation energies (peaks' positions), even at moderate values of ω 0 , are very well reproduced by the asymptotic theory. We point out and emphasize that, while the asymptotic limit holds for an arbitrary system, the speed of the convergence is systemdependent. This is confirmed by Fig. 3 with the use of the fictitious system of the hydrogenic atom with the nuclear charge of Z = 0.25. Since the asymptotic theory is expected to be the more accurate the larger is ω 0 compared to the characteristic excitation energies in a system, in the Z = 0.25 case we observe a much faster convergence compared to the Z = 1. For peaks in Fig. 3 to remain resolved, a large width of the pulse σ = 200 a.u. was chosen in the calculation with Z = 0.25. Further particulars of the solution of the TD Schrödinger equation and the issues of the convergence to the asymptotic limit are presented in Appendix C. Excitation probablity 1s → 2s 3s 4s 5s It is highly instructive to follow the excitation process in time, from the pulse beginning to its end, in order to understand how the system reaches its final state. As can be seen from the derivation [e.g., Eq. (S.15) of the Supplemental Material], the linear response does contribute to the pumping during the pulse action, but it passes a cycle from increasing to decreasing the population of excited states, with the zero net result. On the contrary, the quadratic response does not completely reverse itself, which results in the residual occupancies of the excited states upon the pulse's end. In Fig. 4, we plot the timeevolution of the population numbers of the 2s-and 2p-orbitals of H atom under the action of the pulse of Eq. (14). We observe the principal difference between the change of the occupancies of the s-and p-levels: while the latter gets much more (approximately three orders of magnitude) populated in the middle of the pulse duration, it gives the electron away upon the pulse end. At the same time, the former keeps the accepted electron with a finite probability. This type of behaviour is characteristic of spherically symmetric systems in the high-frequency regime, which is in agreement with our asymptotic theory. This is the linear response that dominates the s → p transition at the time of the pulse duration, which is gone upon the pulse's extinction. In particular, we conclude that the usual dipole selection rules do not hold in this process. At this point we note that, while TDDFT of the electronic response in the high-frequency limit was studied in Ref. 22, it is important to emphasize the principal difference between the physical situation considered in that reference and in the present paper. Ref. 22 deals with the response to the monochromatic field, thus considering a continuous wave. In that regime, the linear response persists in the high-frequency limit and it is, usually, prevailing. On the contrary, here we consider the excitation by a pulse of finite duration, the carrier frequency of which is asymptotically high. We focus on the behaviour of a system after the end of the pulse, in which case we find the total suppression of the linear response, while the nonlinear one is describable in terms of the linear response TDDFT.
With the use of Eqs. (16) and (17), in Fig. 5 we compare the excitation and ionization processes' probabilities for the hydrogen atom initially in its ground-state and exposed to the Gaussian pulse. We conclude that the ionization is dominant for short pulses, in which case a sudden impact strips off electron, while, for longer pulses, transitions to excited bound states become preferential. We also note that transitions to the d-states play insignificant role compared to those to the s-states.

B. Jellium slab
We proceed by considering a slab of the thickness d with the positive constant background charge density n + = ( 4 3 πr 3 s ) −1 , where r s is the 3D density parameter. Within the Kohn-Sham (KS) density-functional theory (DFT) 23 and using the local density approximation (LDA), we calculate the ground-state KS band structure and electron density. To this system, we apply the doubly modulated dipole pulse of Eq. (14), and we use our theory to determine the total energy absorption in the slab in the high carrier frequency regime. The problem being one-dimensional, the operator in Eq. (8) reduces to the Laplacian, and we have by virtue of the Poisson law where Θ(x) is the Heaviside's step-function. Resulting absorption spectra, obtained by Eq. (11) with the use of the adiabatic time-dependent LDA (ATDLDA) in the construction of χ(r, r , ω), 5 are presented in Figs. 6 and 7, for r s = 5 and 2, corresponding to the jellium model of the metallic potassium and aluminum, respectively. The following observations are made: (i) Similar to the case of the hydrogen atom, due to the integration with C 2 (ω) in Eq. (11) and due to the form of the pulse (14), spectra in the left panels of Figs. 6 and 7 as functions of ω are governed by SHG and, accordingly, peaks' positions scale to half the frequencies of the corresponding excitations; (ii) In the linear regime (right panels in Figs. 6 and 7), spectra are dominated by the bulk plasmon (BP) peak, the intensity of which crucially depends on the share of the bulk, i.e., the slab thickness d. On the contrary, the nonlinear spectra in the high-frequency regime (left panels in Figs. 6 and 7) weakly depend on d, suggesting that the surface excitations dominate them. The prevalence of the surface response can be understood by noting that χ(r, r , ω)dr = 0 (no reaction to a constant potential) and, therefore, both the deep interior and exterior of the slab, by Eq. (18), do not contribute appreciably to the integral of Eq. (11).
Notably, in the left panel of Fig. 6 we observe a strong peak with the maximum at 2ω ≈ 0.88ω p . The coun-

Linear response
Energy absorption (arb.u.)   terpart of this peak in the linear response regime (right panel of Fig. 6) is positioned at ω ≈ 0.83ω p , and it is known as the multipole surface plasmon (MP). 24 Because of the BP suppression, MP is very prominent in the left panel of this figure, which makes the high-frequency nonlinear technique an ideal tool to study this otherwise subtle type of excitation. It is instructive to note that F 4 (z) of Eq. (8) provides, effectively, the impact mode of the complementary linear response problem, 25 which is known to be favourable for MP excitation. 26 In Fig. 7 (r s = 2), left panel, we also see a prominent broad peak at 2ω below the BP frequency, while MP is not discernible in the linear response spectrum in the right panel. We, therefore, conclude that the corresponding excitation exists at the surface of metallic aluminum, and the highfrequency nonlinear technique provides a unique way to detect it. At the same time, the traditional method of electron energy loss spectroscopy (EELS) does not possess sufficient sensitivity. 24 The oscillating structures at 2ω > ω p in Figs. 6 and on both sides from ω p in Fig. 7 differ for different slab thicknesses, and they can, therefore, be attributed to the interference effect between the two surfaces of the slabs. Finally, the absence of the conventional (dipole) surface plasmon (SP) peak at ω s = ω p / √ 2 is due to the strictly normal to the surface direction of the exciting field (q = 0), in which case the amplitude of the SP vanishes.
To quantitatively verify the above picture, in Fig. 8 we plot the Fourier transform of the density oscillation in the asymptotic regime [Eq. (10)] and compare it with the linear response density oscillation. Clearly, in the former case, the oscillation is mainly confined to the vicinity of the surfaces of the slab being largely suppressed in the interior. On the contrary, in the linear response regime, oscillations predominantly occur in the bulk of the slab.

C. Jellium sphere
In contrast to a slab, for a sphere, the second derivative in the RHS of Eq. (8) does not reduce to Laplacian and, consequently, F 4 (r) is not given by the positive back-ground density only. Instead, we have where R is the radius of the rigid positive-charge background. Due to the symmetry, the density-response function χ(r, r , ω) splits in angular momentum into χ lm (r, r , ω), the latter acting separately on each harmonic of the externally applied potential. The problem becoming one-dimensional again, we calculate χ 00 and χ 20 , apply them to Eq. (19), and plug the result into Eq. (11). We consider the same form of the doubly modulated pulse of Eqs. (14) as previously. In Fig. 9, results of calculations for two spheres, with radii R = 30 and 40 a.u., and the density parameter r s = 5, are presented, for the nonlinear ω 0 → ∞ and the linear-response regimes, in the left and right panels, respectively. Within the classical electrodynamics, a sphere of the Drude metal supports an infinite series of Mie plasmons ω l = l/(2l + 1)ω p , l = 1, 2, . . . . 27 In the monochromatic linear-response (right panel of Fig. 9), we observe the p-mode only of this series, red-shifted by the quantum size effect.
According to Eq. (19), energy absorption in the nonlinear ω 0 → ∞ regime (left panel of Fig. 9) originates from the superposition of the s-and d-modes. As plotted versus the second modulation frequency ω, it reveals a rich spectrum of the underlying excitations. The leftmost feature near 0.57ω p comes from the d-mode Mie plasmon ω 2 , red-shifted in the quantum calculation. The broad dominating peak with the maximum near 0.80ω p does not have an analog within the classical electrodynamics, and, similar to the multipole plasmon modes in the case of a slab, it becomes accessible with the use of the high-ω 0 nonlinear regime. A signature of the bulk plasmon on the right shoulder of this peak can also be observed, indicating the possibility of the direct recognition of the constituents of nano-particles by their bulk plasmon frequencies ω p with the use of laser pulses. The latter is, obviously, impossible in the linear-response regime. We also note structures above ω p , which are due to the (dressed) single-particle excitations affected by the quantum interference.
Finally, we consider molecular electronic spectroscopy. Referring back to the above-discussed very short pulse spectroscopy, in the dipole interaction case, we used timedependent local density approximation calculations to produce linear response estimates of the high-frequency energy absorption (using Eq. 11) and, in Fig. 10, compare to standard low-frequency energy absorption for the ethylene molecule. The spectra' differences in the two regimes are due to the dipole versus F 4 selection rules, emphasizing the high-frequency spectroscopy's aptitude to probe the excitations forbidden in the linear regime. See Appendix D for details concerning this calculation.

IV. DISCUSSION AND CONCLUSIONS
We have considered excitation of a quantummechanical system by an externally applied electric field of high-frequency ω 0 and finite duration in time. After the end of the pulse, the state of the system being a superposition of the eigenstates of the unperturbed Hamiltonian, the expansion of the corresponding transition amplitudes in the power series in ω −1 0 has been performed, with the leading terms found of the order ω −4 0 for the uniform applied field (dipole case) and of ω −2 0 , otherwise. We have demonstrated that, to the leading order in the inverse frequency, the quadratic, rather than the linear, response determines the excitation process. Nonetheless, we have also shown that all the information necessary to describe this nonlinear excitation regime is contained in the linear density response function of the system under consideration. The problem has been thus reduced to that of the linear response time-dependent density functional theory, for which practical methods of solution, at various levels of accuracy and sophistication, are well established.
Further, we have found that a specific pulse shape, modulation by the second (low) frequency can be advantageous as a probe, delivering spectra of excitations in the nonlinear response regime. In our illustrative applications, to the jellium model nano-films and nanodots, plasmonic modes undetectable or challenging for the detection by the linear optical spectroscopy or electron energy-loss spectroscopy have been discerned. We point out that the high carrier frequency is out of resonance, and its only role is to set the scene for probing the system with the second frequency, the twice of the latter being in resonance with the system's excitations.
Based on our findings, we propose a spectroscopic technique, which we provisionally name the Nonlinear High-Frequency Pulsed Spectroscopy. Our results show that NLHFPS, i.e., exposing an explored system to a finiteduration high-frequency electric field with low-frequency modulation, allows for an efficient nonlinear spectroscopic probe of modes inaccessible or hardly accessible by other techniques. A significant asset of the novel method is its ease of interpretation, enabling a detailed comparison between experiment and theory. This benefit stems from the results' direct dependence on the target material's density-density response function. As demonstrated here, NLHFPS can uncover rich and profound physical phenomena hidden from more conventional methods. The data that support the findings of this study are available from the authors upon reasonable request.
Appendix A: Derivation of Eqs. (5)- (11) In the interaction representatioñ the problem of the solution of Eq. (1) turns into that for the equation or for the equivalent integral equatioñ where by Ψ α we denote the set of eigenfunctions of the Hamiltonian (2), we assume that W (r, −∞) = 0, and the system is initially in its ground-state Ψ 0 . Performing several consecutive integrations by parts in Eq. (A4), assuming the pulse to be of finite duration [Ŵ (r, +∞) = 0] and ω 0 to be large, we obtain, after keeping the terms up to ω −4 0 onlỹ A lengthy derivation of Eq. (A5) is given in full in the Supplementary Material. 28 The commutators in Eq. (A5) can be expanded as (A9)

Non-uniform field case
We evaluate the commutator (A6) to If, furthermore, the factorization of Eq. (4) holds, then we arrive at Eq. (5) with n = 2, where an extra exponent e −iEαt appears in the Schrödinger representation. Equation (A11) gives the transition amplitude to the leading order in ω −1 0 unless the term in the square brackets under the integral is independent on r. However, in the latter case the integration ofn(r) over r produces a constant N , and then the RHS becomes zero because of the zero the matrix element. This is, exactly, what happens if the field is uniform, as can be seen from Eqs. (7) and, therefore, this case requires a separate consideration.

Uniform field case
With the use of Eqs. (A5), (A7), and with the commutator relations and noting that in Eq. (A5) the sum of the 4th, 5th, and 6th terms on the RHS evaluates to zero, as it can be directly verified, we immediately arrive at Eq. (5) with n = 4.

Density oscillations and energy absorbed
The time-dependent density is given by where the last equality is due to the normalization of Ψ(t). Therefore, With account of Eq. (5), we conclude that the leading term in ω −1 0 on RHS of Eq. (A18) is the second one, while, for the same reason, Ψ(t)|Ψ 0 = e iE0t must be set in the latter. Then δn(r, t > T ) = 2 Re e iE0t α =0 Ψ 0 |n(r)|Ψ α Ψ α |Ψ(t) .
(A19) Combining Eqs. (5) and (A19), we have or Recalling the spectral representation of the many-body interacting density response function where η is a positive infinitesimal, we can rewrite Eq. (A21) as Finally, the separation of the real part on the RHS of Eq. (A21) can be dropped since the remaining expression is real already (see the footnote 6). For the total energy absorbed by the system from the pulse, we can write which, with the use of the completeness of the basis set, can be rewritten as and then, by Eq. (5), finally written in the form of Eq. (11).
Appendix B: Derivation in the Kramers-Henneberger's acceleration frame For an arbitrary u(t), if a function Ψ KH ({r}, t) satisfies the equation then the function where satisfies the equation and noting that u (t) = E 0 C(t) cos ω 0 t, we turn Eq. (B4) into Eq. (1) in the case of the dipole applied potential. Expanding in Eq. (B1) up to ω −4 0 , we have with the use of Eq. (B5) which in the interaction picture is written as and, therefore, In Eq. (A3), we expandΨ(r, t) as Ψ(r, t) = lmax l=0 nmax n=0 a n,l (t)F n (r)Y l0 (θ, φ), where F n (r) = λ 3/2 f n (λr), (C2) n (x) are the generalized Laguerre polynomials, and α and λ are positive parameters. The basis set in Eq. (C1) is orthonormal and complete with any α and λ. Although we have been using α = 2 and λ = 1, the convergence of the method has been verified by comparing results with those obtained with other values of these parameters.
Matrix elements of the unperturbed HamiltonianĤ 0 and the time-dependent partŴ (t) were obtained exactly with the use of the recurrence relations for the generalized Laguerre polynomials. 29 The problem was thus reduced to that of the propagation in time of the system of the linear ordinary differential equations for a n,l (t), which was carried out by means of the Magnus expansion. 30 For the hydrogen atom, the Schrödinger equation (1) reads (C4) By scaling the variables r = Zr, t = Z 2 t, we see that Ψ Z (r, t) = Z 3/2 Ψ(Zr, Z 2 t) is the solution to the complementary problem for the hydrogenic atom of the nuclear charge Z From Eq. (C5) we conclude that the frequency ω 0 scales as ω 0 → Z 2 ω 0 , which explains the faster convergence of the solutions to its ω 0 → ∞ limit we have observed in Fig. 3 for Z < 1.

SUPPLEMENTARY MATERIAL TO THE ARTICLE 'HIGH-FREQUENCY LIMIT OF SPECTROSCOPY' BY VLADIMIR U. NAZAROV AND ROI BAER
Derivation of Eq. (A5).
From Eq. (A4), by the integration by parts, we can writẽ and, with the use of Eq. (A3), Continuing in the same waỹ we can rewrite Eq. (S.7) as  As all the previous equations starting from Eq. (S.1), Eq. (S.15) is exact at any time t. Upon the end of the pulse, t → +∞,Ŵ (t) → 0, as do all its time derivatives. Therefore, all the out-of-integrals terms on RHS of Eq. (S.15), except for Ψ 0 , (terms from 2nd to 16th) become zero. We, therefore, can writẽ where all the terms of the order ω −n 0 , n > 4, have been neglected, which allowed us to replace Ψ with Ψ 0 everywhere but in the 2nd term. The last step is to expand the 2nd term to the same order, which is done by using Eq. (S.15) agaiñ Ψ(+∞) = Ψ 0 + 1