Properties of Stark Resonant States in Exactly Solvable Systems

Properties of Stark resonant states are studied in two exactly solvable systems. These resonances are shown to form a biorthogonal system with respect to a pairing defined by a contour integral that selects states with outgoing wave boundary conditions. Analytic expressions are derived for the pseudonorm, dipole moment, and coupling matrix elements which relate systems with different strengths of the external field. All results are based on explicit calculations made possible by a newly designed integration method for combinations of Airy functions representing resonant eigenstates. Generalizations for one-dimensional systems with shortrange potentials are presented, and relations are identified which are likely to hold in systems with three spatial dimensions.


Introduction
Resonance states have been used to solve a wide range of problems in the fields of nuclear physics [1,2], quantum chemistry [3], nonlinear optics [4][5][6], and semiconductor physics [7,8].Despite their widely recognized utility, relatively little is known about their general properties since they do not live in the familiar Hilbert space associated with Hermitian quantum mechanics [9].The properties and issues that are less well-understood than in Hermitian quantum mechanics include inner products, normalization and completeness [10][11][12][13][14][15][16][17], complex expectation values [18,19], and their physical interpretation [20].Despite the mathematical difficulties related to their applications, resonance states do contain valuable physical information and it is important to investigate systems that could provide some guidance.
Here we add to the present understanding of resonance systems by analytically calculating a number of useful quantities for two exactly solvable quantum systems: the 1D Diracdelta potential and 1D square-well models in the presence of a homogeneous field.Despite the latter model being a textbook example, and the former being studied and used in applications for decades (e.g., [21][22][23][24]), their resonances have so far been studied mainly with numerical tools [8,25].The quantities of interest, explicitly evaluated for the first time in this work, are the normalization factors, eigenvalue equations, dipole matrix elements, and off-diagonal transition elements that characterize the dependence of resonant basis states on the external field.
We also generalize our results to more complex models with piecewise constant potentials, and for general onedimensional systems with finite-range potentials.We identify relations between the generalized dipole moments and the gradient of the atomic potentials, which resemble similar properties in systems with self-adjoint Hamiltonians.
Last but not least, all of our new results are based on direct evaluation of integral expressions, for which we have developed a new integration technique that is applicable to functions representing Stark resonances in one dimension with a piecewise constant potential.
Beyond developing a deeper understanding of exactly solvable systems, the additional motivation for this work is in the use of resonance states as a basis for time-dependent Schrödinger evolution, with applications in modeling electron ionization and nonlinear polarization due to a time varying optical pulse field [5].Detailed study of exactly solvable systems with Stark resonant states brings multiple benefits.First, having explicit expressions for complex-valued observables and the ability to study their field and time dependence gives intuition of how one maps these complex values and open-system dynamics back to the real expectation values and the norm-preserving evolution found in 2 Advances in Mathematical Physics Hermitian quantum mechanics of a closed system.Such a connection is crucial for applications in nonlinear optics (e.g., [5,26,27]).Secondly, the ability to compare different resonance systems may indicate which properties or relations are universally valid or common to all resonance systems.For example, we have witnessed in numerical simulations a field-dependent relation connecting the expectation values of the gradient of the atomic potential to the resonant state pseudonorm for one-and three-dimensional systems.

Non-Hermitian Hamiltonians
In this section, we give some background of the class of Hamiltonians that we want to investigate.We begin with a 1D Hamiltonian that is parameterized by the strength of the external field : where the function () represents the atomic potential.To study the Stark resonances, one usually assumes outgoing wave boundary conditions at  → ∞ and seeks solutions of   =     where   is the eigenvalue.With outgoing boundary conditions, the system is open as the particle can escape toward  → ∞ and the operator  is non-Hermitian.Therefore, energy   , along with many other observables, is complex-valued.Without Hermiticity, we lose many of the guarantees of Hermitian quantum mechanics, such as conservation of the number of particles, real-valued observables, and square integrable wave functions.There is not yet a full consensus on how to handle and interpret many of these quantities, including normalization and inner products.
Due to the non-square-integrable character of wave functions   , the standard inner product and normalization prescriptions do not apply, since the integrals normally used to calculate them are divergent.Some regularization method must be used, and a number of approaches can be found in the literature [10,13,28].However, it is important to appreciate that there may not be as much choice as it may seem in how the Stark resonant states should be normalized.For example, if a resonant state expansion of Green's operator exists, the eigenstates appear in it with a definite "norm" [29].
In what follows we utilize biorthogonality of the Stark resonant system and obtain the eigenstates with such preferred normalization factors.
We consider the Hamiltonian  to act on functions living on a complex contour C, where the contour follows the real axis in the vicinity of the atom and then deviates from the real axis far from the origin.To select the space of outgoing wave functions, the contour departs into the upper complex plane as  → ∞.The shape of this contour is inconsequential, except its property that it approaches infinity in the sector of the complex plain in which all outgoing waves, and in particular the resonance states, decay exponentially.One possible example utilizing a piecewise linear path is shown in Figure 1.At the far end of the contour, outgoing wave functions that behave as ∼  →  −Im{} (with positive ) a I{z} R{z} Θ Figure 1: An example contour in the complex plane that serves as a "complexified" spatial axis in a model of an open quantum system in which a non-Hermitian Hamiltonian (1) acts in the space of functions defined along the contour.Both smooth and piecewise linear contours are admissible for our purposes.In this example, the domain of the Hamiltonian would be specified by requiring that ( − ) = ( + ) and that an analogue of Cauchy-Riemann condition,   ( − ) =  −Θ   ( + ), is satisfied for derivatives along the contour for all  ∈ ().We also assume that the potential () has a compact support with a "radius" smaller than , so that nonanalytic potentials can be considered.decay exponentially.Thus, the introduction of the contour is in the spirit of the external complex scaling.
The differential expression  acting in the space of functions defined along the contour results in a non-Hermitian operator that represents an open system.A -product defined as a contour integral is the tool that replaces the standard scalar product in working with non-self-adjoint operators.A formal argument can be made that if two resonant states belong to two different eigenvalues, then they are orthogonal with respect to the above -product [30], and the latter can serve as a definition for a pseudonorm in resonant states.
In this work we aim to avoid any reliance on formal operator properties.Instead, we show by explicit calculation of the underlying contour integrals that the following orthogonality relation holds for outgoing Stark resonances: In particular, we verify the orthogonality of different functions ( ̸ = ), and we evaluate the normalization factor   () explicitly for two exactly solvable models.

Motivation
We have recently presented a proof of principle for application of metastable electronic states to calculate nonlinear response in a time-dependent field of an optical pulse [5].We have also demonstrated that the resulting description is extremely efficient in that a single Stark resonance is sufficient to obtain quite accurate nonlinear polarization for realistic models of atoms [31].Here we recap some of the ideas behind the Metastable Electronic State Approach (MESA) in order to identify the relevant properties of resonant states.
If we represent a particular quantum state Ψ as a sum over resonance states, plus a "background," and ask what is its evolution due to the time-dependent field (), then we can find a system of equations describing the evolution of coefficients   with the help of orthogonality relation (3): where   () is the time derivative of the electric field intensity.
Here we have assumed that the expansion states   are slaved to the time-dependent field and are normalized to unity, ⟨  (()) |   (())⟩ =   , at all times.Here we neglect the coupling to the "background"   which originates from the continuum contribution contained in various resonant state expansions.MESA works with a physically motivated assumption (see [5] for details) that, for systems initially in the ground state, the temporal decay of resonant states represents ionization and that the flow of probability out of the space spanned by {  } manifests as the increase of the norm ‖  ‖.
For the purposes of this work, ( 5) identifies the quantities that we aim to calculate.First, we need the complexvalued energies   .Second, we require normalization and orthogonality relation (3) to be satisfied by all resonance states.To evaluate the induced polarization, we also need the generalized dipole moments (6) and lastly we must calculate the coupling terms ⟨    |   ⟩, which describe the change of the resonant state basis as it evolves slaved to the external field.
The coupling terms can be related to the dipole moment matrix elements with the help of the following argument utilizing the parametric dependence of the Hamiltonian on  [32,33]: Moreover, for the normalized resonances the coupling terms are antisymmetric in indices , , since This means that we have ⟨    |   ⟩ = −⟨  |     ⟩ and for  =  the self-coupling vanishes ⟨    |   ⟩ = 0 as a consequence of the -product symmetry.Thus, the evolution system (5) can be alternatively written with the substitution In the appendix, we present a new integration technique that allows calculation of the coupling terms directly, without reliance on the formal derivation underlying relation (8).

Stark Resonances and Their Properties
In this section we outline the properties of resonant wave functions for two different systems, and we find associated resonant eigenvalue equations.We also calculate explicit normalization factors for all Stark resonances using the orthogonality relation (3).Having established these tools, we continue to calculate dipole moment matrix elements and coupling factors.
In this work we assume that the potential () has a compact support contained in (−, +).Thus, the asymptotic form of both the conventional and resonant wave functions can be obtained as a combination of Airy functions,  1  +  2 .The requirement that the solutions are regular for  → −∞ dictates  2 = 0.For  > 0 one can use combinations  ± () = () ± () which behave as outgoing (+) and incoming (−) waves at  → ∞.Representation of the eigenstates of the original  (i.e., the one that acts on the real axis) which is particularly suitable for our purposes can be written as where  = −(2) 1/3 and  ± () are sought expressions representing the eigenvalue equations for outgoing (+) and incoming (−) wave functions.The fact that the original operator (i.e., the one acting on the real axis) is self-adjoint guarantees that the above states can be normalized to a delta function in energy: It is sufficient to examine the asymptotic behavior of these states to verify that this normalization is obtained with  = 2 2/3  1/6 .While a particular normalization is not crucial for us, the above form of eigenstates allows us to infer the form of the resonant functions sought below.
To obtain the remaining portion(s) of energy eigenstates, one has to "fill in" the wave function in the central region of − to  and in doing so satisfy whatever conditions a given potential imposes on them.In both cases treated in this work, this means to find functions that are continuous and to require continuity of derivatives in the square-well case and a "cusp condition" (13) in the Dirac-delta case.This procedure reveals the concrete form of expressions  ± () for a given ().
The asymptotic form of the energy eigenstates as shown in (9) indicates the form of the resonance wave functions.For example, if we find a complex-energy root of  + () = 0 the incoming part of the wave function  − will be eliminated.At the same time a pole will appear in the projection onto the outgoing wave function.This tells us that the resonance behaves as  and  + for large negative and positive , respectively.
We present our final results, including those on resonant state normalization, in the form that is independent of how one chooses to parameterize the eigenfunctions.Whenever we show intermediate results, it is for the wave functions written as follows.For the Dirac-delta model, we take the unnormalized ansatz for the outgoing resonance in the form For the system with a square-well potential of width 2 and depth  0 we take with coefficients   to be fixed to ensure continuity of wave function value and derivative across well boundaries.For these functions to become resonant states the energies in  = / and   = ( −  0 )/ must be solutions to the eigenvalue equation(s) we present next.

Eigenvalue Equations.
For the Dirac-delta potential () = −(), where  is the depth of the potential, and the eigenstate representations (9) are valid with  = 0.The delta function potential imposes a boundary condition on the wave function's value and derivative: This "cusp condition," when applied to the above eigenstate parameterization, leads directly to the well-known [4] expression for the eigenvalue equation for resonant energies: Complex-valued solutions to  + () = 0 determine the resonant energies of the outgoing Stark functions.Longer calculations are required to obtain the analogous equation [25,34] for the square-well potential.One needs to connect the outer regions with a linear combination of Airy functions (12) and eliminate the unknown coefficients.The result reads where we utilized shorthand notations to compress the otherwise long expression., , and  stand for the corresponding Airy functions, and primes denote derivatives.For a well with depth  0 and half-width , the subscripts indicate on which sides of the well walls the arguments of the functions are evaluated, with 0, 1, 2, 3 representing ( + ) at  = − − , −+, +−, and  = ++, respectively.The value of  is also dependent on where the functions are evaluated.For regions outside of the well  = / and for regions inside the well  = ( −  0 )/.Thus, the subscripts 0 and 1 represent Airy functions evaluated just outside and inside the left boundary of the well, while subscripts 2 and 3 represent arguments at the right boundary of the well.
It is helpful to visualize the resonance energy "landscapes" illustrated in Figure 2.Although the eigenvalue equation for the square-well ( 15) is more complicated, it is apparent that the two systems share some common properties.The main differing feature is the possibility of multiple bound states (highlighted in red) in the square-well potential, while the Dirac-delta system only supports a single bound state.However, there exist two other infinite families of resonances.The "right" family has eigenvalues located along the real axis and corresponds to longer living states, while the "left" family, with energies along the ray arg() = −2/3, are fast decaying, short-lived states.This resonance structure is most likely a generic feature at least in case of short-ranged attractive potentials.

Stark State Normalization.
To establish formulas for normalization and to verify the orthogonality relation (3) for both systems, we make use of the formula (VS 3.50) in Valle and Soares [35] giving a primitive (antiderivative) function for a square of arbitrary combination of Airy functions.
To calculate the contour integral(s), the corresponding primitive functions are evaluated at points of discontinuities of the potential and at both ends of the contour, and this is where the choice of the contour is important.The resonant wave functions decay exponentially for  → −∞, and the corresponding boundary terms vanish.On the other side of the contour, the asymptotic behavior of the primitive function is dominated by where  = −(2) 1/3 ,  =   /, and   is a root of the eigenvalue equation.It is straightforward to verify that asymptotically along the contour, where  ∼  Θ , this function decays for arbitrary fixed  as  → ∞.Thus, the boundary terms brought by this end of the contour also vanish.Moreover, since the integrands are in all cases entire functions (containing no singular points), the precise shape of the integration path does not affect the outcome.As a result, in any piecewise constant atomic potential () it is only the special points of () that give rise to nonzero contribution(s).Thus, direct integration along the contour, followed by simplifications making use of the eigenvalue equation and the Wronskian for Airy functions, yields the following normalization factor for the Stark resonance in the Dirac-delta model: To calculate the normalization factor for the square-well system, a similar but more complicated procedure to evaluate (3) using formula (VS 3.50) results in a surprisingly simple expression for the normalization factor: To the authors' knowledge, this is the first time these results have been presented in the explicit form.
To verify the mutual orthogonality with respect to (3) for different resonant states, we use formula (VS 3.53), together with the fact that the complex energies satisfy the eigenvalue equation(s).
To conclude this subsection, we note that our direct verification of the orthogonality relation (3) and the explicit calculation of the corresponding normalization factors means that there exists an infinite dimensional space of functions that can be expressed as superpositions of Stark resonances.In the spirit of [29], one should ask if this gives us the preferred expansion.Indeed, one can alternatively use the self-adjoint Hamiltonian eigenstates ( 9) and find the projector onto a given resonant state as the residue of   ()  () at the pole in the complex plane that corresponds to its energy  =   .While it is beyond the scope of the present paper, we note that the two approaches in fact lead to the same expansion coefficients.

Dipole Matrix Elements.
Next we calculate the generalized dipole matrix elements, both diagonal and offdiagonal.Expressed with the help of unnormalized resonance eigenfunctions [( +   )], the contour integrals we need to evaluate read Note that these quantities differ from their Hermitian counterparts as there is no complex conjugation in the integrand, and the result is complex-valued [19].Here we assume that the contour C is chosen such that it only starts to deviate from the real axis for  > , that is, outside of the potential support.For both the delta potential model and the square-well potential system, we integrate over each distinct interval of constant potential  making use of the formula (VS 3.51).As the shape of the contour ensures vanishing contributions from its ends, one only needs to evaluate the primitive functions at  = 0 for Dirac-delta and  = ± for the square-well system.
For the diagonal matrix element in the Dirac-delta potential model, we obtain the following expression in terms of Airy functions: These results can be further simplified using the formula for the normalization factor  2 together with the eigenvalue equation ( 14) and are found to be related to the change in the eigenvalue with respect to the field  as follows: where the second equality can be verified by differentiating the eigenvalue equation ( 14) with respect to .
For the off-diagonal dipole matrix elements, we use (VS 3.54).Utilizing  , = ( , ) and  , = ( , ) to shorten the notation, we find that Grouping terms in order to identify Wronskians allows us to simplify the expression down to which in turn can be written solely in terms of the wave function properties at  = 0, with   standing for the derivative of the normalized Stark wave function.
Let us proceed with the dipole moment calculations for the square-well potential.The matrix elements can be found by again making use of (VS 3.51) and (VS 3.54).Thanks to the continuity properties of the wave function, many terms that arise in the course of this calculation cancel, and the resulting diagonal terms are which is simplified further using the normalization factor (again as with the delta model) to It can be checked numerically that the diagonal elements are also equal to −    , just as in the delta potential model.Now we calculate the off-diagonal elements.When using (VS 3.54), taking into account that continuity of the wave function allows us to ignore terms that do not involve  and the sum   +   , as a result every other term cancels, and the resulting equation is surprisingly simple: and one should note the similarity with its counterpart formula for the Dirac-delta model.We thus arrive at the conclusion that all dipole matrix elements can be expressed in simple formulas which only depend on the values of the wave functions (and their derivatives) at the special points given that characterize the potential.It will be interesting to see if these results can be generalized for arbitrary systems with piecewise constant potentials.

Coupling Matrix Elements.
We now turn our attention to the terms ⟨    |   ⟩ identified in (5).These quantities mediate the connections between Stark resonances of a given system at different values of the field  and we call them accordingly coupling terms.They are needed to describe the evolution of the system in a time-dependent ().Instead of trusting the formally derived relation (8), we compute these quantities through direct integration.Thus, our result can be also interpreted as a direct verification of (8) for two model systems.Moreover, the procedure and the particular representation of the results lead us to a generalization for arbitrary one-dimensional systems with short-range potentials.
To calculate the coupling terms, we first differentiate the wave function with respect to the field  and then integrate It is clear that the diagonal ( = ) coupling terms between normalized states vanish, and we only need to calculate offdiagonal elements and therefore can disregard the second term, leaving only the integral ⟨    |   ⟩.For the model with Dirac-delta potential, the integral is where Δ = 1/(3 2  3  2 ), and   ,   ,   stand for the corresponding Airy functions or their combinations evaluated at   (i.e., at  = 0).
For the model with the square-well potential we have a combination of integrals of similar types; namely, where  and  are the coefficients that guarantee continuity of   and   , respectively, at the well boundaries.These factors are expressed in terms of Airy functions and depend on  through ,   , and   .Their primed notations represent derivatives with respect to .
The above integrals for both models contain terms of two kinds.The first group can be evaluated making use of known, previously published Airy integrals.These have the form ∫    d, where ,  are a linear combination of Airy functions with arguments ( + ), where  is different in  and .Thus, the first two lines in (29) and first three lines in (30) can be dealt with, although each integral requires lengthy computations, especially in the square-well case.Then there are integrals of the types ∫      d and ∫      d.No known formulas are available for these, and we have developed a new integration technique that we outline in the appendix.As a result, analytic results can be obtained for both systems.
For the Dirac-delta model, the individual integrals on the RHS of (29) simplify pairwise.The first pair is calculated using (VS 3.53), while the last two pairs use the identities derived in the appendix.The intermediate expression reads The first two lines can be simplified using the Wronskian and normalization factors, while the third can be rewritten with the help of the eigenvalue equation.Combining these shows that the first three lines sum up to zero, leaving only the last term which we write in terms of derivatives of the wave functions at the origin: Comparing this expression to (24), it is clear that the relation between coupling term and off-diagonal dipole matrix elements ( 8) is valid.
Coupling terms (30) for the square-well system are also calculated in analytic form but due to excessive number of terms comprising the result they are not listed here.In principle, a similar procedure to simplify the square-well coupling term (30) should work.However, the resulting expression is extremely large and we could not find a practical way to compress it to a manageable length.The main difficulty in simplifying the square-well result is that the eigenvalue equation is much more complicated (cf.( 14) to ( 15)).Nevertheless, having explicit formulas in terms of Airy functions, we verified numerically that the integrated result does relate to the dipole matrix element as suggested by (8): To conclude this section, we have shown that the coupling terms can be directly calculated using a new Airy integral technique detailed in the appendix.While our explicit calculations do not justify the formal steps taken to obtain (8), they do corroborate that the relation between the coupling terms and the dipole moment holds.Pragmatically, one should choose to use the dipole matrix elements, since they are easier to calculate numerically.

Generalization for Arbitrary Piecewise Constant Potentials
Comparing results (18) and (27) we see that they have a similar form, with the dipole moment matrix elements expressed in terms of wave function values evaluated at discontinuities of the potential .While the same result can be easily written explicitly with Airy functions, this particular form indicates that the expressions are in fact sums over atomic potential discontinuities, with weights corresponding to the potential-value jumps.This suggests the following generalization of the dipole matrix element formulas for a system with arbitrary piecewise constant potential: where the sums run over all potential discontinuities.It is in fact not too difficult to realize that the procedures utilized above can be modified for a more general case of piecewise constant potential.From there, one can take a continuum limit, approximating an arbitrary potential as a limit of piecewise constant functions, and arrive at We have assumed that the potential is short-ranged and the integration in these formulas is along the real axis (i.e., the contour C is not necessary for convergence).This is an intriguing result, because an identical formula can be derived for the discrete-energy eigenstates of a selfadjoint Hamiltonian by evaluating its double commutator with the position operator.However, here we have Stark resonances represented by complex-valued functions living on the contour C.So it seems that as long as the normalization and "scalar product" are defined with the help of pairing (3), the Stark resonance states satisfy relations analogous to those obeyed by their self-adjoint counterparts.
We have used numerical simulations (not shown here) to verify that the relation between the Stark resonance pseudonorm and the generalized expectation value of the "atomic" potential gradient could also be valid for threedimensional systems.It is tempting to speculate that the offdiagonal dipole element relation (37) could be generally valid for Stark resonant states in higher dimensions.

Conclusion
We have derived analytic expressions for a number of quantities that characterize the Stark resonance states in two exactly solvable systems.The first model studied in this work is the one-dimensional particle in a Dirac-delta potential with additional homogeneous field, and the second has the squarewell potential.
We have studied these systems as open, non-Hermitian models, and identified a natural choice for the pairing connecting the states in the domain of the Hamiltonian with the states in the domain of its adjoint operator.With respect to this pairing, Stark resonances form an orthogonal system, and many of their properties can be evaluated analytically.
Despite the fact that both models have been studied for years, explicit expressions for their (pseudo-) norms, dipole moment expectation values, and their relations connecting the resonant state wave functions at different field values are new.Our results thus further the understanding of the mathematical properties that underline the Stark effect.
Moreover, we have shown that certain results naturally extend to a wide class of one-dimensional models and we have also identified relations that appear to be candidates for properties generally applicable to three-dimensional Stark systems.In particular, we have found that the generalized dipole moment matrix elements between the nonphysical resonant states can be related to the expectation values of the atomic potential gradient in a way that is completely analogous to relations that hold for real-valued discreteenergy eigenfunctions in Hermitian systems.We speculate that these relations apply generally in three dimensions and could be used in numerical calculations to assess the fidelity of the resonant eigenfunctions.
An important by-product of this study is a new integration technique applicable to combinations of Airy functions that represent Stark resonances in one-dimensional models with piecewise constant potentials.
Results presented in this work have also an immediate practical impact on modeling of light-matter interactions in strong time-dependent optical fields in the framework of the Metastable Electronic State Approach.

Figure 2 :
Figure 2: Outgoing resonance eigenvalue equation landscapes for Dirac-delta ((a)  = 1) and square-well potential ((b)  = 4 and  0 = −0.5)systems in the external field  = 0.03.To visualize the locations of the energy eigenvalues, we evaluate (14) and (15) over a range of  in the complex plane and convert || to a height map via the formula (1 − (1 + || 0.3 ) −1 + ) −1 , so that its roots are represented by poles that are easy to locate ( ≈ 0.1).