Open Access
Add to BibSonomy
Abstract
We show that, for the case of resonant media, the available models for unidirectional propagation of short pulses can face serious challenges with respect to numerical efficiency, accuracy, or numerical artifacts. We propose an alternative approach based on a propagator operator defined in the time domain. This approach enables precise simulations using short time windows even for resonant media and facilitates coupling of the propagation equation with first-principle methods such as the time-dependent Schödinger equation. Additionally, we develop a numerically efficient recipe to construct and apply such a propagator operator.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Short optical pulses provide a universal and important tool for investigation of light-matter interaction as well as of ultrafast physical and chemical processes. Techniques such as time-resolved spectroscopy (see e.g. [1–3]), high harmonic generation [4,5], ultrafast material modification [6], studies of time-resolved molecular dynamics [7], etc. rely on the generation and propagation of ultrashort pulses with sub-100-fs durations and below.
Accurate modelling of pulse propagation in various materials is an essential component of understanding and optimization of light interaction with matter. Numerous approaches were developed for simulation of pulse propagation, including methods suitable for modelling propagation in inhomogeneous media. These methods, such as the finite-difference time-domain approach [8], provide accurate results but are associated with a significant numerical effort. Therefore we concentrate here on methods which ignore back-reflection and treat propagation as unidirectional; such methods are both numerically efficient and accurate in many practical cases with no sharp interfaces and/or inhomogeneities, where back-reflection can be neglected.
For the linear propagation of pulses, the group velocity dispersion (GVD) is, together with loss, the key effect, which is determined by the frequency-dependent wavenumber. Historically, the first approach to modelling the GVD is to represent the wavenumber in a Taylor series as a function of frequency [9]. While suitable for relatively long (100 fs and above) pulses, this approach fails for broad pulses and/or media with complicated wavenumber behaviour, which cannot be sensibly represented in form of Taylor series, such as resonant media. Therefore, as the accessible pulse duration grew shorter, the paradigm has shifted towards including the frequency-dependent wavenumber in an exact way, using equations such as the Forward Maxwell equation [10], the unidirectional propagation equation [11], etc. For a detailed overview and comparison of the available methods and equations, see [12].
In the recent years, the need to accurately simulate even shorter pulses (10 fs and well below) has grown considerably, stimulated by the success in generation of such pulses, e.g. [13,14]. At the same time, simulation of materials with resonant response has become an important topic. This is in part motivated by the development of first-principle methods such as time-dependent Schrödinger equation (TDSE) solvers and their combination with Maxwell equations [15,16]. The polarization yielded by the TDSE for resonant media and sufficiently broad pulses features oscillations (so-called "ringing") after the pulse, which originate from coherent population of the excited states. Extending the numerical time window to include such ringing is not practical due to the absence of decay in the standard TDSE and the corresponding theoretically infinite decay time of ringing. Surprisingly, as will be shown below, there is currently no efficient, accurate, artifact-free, and causal method to simulate unidirectional short pulse propagation in a resonant medium – even in the linear regime.
In this paper we provide such a method in form of the propagator operator approach. We first review the available approaches to identify their issues and drawbacks. After that, we propose an alternative approach based on a pulse propagator operator in the time domain. We derive efficient numerical recipes for both the construction of such an operator and its application for pulse propagation, and discuss possible extensions of our method.
2. Issues of the available methods
In the framework of the unidirectional approach, light propagation is described by the following equation:
where $E(z,t)$ is the electric field, $z$ is the propagation coordinate, and $P(z,t)$ is the polarization. Here we do not use the slowly-varying envelope approximation, so that $E(z,t)$ is the real-valued field, and consider in this section only the linear part of the polarization of a low-density medium. We note that the aim of this investigation of the linear terms is not the – relatively trivial – linear propagation per se; rather, we study here the linear-term numerical techniques which will augment the nonlinear propagation equation.The most general expression for the polarization is
where $G(\tau )$ is the Green function which relates the linear polarization and the electric field.In the frequency domain, the above relation can be expressed as $P(z,\omega )=\epsilon _0\chi (\omega )E(z,\omega )$, with $\chi (\omega )=G(\omega )/\epsilon _0$, and the field after propagation of a distance $dz$ is
In many cases modelling Eq. (3) is sufficient, since for an infinite time window it provides an exact description of the dispersion and loss. In particular for pulses far from resonances, it is enough to take the numerical time window only a few times longer than the pulse itself, since the polarization is temporally localized to the vicinity of the pulse. This is equivalent to the absence of sharp features in the $G(\omega )$ profile, which need to be resolved by using a small step in the frequency domain (or equivalently a long time window, since the step in the frequency domain $\delta \omega$ is generally related to the time window $T$ by $\delta \omega =2\pi /T$). Another possibility is to use a decomposition of the $G(\omega )$ in Pade-type terms [17]; this, however, does not solve the problem of the excessively long time window outlined below.
Sharp features in $G(\omega )$ do occur for many systems with relatively long relaxation times, such as gases at normal conditions, which can have relaxation times of up to nanoseconds. As a simple illustration, it is sufficient to consider a model one-dimensional atom. The 1D TDSE is solved for the soft-core Coulomb potential [18],
where $I_\textrm {at}=27.211$ eV, $x_\textrm {at}=0.529$ Å, and the value of $a$, $a=1.4142$, is chosen to reproduce the ground state of the hydrogen atom. The TDSE is propagated on a time and spatial grid, using a 5 point finite differences representation of the kinetic energy operator and the Crank-Nicolson scheme for the time propagation.In Fig. 1, we show the real and imaginary parts of the susceptibility $\chi (\omega )$ for the atomic hydrogen, which was calculated based on the energies and dipole momenta of the levels obtained numerically. The corresponding Green function $G(\tau )$ is given by
whereby we have multiplied the Green function obtained from TDSE by a relaxation term $e^{-\nu \tau }$. Here $N$ is the atom density, $\mu _{0k}$ are the dipole momenta between the ground state and the excited state $k$, and $E_0$ and $E_k$ are the corresponding energies. The relaxation was not included in the TDSE itself, as well as many-atom effects (due to a relatively low pressure). One can clearly see the presence of sharp resonances; their widths are estimated by the relaxation rate $\nu =\langle v\rangle N\pi (2n^*r)^2$, where $\langle v \rangle \sim 2500$ m/s is the average velocity of hydrogen atoms at room temperature, $N=2.69\times 10^{25}$ m$^{-3}$ is the particle density given by the Loschmidt number at the pressure of 1 bar, $n^*=2$ is the characteristic number of the excited level involved in the process, $r$ is the radius of the atom equal in this case to the Bohr radius of $0.529$ Å.
Fig. 1. Real (thick red curve) and imaginary (thin green curve) parts of the susceptibility for 1D atomic hydrogen as a function of frequency.
For the above conditions, we get $\nu =2\times 10^{11}$ s$^{-1}$ which corresponds to a relaxation time of 5 ps. To incorporate the decay of polarization, one would have to use a time window of roughly 20 ps or above, which is clearly unpractical and leads to unnecessary numerical effort for pulses in the range of 20 fs.
A possible alternative approach would be to use Eq. (3) in conjunction with the frequency-dependent $\chi (\omega )$, while keeping the time window short. The result is shown in Fig. 2, where the exact polarization of a delta-like pulse (green curve, calculated by using an extremely long time window) is compared to the result using a short time window of roughly 48 fs (red curve). One can see that the polarization starts to "wrap around" the computational time window, yielding non-causal values of polarization before the input pulse, which is placed at the center of the time window. This wrap-around originates from the large group delay of near-resonant components, and from the fact that group delays differing by integer numbers of $T$ are indistinguishable for the Fourier transform. In addition, due to multiple wrap-arounds, the values of polarization are much higher than the exact result. This behavior stems from the insufficiently large time window or, equivalently, from the insufficiently small step in the frequency domain compared to the width of the resonance lines in $G(\omega )$ given by $\nu$.
Fig. 2. Exact calculation of the polarization for a delta-like pulse in the center of the time window (thick green curve) compared to the calculation using Eq. (3) and a short time window (thin red curve).
Still another possibility is to solve the propagation equation directly in the time domain:
For infinitesimally small $dz$, this relation is exact. If we write
where the linear operator $\hat {A}$ is defined by an integrodifferential expressionThe solution based on Eq. (6) is then usually a result of expansion, assuming that $\hat {A}dz$ is small, and $\exp [\hat {A}dz] \approx 1+\hat {A}dz$.
Unfortunately, simulations according to Eq. (6) fail even for relatively short $dz$. As an example, we consider a step of 5 $\mu$m of propagation in atomic hydrogen. Note that typical experimental propagation lengths in gases (e.g. in hollow fibers) are several tens of centimeters or above, which would already yield $\sim$ 10000 grid points for the 5 $\mu$m step, meaning that the step cannot be easily decreased by several orders of magnitude. The comparison of the relative intensity after the propagation to the input intensity, for the delta-like input pulse, is given in Fig. 3. One can see that around resonances, the intensity increases strongly, by a factor significantly above unity, which is clearly an artifact, since an atom in the ground state cannot provide gain in the linear regime.
Fig. 3. Ratio of intensities after and before the propagation of 5 $\mu$m in atomic hydrogen as a function of frequency. A delta-like input pulse is considered.
3. Propagator operator
In view of the above, it would be useful to develop a numerical method which is free of such challenges. Fortunately, Eq. (10) offers a possibility to formally define the propagation operator
This operator connects the fields at positions $z$ and $z+dz$, and since operator $\hat {A}$ is a linear integral operator, see Eq. (9), operator $\hat {B}_{dz}$ is also a linear integral operator:
For a numerical scheme with a fixed $dz$ step, this operator needs to be calculated only once and then applied at each step.
The viability of this concept relies upon an accurate and efficient algorithm for calculating the function $H_{dz}$. We can start by finding $H_{\delta z}$ for some extremely small $\delta z$, for which one can still use the expansion $\exp [\hat {A}\delta z] \approx 1+\hat {A}\delta z$ with
After that, a naïve approach would be to subsequently multiply $\hat {B}_{\delta z}$ $dz/\delta z$ times: $\hat {B}_{2\delta z}=\hat {B}_{\delta z}\times \hat {B}_{\delta z}$, $\hat {B}_{3\delta z}=\hat {B}_{\delta z}\times \hat {B}_{2\delta z}$, $\hat {B}_{4\delta z}=\hat {B}_{\delta z}\times \hat {B}_{3\delta z}$ etc. However, our preliminary simulations for typical parameters show that there is a balance between the accuracy of the expansion of the $\exp [\hat {A}\delta z]$, which improves for small $\delta z$, and the accumulation of arithmetic errors, which accelerates with $dz/\delta z$. The optimum is reached for $dz/\delta z \sim 10^9$, which corresponds to a $\delta z$ on the order of several femtometers and makes the above sequential multiplication numerically impractical. Fortunately, if one chooses $dz/\delta z$ to be a power of 2, i.e. $dz/\delta z=2^k$, one can use a much faster approach:
This reduces the number of operator multiplications from $dz/\delta z$ to $\log _2dz/\delta z$, i.e. from 2$^{30}$=1073741824 to 30 for the above example. Note that, except for the arithmetic (rounding and truncation) errors, the procedure described by Eq. (14) is exact. However, when linear step is used in conjunction with a nonlinear step, as is the case in most simulations, the total error will be determined by the interaction of two terms and need to be analyzed on a case-by-case basis.
The square of operators is given by
which, if implemented straightforwardly, would require $N_t^2$ number of operations, where $N_t$ is the number of points in the time grid. However, here we also propose a faster approach. First, we map $H_{\delta z}$ from a grid with $N_t$ points to a grid with $2N_t$ points, padding the grid points from $N_t+1$ to $2N_t$ with zeros. Second, we perform the multiplication in the Fourier domain, using where $\hat {F}$ and $\hat {F}^{-1}$ are the forward and backward Fourier transforms, respectively. The obtained $H_{2\delta z}$ is then mapped back on the $N_t$-points grid by discarding the values from $N_t+1$ to $2N_t$. Since $\hat {B}_{\delta z}$ produces a maximum time delay of $T$, the mapping on the $2N_t$-points grid eliminates the wrap-around in the time domain and the associated non-causal artifacts.Similarly, the application of the operator $\hat {B}_{dz}$ can be implemented in the Fourier domain. First, we pad $H_{dz}$ and $E(z)$ by zeros to $2N_t$-point grid, then we calculate
followed by discarding the values of $E(z+dz)$ at time points from $N_t+1$ to $2N_t$.We summarize the proposed numerical recipe as follows: i) obtain the Green function $G(\tau )$ from the numerical simulation or from the data on the energy level structure, e.g. given by Eq. (5); ii) calculate $H_{\delta z}(t)$ using Eq. (13); iii) calculate $H_{\delta z}(t)$, $H_{2\delta z}(t)$, $H_{4\delta z}(t)$, …, $H_{dz}(t)$ using Eq. (16); iv) finally, use Eq. (17) to propagate the electric field in space.
Equations (11), (13), (14), and (16) constitute the core result of this work. They establish the formalism of the propagator operator and allow one to calculate the corresponding Green function accurately and efficiently, using $N_t\log (N_t)\log (dz/\delta z)$ operations instead of $N_t^2dz/\delta z$ operations. Also, the use of the propagator operator does not require larger numerical effort in comparison with the conventional Fourier-domain method given by Eq. (3), since both require $N_t\log (N_t)$ operations.
4. Example of simulation using the propagator operator
To illustrate the use and advantages of the propagator operator, we have simulated the propagation of a 50 TW/cm$^2$, 5-fs (FWHM) pulse with a central wavelength of 517 nm through 16 cm of 1-atmosphere atomic hydrogen. For consistency with the previous sections, we have used the linear dispersion as calculated by the 1D TDSE and shown in Fig. 1, in conjunction with a Kerr-type nonlinear term characterized by n$_2$=1.7x10$^{-7}$ cm$^2$/W as obtained from the TDSE. We investigate three approaches to the simulation of linear dispersion: using the conventional approach in the Fourier domain described by Eq. (3), the first-order expansion of the $\hat {B}$ operator as given by Eq. (6), and the propagator operator approach of Eq. (11).
The results are presented in Fig. 4. When the propagator operator is used, one can see a clear spectral broadening and the generation of the third harmonic in the spectrum [Fig. 4(b)] accompanied by loss around 3$\omega _0$. In the temporal domain [Fig. 4(a)] we see, as expected, a shift of the pulse in the coordinate frame moving with the velocity of $c$ accompanied by pulse elongation.
Fig. 4. Simulation of propagation using (a),(b) the propagator operator approach of Eq. (11), as well as (c),(d) the first-order expansion of the $\hat {B}$ operator as given by Eq. (6), and (e),(f) the conventional approach in the Fourier domain described by Eq. (3). In panels (a),(c),(e) the temporal profiles are shown, while in (b),(d),(f) the spectra are presented. The propagation distance is 1.6 cm (thick solid red curves), 4.8 cm (dashed green curves), and 16 cm (thin solid blue curves). The parameters are given in the text. In addition to the above, in (g) the comparison of the field temporal profiles after 16 cm for the conventional approach in the Fourier domain (thin red curve) and the propagator operator approach (thick green curve) is shown.
The case of the first-order expansion of the $\hat {B}$ operator, as given by Eq. (6), is shown in Fig. 4(c),(d). In the temporal domain the total energy increases, and also in the spectral domain one can see significant distortions of the spectrum in comparison to the exact case of the propagator operator. Note that we had to reduce the step in the $z$ coordinate by a factor of two to propagate up to 16 cm without numerical overflow; even with the reduced step, overflow still occurs with further propagation.
The results for the conventional approach in the Fourier domain, as described by Eq. (3), are presented in Fig. 4(e),(f) and (g). At first glance, the results in the temporal domain [Fig. 4(e)] are similar to those of the propagator operator method [Fig. 4(a)]. However, a closer examination of the temporal range before the pulse, shown in Fig. 4(g), shows that in contrast to the propagator-operator method, the conventional method results in a wrap-around in the temporal domain, in particular, of the spectral components around 3.6 fs$^{-1}$, which are close to a resonance and are characterized by a high velocity (delay) relative to the moving frame. This wrap-around results in the appearance of the nonphysical field before the pulse, and also has implications in the spectral domain: the maximum around 3.6 fs$^{-1}$, since these components are nonphysically preserved in the numerical domain, as well as a general deviation from the exact spectrum shape shown in Fig. 4(b).
5. Extensions for an arbitrary velocity of the coordinate frame, slowly-varying-envelope approach, and propagation accompanied by TDSE
For the propagation of light in a material with refractive index significantly different from 1, the coordinate frame moving with a velocity of light $c$ is not useful. The reason is that the group velocity strongly deviates from $c$ and a pulse would shift out of the time frame. Therefore, one typically uses the coordinate frame which moves with the group velocity of light at the central frequency $\omega _0$, and the propagation equation in the frequency domain becomes:
In order to construct the propagator operator in this case, for an extremely small propagation distance $\delta z$, we obtain
withThe difference between Eqs. (19) and (12) are the integration limits; in the case of Eq. (19), the lower integration limit can be below zero. This reflects the fact that for the coordinate frame moving with velocity $c/n_g(\omega _0)$, where $n_g(\omega _0)>1$, some frequency components of the radiation can move faster than the coordinate frame, resulting in a seeming breach of causality. In contrast, for the case of Eq. (12), the coordinate frame moves with velocity of $c$, which leads to $H_{\delta z}(\tau )$ being equal to zero for negative $\tau$.
The propagator operator for the non-infinitesimal propagation distance can be constructed using Eq. (14). Compared to Section 3, here, the mapping needs to be modified due to $H_{\delta z}(\tau )$ being nonzero for negative arguments. Namely, starting from a grid with $N_t$ points, we map $H_{\delta z}$ onto the middle section of a grid with $3N_t$ points, i.e. from point $N_t+1$ to $2N_t$, padding the points from 1 to $N_t$ and from $2N_t+1$ to $3N_t$ with zeros. Second, we perform the multiplication in the Fourier domain: $H_{2\delta z}=\hat {F}^{-1}[(\hat {F}H_{\delta z})^2]$. The obtained $H_{2\delta z}$ is then mapped back onto the $N_t$-points grid by discarding the values from $1$ to $N_t$ and from $2N_t+1$ to $3N_t$. The mapping onto the $3N_t$-points grid is required to avoid wrap-arounds due to both negative and positive group delays.
Note that Eq. (18) formally coincides with the propagation equation for the envelope in the slowly-varying-envelope approach, if one redefines the real-valued electric field $E(z,t)$ as a complex-valued envelope. Therefore, the above formulation of the propagator operator, given by Eqs. (20) and (14), as well as the numerical recipes to calculate such an operator, can be applied also in this case – with the only difference that all the Green’s functions become complex.
In the recent past, first-principle simulations, which combine propagation equations with the TDSE for the calculation of the polarization, have attracted increasing attention. The propagator operator approach is particularly suitable for such simulations. In order to apply it, one needs to calculate the Green response function, $G(\tau )$, by using a delta-function-shaped input pulse, i.e. with a non-zero electric field only at one point of the numerical grid. This simulation has to be done only once. After this, one can easily split the propagation step into two, similar to the idea of the split-step Fourier approach: the linear step and the nonlinear step. For the nonlinear step, one subtracts the linear part from the total dipole using the Green response function $G(\tau )$, and solves Eq. (1). For the linear step, the propagator operator constructed according to Eq. (14) is used. We have successfully implemented this approach for a highly nonlinear (2+1)D simulation of pulsed beam propagation in atomic hydrogen (results to be published elsewhere).
6. Conclusion
We have presented a novel concept for calculating the linear propagation of light pulses in resonant media. We have shown that our method is exact and preferable to existing standard methods, as it is free of nonphysical, non-causal artifacts and numerical instabilities. To demonstrate its advantages, we have performed benchmark calculations with 1D atomic hydrogen, and compared the results with those from standard methods. Extensions of the propagation operator method to other physically and numerically relevant situations have been discussed.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. F. Krausz and M. Ivanov, “Attosecond physics,” Rev. Mod. Phys. 81(1), 163–234 (2009). [CrossRef]
2. E. Lindroth, F. Calegari, L. Young, M. Harmand, N. Dudovich, N. Berrah, and O. Smirnova, “Challenges and opportunities in attosecond and XFEL science,” Nat. Rev. Phys. 1(2), 107–111 (2019). [CrossRef]
3. N. Berrah, A. Sanchez-Gonzalez, Z. Jurek, R. Obaid, H. Xiong, R. J. Squibb, T. Osipov, A. Lutman, L. Fang, T. Barillot, J. D. Bozek, J. Cryan, T. J. A. Wolf, D. Rolles, R. Coffee, K. Schnorr, S. Augustin, H. Fukuzawa, K. Motomura, N. Niebuhr, L. J. Frasinski, R. Feifel, C. P. Schulz, K. Toyota, S.-K. Son, K. Ueda, T. Pfeifer, J. P. Marangos, and R. Santra, “Femtosecond-resolved observation of the fragmentation of buckminsterfullerene following X-ray multiphoton ionization,” Nat. Phys. 15(12), 1279–1283 (2019). [CrossRef]
4. M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’Huillier, and P. B. Corkum, “Theory of high-harmonic generation by low-frequency laser fields,” Phys. Rev. A 49(3), 2117–2132 (1994). [CrossRef]
5. J. L. Krause, K. J. Schafer, and K. C. Kulander, “High-order harmonic generation from atoms and ions in the high intensity regime,” Phys. Rev. Lett. 68(24), 3535–3538 (1992). [CrossRef]
6. P. Jürgens, B. Liewehr, B. Kruse, C. Peltz, D. Engel, A. Husakou, T. Witting, M. Ivanov, M. J. J. Vrakking, T. Fennel, and A. Mermillod-Blondin, “Origin of strong-field-induced low-order harmonic generation in amorphous quartz,” Nat. Phys. 16(10), 1035–1039 (2020). [CrossRef]
7. S. Baker, J. S. Robinson, C. A. Haworth, H. Teng, R. A. Smith, C. C. Chirila, M. Lein, J. W. G. Tisch, and J. P. Marangos, “Probing proton dynamics in molecules on an attosecond time scale,” Science 312(5772), 424–427 (2006). [CrossRef]
8. K. S. Yee, “Numerical solution of initial value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14(3), 302–307 (1966). [CrossRef]
9. R. K. Bullough, P. M. Jack, P. W. Kitchenside, and R. Saunders, “Solitons in Laser Physics,” Phys. Scr. 20(3-4), 364–381 (1979). [CrossRef]
10. A. V. Husakou and J. Herrmann, “Supercontinuum generation of higher-order solitons by fission in photonic crystal fibers,” Phys. Rev. Lett. 87(20), 203901 (2001). [CrossRef]
11. M. Kolesik, J. V. Moloney, and M. Mlejnek, “Unidirectional Optical Pulse Propagation Equation,” Phys. Rev. Lett. 89(28), 283902 (2002). [CrossRef]
12. M. Kolesik and J. V. Moloney, “Nonlinear optical pulse propagation simulation: From Maxwell’s to unidirectional equations,” Phys. Rev. E 70(3), 036604 (2004). [CrossRef]
13. T. Brabec and F. Krausz, “Intense few-cycle laser fields: Frontiers of nonlinear optics,” Rev. Mod. Phys. 72(2), 545–591 (2000). [CrossRef]
14. H. Liang, P. Krogen, Z. Wang, H. Park, T. Kroh, K. Zawilski, P. Schunemann, J. Moses, L. F. DiMauro, F. X. Kärtner, and K.-H. Hong, “High-energy mid-infrared sub-cycle pulse synthesis from a parametric amplifier,” Nat. Commun. 8(1), 141 (2017). [CrossRef]
15. E. Lorin, S. Chelkowski, and A. Bandrauk, “A numerical Maxwell–Schrödinger model for intense laser–matter interaction and propagation,” Comput. Phys. Commun. 177(12), 908–932 (2007). [CrossRef]
16. M. B. Gaarde, J. L. Tate, and K. J. Schafer, “Macroscopic aspects of attosecond pulse generation,” J. Phys. B: At., Mol. Opt. Phys. 41(13), 132001 (2008). [CrossRef]
17. Sh. Amiranashvili, U. Bandelow, and A. Mielke, “Padé approximant for refractive index and nonlocal envelope equations,” Opt. Commun. 283(3), 480–485 (2010). [CrossRef]
18. J. Javanainen, J. H. Eberly, and Q. Su, “Numerical simulations of multiphoton ionization and above-threshold electron spectra,” Phys. Rev. A 38(7), 3430–3446 (1988). [CrossRef]
