{'fraction_non_alphanumeric': 0.04137274844975229, 'fraction_numerical': 0.018898257816857507, 'mean_word_length': 4.302713987473904, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fictitious time wave packet dynamics: I. Nondispersive wave packets in the quantum Coulomb problem

Nondispersive wave packets in a fictitious time variable are calculated analytically for the field-free hydrogen atom. As is well known by means of the Kustaanheimo-Stiefel transformation the Coulomb problem can be converted into that of a four-dimensional harmonic oscillator, subject to a constraint. This regularization makes use of a fictitious time variable, but arbitrary Gaussian wave packets in that time variable in general violate that constraint. The set of"restricted Gaussian wave packets"consistent with the constraint is constructed and shown to provide a complete basis for the expansion of states in the original three-dimensional coordinate space. Using that expansion arbitrary localized Gaussian wave packets of the hydrogen atom can be propagated analytically, and exhibit a nondispersive periodic behavior as functions of the fictitious time. Restricted wave packets with and without well defined angular momentum quantum n umbers are constructed. They will be used as trial functions in time-dependent variational computations for the hydrogen atom in static external fields in the subsequent paper [T. Fab\v{c}i\v{c} et al., submitted].


I. INTRODUCTION
Wave packets play an important role in the description of atoms, e.g., for the understanding of ionization processes in microwave experiments [1,2] or in experiments with short laser pulses [3,4]. Theoretically, the wave packet propagation can be calculated by exact quantum computations [5] or approximately with, e.g., semiclassical [6] or variational [7] techniques. Contrary to the harmonic oscillator, where Gaussian wave packets and coherent states [8] can easily be described analytically, the evolution of arbitrary Rydberg wave packets is nontrivial already in the pure Coulomb problem where spreading and revival phenomena are observed in the long-time computer simulation of a quantum wave packet [9]. Coherent states for the hydrogen atom have been constructed by Klauder [10] and by Majumdar and Sharatchandra [11], however, Bellomo and Stroud [12,13] have shown that these states do not move quasiclassically but spread rapidly over the Keplerian orbit. Dispersion is a general property of Rydberg wave packets with the exception of nondispersive electronic wave packets existing in periodically driven atoms such as the hydrogen atom in microwave fields [14,15].
The success of applying variational methods to wave packet propagation crucially depends on the choice of the trial function. Gaussian wave packets (GWPs) are certainly well suited for smooth and nearly harmonic potentials [16,17]. The Coulomb potential is not a promising candidate for successfully propagating GWPs directly. Nevertheless, the GWP method based on the local harmonic approximation has been applied in one dimension to the singular Coulomb potential [18,19,20]. In the three-dimensional space a regularization in Kustaanheimo-Stiefel (KS) coordinates [21,22] originally introduced for the Kepler problem in classical celestial mechanics, but also adapted to the hydrogen atom [23], transforms the Coulomb potential to a harmonic potential with a constraint. In the regularized hydrogen atom the application of the GWP method should therefore be capable of yielding exact results when the constraint can be handled. The regularization implies a fictitious time variable which has been shown to be the eccentric anomaly of the corresponding classical orbit [24]. Various approaches have been made to construct coherent states for the hydrogen atom in the fictitious time [25,26,27,28,29,30,31] in analogy with the coherent states of the harmonic oscillator. These approaches construct the coherent states as the eigenstates of the lowering operators associated with the harmonic potential.
In this paper we consider the field-free hydrogen atom and show that contrary to the dynamics in the real physical time the exact propagation of arbitrary initial Gaussian wave packets in the fictitious time can be described analytically and exhibit a nondispersive periodic time dependence. In Sec. II the Coulomb problem is transformed to the problem of the four-dimensional (4D) harmonic oscillator in Kustaanheimo-Stiefel coordinates subject to a constraint. The consequences of the constraint for Gaussian wave packets are discussed in Sec. III and the physically allowed set of "restricted Gaussian wave packets" is constructed. In Sec. IV the analytical time evolution is derived for initially three-dimensional (3D) Gaussian wave packets in the physical space, and also for twodimensional (2D) and one-dimensional (1D) wave packets with cylindrical and spherical symmetry, respectively.
In the subsequent paper [32] the investigations are extended from the pure Coulomb problem to the hydrogen atom in static external electric and magnetic fields. Wave packets are propagated by application of the timedependent variational principle in such a way that the dynamics is exact for the Coulomb problem and approximations in the variational approach are only induced by the external fields. Quantum spectra of the nonintegrable systems are then obtained by the frequency analysis of the time autocorrelation function of the propagated wave function.

II. REGULARIZATION OF THE HYDROGEN ATOM
To make our presentation self-contained we briefly review the Kustaanheimo-Stiefel transformation for the Coulomb problem. The time-independent Schrödinger equation for the hydrogen atom reads with ∆ 3 the Cartesian form of the Laplacian. A regularization of the singular Coulomb potential is obtained by using Kustaanheimo-Stiefel (KS) coordinates u = (u 1 , u 2 , u 3 , u 4 ) [21,22] which are introduced here, according to Ref. [25], differing by a factor of two from the original definition, By adding a fourth component with the constant value zero to the physical position vector, i.e. x = (x, y, z, 0), the transformation can be written in matrix notation, with The introduction of the auxiliary degree of freedom, which renders the originally three-dimensional problem four-dimensional, entails a constraint on physically allowed wave functions ψ, i.e.
where ∆ 4 denotes the 4D Cartesian form of the Laplacian. Multiplication with u 2 and reordering of the terms yields Eq. (7) is not a standard linear eigenvalue problem. Scaling the coordinates and setting leads to the time-independent Schrödinger equation which represents the Schrödinger equation of the 4D harmonic oscillator subject to the constraint (5). The scaling parameter, which takes only integer values n = 1, 2, 3, . . . , turns out to be the principal quantum number of the hydrogen atom. The 4D isotropic harmonic oscillator is invariant under the unitary group U (4) and thus the eigenstates of Eq. (10) are not unique. The simple product of four eigenstates of the 1D harmonic oscillators, i.e., the separation of Eq. (10) in the Cartesian coordinates (u 1 , u 2 , u 3 , u 4 ), in general violates the constraint (5) and therefore represents unphysical solutions. The constraint (5) rather suggests the introduction of two sets of polar coordinates, viz. the semiparabolic coordinates with the associated angular momenta (p j = 1 i ∂ ∂uj ) Eq. (12) yields the constraint (5) in the form The relation between the physical Cartesian coordinates and the semiparabolic coordinates is obtained using the definitions (2) and (11) x = µν(cos ϕ µ cos ϕ ν − sin ϕ µ sin ϕ ν ) = µν cos(ϕ µ + ϕ ν ) = µν cos ϕ , y = µν(cos ϕ µ sin ϕ ν + sin ϕ µ cos ϕ ν ) = µν sin(ϕ µ + ϕ ν ) = µν sin ϕ , with the physical azimuthal angle ϕ = ϕ µ + ϕ ν . In the semiparabolic coordinates the Schrödinger equation reads with Eq. (15) is separated in two uncoupled 2D harmonic oscillators in the coordinates µ, ϕ µ and ν, ϕ ν , respectively. The solution can be taken in the product form where − 1 2ρ with ρ = µ, ν and N ρ = 0, 1, 2, . . . . The coordinate representation of the eigenstates is with the associated Laguerre polynomials L |m| Nρ . For the principal quantum number n introduced above we obtain the relation and therefore via Eq. (9) the correct Rydberg spectrum.
To perform time-dependent computations it is necessary to formulate the time-dependent version of the Schrödinger equation (10). By analogy with the usual identification E → i ∂ ∂t where t is the physical time, the "fictitious time" variable τ is introduced, (see, e.g., Refs. [25,27]) as the conjugate variable to the principal quantum number n, The regularized Schrödinger equation for the hydrogen atom in the fictitious time then reads In the following the "fictitious time" τ will simply be denoted as "time" for brevity, whereas t will be named "physical time". The form of Eq. (22) suggests that it can simply be solved by Gaussian wave packets in Kustaanheimo-Stiefel coordinates u, i.e.
where A is a complex symmetric 4 × 4 matrix with positive definite imaginary part and the momentum π and the center q are real, 4D vectors in the Kustaanheimo-Stiefel coordinates. Those vectors represent the expectation values of the position and the momentum operator, respectively, i.e. q = g|u|g / g|g and π = g| 1 i ∇ 4 |g / g|g . The phase and normalization is given by the complex scalar γ. Collectively the parameters are denoted by y = (A, π, q, γ), which is a set of 4 × (4 + 1)/2 + 4 + 1 = 15 complex parameters when the two real vectors q and π are counted as a single complex vector. Inserting (23) into the Schrödinger equation (22) yields a set of ordinary differential equations for the parameters y, which can be solved analytically for any given initial GWP. As the potential is harmonic the timedependent GWPs are exact solutions of Eq. (22).
The problem, however, is that in Kustaanheimo-Stiefel coordinates the constraint (5) on physically allowed wave functions must be taken into account. The GWPs (23) in general violate that constraint. The question is whether the constraint can be fulfilled exactly by a Gaussian at all, and, if so, whether GWPs fulfilling the constraint are reasonable trial functions in the sense that they still present a complete basis set.
The advantage of the formulation of the Hamiltonian in semiparabolic coordinates (15) is that the constraint (5) is already incorporated. However, in semiparabolic coordinates it is nontrivial to find Gaussian type trial functions for the exact solution of the time-dependent Schrödinger equation because the Laplacian is not of Cartesian form, see Eq. (16), and includes centrifugal barriers. In the following we will therefore use the Cartesian type KS coordinates rather than semiparabolic coordinates to investigate the impact of the constraint (5) on a Gaussian trial function and to discuss the properties of the resulting "restricted Gaussian wave packets".

III. RESTRICTED GAUSSIAN WAVE PACKETS
The regularization of the hydrogen atom in Sec. II has transformed the Coulomb potential to a harmonic potential in the Schrödinger equation (10). The goal now is to perform exact wave packet propagation in the fictitious time τ for the hydrogen atom. The constraint can be fulfilled by a 4D GWP if the space of admissible configurations of the parameters y is restricted. The structure matrix allows for the compact notation of the constraint (5), Letting the constraint operator X act on the trial function (23) yields Thus the result is a quadratic polynomial in the coordinates u, multiplied by the GWP itself. The constraint (26) has to be satisfied pointwise for all u ∈ R 4 . For nontrivial wave packets the polynomial in (26) must vanish, and an algebraic equation remains, which is only possible if all coefficients of the second order polynomial in (27) are zero. Let us first investigate the term linear in u whose coefficients must vanish, i.e. π = 2Aq. Inserting this condition into the wave func- which means that without loss of the variational freedom we can set because any nonzero q and π vectors only change the scalar γ. The bilinear form in (27) must also be zero, i.e. 2uJ T Au = 0. This requires the matrix of the bilinear form to be skew-symmetric, i.e., the diagonal elements of J T A must vanish, a 12 = a 34 = 0, and from the off diagonal elements we obtain a 11 = a 22 , a 33 = a 44 , a 24 = −a 13 , and a 23 = a 14 . With the definitions a 11 = a µ , a 33 = a ν , a 13 = a x , a 14 = a y , the matrix A must be of the form Eqs. (28) and (29) imply that the set of 15 complex parameters y of the general GWP (23) is reduced to only five parameters y = (a µ , a ν , a x , a y , γ) for a restricted GWP satisfying the constraint (5)! The question arises whether the restricted GWPs form a complete basis set, such that any physically allowed state can be expanded in this basis. It is not evident that a superposition of restricted GWPs whose centers are all located at the origin and which only differ by their complex widths, is flexible enough to represent arbitrary quantum states. The usual form of the resolution of the identity [33] for a continuous basis set of normalized, unrestricted GWPs of the form (23) is where the width of each GWP basis state is kept fixed. Obviously, Eq. (30) cannot be applied to the restricted GWPs, since both parameters π and q are set to zero (see Eq. (28)). However, it is sufficient to require the restricted GWPs to be complete in the 3D physical space only. To verify the completeness in the 3D space we transform the restricted GWP in KS coordinates back into the original 3D Cartesian coordinates, where we have exploited that in semiparabolic coordinates µ 2 = r + z, ν 2 = r − z. In Eq. (31e) the set of parameters (a µ , a ν , a x , a y ) has been replaced by an equivalent set of complex parameters, defined by For p r = 0 and real-valued parameters p x , p y , p z the restricted GWP in Cartesian coordinates (31e) reduces to a plane wave e ip·x . (Note that in that case the imaginary part of the matrix A is not positive definite and the wave function cannot be normalized.) Since plane waves are known to form a complete basis it is proved that the restricted GWPs, forming a superset of plane waves, are also complete, or even over-complete. However, they do not form a complete basis set of the 4D harmonic oscillator (10). The GWPs (31) satisfy the constraint (5), but "for the price" of the condition (28), i.e., they are localized around the origin with zero mean velocity. It seems impossible that the time propagation of a single restricted GWP exhibits any meaningful dynamics with, in particular, a classical limit in the sense of the correspondence principle. However, a wave packet localized in the physical coordinate and momentum space of the hydrogen atom can be constructed as a superposition of the restricted GWPs. The expansion and exact time evolution of wave functions in the basis (31), and in addition, in modified bases sets for certain symmetry subspaces of the hydrogen atom are the subjects of the next section.

IV. ANALYTICAL WAVE PACKET DYNAMICS IN THE HYDROGEN ATOM
An arbitrary wave packet of the hydrogen atom can be propagated analytically in the fictitious time. This is achieved by expanding the initial wave packet in terms of the restricted GWPs, whose time-dependence is derived and shown to be given by simple analytical formulae. We consider three different cases. In Sec. IV A the time propagation of Gaussian wave packets in the 3D physical space is discussed. These wave packets are not eigenstates of the angular momentum operator. We then introduce the time propagation of wave packets with symmetries, viz. in Sec. IV B cylindrically symmetric 2D wave packets with well defined magnetic quantum number m, and in Sec. IV C spherically symmetric 1D wave packets with well defined angular momentum quantum numbers l and m.

A. Propagation of 3D Gaussian wave packets
The aim is an exact time propagation of arbitrary wave functions in the hydrogen atom. The initial wave function ψ(0) is taken as a superposition of the restricted GWPs (31). These basis states are then propagated analytically in time. We start with the derivation of the time evolution of the basis states. Inserting the ansatz (31a) in the time-dependent Schrödinger equation (22) yields The equations of motion for the 4 × 4 width matrix A given in Eq. (29) and the complex phase factor γ are now directly obtained from Eq. (34) aṡ It is important to note that although the parameter set of the GWP (31a) is restricted to only five complex parameters (as compared to 15 parameters in (23)) the differential equations (34) still describe the exact dynamics of the wave packet without any approximation, i.e., with the initial condition (28) the wave packet stays centered around the origin for all times, and the width matrix A keeps the form (29) because the matrixȦ in (34a) has the same structure (29) as A itself. Both equations (34) can be solved analytically. Equation (34a) is solved most easily when two auxiliary complex matrices B and C are introduced with A = 1 2 BC −1 , and the initial conditions B(0) = 2A(0) and C(0) = 1 [34]. Then (34a) is replaced with the two equationṡ C = B andḂ = −C, or equivalentlyB = −B. The matrices B and C have the same structure as the width matrix A (29), with the solution for the matrix B and for the matrix C Here and in the following the superscript 0 indicates parameters of the initial state at time τ = 0. The matrix A is obtained from the above definition A = 1 2 BC −1 and the four elements in the form of Eq. (32) read where Z(τ ) abbreviates the expression With the matrix A at hand it is possible to integrate Eq. (34b) to obtain γ. The quantity e −iγ , i.e., the phase and normalization of the wave function, then reads, with the initial value γ(0) = 0 The analytical time evolution of the wave function is obtained by inserting the time-dependent parameters in the wave function (31e) which finally yields Now that the time evolution of the basis states (31e) is known we can expand an arbitrary initial state in this basis. The time evolution of that wave function is then analytically given by the superposition of the timedependent restricted GWPs (40). In general, the expansion of an arbitrary state ψ(x) in an over-complete set of Gaussian wave functions g(y, x) is a nontrivial task. A procedure for finding the "optimal" expansion for a given number N of basis states is to minimize the deviation ∆ = ||ψ(x) − N k=1 g(y k , x)|| 2 , e.g., by searching for stationary points ∂∆ ∂y k = 0 with respect to the parameters y k , k = 1, . . . , N [36]. This procedure presents a highly nonlinear minimization problem. Many stationary points, i.e., local minima may exist, and the difficulty is to find the true global minimum.
Here, we concentrate on the propagation of 3D Gaussian wave packets which are localized around a point x 0 with width σ in coordinate space and around p 0 in momentum space, and present a direct approach for the expansion of the 3D GWPs given as Using the approximation r ≈ x · x 0 /|x 0 |, which is valid in the vicinity of x 0 where ψ(x) is localized Eq. (42) can be written as where the g(y, x) are the restricted GWPs (31e) for the set of parameters y given by (p r = const., p − p r x 0 /|x 0 |, γ = 0). As can be easily shown the restricted GWPs (31e) for constant p r and γ = 0 are a complete continuous basis as functions of the momentum p, i.e.
Note that the completeness relation (30) is valid for frozen Gaussians with arbitrary localization in coordinate and momentum space while in Eq. (44) the restricted GWPs in KS coordinates are all located around the origin but the width matrix A is varied. Note also that Eq. (44) is exact for p r = 0, i.e., for the restricted GWPs given as plane waves, and becomes an approximation for p r = 0. This means that p r should be chosen as zero or close to zero in practical applications. In numerical computations it is convenient to approximate the 3D Gaussian wave packet (41) by a finite number of restricted GWPs rather than using the integral representation (44). This is most efficiently achieved by evaluating the integral in (44) with a Monte Carlo method using importance sampling of the momenta with a normalized Gaussian weight function The initial wave packet then reads with y k = (p r , p k − p r x 0 /|x 0 |, 0), and the p k , k = 1, . . . , N distributed randomly according to the normalized Gaussian weight function (45). The wave function ψ(x) in Eq. (46) is an approximation to the 3D Gaussian wave packet (41), and the accuracy depends on how many restricted GWPs are included. However, it is important to note that a wave packet which is strongly localized around x 0 in coordinate space and around p 0 in momentum space can be described even with a rather low number N (of the order of 50-100) restricted GWPs. The time propagation of an initial state (46) in the fictitious time τ is now obtained exactly and fully analytically in a simple way. In the unperturbed Coulomb problem a superposition of time-dependent restricted GWPs is uncoupled, i.e., all basis functions propagate independently. Thus, the time propagation of the initial state (46) simply reads with the same parameters y k as above and the time propagation of the restricted GWPs given by Eq. (40). Eq. (47) is the final result of this section, and is illustrated in Fig. 1 for an initially Gaussian wave packet (41) whose center is moving in the z = 0 plane. The wave packet is expanded and propagated analytically in a basis of N = 10000 restricted GWPs. This is possible since the analytical approach allows for a large number of basis states. We note, however, that results of similar quality can be obtained using N = 50 ∼ 100 GWPs only. The probability density in the z = 0 plane is plotted at equidistant times with step size ∆τ = π/5. As mentioned above Eq. (44) is exact for p r = 0, however, in our numerical calculations we choose p r = iǫ with a small ǫ > 0. The damping enables normalization of the restricted GWPs, and improves the convergence of the Monte Carlo integral. The initial GWP presented in Fig.  1 for τ = 0 is centered at x 0 = (8, 0, 0) with the mean momentum p 0 = (1, 2, 0), and p r is set to p r = iǫ with ǫ = 0.01. Classically the electron with these initial conditions is running on the Kepler ellipse plotted by dots on the bottom of each panel in Fig. 1. For every time step, that part of the ellipse that has been passed by the electron so far is shown by a black solid line for time resolved comparison. The position of the maximum of the probability density agrees well with the classical position of the electron on the ellipse for all times. The π-periodicity of the motion is reflected by the coincidence of the wave packet after one period at τ = π with the initial GWP at τ = 0.

B. Propagation of 2D cylindrically symmetric wave packets
In this section basis functions based on the restricted GWP (31c) with a well defined angular momentum component l z = m are derived. This case is especially important when a cylindrically symmetric external field, e.g. a magnetic field, is applied to the hydrogen atom, as discussed in the following paper [32]. First wave packets with definite l z are constructed and their exact, analytical dynamics in the hydrogen atom is discussed. Then we introduce a procedure to expand quantum states of defined l z in terms of the basis states.
The solution of Eq. (52b) is the solution of Eq. (34b) multiplied by the factor (1+|m|), and the phase and normalization factor of the wave function reads, with γ(0) = 0 The time evolution of the wave packet (49) then is given by The time-dependent basis states (56) are the analogue of the restricted GWPs (40) for the propagation of wave packets with constant magnetic quantum number m. They obey the constraint (5), but they are not sufficiently general to describe the dynamics of 2D wave packets localized around a given point (ξ 0 , η 0 , p ξ0 , p η0 ) in the parabolic coordinate phase space. Such localized states can now be constructed in a similar way as described in Sec. IV A. We use the formal plane wave expansion of a Gaussian wave packet in parabolic coordinates, viz.
In Fig. 2 the expansion of a GWP (57) is shown with the center at ξ 0 = η 0 = 25.0 and the momenta p ξ0 = 0.535 and p η0 = −0.117 or in terms of cylindrical coordinates ρ 0 = 25.0, z 0 = 0.0 and p ρ0 = 0.419, p z0 = 0.652 and the width σ = 4.472. Note that a Gaussian shape of a wave packet in parabolic coordinates is nearly Gaussian also in cylindrical coordinates (see e.g. the wave packet at τ = 0 in Fig. 2). The wave function shown has zero angular momentum component l z = 0. For reasons of presentation the originally positive radial coordinate ρ is extended to negative values, and the symmetry ψ(−ρ) = ψ(ρ) is used. The probability density ρ|ψ(ρ, z)| 2 is plotted in the (ρ, z) plane. The Kepler ellipses plotted on the bottom in each panel of Fig. 2 show the corresponding classical motion of the particle with the initial conditions given above. The second ellipse again is obtained by reflection symmetry as the intersection of the torus, which is obtained from rotating the ellipse around the z axis. At times τ ≈ 2π/5 and τ ≈ 4π/5 the high probability density close to the z axis leads to interference patterns. After one period τ = π the initial wave function τ = 0 is recovered. A number of N = 5000 modified basis states (49) with ǫ = 0.05 are employed. Results for the case m = 0 are not shown since they differ only qualitatively by avoiding the crossing of the z axis due to the rotational barrier.

C. Propagation of 1D spherically symmetric wave packets
The procedure of the two previous subsections is applied to quantum states with conserved angular momentum. First an extension of the basis states (31e) to basis states with well defined angular momentum quantum numbers lm is presented, and they are shown to be exact solutions of the time-dependent Schrödinger equation of the regularized hydrogen atom. Then the procedure of expanding states with definite lm in the constructed basis states together with an example are presented. For radial symmetry the complex width matrix A (29) of the restricted GWP must be a multiple of the identity matrix A = a1, i.e. a x = a y = 0, a µ = a ν ≡ a = p r /2. The restricted GWP reduces to g 00 (r) = e i(2ar+γ) = e i(pr r+γ) .
Wave packets with well defined angular momentum quantum numbers l and m can now be expanded in the basis (62). For the radial coordinate r the same procedure as introduced in Sec. IV B for the parabolic coordinates ξ and η can be applied. The plane wave expansion of a Gaussian wave packet localized around (r 0 , p r0 ) reads = Aσ √ π dp r e −σ 2 (pr +iǫ−pr 0 ) 2 +ǫr0−ipr r0 e i(pr+iǫ)r = Aσ √ π dp r e −σ 2 (pr +iǫ−pr 0 ) 2 +ǫr0−ipr r0 g 00 (y, r) , where g 00 (y, r) is the spherically symmetric restricted GWP (61) for the set of parameters y = (p r + iǫ, γ = 0). An initial state with given angular momentum quantum numbers l and m is obtained as Using the Monte Carlo evaluation of the integral in (69) and the replacement of the initial basis states (61) with the corresponding time-dependent solutions g 00 (τ, r) = exp{i[2a(τ )r + γ(τ )]} (with the time-dependent parameters given in Eqs. (66) and (68)) we finally obtain with sampling points p k r randomly distributed around p r0 according to the weight function w(p) = (σ/ √ π)e −σ 2 (p−p0) 2 , and the parameter sets y k = (p k r + iǫ, γ = 0).
The results of the propagation of the wave function ψ lm (r, θ, ϕ) = r l ψ 00 (r)Y lm (θ, ϕ) with ψ 00 (r) given in Eq. (69) and the initial values r 0 = 10, p r0 = −0.5 and the width σ = 3 is presented in Fig. 3 for different times 0 ≤ τ ≤ π. The imaginary parts of p k r = 2a k are set to ǫ = 0.2 and the number of basis states (62) is N = 10000. In Fig. 3(a) the angular momentum is set to m = l = 0 and in Fig. 3(b) the components of the angular momen-tum are l = 5 and m = 0. Due to the negative initial value of the radial momentum p r0 the wave is initially running towards the nucleus located at the origin. The wave function with zero angular momentum in Fig. 3(a) comes close to the origin r = 0. Similar to the radial symmetric case in Sec. IV B there appears to occur some interference pattern due to the overlapping parts of the incoming wave function at the inner turning point (see τ = 2π/5). In the nonvanishing angular momentum case ( Fig. 3(b)) the barrier of rotational energy prevents the wave function from reaching the nucleus. Instead there is a turning point whose distance from the nucleus increases with growing angular momentum. Again at τ = 2π/5 an interference pattern is observed close to the inner turning point. In both panels the maximum of the probability density overshoots the position of the initial maximum r 0 = 10 at τ = 4π/5 due to the initial kinetic energy and returns to the initial wave packet after the period τ = π, indicating the periodicity of the wave function.

V. CONCLUSION
In this paper we have derived the wave packet dynamics for the field-free hydrogen atom in a fictitious time variable. The Coulomb problem has been transformed to the four-dimensional harmonic oscillator in Kustaanheimo-Stiefel coordinates with a constraint. The "restricted Gaussian wave packets" obeying that condition have been constructed and their exact time dependence is calculated analytically. The wave packets with and without symmetries exhibit a nondispersive periodic behavior in the fictitious time.
It should be noted that the wave packet propagation in the fictitious time substantially differs from the physical time dynamics, and thus cannot provide the analytical propagation of dispersive wave packets in the physical time [18,19]. Nevertheless, the fictitious time dynamics can be used to solve the Schrödinger equation for Coulomb systems with strong time-independent perturbations, e.g., the hydrogen atom in static external electric and magnetic fields. The restricted Gaussian wave packets are the basis for the application of the time-dependent variational principle to the hydrogen atom in external fields and the computation of quantum spectra by frequency analysis of the time autocorrelation function in the following paper [32]. As a consequence of using the fictitious time variable the method is exact for the fieldfree hydrogen atom and approximations in the variational approach are only induced by the external fields.