{'fraction_non_alphanumeric': 0.03259909949677761, 'fraction_numerical': 0.008629822547894411, 'mean_word_length': 4.464182344428365, 'pattern_counts': {'":': 1, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exact energy-eigenstates of the Coulomb-Stark Hamiltonian

An approximation-free, numerically efficient algorithm is presented for the Hamiltonian eigen-states of the Stark-Hydrogen problem describing a quantum particle exposed to the central Coulomb force and a homogeneous external field. As an example of application in a state-expansion with continuous energy, we calculate the time-dependent wavefunction of an electron tunneling from a hydrogen atom suddenly exposed to an external electric field.


I. INTRODUCTION
In calculating the time evolution of a pure quantum state, the most fundamental approach is to expand the wave function in terms of the eigenstates of the Hamiltonian. While this is relatively straightforward in finitedimensional Hilbert spaces, infinite-dimensional cases are more difficult. Calculating the eigenstates of the system and using them as a basis for the space of quantum states is particularly challenging when the Hamiltonian has a continuous spectrum. In such a situation, at least some of the eigenstates are not normalizable in the usual sense, and normalization to the Dirac-delta function in energy is rather subtle from the numerical standpoint. For this reason, in most practical situations one resorts to approximations, for example by expanding solutions into approximate, "discretized-continuum" basis.
An important example of a system with a continuous spectrum is the Stark-Coulomb Hamiltonian [1], i.e. the problem of a particle subject to a constant external electric field and a central Coulomb force. Even an arbitrarily weak external field causes a qualitative transformation of the original hydrogen-atom spectrum into infinitely degenerate continuum encompassing the whole real axis. The zero-field discrete-energy states mutate into Stark resonances [2] which further complicate exact treatments. Perhaps not surprisingly, numerically exact calculation of the energy-eigenstates for this system has not yet been reported. In this paper we present a robust, approximation-free method to calculate all the eigenstates in a way that allows practical applications, including expansion of time-dependent wave functions.
Obviously, the Stark-Coulomb Hamiltonian is important for a number of applications. Recently, experiments on atomic hydrogen [3] exposed to strong electric fields provide a vital testing ground for theory and simulations [3][4][5][6][7][8]. Moreover, in ionization of atoms and molecules [9], electrons liberated from the neutral systems experience the Stark-Coulomb potential at larger distances from the nucleus [10][11][12]. This composite potential affects the properties of the electronic wave functions contributing to strong-field ionization [13][14][15][16], and high-harmonic generation [11]. Thus, the numerically exact wave functions of such a problem can be of great practical use.
The rest of this paper is organized as follows: We start by laying out the theoretical foundations, by reviewing the Stark-Coulomb Hamiltonian in parabolic coordinate system. We proceed by calculating the exact eigenstates including their proper normalization to the Dirac-delta function in energy. We then demonstrate the power of this method by using it to describe the time evolution of an electron tunneling from the ground state of hydrogen after a sudden exposure to an external electric field.

II. THE COULOMB-STARK PROBLEM
We start with the formulation of the problem of a quantum particle (electron) in a Coulomb potential and subject to an external electric field. It has been recognized that for the calculation of Stark resonances the parabolic coordinate system is the most suitable [17][18][19][20][21]. Naturally, this is also the case for the continuum-energy eigenstates we are interested in, and we shall use this frame of reference as well. In this section we recall the equations which serve as our point of departure.
The coordinate relations between the parabolic coordinates (u, v, φ) and Cartesian coordinates (x, y, z) are: y = (uv) 1/2 sin φ, φ = tan −1 y x u, v ∈ [0, ∞), φ ∈ [0, 2π) , and the volume element reads: The main advantage of this coordinate system is that the Schrödinger equation of an electron in a Coulomb potential of a singly-charged nucleus (hydrogen) remains separable even in the presence of an external electric field. The time independent Schrödinger equation (TISE) becomes: where F is the strength of the homogeneous external field. Without any loss of generality we will assume that this quantity is positive.
As we will see in a moment, the energy eigenstates Ψ qW for this Hamiltonian can be labeled by two kinds of quantum numbers; The first labels a continuum of energies, W ∈ R, encompassing the whole real axis, and the second is a composite discrete pair q = {n v , m} where m is an integer standing for the usual "magnetic" quantum number, and n v is a whole number counting the zerocrossing along the parabolic axis v.
An arbitrary time-dependent wave function can be expanded in this eigenenergy basis as: where A q (W ) is the wave function in the energy representation, and it can be understood as the overlap integral between the given quantum state (at time t = 0) and the corresponding eigenstate. Note that Ψ qW depend parameterically on F , so there is a continuum of different bases, distinguished by the value of the external field. For each fixed F , we have an infinitely degenerate continuum of energies. Our task is to design an accurate and efficient numerical algorithm to evaluate these states without any approximations.
Because of the symmetry, the sought after eigenfunctions depend on the azimuthal quantum number as Ψ ∼ e imφ . The methods required for different m are completely analogous, so for the sake of simplicity we will restrict ourselves to the case m = 0. In order to further simplify our notation, we will use a shorthand, n = {n, m = 0}, for the discrete part of the eigenstate label.
Using the method of separation of variables we can turn the partial differential equation into a pair of ordinary differential equations. With this in mind, we rewrite the TISE as and note that it has the form We use the following ansatz leading to these separated equations for U nW (v) and with the two separation constants tied by the constraint Note that the above equations are well-known from the Stark-resonance problem in atomic hydrogen [17]. Unlike in the discrete non-Hermitian resonance calculation, we look for a continuum of standard, i.e. Hermitian, real-valued energy eigenstates. The distinction between these two kinds of eigen-value problems boils down to the boundary conditions imposed on the eigenfunction at infinity. While the frequently-studied Stark resonances with complex-valued energies must behave as outgoing waves, the Hermitian eigenfunctions possessing real-valued energies behave as standing waves at infinity. In other words, the main challenge here is not in finding the eigenvalues, but in designing an approach to obtain the properly normalized wave function, and it is the continuum nature of the spectrum that makes the problem difficult.
An important aspect of this work is that our approximation-free numerical calculation is sufficiently efficient and accurate at the same time, so that it makes it possible to use the resulting continuum-energy basis for state expansion of arbitrary wave functions.

III. CALCULATION OF EIGENSTATES
The most challenging step in the calculation of the wave functions corresponding to a continuum of energies is to ensure the correct normalization for U nW . For applications in which the set of eigenstates is used as a basis, it is necessary that the resolution of unity in the energy space is achieved. This will require construction of the inner and outer solutions which smoothly connect to each other, each supplying a crucial piece of information. We take the following steps to calculate U nW with the correct normalization factor.
Inner solution via "analytic continuation": Both V nW (u) and U nW (u) can be easily obtained in the form of a series in the vicinity of an arbitrary point, say u = a, provided that the function value and its derivative are known at this point. Utilizing a large number of terms in such a series, such initial data can be obtained at u = a+ δa, and the series expansion can be re-done around this new "center." This is akin to the analytic continuation of holomorphic functions of a complex variable. The process can be extended to an arbitrary distance from the origin, but it leaves the normalization of U nW (u) undetermined.
Outer solution via carrier-envelope method: In order to fix the normalization, we examine the large-u behavior of the differential equation for U and split the function into a "carrier-wave" and its slowly changing "envelope." Such a form allows us to ensure the correct normalization to a Dirac-delta function in energy. However, it will not fix the relative phase-shift between the incoming and outgoing waves.
Joining the inner and outer solutions: Because both representations solve the same differential equation, the remaining degrees of freedom, namely the normalization of the inner part and the phase shift of the outer part can be found by requiring that the two agree at two arbitrarily chosen points.
Before we execute this plan, it should be useful to point out some qualitative differences between V nW (v) and U nW (u), so we start with their spectral properties next.
A. Discrete and continuous spectra and the eigenstate normalization One can see that for large values of u and v the equations for U and V will turn into: and − 1 2 These equations resemble the TISE of a particle living on a half-line with a potential that is pulling it from the origin (in the equation for U ) and/or pushing it towards the origin (in the equation for V , which is sometimes called the "quantum bouncer" [22]). Hence the spectrum of U is continuous while the spectrum of V is discrete in the following sense. For any given real value of W (with a fixed F ) there is a discrete (infinite) set of eigenvalues z v (n, F, W ). In contrast, for a given z u = 1 − z v , the equation for U has a deltanormalizable solution for any W ∈ R. In other words, the normalization for V nW s is to "Kronecker delta," while the normalization of U nW s is to "Dirac delta." The nature of the spectrum is reflected in the following normalization relation which the eigenstates must obey, This in turn implies the orthogonality relations between the U and V functions: and It is the second, Dirac-delta normalization condition that requires more attention if one aims for a practical method to utilize these states in expansion of arbitrary wave functions.
The third condition holds in the distribution sense and is actually independent of the normalization. While it may seem that the term proportional to v in the normalization integral (15) was left out to obtain (17), one must evaluate all wave function overlaps with the whole volume element proportional to (u + v)/4. The detailed derivation of the above normalization conditions is given in the Appendix.
For readers who may prefer a more intuitive normalization argument, we note that these conditions may be obtained by following the recipe given in Refs. [23,24]. This approach splits the wave function into the incoming and outgoing waves, and requires that the total probability outflow in the latter equals 1/2π for all continuum energies. It is straightforward to demonstrate that the above normalization satisfies this requirement.

B. VnW eigenstates
While for F > 0 closed-form expressions for the V nW solutions to (11) are not known, they are normalizable and can be easily calculated using the method of series expansion. We will use the following representation in which M is a series in v normalized such that M (W, z, F |v = 0) = 1 and v n0 represents the first coefficient in the series expansion for a normalized eigenstate. (note that argument z of M here and in what follows stands for the separation constant z v or z u , not a coordinate). The series expansion can be written as: and substituting this in the differential equation we obtain Using this recurrence relation, hundreds of terms can be calculated efficiently. Note that this series expansion does not represent the solution to our problem until the separation constant z = z v is determined. For a general z, the above function diverges at infinity, while we seek a normalizable V nW (v). We can numerically calculate z v by demanding that the series expansion approximation remains reasonably close to zero at a sample v 0 far enough from the origin and adjusting z v gradually until the function tends to zero for large arguments v. In this process, one obtains a function which can have none or several zeros. The number n of these nodes is the first part of the discrete label of the energy eigenstate. At this point we have determined for any fixed W and positive F . The behavior of z v as a function of energy W is illustrated in Fig. 1. Having fixed z v (n, W, F ), formulas (20) and (21) provide the shape of the function V nW (v). Figure 2 shows a few as an example. Finally, we determine the value v n0 by imposing the normalization condition (16). The integrals required for this calculation can be obtained analytically in terms of the series-expansion coefficients, so v n0 can be highly accurate.
From the numerical standpoint, calculation of the Vpart of the energy eigenstate is straightforward, and various algorithms can be utilized to find z v to a high accuracy. For the illustrations in this work, we utilized between hundred and two hundred terms for the series expansion around v = 0. This is sufficient to reach v ≈ 40 − 50 with sufficient accuracy. Should V nW (v) be needed for even larger v, the technique described next can also be used.

C. UnW eigenstates
The more challenging part of the problem is to calculate the U eigenfunctions which are not normalizable and have a continuous spectrum. After we determined the z v using the properties of the V -functions, we can use the relation between z v and z u eigenvalues and calculate The behavior of U at small u can be calculated using the same method as the V eigenfunctions; We use a series expansion again, and since their differential equations are related by the transformation we can use the same series by simply applying the above modification. The small u behavior of U nW (u) can thus be written as: As we know the U eigenfunctions are not normalizable, so u n0 cannot be determined as easily as v n0 . This is in fact the most subtle part of the whole procedure and we will address it in its due course. However, it will be necessary to calculate U nW (u) for a large u, typically in the range of u ∼ 10 3 . This is not possible with the series for the vicinity of the origin -with hundred to two hundred terms in the series one can obtain the function with a good accuracy for up to u ∼ 40. However, the series can be easily analytically continued to reach even very large values of u as follows.
Let us assume that we have used the above series to calculate f = M (a) and p = M ′ (a) for some a > 0, utilizing a large number of terms (e.g. 150, to ensure sufficient accuracy). As a next step we can obtain a series expansion of M centered at a: The recursion relations for the coefficients can be obtained by inserting into (11) (note the we use the same master function M for both V and U functions), giving This specifies a series expansion valid around u = a, for the function value and its derivative at this point set to f , and p, respectively. Note that this recursion scheme applies equally to both U and V functions since both are expressed in terms of the same master function M ; the only difference is in the parameters passed to M . Obviously, this is a more complicated recursion than (21), because it solves a more general problem characterized by three additional parameters (f, p, a). However, this is the most important piece of the whole scheme because it provides locally exact solutions which can be evaluated to an arbitrary precision. Next, one can repeat the same procedure, obtaining a new series centered around u = b > a, and keep repeating the re-expansion until reaching the uregion of interest. Typically, with long-double precision and utilizing between hundred and two hundred terms, one can create a set of power-series expansions, center of each shifted further beyond the previous. Such a representation of the function is fast to evaluate and remains accurate for u as large as several thousands. One could be worried that errors can accumulate by a repeated reexpansion, and they indeed could. However, because one can easily increase the number of terms in the expansion, the accuracy is only limited by the number of bytes used in the floating-point type. So the error propagation and accumulation turns out to be a non-issue with the long double type (simple double type can also be used with smaller distance between patch-centers b − a) when the two values (f and p) passed from a patch to the next are calculated to fifteen to eighteen digits. The functional shape of U nW (u) is distinct from that of V nW (v), as should be expected for the continuous energy spectrum. Figure 3  figure, the oscillation frequency increases toward infinity, reflecting the fact that the particle is accelerated by the external field.
At this point, we have a practically usable implementation of the eigenstate up to a multiplicative constant u n0 . This remaining piece is a function of W and F , and it carries information crucial for the set of functions to be used as a continuum basis. It will be obtained next.

D. Large-u behavior and normalization of UnW
In order to fix u n0 , we need to turn to the Dirac-delta normalization requirement (17). Obviously, an asymptotic solution for large u is needed, and the strategy is to use an ansatz to separate a fast changing "carrier" of the wave function from its slowly changing "envelope". Looking back at the asymptotic differential equation (13), one can see that the relevant solutions should behave essentially as Airy functions [25]. The following linear combinations of Airy functions, is suitable to serve as the carrier. More concretely, our ansatz for U can be written as: (29) in which N U is some normalization constant, α = −(2F ) 1/3 , δ(W, F ) is as yet undetermined phase shift, and H(u) is the slowly changing envelope of the wave function. The rationale behind this is that the envelope can be calculated with very modest numerical effort. In fact, its analytic asymptotics will be sufficient for many purposes. Substituting this ansatz in the differential equation for U we get: stands for the logarithmic derivative of Ci + . Unlike the function itself, R changes slowly with u and as a consequence, the differential equation for the envelope is easy to solve numerically. Indeed, such a calculation can yield H(u) which can be used in practice. Nevertheless, here we will only utilize an asymptotic solution for H valid for large u, so that we will avoid numerical treatment of H. One reason to avoid numerical ODE solution and favor an analytic approach is that any numerical ODE solver is designed for universal usage, and as such it can not take advantage of the concrete equation solved. Another reason to avoid numerical ODE solution is the superior accuracy and speed of the algorithm given next.
Using the asymptotic behavior of the Airy functions, lengthy but straightforward calculations give, for u → ∞, where more terms can be calculated with some effort.
We only need the first term in (32) to obtain N U , which is fixed such that the normalization condition is satisfied. This requires a calculation involving the asymptotic behavior of U that is determined by the carrier wave Ci + alone, and leads to the normalization factor Having properly normalized the eigenstates to the Diracdelta in energy, we can proceed to calculate the last remaining unknowns namely the phase shift δ n (W, F ) and u 0n (W, F ). The inner solution in terms of the series expansion (24) and the outer solution expressed in (29) satisfy the same second-order differential equation, so in order to obtain these two unknown parameters it suffices to require that the two representations agree at two arbitrary locations. Taking two arbitrary points u 1,2 we have a system of two equations (i = 1, 2): from which u n0 (W, F ) together with δ n (W, F ) are calculated numerically. In practice, the left hand side of these equations is obtained by the analytic continuation with the series-center a in the vicinity of the u 1 ∼ u 2 . The advantage of the "analytically continued" representation of U is that the inner-outer join points u 1,2 can be taken to a very large distance form the origin where the asymptotics of the envelope (32) can be used. Thus, no numerical-ODE solution of (30) is required. Quantities u n0 (W ) and e iδn(W,F ) are of central importance here. They are illustrated in Fig. 4. The "jump" around W ≈ −0.5 is due to the Stark resonance born from the ground state. Right: u00(W ) exhibits a sharp peak due to the same resonance (it is cut off in this graphics).
While the phase shift is a useful quantity when one aims to calculate, for example, the quasi-classical approximation of the wave function at large distances from the nucleus, u n0 (W ) is necessary to finalize the calculation of the inner solution for U nW (u). A most prominent feature in u n0 (W ) is the peak corresponding to the Stark resonance. This requires that u n0 (W ) is sampled on a fine grid in the vicinity of W ≈ −1/2. A fit with a simplepole function can yield the value of the complex-valued resonance energy accurate to a part in a million [26,27].

IV. SUMMARY OF THE ALGORITHM
We have thus arrived at the central result of this work, which is the algorithm for calculating the eigenstates of the Stark Coulomb problem. The method can be summarized as follows: 1. Fix W , F and n: We start by choosing a fixed value of the external field F , a real energy eigenvalue W ∈ (−∞, +∞), and an integer n which is the desired number of zero-crossing the wave function has along the v-axis. For many experiments utilizing femtosecond optical pulses, F is around 0.01 − 0.05 in atomic units. The relevant interval for W is from about -0.6 to +0.6; 2. Calculate z v (n, W, F ): Using (20) and (21) find numerically a value of z = z v (n, W, F ) which ensures that the function M (W, z v , F |v → ∞) → 0 converges to zero at large values of v, and that there are n zeros of M . This gives the solution to the eigenvalue problem (11).
3. Construct V nW : Having found z v (n, W, F ), we have all necessary ingredients to calculate V nW (v). The remaining piece is the normalization factor v n0 (W, F ) in (19) which we can calculate by integrating over |V nW (v)| 2 to satisfy the normalization condition (16). This completes calculation of the normalized V nW (v) function.
To fix the normalization, we calculate analytically continued series expansions (25) for u up to several thousand. Finally, u n0 (W ) is calculated by solving the system of two equations given in (35) using asymptotic representation (32) for the envelope. As a by-product in this step, the phase-shift δ n (W, F ) is also obtained. This completes the calculation of U nW (u) satisfying the Dirac-delta normalization condition (33).
Note that the same "continuation method" employing formulas (27,26) can be also applied to functions V as the only difference is in the input, where the sign of F needs to be reversed and the appropriate separation constant z v or z u must be used.
5. Obtain the energy eigenstate Ψ nW : The eigenstate for the chosen W and n is given by the product in (9), expressed as

IMPLEMENTATION AND USAGE CONSIDERATIONS
The method presented in this work deviates from the numerical solutions implemented on grids (e.g. [24]) in that the resulting algorithm is essentially a formula which is both simple to implement and very fast to execute. The code can be written in less than three hundred lines of c++, and its core for the series expansion is much smaller yet. The resulting speed of the algorithm is such that the compute time is not an issue. The computational complexity is comparable to that of some special functions.
The implementation does utilize a few "hyperparameters" which control the evaluation. These include the order of the series expansion (o v and o u ), the size of the patch p v and p u which is the distance between the centers of the series-expansion components (it is b − a in formula (27)), and the two arbitrary points u 1 and u 2 selected for joining the outer and inner solutions.
Of course, the results must be insensitive to the choice of these values. Thanks to the fast evaluation speed, it is straightforward to establish that they indeed are by repeating the calculation after variation of these hyperparameters.
To give an example of appropriate values, for our illustrations we have used o v = 100 and p v = 20 for the V -solutions and o u = 25 with p u between one and three. Increasing the expansion orders by fifty percent did not change the results by more than one part in 10 10 indicating that these values represent a very safe choice.
The least trivial choice of the hyper-parameters is that of u 1,2 , but this is not because they would need to be finetuned in any way. They merely need to be far from the origin to make the asymptotic envelope representation (32) accurate. So we used u 1,2 of the order of several thousands to ensure accuracy of the U -function amplitude u 00 (W ) to a part in a million or better. By moving u 1,2 further, to u i ≈ 10 4 , the accuracy can be further increased to a relative level of 10 −10 − 10 −12 .
The distance between the join-points u 2 −u 1 was varied between unity and several hundred without a significant impact on accuracy (change in u 00 (W ) less than 10 −10 ). In order to demonstrate the robustness and the accuracy of the algorithm, Fig.5 depicts data obtained for the amplitude u 00 (W ) for the energy in a close neighborhood of the Stark resonance, which is located at S R = −0.500225560455... for F = 1/100. We choose this quantity here because it is the most sensitive as it varies over many orders of magnitude in an extremely narrow interval of energies as shown in the left panel of the figure. The panel on the right depicts the difference between u 00 -values obtained for two very different choices of u 1,2 , i.e. the locations used to join the outer and in-ner solutions. One can see that despite moving u 1,2 by a thousand atomic units, the resulting numerical variation in u 00 is extremely small. Even as the value changes by five orders of magnitude, the relative error (as estimated by the difference shown) remains on the level of 10 −6 . It should be emphasized that the robustness of the u 00 -value reflects the accuracy of the wave function everywhere as this quantity "connects" the asymptotic region to the origin.
An additional evidence for the fidelity of our algorithm is in the high accuracy of the Stark resonances. For the example shown in Fig. 5, we have estimated the resonance location S R by finding the point in the middle of the peak for higher and higher value(s) of u 00 . We have obtained an estimate which agrees with the value from the literature [26,27] to twelve digits.

VI. RESONANT CONTRIBUTION
We have seen that there is a potentially non-trivial W -dependent feature in the u n0 (W ) caused by the Stark resonance born from the originally stable ground state. Depending on the field it can be so narrow that it could remain unresolved, or even completely missed when W is sampled on a coarse grid. So a question arises how to deal with it, especially for applications in time-dependent problems. Here we sketch how the resonant contribution can be extracted and treated separately from the continuum background.
Our previous study of toy models in one spatial dimensions suggested the functional shape of u n0 (W ) in the vicinity of a resonance [28]. More specifically, every resonance causes a "step" in the phase shift δ, and in its neighborhood we can use where the phase is controlled by the Stark resonance located at S(F ) = S R (F )+iS I (F ). The effect of the square root ratio is a very fast phase change with W and this is because the imaginary part S I of S(F ) is very small (see Fig. 4). Incorporating this into the carrier-envelope representation of the wave function, we have For a sharp resonance (i.e. a weak field), the Airy combination Ci + is dominated by Bi, and we can further approximate where H(0) is evaluated for W = S R . Now we want to use an observation based on our numerical data. Namely, we have found that H(W = S R , u = 0) = iH I (0) is purely imaginary, which leads us to the following representation for the resonant part contribution to the probability density (in energy): . (40) Here, the last term tends to a delta function in a weak field when S I → 0. The fit shown as full line in Fig. 5 demonstrates that this functional shape indeed dominates the resonant portion of u 00 (W ). Symbols represent the calculated values, and the solid lines are fits based on the above expression. The fit was obtained with merely fifteen datums closest to the resonance. Nevertheless, the agreement with the data remains extremely good over the whole interval shown in the figure.
We have thus shown that in case the field F is weak the resonance contribution to the wave function can be accurately extracted with the help of the analytic structure built into the carrier-envelope representation. This can present an advantage in applications, as it offers the option to represent the spectral amplitude A n (W ) as consisting of a smooth background plus a resonant pole contribution.

VII. ILLUSTRATION: TUNNELING DYNAMICS
In general, using a continuum-energy basis to expand time-dependent wave functions is a highly nontrivial task from the numerical point of view. For instance, one of the challenges in methods utilizing discretized continuum basis sets, is that the wave function can only be calculated within a relatively small computational box surrounding the system. Our illustration aims to emphasize that with an accurate algorithm to evaluate all energy eigenstates, we are free of this problem and wave functions can be obtained even for very large distances from the origin.
We choose to consider an electron tunneling from the hydrogen atom, and calculate its time-dependent wave function. We assume that the initial state is the zerofield ground state wave function, and that the field is suddenly set to a constant value F . Motivated by its simplicity, this is an idealized scenario previously studied in exactly solvable one-dimensional models [28]. Here we investigate the same dynamics in a realistic threedimensional system.
The time-dependent solution can be expressed as a superposition of Hamiltonian eigenstates calculated for a fixed F , in which the energy-representation A n (W ) must be set such that the initial condition is satisfied. Using the resolution of unity in the energy space, Ψ nW |Ψ mE = δ mn δ(W − E), spectral amplitude A n (W ) is obtained as an overlap integral with the initial wave function, Figure 6 illustrates the behavior of A n (W ) for several lowest ns. It reveals that the sudden application of the field excited multiple continua of higher-energy states. One can see that the contribution of higher states diminishes quickly with their energy, yet the resulting spectrum is very broad and it implies very fast evolution which we will see shortly.
Because A n (W ) is nothing but the energyrepresentation of the initial-time wave function, |A n (W )| 2 gives the probability distribution that the given energy is excited. From Fig. 6 one can see that the excitation of energies above −0.4 is small, of the order of 10 −4 . This means that almost all of the wave function remains "concentrated" around the Stark-resonance peak, which is the feature in A 0 (W ) in the vicinity of the original ground-state energy. This peak is so narrow that it is impossible to resolve properly in Fig. 6. Figure 7 shows a logarithmic-scale view, and allows one to appreciate the presence of the resonance contribution. It also indicates that two very different time scales govern the evolution of the wave function after the field is applied.
The time-energy uncertainty relation suggests that the narrow resonance peak gives rise to a slowly evolving part of the wave function, and the wide peak results in a fast changing part of the wave function. The former remains localized and very similar to a bound state, only slightly deformed by the external field. The probability that the electron occupies this state can be obtained by fitting the resonance functional shape (40) to the |A 0 (W )| 2 finely Spectral amplitude A0(W ) as a function of the energy W exhibits a narrow, ground-state Stark resonance. Note that on the scale of this figure the resonance peaks remain un-resolved and are cut-off by the plot-range. In contrast, the higher-energy resonances are blurred as a consequence of the suddenly turned-on field. The sharp feature in the curve for F = 1/100 is an example of the resonance "remnant." sampled around W ∼ S R . This "amplitude" is the survival probability of the system's state after sudden exposure to the field. For our illustration case shown here, one obtains the survival probability of 99.9%. This is the part of the wave function which remains locked in the Stark resonance. It gives rise to an exceedingly slow but steady leakage of the probability density from the vicinity of the nucleus toward infinity. But the survival rate of 99.9% means that with probability of about 0.1%, the electron escapes from the atom, and it occurs very fast. This is a consequence of the non-adiabatic change in the field, and it manifests itself as a "pulse" in which the electron escapes from the atom. Figure 8 illustrates this in a "time-of-arrival" picture, where the evolution of the wave function is observed at a fixed point z = 250, x = y = 0. Because the localized portion of the wave function is exceedingly small at this distance, one can only see a pulse for each component n = 0, 1, 2, 3, . . . as it propagates away from the nucleus.
For a snapshot-view, Fig. 9 shows the n = 0 component depicted at time t = 60 after the external field was turned on. Figure 10 shows the n = 1 component which exhibits spatial transverse structure reflecting directional distribution of the tunneling particles. The higher n components of the escaping wave function have correspondingly richer structure.
These figures illustrate the part of the tunneling wave function which can be called non-adiabatic, because it is caused by the sudden turn on of the field. Given the distance from the nucleus, and the accuracy required, this would be extremely difficult to evaluate with standard numerical methods.
It is worth noting that our wave function, as shown in Fig. 9 contains both the eventually freed part and the part which is bound to the nucleus for a very long time The cut-off part in the center (the "stump") is the (slightly deformed by the field) wave function remaining in the metastable ground state which will "survive" for ∼ 10 10 atomic units of time. In contrast, the waveform will propagate away from the atom very quickly, and thus contributes to the non-adiabatic ionization. and we do not explicitly distinguish between the two in the formulas presented. However, at larger distances the part that is bound to the nucleus is exponentially small in comparison to the "freed" part, and what one obtains far from the nucleus is for all practical purposes the part of the wave function which represents the "freed" electron.
Our example is therefore relevant for the ongoing debate concerning the question of the tunneling time. More specifically, there is no universal agreement about the time an electron spends while tunneling from an atom under the influence of an external electric field (a brief The spatial structure indicates that one can't assign a single classical trajectory to this tunneling particle. discussion of different schools of thought can be found in [28]).
In this context, our example is the first exact solution in a three-dimensional system that suggests the tunneling time and velocity might not be zero. In this model there is no single well-defined tunnel exit, and therefore the notion of "tunneling time" is also ambiguous. As it can be seen from Fig. 9 and 10, we can also conclude that in this three-dimensional case we cannot ascribe a single classical trajectory to the tunneling portion of the wave function. The tunneling dynamics continues to be quantum in nature even beyond the emergence of the particle from the classically forbidden region. Consequently, proper quantum description should be preferred also in the outer regions, and our representation of the Stark-Hydrogen solutions represent an ideal tool for such an application.
For example, the same kind of calculations as above only at larger distances can be used to create data for the back-propagation method to find a family of classical trajectories (see e.g. Ref. [29] for various classifications of the exit and initial velocity).
In relation to the debate about the nature of the tunneling time, we have to admit that the present model does not yet completely reflect a typical experiment since such experiments are done using a circularly polarized (or "almost" circularly polarized) laser beam. [3,4] Moreover, our sudden excitation drives the system into a more nonadiabatic regime than an optical pulse with a smoothly varying intensity envelope would. Thus, it should be interesting to generalize our model to one which reflects the circular polarization of the driving field. The Stark-Hydrogen states could serve as a suitable point of departure to describe the ionization and tunneling dynamics in a reference frame co-rotating with the electric field.
Another currently relevant physics question concerns strong-field ionization in a long-wavelength (i.e. slowly evolving) optical field [30]. We have seen that there are two qualitatively different channels for ionization. Besides the well-known adiabatic tunneling ionization which is mediated by the Stark resonance [9,[31][32][33], there is also ionization current in the form of a short-duration "pulse." This is the reaction of the system to the increment of the external field, so it is non-adiabatic in nature. Post-adiabatic corrections are universally believed to become negligible for a very slowly evolving external fields, but accurate evaluations are not available at this time.
Having seen here how strong the non-adiabatic contribution is in our case, a detailed investigation should be of great interest. The technique used in our example can be generalized to evalaute the relative strength of the adiabatic and post-adiabatic ionization yields in slowly evolving external fields.

VIII. CONCLUSION
We have presented an accurate, approximation-free algorithm to calculate the continuum-energy eigenstates of the Stark-Hydrogen problem for a particle subject to Coulomb and a homogeneous electric fields.
Importantly for future applications, the algorithm presented in this work takes advantage of the analytic properties of the problem. In particular, the wave functions are represented with the help of a "carrier wave" which captures the tunneling portion of the wave function. This is then modified by an "envelope," or slowly changing function, which is relatively easy to find numerically, or analytically in the form of an asymptotic expansion. This representation of the wave functions can be useful in constructing quasi-classical solutions for electrons escaping from single-charged quantum systems.
Complementary to the carrier-envelope representation, we have also designed an algorithm suitable for small and intermediate (up to a few thousands of atomic units) distances from the nucleus. This approach is based on a series expansion "analytically continued" away from the origin, and it can be used to study the properties and eventually the temporal dynamics of the electrons driven by strong long-wavelength optical fields.
We have illustrated the usage of the continuum-energy basis to expand a time-dependent wave function of a tunneling particle after an excitation due to a sudden turn-on of an external field. Application to a general case with time-dependent F (t) requires evaluation of the spectral amplitude A nW (t) which, too, depends on time. This generalization will be addressed elsewhere.
To the best of our knowledge, ours is the first practically usable state-expansion method utilizing a continuous energy basis for a realistic three-dimensional quantum system. Given that the superposition of the Coulomb and homogeneous fields appears in many situations, we trust that the capability to treat these nontrivial quantum states without approximations will prove useful in various situations, including for example strong-field ionization of atoms and molecules, and highharmonic generation.

Acknowledgments
This material is based upon work supported by the Air Force Office of Scientific Research under award number FA9550-18-1-0183.

IX. APPENDIX: NORMALIZATION
The normalization condition we seek to satisfy reads and expands into two contributions which we will refer to as u-and v-terms: together giving rise to the delta function in energy like so Let us look first at the integrals over v. Because the V -functions are integrable, they can be normalized to unity. It can also be shown in a standard way that they are orthogonal to each other, so we can assume For the other integral, it turns out that the non-diagonal part will not be needed, and one can show that the diagonal portion can be simplified as where prime denotes partial derivative with respect to W . To obtain this result, one takes the "expectation value" equation for h v , differentiates on both sides w.r.t. W , and uses the fact that the normalization is fixed as in (48). Now we turn to the integrals over u. With the normalization of V accounted for, the I u term turns into The relevant contribution to this integral is from the region of large u, where we use the asymptotic behavior of the U -functions, where H is the envelope which behaves as H ∼ u −1/2 at infinity and thus cancels u in the integrand. The integral over large values of u tends to (we leave out index n in the integrand) Eliminating terms which oscillate even for E = W and do not contribute to the delta-function one gets Using the leading term in the asymptotic expansion of the Airy functions, one obtains Recall that the phase-shift terms δ(., F ) belong to the same m = n, and because they cancel for W → E they can be dropped. After substitution k = u 1/2 F −1/2 and representing the complex conjugate term as an integral over −k, we obtain (56) or Next we turn to the v-term given by the integral We will see that only the diagonal portion of V (1) mn (W, E) plays a role, because the second term turns out propor-tional to δ mn . It requires to evaluate which is a distribution because we deal with a continuous spectrum. Next we show that for m = n the above equals zero in the distribution sense. We start by rewriting it equivalently like so One can integrate by parts twice in the first term on the right-hand side to get where the "surface term" reads and the "volume contribution" is This is the point where one needs to use the specific properties of the U -functions. When W = E both terms on the right-hand size oscillate themselves to a distribution-zero. As a consequence, I m =n (W, E) is zero in the distribution sense, meaning that I m =n acting on an arbitrary smooth test function T (E) yields zero: However, this argument does not hold for m = n and W → E. We have to investigate such a case more closely, starting with the "volume" term. We know from (57) that the integral part is proportional to the Dirac-delta in energy, so the prefactor can be replaced by its limit for W → E, As a consequence of (49) and because z ′ v = −z ′ u , the contribution to I v simplifies to which, interestingly, exactly cancels I u contribution to the normalization integral. The surviving contribution is therefore given by the "surface terms" V The relevant behavior is of W ∼ E, so expanding accordingly we get Once again, the derivative of the separation constant z ′ u gets canceled by the integral over v and the contribution to the whole I v is: In summary, the normalization integral gives three contributions which are equal in their absolute values while one of them is negative. In effect the u-term remains, while the v-term vanishes (as a distribution). As a result the normalization to the Dirac-delta in energy requires to set N U such that or With the normalization factor set like this, it can be shown that the total probability outflow in the outgoing part of the wave function is equal to 1/2π independently of the energy, as it is expected to be [23] ( §21.)