text:
# Wannier-Stark states of a quantum particle in 2D lattices

## I Introduction

The quantum states of a particle in a periodic potential plus homogeneous field (known nowadays as the Wannier-Stark states, WS-states in what follows) are one of the long-standing problems of single-particle quantum mechanics. The beginning of the study of this problem dates back to the paper by Bloch of 1929, followed by contributions of Zener, Landau, Wannier, Zak and many others history . In the late eighties the problem got a new impact by the invention of semiconductor superlattices. The unambiguous observation of the WS-spectrum in a semiconductor superlattice mendez  ended a long theoretical debate about the nature of WS-states, and now it is commonly accepted that they are the resonance states of the system. Besides, WS-states were recently studied in a system of cold atoms in an optical lattice atoms  and some other (quasi) one-dimensional systems.

Although WS-states are resonances, i.e. metastable states, in the theoretical analysis of related problems they were usually approximated by stationary states (one-band, tight-binding, and similar approximations). Beyond the one-band approximation, WS-states in the semiconductor and optical lattices were studied in recent papers preprint  and papers  by using the scattering matrix approach of Ref. PRL1  (see also Ref. JOB2  for details). This approach actually solves the one-dimensional Wannier-Stark problem and supplies exhaustive information about 1D WS-states. In the present letter we extend the method of Ref. PRL1 ; JOB2  to the case of two-dimensional lattices. For the first time we find the complex spectrum of 2D WS-states and analyze its general structure.

To be concrete, we choose the following system:

$$H=𝐩^2/2+V(𝐫)+𝐅𝐫,𝐫=(x,y),$$
(1)
$$V(𝐫)=\mathrm{cos}x+\mathrm{cos}yϵ\mathrm{cos}x\mathrm{cos}y,$$
(2)
where $`0ϵ1`$ hbar . Two limiting cases $`ϵ=0`$ and $`ϵ=1`$ correspond to an ‘egg crate’ potential, for which the system is separable, and a ‘quantum well’ potential, where the coupling between two degrees of freedom is maximal (see Fig. 1). Let us also note that the choice $`ϵ=1`$ corresponds to a 2D optical potential created by two standing laser waves crossing at right angle. Thus the results presented below can be directly applied to the system of cold atoms in a 2D optical lattice.

## II 2D Wannier-Bloch spectrum

We briefly recall the key points of the 1D theory. The spectrum of the Bloch particle in the presence of a static field consists of several sets of equidistant levels

$$_{\alpha ,l}=E_\alpha +2\pi Fl\mathrm{i}\mathrm{\Gamma }_\alpha /2,$$
(3)
known as Wannier-Stark ladders of resonances. In Eq. (3), $`2\pi `$ stands for the lattice period, $`F`$ is the amplitude of the static force, $`l=0,\pm 1,\mathrm{}`$ is the site index and the index $`\alpha =0,1,\mathrm{}`$ labels different ladders. The lifetime of WS-states $`\mathrm{\Psi }_{\alpha ,l}(x)`$ is defined by the resonance width $`\mathrm{\Gamma }_\alpha `$ as $`\tau _\alpha =\mathrm{}/\mathrm{\Gamma }_\alpha `$. Typically, the lifetime $`\tau _\alpha `$ rapidly decreases with increasing index $`\alpha `$. Because of this only the first few WS-ladders are of physical importance.

Along with the WS-states $`\mathrm{\Psi }_{\alpha ,l}(x)`$, one can also introduce Wannier-Bloch states (WB-states) by

$$\psi _{\alpha ,k}(x)=\underset{l}{}\mathrm{\Psi }_{\alpha ,l}(x)\mathrm{exp}(\mathrm{i2}\pi kl).$$
(4)
As follows from the definition (4), the continuous evolution of WB-states obeys the equation $`\psi _{\alpha ,k}(x,t)=\mathrm{exp}(\mathrm{i}_\alpha t/\mathrm{})\psi _{\alpha ,kFt/\mathrm{}}(x)`$, where $`_\alpha =E_\alpha \mathrm{i}\mathrm{\Gamma }_\alpha /2`$. Thus, WB-states can be alternatively defined as the eigenfunction of the evolution operator over the Bloch period $`T_B=\mathrm{}/F`$ niu2 . (Note that the eigenvalues of the evolution operator form degenerate bands $`_\alpha (k)=_\alpha `$). Additionally, to ensure that $`\psi _{\alpha ,k}(x)`$ are resonance states of the system, the eigenvalue equation for the evolution operator should be accomplished by the specific non-hermitian boundary condition. It was proven in Ref. JOB2  that the required boundary conditions are imposed by the truncation of the evolution operator matrix in the momentum representation.

We proceed with the two-dimensional case. As mentioned above, WB-states in a 1D lattice can be defined as the non-hermitian eigenstates of the evolution operator over one Bloch period. In the 2D problem there are two different Bloch periods associated with the two components of the static field. Therefore the notion of the WB-states can be introduced only in the case of commensurate periods, i.e., in the case of ‘rational’ direction of the field ($`q,r`$ are coprime integers):

$$F_x=\frac{qF}{(r^2+q^2)^{1/2}},F_y=\frac{rF}{(r^2+q^2)^{1/2}}.$$
(5)
Provided condition (5) is satisfied, we define 2D WB-states as the non-hermitian eigenfunctions of the system evolution operator over the common Bloch period $`T_B=(r^2+q^2)^{1/2}\mathrm{}/F`$. Using the Kramers-Henneberger transformation, which is just the gauge which transforms the static term into the vector potential, the evolution operator can be presented in the form

$$\widehat{U}(T_B)=\mathrm{e}^{\mathrm{i}qx}\mathrm{e}^{\mathrm{i}ry}\widehat{\mathrm{exp}}\left(\frac{\mathrm{i}}{\mathrm{}}_0^{T_B}dt\stackrel{~}{H}(t)\right),$$
(6)
$$\stackrel{~}{H}(t)=\frac{(\widehat{p}_xF_xt)^2}{2}+\frac{(\widehat{p}_yF_yt)^2}{2}+V(x,y),$$
(7)
which reveals its translational invariance (the hat over the exponent sign denotes time ordering). Alternatively, we can rotate the coordinates so that the direction of the field coincides with the $`x^{}`$-axis:

$$x^{}=\frac{qx+ry}{(r^2+q^2)^{1/2}},y^{}=\frac{qyrx}{(r^2+q^2)^{1/2}}.$$
(8)
Transformation (8) introduces a new lattice period $`a=2\pi (r^2+q^2)^{1/2}`$ and reduces the size of the original Brillouin zone $`s=r^2+q^2`$ times. Associated with the new lattice period is a new Bloch time $`T_a=(r^2+q^2)^{1/2}\mathrm{}/F`$, which is $`s`$ times shorter than the original Bloch time $`T_B`$. Using $`\widehat{p}_x^{}=\mathrm{i}\mathrm{}/x^{}`$ and $`\widehat{p}_y^{}=\mathrm{i}\mathrm{}/y^{}`$, the time evolution operator over the new Bloch time $`T_a`$ in the rotated coordinates has the form

$$\widehat{U}^{}(T_a)=\mathrm{e}^{\mathrm{i2}\pi x^{}/a}\widehat{\mathrm{exp}}\left(\frac{\mathrm{i}}{\mathrm{}}_0^{T_a}dt\stackrel{~}{H}^{}(t)\right),$$
(9)
$$\stackrel{~}{H}^{}(t)=\frac{(\widehat{p}_x^{}Ft)^2}{2}+\frac{\widehat{p}_y^{}}{2}+V(x^{},y^{}).$$
(10)
Then, presenting the wave function as

$$\psi (𝐫^{})=\mathrm{e}^{i𝐤^{}𝐫^{}}\underset{𝐧^{}}{}c_𝐧^{}𝐫^{}|𝐧^{},𝐫^{}|𝐧^{}=\frac{1}{a}\mathrm{e}^{\mathrm{i2}\pi 𝐧^{}𝐫^{}/a},$$
(11)
we get the matrix equation

$$\underset{𝐦^{}}{}U_{𝐧^{}𝐦^{}}^{(𝐤^{})}c_𝐦^{}=\mathrm{e}^{\mathrm{i}T_a/\mathrm{}}c_𝐧^{},$$
(12)
where $`U_{𝐧^{}𝐦^{}}^{(𝐤^{})}`$ denotes the $`k^{}`$-dependent matrix elements of the operator (9):

$$U_{𝐧^{}𝐦^{}}^{(𝐤^{})}=𝐧^{}|\mathrm{e}^{\mathrm{i}𝐤^{}𝐫^{}}\widehat{U}^{}(T_a)\mathrm{e}^{\mathrm{i}𝐤^{}𝐫^{}}|𝐦^{}.$$
(13)
Similar to the 1D case, the truncation of the infinite unitary matrix (13),

$$|n_x^{}|,|m_x^{}|N\mathrm{},|n_y^{}|,|m_y^{}|M\mathrm{},$$
(14)
which is presumed in the numerical calculations, automatically imposes the non-hermitian boundary condition along the $`x^{}`$-direction. (Truncation of the matrix over the index $`n_y^{}`$, $`m_y^{}`$ does not change the hermitian boundary condition along the $`y^{}`$-direction.) Then the eigenvalues $``$ obtained by numerical diagonalization of the truncated matrix correspond to the quantum resonances.

In the transformed coordinates, the unit cell with area $`a^2=(2\pi )^2s`$ contains $`s`$ different sublattices, and each of them supports its own WB-states. The sublattices are related by primitive translations of the unrotated lattice, and correspondingly the energies of their WB-states differ by multiples of $`aF/s`$. Furthermore, as function of the quasimomentum, the energies $`=_\beta ^{(i)}(k_x^{},k_y^{})`$ (here $`\beta =0,1,\mathrm{}`$ is the ‘Bloch band’ index and $`i=1,\mathrm{},s`$ is the sublattice index) do not depend on $`k_x^{}`$. This follows from the fact that a change of $`k_x^{}`$ in Eq. (13) can be compensated by shifting the time origin in Eq. (9). For the $`y^{}`$-degree of freedom the Bloch theorem can be applied, and therefore $`_\beta ^{(i)}(k_x^{},k_y^{})`$ is a periodic function of $`k_y^{}`$ with generally nonzero amplitude $`\mathrm{\Delta }_\beta `$. Thus, assuming a rational direction of the field, in each fundamental energy interval $`aF`$, the static field induces $`s=r^2+q^2`$ identical sub-bands, separated by the energy interval $`aF/s`$. Simultaneously, the size of the Brillouin zone is reduced by a factor $`s`$. This result resembles the one obtained for a 1D lattice affected by a time-periodic perturbation niu  or that for a 2D lattice in a magnetic field azbel . In these cases – provided the condition of comensurability between the Bloch period and the period of the driving force or the condition of ‘rationality’ for the magnetic flux through a unit cell, respectively, is fulfilled – the (quasi)energy spectrum of the system has a similar structure.

We conclude this section with a remark concerning the numerical procedure. Although the reduced Brillouin zone approach described above is the most consistent, we found it more convenient to diagonalize the evolution operator without preliminary rotation of the coordinate. In other words, in order to find the WB-spectrum, we solve the eigenvalue equation (12) with the truncated matrix constructed on the basis of the operator (6). As a result of the diagonalization, one obtains eigenvalues $`_\beta (k_x,k_y)`$ with quasimomentum $`𝐤=(k_x,k_y)`$ defined in the original Brillouin zone. Because the WB-bands are uniform along the direction of the field, $`_\beta (k_x,k_y)`$ is a periodic function of both $`k_x`$ and $`k_y`$ with periods $`1/r`$ and $`1/q`$ respectively. The energies obtained in this way can then be used to construct the complete WB-spectrum $`_\beta ^{(i)}(k_x^{},k_y^{})`$, $`i=1,\mathrm{},s`$. In the next section we present results of a numerical calculation of the dispersion relation $`_\beta (k_x,k_y)`$ for the periodic potential (2) and moderate values of the static field $`𝐅=(F_x,F_y)`$, $`|𝐅|=F=const`$.

## III Numerical results

It is instructive to begin with the separable case $`ϵ=0`$. In this case, 2D WB-states are given by the product of 1D states and 2D WB-energies are just the sum of 1D energies. In what follows we restrict ourselves to analyzing only the ground band. First we consider the real part of the spectrum $`E_0=\mathrm{Re}(_0)`$.

It was shown in the previous section that for rational directions of the field the ground WB-subbands repeat with energy splitting $`aF/s`$. As an example, Fig. 2 shows the relative positions of these subbands as a function of the angle $`\theta =\mathrm{arctan}(r/q)`$ for $`\mathrm{}=2`$ and $`F=0.08\sqrt{2}`$. We recall that in the considered case of a separable potential the bands have zero width for any $`\theta 0,\pi /2`$.

The main difference between separable and non-separable potentials is that the subbands $`E_0^{(i)}(𝐤)`$ have a finite width in the latter case. This is illustrated by Fig. 3(a) which shows the dispersion relation $`E_0(k_x,k_y=0)`$ for the potential (2) with (from top to bottom) $`ϵ=0`$, $`0.1`$, $`0.5`$, and $`1`$. The direction of the field is $`\theta =\pi /4`$, i.e. $`r=q=1`$. The amplitude of the static field and the value of the scaled Planck constant are the same as in Fig. 2. It is seen in Fig. 3(a) that the WB-bands gain a finite width as $`ϵ`$ is increased. We also calculated the dispersion relation $`E_0(k_x,k_y=0)`$ for different angles $`\theta =\mathrm{arctan}(r/q)`$, with $`r,q6`$. It was found that the band widths $`\mathrm{\Delta }E_0=\mathrm{\Delta }E_0(r,q)`$ are typically much smaller than the mean energy separation between the subbands. Thus, for practical purpose, one can neglect the band width for the real part of the spectrum. (An exception is the case $`\theta =0,\pi /2`$ where the width of the WB-bands approximately coincides with the width of the Bloch band in the absence of the static field.) Neglecting the width of the bands they were found to form a structure similar to that shown in Fig. 2.

We proceed with the analysis of the decay rate of the WB-states, which is determined by the imaginary part of the complex energy, $`\mathrm{\Gamma }_0=2\mathrm{I}\mathrm{m}(_0)`$. In the case of a separable potential the dependence $`\mathrm{\Gamma }_0=\mathrm{\Gamma }_0(F,\theta )`$ is obviously given by the equation

$$\mathrm{\Gamma }_0(F,\theta )=\mathrm{\Gamma }_0^{}(F\mathrm{cos}\theta )+\mathrm{\Gamma }_0^{}(F\mathrm{sin}\theta ),$$
(15)
where $`\mathrm{\Gamma }_0^{}(F^{})`$ stands for the width of 1D WS-resonances. For the parameters used ($`\mathrm{}=2`$ and $`F=0.08\sqrt{2}`$) the dependence (15) is shown in Fig. 4 by a solid line. The maximum around $`\theta =\pi /2`$ originates from a peak-like behavior of $`\mathrm{\Gamma }_0^{}(F^{})`$ and is explained by the phenomenon of 1D resonant tunneling JOB2 .

For a non-separable potential and rational direction of the field the decay rate depends on the quasimomentum. For the particular case $`\theta =\pi /4`$ this dependence is depicted in Fig. 3(b). We would like to note the complicated behavior of $`\mathrm{\Gamma }_0(𝐤)`$. The oscillating character of the decay rate is an open problem for the present day. Because the decay rate depends on the quasimomentum it might be convenient to introduce the notion of $`\overline{\mathrm{\Gamma }}_0`$, where the average is taken over the reduced Brillouin zone. The dots in Fig. 4 show the values of $`\overline{\mathrm{\Gamma }}_0`$ for some rational direction of the field and two different values of $`ϵ`$. It is seen that for a small $`ϵ=0.1`$ the ratio $`\mathrm{\Delta }\mathrm{\Gamma }_0/\overline{\mathrm{\Gamma }}_0`$ is small and the obtained dependence $`\overline{\mathrm{\Gamma }}_0=\overline{\mathrm{\Gamma }}_0(r,q)`$ essentially reproduces that of the separable case. However, this is not valid for $`ϵ=1`$, where the decay rate varies wildly. Thus, in the case of strong coupling between two degrees of freedom the description of WS-state by a mean decay rate is insufficient.

## IV Conclusion

We studied Wannier resonances in a 2D system, mainly discussing the complex energy spectrum of the Wannier-Bloch states. However, because the latter are related to the Wannier-Stark states by a Fourier transformation, the obtained results can be easily reformulated in terms of the Wannier-Stark resonances. Then the following is valid. (i) Neglecting the asymptotic tail, WS-states are localized functions along the direction of the field. (This follows from the degeneracy of WB-bands along the field direction.) (ii) For any rational direction of the field \[see Eq. (5)\] WS-states are Bloch waves in the transverse direction. (iii) For a non-separable potential the corresponding energy bands have a finite width. (iv) For the real part of the spectrum, the band widths are small and can be well neglected for $`r,q>1`$.

We also found a nontrivial dependence of the resonance width (inverse lifetime of WS-states) on the direction of the field. Because the value of the resonance width defines the decay of the probability, a complicated behavior of the survival probability is expected when the direction of the field is varied. The detailed study of the probability dynamics is reserved for future publication.


format: text

text:
# Critical exponents of the random-field O(𝑁) model

## Abstract

The critical behavior of the random-field Ising model has been a puzzle for a long time. Different theoretical methods predict that the critical exponents of the random-field ferromagnet in $`D`$ dimensions are the same as in the pure $`(D2)`$-dimensional ferromagnet with the same number of the magnetization components. This result contradicts the experiments and simulations. We calculate the critical exponents of the random-field O$`(N)`$ model with the $`(4+ϵ)`$-expansion and obtain values different from the critical exponents of the pure ferromagnet in $`2+ϵ`$ dimensions. In contrast to the previous approaches we take into account an infinite set of relevant operators emerging in the problem. We demonstrate how these previously missed relevant operators lead to the breakdown of the $`(6ϵ)`$-expansion for the random-field Ising model.

The role of quenched disorder in condensed matter depends on its strength but even weak disorder can strongly modify the nature of the phase transition. This effect is most prominent in the case of random-field (RF) disorder which breaks both the translational symmetry and the symmetry with respect to transformations of the order parameter. In particular, in some systems the arbitrarily weak disorder of this type destroys long-range order . The strong effect of weak disorder makes it difficult to apply the standard perturbative methods of phase transition theory to RF problems. Besides, it is much more difficult to solve exactly a random model than a pure one. As a result, a theory of phase transitions in the presence of random fields is still absent.

The large number of systems with RF disorder provides a strong motivation to develop such theory. Some of these systems are known for a long time. Examples are disordered antiferromagnets in the external magnetic field , binary liquids in random porous media  and vortices in disordered superconductors . Recently a lot of attention was devoted to disordered liquid crystals  and liquid He-3 in aerogels . In contrast to the pure systems even the question of ordering at low temperatures in the presence of random fields is nontrivial. The more difficult problem of the critical behavior is still open in spite of two decades of investigations.

After a lot of controversy the lower critical dimension of the RF Ising model was found  but the structure of the phase diagram and the nature of the transition to the ferromagnetic state are still unclear. Different theoretical approaches to the paramagnet-ferromagnet transition predict the dimensional reduction : the critical exponents of the RF O$`(N)`$ model in $`D`$ dimensions should be the same as in the pure O$`(N)`$ model in $`D2`$ dimensions. This prediction relates the three-dimensional RF Ising model to the one-dimensional pure system in contradiction to the presence of long-range order in the former and its absence in the later. Moreover, the high-temperature expansion  shows that the dimensional reduction rule is invalid in any dimension less than the upper critical dimension 6. The most elegant derivation of the dimensional reduction  provides an (incorrect) exact solution of the RF O$`(N)`$ model at zero temperature. The failure of that approach was explained by Parisi on the basis of the complicated energy landscape of the model . However, such explanation is insufficient for the renormalization group (RG) and $`1/N`$-expansion. It was conjectured  that RG fails because of non-perturbative corrections to the expansion in $`ϵ`$. However, such corrections are not found and recent numerical results  suggest that the existing predictions for the coefficients of the series in $`ϵ`$ are wrong .

A possible reason why the perturbation theory does not provide satisfactory results is the appearance of some additional relevant operators missed by the existing approaches. This phenomenon is responsible for the failure of the standard RG theory to predict the order of the phase transition in some systems . The possibility of a similar phenomenon in the RF O$`(N)`$ model was first suggested in Ref.  but no new values of the critical exponents were found. Recently a similar idea was used in Ref.  but it also did not allow to calculate the exponents. Besides, the approach  involves delicate manipulations with the replica limit and its validity is not clear . In the present paper we demonstrate that there are indeed some additional relevant operators in the problem although not those found in Ref. . We calculate the critical exponents of the RF O$`(N)`$ model in $`4+ϵ`$ dimensions and show that they do not satisfy the dimensional reduction. We also demonstrate how the additional relevant operators lead to the breakdown of RG in $`6ϵ`$ dimensions. Some phenomenological approaches  allowed one to obtain the critical exponents different from the dimensional reduction prediction. However, their results based on different unproven assumptions contradict each other. An important breakthrough was made by Mezard and Young  who considered the possibility of the replica symmetry breaking in the RF O$`(N)`$ model at large $`N`$. Unfortunately, the approach  did not allow the calculation of the exponents. The present paper contains the first systematic method that explains the failure of RG and allows us to find the critical exponents of the RF systems.

In the critical point the RF O$`(N)`$ model can be described by the connected and disconnected correlation functions :

$$G_{\mathrm{con}}(𝐪)=[𝐧(𝐪)𝐧(𝐪)][𝐧(𝐪)][𝐧(𝐪)]q^{2+\eta };G_{\mathrm{dis}}(𝐪)=[𝐧(𝐪)𝐧(𝐪)]q^{4+\overline{\eta }},$$
(1)
where $`𝐧(𝐪)`$ is the Fourier component of the magnetization, the square brackets denote the thermal average and the angular brackets denote the disorder average. To calculate the critical exponents $`\eta `$ and $`\overline{\eta }`$ we develop the $`ϵ`$-expansion near the lower critical dimension 4 following the line of Ref. . Our starting point is the Hamiltonian of the RF O$`(N)`$ model

$$H=d^Dx[J\underset{\mu }{}_\mu 𝐧(𝐱)_\mu 𝐧(𝐱)\underset{k}{}(𝐡_k(𝐱)𝐧(𝐱))^k],$$
(2)
where the unit vector of the magnetization $`𝐧(𝐱)`$ has $`N`$ components and $`𝐡_k`$ are random fields with zero average, $`h_{k,\alpha }(𝐱)h_{q,\beta }(𝐲)=H_k\delta _{\alpha \beta }\delta _{kq}\delta (𝐱𝐲)`$. The Hamiltonian includes an infinite set of random anisotropies of different ranks. These contributions are allowed by symmetry and turn out to be relevant in the RG sense. The replica Hamiltonian has the form

$$H_R=d^Dx[\underset{a}{}\frac{1}{2T}\underset{\mu }{}_\mu 𝐧_a_\mu 𝐧_a\underset{ab}{}\frac{R(𝐧_a𝐧_b)}{T^2}+\mathrm{}],$$
(3)
where $`a,b`$ are replica indices, $`R(z)`$ is some function, $`T`$ the temperature and the dots denote the irrelevant terms. Near the zero-temperature fixed point the whole function $`R(z)`$ is relevant. Indeed, to ensure the fixed length condition $`𝐧_a^2=1`$ at each RG step, we ascribe the scaling dimension 0 to the magnetization $`𝐧`$. The scaling dimension of the temperature is $`2+O(ϵ)`$. The relevance of any operator is determined by the number of the derivatives in it and the power in which it contains the temperature. This shows that all operators $`R_k=_{ab}(𝐧_a𝐧_b)^k/T^2`$ are relevant in the same space dimensions. The functional RG analysis in $`4+ϵ`$ dimensions repeats the same steps as the $`(4ϵ)`$-expansion developed in Ref.  and results in the following RG equation for the disorder strength

$`0={\displaystyle \frac{dR(\varphi )}{d\mathrm{ln}L}}=ϵR(\varphi )+(R^{\prime \prime }(\varphi ))^22R^{\prime \prime }(\varphi )R^{\prime \prime }(0)`$ (4)
$`(N2)[4R(\varphi )R^{\prime \prime }(0)+2\mathrm{c}\mathrm{t}\mathrm{g}\varphi R^{}(\varphi )R^{\prime \prime }(0)\left({\displaystyle \frac{R^{}(\varphi )}{\mathrm{sin}\varphi }}\right)^2]+O(R^3,T),`$ (5)
where we define $`\mathrm{cos}\varphi =𝐧_a𝐧_b`$ to make the equation more compact. The RG equation is valid in the zero-temperature fixed point that describes the phase transition in the RF ferromagnet. The critical exponents (1) can be expressed  via the RG charge $`R(\varphi )`$:

$$\eta =2R^{\prime \prime }(0);\overline{\eta }=2(N1)R^{\prime \prime }(0)ϵ.$$
(6)
Note that $`R^{\prime \prime }(0)<0`$ since $`G_{\mathrm{con}}(𝐫)r^{2(N1)R^{\prime \prime }(0)}`$ must be limited. The solution of Eq. (5) can be found numerically with shooting. Equation (5) has only a discrete set of the solutions which are everywhere finite. The region of possible initial conditions $`R^{\prime \prime }(0)`$ is limited by the restrictions

$$ϵ/2(N3)R^{\prime \prime }(0)=\eta /2ϵ/2(N2).$$
(7)
The first inequality follows from the Schwartz-Soffer inequality  and is the stability condition for the fixed point . The second one can be derived from the third term of the expansion of the function $`R(\varphi )`$ at small $`\varphi `$:

$$R(\varphi )=\frac{(N1)R^{\prime \prime }(0)^2}{ϵ+4(N2)R^{\prime \prime }(0)}+\frac{R^{\prime \prime }(0)}{2}\varphi ^2\pm \sqrt{\frac{R^{\prime \prime }(0)ϵ+2(N2)R^{\prime \prime }(0)^2}{18(N+2)}}|\varphi |^3+\mathrm{}$$
(8)
Since coefficients of Eq. (5) are singular at $`\varphi =0,\pi `$, one has to use the expansions of $`R(\varphi )`$ in powers of $`|\varphi |,(\pi \varphi )`$ near those points at solving numerically the RG equation. Note that the function $`R(\varphi )`$ is non-analytical at $`\varphi =0`$. We have solved the RG equation at $`N5`$. It turns out that for $`5N3`$ Eq. (5) has only one solution compatible with the condition (7). The critical exponents calculated for those $`N`$ are given in Table I. One can see that the critical exponents are different from the dimensional reduction values $`\eta (4+ϵ)=\overline{\eta }(4+ϵ)=\eta _{\mathrm{pure}}(2+ϵ)=ϵ/(N2)`$, where $`\eta _{\mathrm{pure}}`$ is taken from Ref. .

The second of the inequalities (7) shows that for the XY model $`(N=2)`$ there is no appropriate solution of the RG equation. This can be explained by the fact  that the RF XY model possesses quasi-long-range order in dimension 3. Hence, it has a phase transition in 3 dimensions and thus its lower critical dimension near which the $`ϵ`$-expansion should be developed is less than 4.

It was claimed  that in the first order in $`ϵ`$ the critical exponents of the RF ferromagnet in $`4+ϵ`$ dimensions are the same as in the pure O$`(N)`$ model in $`2+ϵ`$ dimensions. That prediction could be obtained from Eq. (5) with the wrong assumption that there is an analytical solution $`R(\varphi )`$ controlling the critical point.

We have demonstrated the existence of an infinite set of relevant operators near 4 dimensions in the RF O$`(N)`$ model. A possible explanation of the failure of the ($`6ϵ`$)-expansion is thus the appearance of additional relevant operators below some dimension $`D<6`$. However, we show that another scenario takes place: an infinite set of relevant operators emerges in any dimension less than the upper critical dimension 6. Since the upper critical dimension is the same for the Ising and O$`(N)`$ models our approach allows us to consider both cases in the same way.

Near 6 dimensions it is convenient to use the Ginzburg-Landau model with the random field $`𝐡(𝐱)`$

$$H=d^Dx[(𝐦)^2+g(𝐦^2)^2𝐡(𝐱)𝐦]$$
(9)
that can be described by the replica Hamiltonian

$$H_R=d^Dx[\underset{a}{}(𝐦_a)^2+\underset{a}{}g(𝐦_a^2)^2\mathrm{\Delta }\underset{ab}{}𝐦_a𝐦_b],$$
(10)
where $`𝐦_a`$ are the replicas of the $`N`$-component magnetization. The standard power counting suggests that all operators which are relevant in $`6ϵ`$ dimensions are included into Eq. (10). However, if it were so the dimensional reduction would be correct. The hint why the power counting fails is given by the theory of the metal-dielectric transition. It was argued that in the nonlinear sigma-model for the metal-dielectric transition one should include an infinite set of relevant operators  missed by the power counting. For any eigen-operator $`A_k`$ of the RG transformation the RG equation has the structure

$$\frac{dA_k}{d\mathrm{ln}L}=(C_0^k+C_1^kϵ^s+o(ϵ^s))A_k+O(A^2).$$
(11)
The power counting is based on the sign of the constant $`C_0`$ and predicts that the operator $`A_k`$ is irrelevant as $`C_0^k<0`$. However, if the ratio $`C_1^k/C_0^k`$ grows up to infinity at large $`k`$ then for any fixed $`ϵ`$ there is such $`k`$ that the correction $`C_1^kϵ^s`$ is greater than $`C_0^k`$. Hence, the role of the operators $`A_k`$ depends on the signs of $`C_1^k`$. In the $`\varphi ^4`$-theory without the random field this sign is negative for any operator which is expected to be irrelevant from the power counting . This agrees with the success of the RG approach in that problem. On the other hand, we shall see that the signs of these constants are positive for some operators in the RF problem. This signals that additional relevant operators emerge.

For the pure $`\varphi ^4`$-model one can define the rank of an operator with $`q`$ derivatives and the $`p`$th power of the order parameter as $`r=p+q4`$. Since only the operators of the same rank can mix , for the calculation of the anomalous dimension $`C_1^kϵ^s`$ of an operator with rank $`r_k`$ one can ignore the diagrams which produce the operators of the higher ranks. It is easy to check that the definition of the rank has to be modified in the RF problem to ensure that the operators of different ranks do not mix: $`r=p+q+2t6`$, where $`t`$ is the number of the different replica indices in the operator.

To avoid the difficult problem of operator mixing we guess a set of relevant operators which do not mix with the other operators of the same rank up to the second order in $`ϵ`$. It turns out that the following operators are relevant and do not mix with the other: $`A_k=_{ab}[(𝐦_a𝐦_b)^2]^k`$. At $`k>1`$ their canonical dimensions

$$d_c^k=(4ϵ)k6+ϵ$$
(12)
are positive and proportional to $`k`$. We shall see that the anomalous dimensions are negative and proportional to $`(kϵ)^2`$. If one imposes the fixed length restriction $`𝐦^2=1`$ so that the Ginzburg-Landau model reduces to the O$`(N)`$ model then the operators $`A_k`$ can be represented as the linear combinations of the random anisotropies $`B_k=_{ab}(𝐦_a𝐦_b)^k`$. This is natural since the operators $`B_k`$ are relevant in $`4+ϵ`$ dimensions as one could see above.

To the first order in the quartic vertex $`g`$ the anomalous dimensions of the operators $`A_k`$ are zero. To prove this we demonstrate that all diagrams with one vertex $`gm_a^4`$ and one vertex $`A_k^{ab}=[(𝐦_a𝐦_b)^2]^k`$ either are equal to zero or produce operators of higher ranks. Indeed, any line of the Feynman diagram corresponds to the product of a momentum-dependent factor and two operators which differentiate the vertices at the ends of the line with respect to replicas $`𝐦_a`$ of the magnetization. For example, line 2 in Fig. 1 acts on the vertex $`A_k^{ab}`$ as the differential operator $`_2=/𝐦_a`$. Line 1 in which the vertex $`\mathrm{\Delta }`$ is inserted acts on the vertex $`A_k^{ab}`$ as the operator $`_1=/𝐦_a+/𝐦_b`$. Obviously, $`_1A_k^{ab}=0`$. Hence, any diagram including line 1 is equal to zero. On the other hand, any diagram with one vertex $`g`$, one vertex $`A_k`$, and without lines in which $`\mathrm{\Delta }`$ is inserted produces an operator of a higher rank. Fig. 2 shows the only nonzero diagram of the order $`g^2`$ that does not increase the rank of the operators $`A_k`$. Calculating this diagram one obtains the anomalous dimension

$$d_{\mathrm{an}}^{k,N}=k(N+2)ϵ^2/2(N+8)^2ϵ^2[Nk(2k+1)+16k^210k]/2(N+8)^2$$
(13)
The anomalous dimension (13) is negative and proportional to $`k^2`$. Thus, at any fixed $`ϵ`$ one expects that the operators $`A_k`$ with $`1/ϵ^2k`$ are relevant but missed by the existing theoretical methods. Although the rigorous analysis requires consideration of all orders in $`ϵ`$, the appearance of an infinite set of relevant operators in $`6ϵ`$ dimensions is plausible since such set exists near 4 dimensions and since the alternative explanation  of the failure of RG due to non-perturbative corrections is hardly compatible with the existing numerical results .

In conclusion, we have calculated critical exponents of the O$`(N)`$ model in $`4+ϵ`$ dimensions and demonstrated that they do not obey the dimensional reduction. The failure of the dimensional reduction is related to the appearance of an infinite set of relevant operators.

I thank Y. Gefen, I.V. Lerner, V.L. Pokrovsky, and M.G. Stepanov for useful discussions and R. Whitney for the critical reading of the manuscript. This work was supported by the Koshland fellowship and RFBR grant No. 00-02-17763.


text:
# Numerical Evidence that the Perturbation Expansion for a Non-Hermitian 𝒫𝒯-Symmetric Hamiltonian is Stieltjes

## 1 Introduction

Hamiltonians describing fundamental interactions traditionally possess two symmetries, the continuous symmetry of the proper Lorentz group and the discrete symmetry of Hermiticity. Lorentz invariance is a physical requirement. Hermiticity is a useful mathematical constraint that guarantees that the spectrum is real, although recent work shows that Hermiticity is only a sufficient condition and is not necessary for the reality of eigenvalues. From the assumptions of Lorentz invariance and positivity of the spectrum of the Hamiltonian one can prove the $`𝒫𝒞𝒯`$ theorem and thereby establish the physical symmetry of $`𝒫𝒞𝒯`$ invariance.

What happens if we impose only the more physical symmetries of Lorentz invariance and $`𝒫𝒞𝒯`$ invariance when we construct a Hamiltonian? The constraint of $`𝒫𝒞𝒯`$ invariance is weaker than Hermiticity, so Hamiltonians having this property need not be Hermitian. While it has not been proved, there is compelling analytical and numerical evidence supporting the conjecture that, except when $`𝒫𝒞𝒯`$ symmetry is spontaneously broken, the energy levels of such Hamiltonians are all real and positive. The reality and positivity of the spectrum is apparently a consequence of the $`𝒫𝒞𝒯`$ symmetry of $`H`$.

Many examples of $`𝒫𝒞𝒯`$-symmetric Hamiltonians in quantum field theory have been studied. In quantum mechanics, where the $`𝒞`$ operator is unity, many examples of $`𝒫𝒯`$-symmetric Hamiltonians have also been studied. A simple example of such a quantum-mechanical Hamiltonian is $`H=p^2+ix^3`$. Hamiltonians such as this may be regarded as complex deformations of conventional Hermitian Hamiltonians. To understand this deformation we consider the Hamiltonian

$$H=p^2(ix)^{2+ϵ},$$

where $`ϵ0`$. When $`ϵ=0`$, we have the conventional harmonic oscillator Hamiltonian, whose spectrum is real and positive. As $`ϵ`$ increases from $`0`$, the entire spectrum of the Hamiltonian smoothly deforms as a function of $`ϵ`$ and remains real and positive for all positive values of $`ϵ`$. Thus, these theories are in effect the analytic continuation of conventional quantum mechanics into the complex plane.

These non-Hermitian theories exhibit some remarkable properties. Most interesting is that the expectation value of the operator $`x`$ in quantum mechanics (and of the field $`\varphi `$ in quantum field theory) is nonzero when $`ϵ>0`$. This is true even for the $`p^2x^4`$ Hamiltonian that one obtains at $`ϵ=2`$ and it is also true for the $`g\varphi ^4`$ scalar quantum field theory. The $`g\varphi ^4`$ quantum field theory is particularly surprising because it has a positive real spectrum and exhibits a nonzero value of $`\varphi `$. In four-dimensional space-time it has a dimensionless coupling constant, is renormalizable, and is asymptotically free (and thus nontrivial). It may thus provide a useful setting to describe the Higgs particle.

We are struck by the close similarity between the properties of non-Hermitian $`𝒫𝒯`$-symmetric quantum-mechanical Hamiltonians and conventional Hermitian Hamiltonians. Moreover, in mathematical terms, we are struck by the strong resemblance between self-adjoint Sturm-Liouville problems and these new complex Sturm-Liouville problems. The purpose of this paper is to present further evidence of this strong similarity by investigating various aspects of Padé summation and Padé prediction of the Rayleigh-Schrödinger perturbation series for the ground-state energy of the complex $`𝒫𝒯`$-symmetric Hamiltonian

$$H(\lambda )=p^2+\frac{1}{4}x^2+i\lambda x^3.$$
(1)
Note that this Hamiltonian is $`𝒫𝒯`$ symmetric because under parity reflection $`𝒫:pp`$ and $`𝒫:xx`$ and under time reversal, which is an antiunitary operation, $`𝒯:pp`$, $`𝒯:xx`$, and $`𝒯:ii`$.

The large-order behavior of the divergent Rayleigh-Schrödinger perturbation series

$$E_0(\lambda )\frac{1}{2}+\underset{n=1}{\overset{\mathrm{}}{}}b_n\lambda ^{2n}(\lambda 0^+),$$
(2)
for the ground-state energy eigenvalue of the Hamiltonian in (1) has already been examined in , where the first 46 terms of the perturbation expansion had been generated using recursion formulas. It was observed there that the coefficients $`b_n`$ are all integers, that they alternate in sign, and that their magnitude grows rapidly with $`n`$. The first 10 coefficients are listed in Table 1.

To calculate the coefficients $`b_n`$ we make the ansatz that the wave function is a formal series in powers of the coupling constant $`\lambda `$ and that the coefficient of $`\lambda ^n`$ has the form of a Gaussian $`\mathrm{exp}(x^2/4)`$ times a polynomial of degree $`3n`$ in the variable $`x`$. The eigenvalue $`E_0(\lambda )`$ automatically appears as a series in powers of $`\lambda ^2`$. Thus, for each additional coefficient in the series for $`E_0`$ it is necessary to calculate two orders in powers of $`\lambda `$ for the wave function.

In  it was pointed out that the Hamiltonian (1) describes a $`0+1`$ dimensional $`\varphi ^3`$ field theory and that $`\varphi ^3`$ theories were the first quantum field theories in which the divergences of perturbation theory were studied . Using the standard methods for determining the large-order behavior of perturbation theory  it can be shown that the leading large-$`n`$ behavior of the coefficients $`b_n`$ is given by the formula

$$b_n=(1)^{n+1}\mathrm{\hspace{0.17em}60}^{n+1/2}(2\pi )^{3/2}\mathrm{\Gamma }(n+1/2)\left[1+\mathrm{O}\left(1/n\right)\right](n\mathrm{}).$$
(3)
This asymptotic behavior was verified numerically in . There, it was also shown that the first correction term to this leading asymptotic behavior, which is proportional to $`1/n`$, is negative. Although divergent, the series in (2) is Borel summable . If the factor of $`i`$ were absent from the Hamiltonian (1), then the perturbation coefficients $`b_n`$ would not alternate in sign and the perturbation series would not be Borel summable.

It is interesting that to derive the asymptotic formula in (3) one must use dispersion-relation techniques that rely on crucial assumptions about the analyticity of the function $`E_0(\lambda )`$. These assumptions are justified for the Hermitian Hamiltonian of the anharmonic oscillator ,

$$(\beta )=p^2+x^2+\beta x^4.$$
(4)
However, the validity of these assumptions is unproved for the non-$`𝒫𝒯`$-symmetric $`\lambda x^3`$ oscillator.

The eigenvalues $`(\beta )`$ of the quartic anharmonic oscillator possess several other properties which closely resemble those of the eigenvalues $`E(\lambda )`$ of the $`𝒫𝒯`$-symmetric Hamiltonian (1). For example, the ground-state energy eigenvalue $`_0(\beta )`$ of the quartic anharmonic oscillator possesses a divergent weak-coupling perturbation expansion, which also diverges factorially :

$`_0(\beta )`$ $``$ $`1+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}_n\beta ^n(\beta 0^+),`$ (5)
$`_n`$ $`=`$ $`(1)^{n+1}{\displaystyle \frac{4}{\pi ^{3/2}}}\left({\displaystyle \frac{3}{2}}\right)^{n+1/2}\mathrm{\Gamma }(n+1/2)\left[1+\mathrm{O}\left(1/n\right)\right](n\mathrm{}).`$ (6)
A comparison of the large-order asymptotics (3) and (6) shows that the two perturbation expansions (2) and (5) possess the same rate of divergence if we choose

$$\beta =\mathrm{\hspace{0.33em}40}\lambda ^2.$$
(7)
In view of these striking similarities between the ground-state eigenvalues $`E_0(\lambda )`$ and $`_0(\beta )`$ it should be interesting to investigate what other similarities do exist. In particular, we are interested in similarities that could provide evidence that the Padé summation of the divergent perturbation series (2) for the ground-state energy of the $`𝒫𝒯`$-symmetric Hamiltonian (1) converges.

In the case of the ground-state energy shift $`\mathrm{\Delta }_0(\beta )`$ of the quartic anharmonic oscillator, which is defined by

$$_0(\beta )=\mathrm{\hspace{0.33em}1}+\beta \mathrm{\Delta }_0(\beta ),$$
(8)
it was shown rigorously by Simon (Theorem IV.2.1 of ) that the corresponding perturbation series is a Stieltjes series. This has some far-reaching consequences. In the case of Stieltjes series, Padé approximants possess a highly developed convergence theory, as we discuss in detail in Section 2. In particular, the Stieltjes nature of the perturbation series (5) guarantees that certain subsequences of the Padé table converge to a uniquely determined Stieltjes function.

Although we cannot prove it rigorously, we believe that for the $`𝒫𝒯`$-symmetric Hamiltonian (1) the ground-state energy shift $`\mathrm{\Delta }E_0(\lambda ^2)`$ defined by

$$E_0(\lambda )=\frac{1}{2}+\lambda ^2\mathrm{\Delta }E_0(\lambda ^2),$$
(9)
in which $`E_0`$ is considered as a function of $`\lambda ^2`$, is also a Stieltjes function. This implies that the corresponding perturbation series is a Stieltjes series.

It is the intention of this article to provide numerical evidence supporting this conjecture. We do this by comparing Padé summations and Padé predictions (see  and references therein) of the perturbation expansions for the ground-state energy shift $`\mathrm{\Delta }E_0(\lambda ^2)`$ and for the analogous ground-state energy shift $`\mathrm{\Delta }_0(\beta )`$ of the quartic anharmonic oscillator.

Here, one might argue that one should also investigate the summation of the perturbation expansion for the ground-state energy shift $`\mathrm{\Delta }E_0(\lambda ^2)`$ with the help of the sequence transformations that were described in Section 7 and 8 of  and which produced very good results in the case of the anharmonic oscillators . However, the convergence theory of these sequence transformations, which in the case of power series also produce rational approximants, is still very much in its infancy and no theoretical results concerning the transformation of Stieltjes series are known so far. Consequently, we would only produce numbers without gaining any further mathematical insight.

In Section 2, we present the relevant details of Padé approximants, the computation of Padé approximants by means of Wynn’s recursive epsilon algorithm , and we discuss Stieltjes series and their associated Stieltjes functions. In Section 3, we show that the Padé summation of the perturbation expansions for $`\mathrm{\Delta }E_0(\lambda ^2)`$ and $`\mathrm{\Delta }_0(\beta )`$, respectively, produces results of identical quality if the two coupling constants $`\lambda `$ and $`\beta `$ satisfy (7). In Section 4, we discuss the prediction of unknown perturbation coefficients with the help of Wynn’s epsilon algorithm, and we show that the coefficients of the perturbation expansions for either $`\mathrm{\Delta }E_0(\lambda ^2)`$ or $`\mathrm{\Delta }_0(\beta )`$ can be predicted equally well. Finally, in Section 5 we give a brief summary.

## 2 Padé Approximants, Wynn’s Epsilon Algorithm, and Stieltjes Series

In recent years, Padé approximants have become the standard tool in theoretical physics to overcome problems with slowly convergent or divergent power series. Accordingly, there is a vast literature on the mathematical properties of Padé approximants as well as on their applications in theoretical physics. Any attempt to provide a reasonably complete bibliography would be beyond the scope of this article (see for example the extensive bibliography compiled by Brezinski ). We just mention that the popularity of Padé approximants in theoretical physics can be traced back to a review by Baker , that the first specialized monograph on Padé approximants is due to Baker , and that currently the most complete source of information is the second edition of the monograph by Baker and Graves-Morris . In addition to treatments in more mathematically oriented books on continued fractions and related topics , Padé approximants are also discussed in book on mathematical and theoretical physics, for example in Section 8 of the book by Bender and Orszag , or in Part III of a book by Baker on critical phenomena . Then, there is a book by Pozzi  on the use of Padé approximants in fluid dynamics. Finally, there is even a monograph  and two articles  on the history of Padé approximants and related topics.

A Padé approximant $`P_m^l(z)`$ to a function $`f`$ possessing a (formal) power series expansion

$$f(z)=\underset{\nu =0}{\overset{\mathrm{}}{}}\gamma _\nu z^\nu ,$$
(10)
which may converge or diverge, is the ratio of two polynomials $`A_l(z)`$ and $`B_m(z)`$ of degrees $`l`$ and $`m`$ in $`z`$ (p. 383 of ):

$$P_m^l(z)=\frac{A_l(z)}{B_m(z)}=\frac{a_0+a_1z+a_2z^2+\mathrm{}a_lz^l}{1+b_1z+b_2z^2+\mathrm{}b_mz^m}.$$
(11)
An alternative notation for Padé approximants, which is used in the books by Baker and Graves-Morris , is $`P_m^l(z)=[l/m]_f(z)`$. This notation is usually simplified further to $`P_m^l(z)=[l/m]`$ if explicit references to $`f`$ and $`z`$ are not necessary.

The coefficients $`a_\lambda `$ and $`b_\mu `$ of the polynomials $`A_l(z)`$ and $`B_m(z)`$ in (11) are chosen in such a way that the Taylor expansion of $`f(z)`$ and of its Padé approximant agree as far as possible:

$$f(z)P_m^l(z)=\mathrm{O}\left(z^{l+m+1}\right)(z0).$$
(12)
This asymptotic error estimate leads to a system of linear equations by means of which the coefficients $`a_0`$, $`a_1`$, …, $`a_l`$ and $`b_1`$, $`b_2`$, …, $`b_m`$ in (11) can be computed . Moreover, several algorithms are known which permit a recursive computation of Padé approximants. A discussion of the merits and weaknesses of the various computational schemes can for instance be found in Section II.3 of the book by Cuyt and Wuytack .

Probably, the best known recursive algorithm for Padé approximants is Wynn’s epsilon algorithm :

$`\epsilon _1^{(n)}`$ $`=`$ $`0,\epsilon _0^{(n)}=s_n,n_0,`$ (13a)
$`\epsilon _{k+1}^{(n)}`$ $`=`$ $`\epsilon _{k1}^{(n+1)}+\mathrm{\hspace{0.17em}1}/[\epsilon _k^{(n+1)}\epsilon _k^{(n)}],k,n_0.`$ (13b)
A compact FORTRAN program for the epsilon algorithm as well as the underlying computational algorithm is described in Section 4.3 of .

If the input data $`\epsilon _0^{(n)}=s_n`$ of Wynn’s epsilon algorithm are the partial sums

$$f_n(z)=\underset{\nu =0}{\overset{n}{}}\gamma _\nu z^\nu $$
(14)
of the formal power series (10) according to $`\epsilon _0^{(n)}=f_n(z)`$, then the elements $`\epsilon _{2k}^{(n)}`$ with *even* subscripts are Padé approximants to $`f(z)`$ according to

$$\epsilon _{2k}^{(n)}=P_k^{k+n}(z).$$
(15)
The elements $`\epsilon _{2k+1}^{(n)}`$ with *odd* subscripts are only auxiliary quantities, which diverge if the whole transformation process converges.

The epsilon algorithm is a useful numerical algorithm that is applied successfully in a large variety of different fields. Accordingly, there is an extensive literature dealing with it. A fairly complete coverage of the older literature can be found in a book by Brezinski . It may be interesting to note that the epsilon algorithm is not restricted to scalar sequences but can be generalized to cover vector sequences. A recent review of these developments can be found in .

If one tries to sum a divergent power series by converting its partial sums (14) to Padé approximants, it is usually a good idea to use *diagonal* Padé approximants, whose numerator and denominator polynomials have equal degrees. If this is not possible one should use Padé approximants with numerator and denominator polynomials whose degrees differ as little as possible. If we use the epsilon algorithm for the computation of the Padé approximants, then (15) implies that we should use the elements of following staircase sequence in the Padé table as approximations to $`f(z)`$ (see Eq. (4.3-7) of ):

$$P_0^0(z),P_0^1(z),P_1^1(z),\mathrm{},P_\nu ^\nu (z),P_\nu ^{\nu +1}(z),P_{\nu +1}^{\nu +1}(z),\mathrm{}.$$
(16)
This staircase sequence exploits the available information optimally if the partial sums $`f_n(z)`$ with $`n0`$ are computed successively and if, after the computation of each new partial sum, the element of the epsilon table with the highest possible even subscript is computed. With the help of the notation $`[[x]]`$ for the integral part of $`x`$, this staircase sequence can be written compactly as follows:

$$\epsilon _{2[[n/2]]}^{(n2[[n/2]])}=P_{[[n/2]]}^{n[[n/2]]}(z)(n=0,1,2,\mathrm{}).$$
(17)
As remarked above, Padé approximants are now used almost routinely to overcome problems with slowly convergent or divergent power series. Hence, their practical usefulness is beyond question. However, from a theoretical point of view, the situation is not so good. So far, a completely satisfactory *general* convergence theory of Padé approximants for essentially arbitrary power series does not exist.

Nevertheless, there is a special class of series, the so-called *Stieltjes* series, which possess a highly developed and elegant convergence theory. In this section we will only discuss those properties of Stieltjes series and Stieltjes functions that are needed to provide numerical evidence that the perturbation expansion for the ground-state energy shift $`\mathrm{\Delta }E_0(\lambda ^2)`$ of the complex $`𝒫𝒯`$-symmetric Hamiltonian (1), considered as a function of $`\lambda ^2`$, is a Stieltjes series. Detailed discussions of the properties of Stieltjes series and their special role in the theory of Padé approximants can be found in Section 8.6 of  or in Section 5 of .

A function $`F(z)`$ with $`z`$ is called a *Stieltjes* function if it can be expressed as a Stieltjes integral according to

$$F(z)=_0^{\mathrm{}}\frac{\mathrm{d}\mathrm{\Phi }(t)}{1+zt}[|\mathrm{arg}(z)|<\pi ].$$
(18)
Here, $`\mathrm{\Phi }(t)`$ is a bounded, nondecreasing function taking infinitely many different values on $`0t<\mathrm{}`$. Moreover, the moment integrals

$$\mu _n=_0^{\mathrm{}}t^nd\mathrm{\Phi }(t)(n_0)$$
(19)
must be real and finite for all finite values of $`n`$. A Stieltjes function can be expressed by its corresponding *Stieltjes* series:

$$F(z)=\underset{\nu =0}{\overset{\mathrm{}}{}}(1)^\nu \mu _\nu z^\nu .$$
(20)
Whether this series converges or diverges depends on the behavior of the Stieltjes moments $`\mu _n`$ as $`n\mathrm{}`$.

In a typical Stieltjes summation problem, as it occurs in the context of divergent perturbation expansion, only the numerical values of a finite number of Stieltjes moments $`\mu _n`$ are known. Thus, one has to find a way of constructing an approximation to the unknown Stieltjes function $`F(z)`$ from a finite string of moments.

Of course, one would also like to have some theoretical evidence that $`F(z)`$ exists and is uniquely determined by the Stieltjes moments $`\{\mu _n\}_{n=0}^{\mathrm{}}`$. Many necessary and sufficient conditions that guarantee this are known in literature. Unfortunately, it is by no means easy to apply them. For example, a necessary condition that the series (20) is indeed a Stieltjes series is that the Hankel determinants

$$D(m,n)=\left|\begin{array}{cccc}\mu _m& \mu _{m+1}& \mathrm{}& \mu _{m+n}\\ \mu _{m+1}& \mu _{m+2}& \mathrm{}& \mu _{m+n+1}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mu _{m+n}& \mu _{m+n+1}& \mathrm{}& \mu _{m+2n}\end{array}\right|$$
(21)
are positive for all $`m,n0`$ (see Theorem 5.1.2 on p. 197 of ).

Then, there is a sufficient criterion, the so-called Carleman condition (see p. 410 of  or pp. 239 - 240 of ), which requires that the series $`_{j=1}^{\mathrm{}}(\mu _j)^{1/(2j)}`$ diverges and thus limits the admissible growth of the moments $`\mu _n`$ as $`n\mathrm{}`$ . If the condition $`D(m,n)>0`$ on the determinants as well as the Carleman condition are both satisfied, then the Padé approximants $`P_m^{m+j}(z)`$ constructed from the partial sums of the moment expansion (20) converge for every $`j1`$ to the corresponding Stieltjes function $`F(z)`$ as $`m\mathrm{}`$ (see for example Theorem 5.5.1 on p. 240 of ).

If only a finite number of moments are known, it is impossible to prove that $`D(m,n)>0`$ holds for all $`m,n0`$, and it is also not possible to prove rigorously that the Carleman condition is satisfied, although we would like to emphasize that the large-order formula (3), which was verified numerically in , is in agreement with the Carleman condition.

Therefore, we prefer an indirect approach in order to provide evidence that the perturbation series for the energy shift $`\mathrm{\Delta }E_0(\lambda ^2)`$ of the $`𝒫𝒯`$-symmetric Hamiltonian (1) is indeed a Stieltjes series. For that purpose, let us assume that the moment expansion (20), whose Stieltjes nature we want to establish, is a Stieltjes series. Padé approximants to Stieltjes series possess a highly developed convergence theory, and many conditions and inequalities are known that Padé approximants to a Stieltjes series must satisfy. For example, Padé approximants constructed from the partial sums

$$F_n(z)=\underset{\nu =0}{\overset{n}{}}(1)^\nu \mu _\nu z^\nu $$
(22)
of the moment expansion (20) for a Stieltjes function $`F(z)`$ satisfy for $`z>0`$ the following inequalities (Theorem 15.2 on p. 215 of ):

$`(1)^{j+1}\left\{P_{m+1}^{m+j+1}(z)P_m^{m+j}(z)\right\}`$ $`0,`$ (23)
$`(1)^{j+1}\left\{P_m^{m+j}(z)P_{m1}^{m+j+1}(z)\right\}`$ $`0,`$ (24)
$`P_m^m(z)F(z)P_m^{m1}(z).`$ (25)
It follows from inequality (23) that the Padé sequence $`\left\{P_m^{m+j}(z)\right\}_{m=0}^{\mathrm{}}`$ is *increasing* for $`z>0`$ if $`j`$ is *odd*, and it is *decreasing* if $`j`$ is *even*. Moreover, if we set $`j=1`$ in (24) and replace $`m`$ by $`m+1`$, we obtain the inequality

$$P_{m+1}^m(z)P_m^{m+1}(z)(m_0).$$
(26)
Thus, if we use Wynn’s epsilon algorithm (13) to convert the partial sums (22) to Padé approximants and choose the approximants to $`F(z)`$ according to (17), then it follows from inequalities (25) and (26) that these Padé approximants satisfy the following inequality if the moment expansion (20) for $`F(z)`$ is a Stieltjes series:

$$P_m^{m+1}(z)F(z)P_{m+1}^{m+1}(z)(z>0,m_0).$$
(27)
Thus, the approximants (17) produced by Wynn’s epsilon algorithm yield for $`z>0`$ two nesting sequences $`P_m^{m+1}(z)=\epsilon _{2m}^{(1)}`$ and $`P_{m+1}^{m+1}(z)=\epsilon _{2m+2}^{(0)}`$ of lower and upper bounds to $`F(z)`$ if the the moment expansion (20) is a Stieltjes series.

If only the numerical values of a finite number of Stieltjes moments are available, then it is of course not possible to prove rigorously that the series under consideration is a Stieltjes series. Nevertheless, we can provide considerable evidence that this hypothesis is true if inequality (27) is valid in all cases that can be checked.

## 3 Summation Results

In this section we want to show that the Padé summation of the perturbation expansion

$$\mathrm{\Delta }E_0(\lambda ^2)=\underset{\nu =0}{\overset{\mathrm{}}{}}b_{\nu +1}\lambda ^{2\nu }$$
(28)
for the ground-state energy shift of the of the $`𝒫𝒯`$-symmetric Hamiltonian (1) and of the perturbation expansion

$$\mathrm{\Delta }_0(\beta )=\underset{\nu =0}{\overset{\mathrm{}}{}}_{\nu +1}\beta ^\nu $$
(29)
for the ground-state energy shift of the quartic anharmonic oscillator yield results of virtually identical quality if the two coupling constants $`\lambda `$ and $`\beta `$ satisfy (7). Moreover, we want to demonstrate numerically that inequality (27), which is satisfied in the case of the Padé summation a Stieltjes series, is apparently also satisfied. In our summation calculations, we used all coefficients $`b_\nu `$ and $`_\nu `$ with $`0\nu 193`$ which had been computed recursively.

In this article, we compute all Padé approximants with the help of Wynn’s epsilon algorithm (13). Thus, the partial sums

$$s_n(\lambda ^2)=\underset{\nu =0}{\overset{n}{}}b_{\nu +1}\lambda ^{2\nu }$$
(30)
and

$$\sigma _n(\beta )=\underset{\nu =0}{\overset{n}{}}_{\nu +1}\beta ^\nu $$
(31)
are used as input data for Wynn’s epsilon algorithm according to $`\epsilon _0^{(n)}=s_n(\lambda ^2)`$ or $`\epsilon _0^{(n)}=\sigma _n(\beta )`$, respectively, and the approximations to the energy shifts are chosen according to (17).

In Table 2 we present illustrative results of the Padé summation of the perturbation expansion (28) for the ground-state energy shift $`\mathrm{\Delta }E_0(\lambda ^2)`$ with $`\lambda =1/7`$ and of the perturbation expansion (29) for the ground-state energy shift $`\mathrm{\Delta }_0(\beta )`$ with $`\beta =40/49`$. Thus, the two coupling constants $`\lambda `$ and $`\beta `$ satisfy (7). This implies that the two perturbation expansions, whose partial sums are displayed in columns 2 and 3, should show the same rate of divergence.

Here, we must remember that the larger-order estimates (3) and (6), respectively, imply that the partial sums (30) for the ground-state energy shift $`\mathrm{\Delta }E_0(\lambda ^2)`$ of the $`𝒫𝒯`$-symmetric Hamiltonian (1) are for all $`n0`$ at least one order of magnitude greater than the partial sums (31) for the ground-state energy shift $`\mathrm{\Delta }_0(\beta )`$ of the quartic anharmonic oscillator. Otherwise, the observed rates of divergence in columns 2 and 3 are virtually identical.

Moreover, the Padé approximants in columns 4 and 5 apparently satisfy the inequality (27), which holds if the series to be transformed is a Stieltjes series. If the index $`n`$ in column 1 is even ($`n=2m`$), then the diagonal Padé approximants $`P_m^m(\lambda ^2)`$ and $`P_m^m(\beta )`$ provide upper bounds that strictly decrease with increasing $`m`$ and, if $`n`$ is odd ($`n=2m+1`$), then the Padé approximants $`P_m^{m+1}(\lambda ^2)`$ and $`P_m^{m+1}(\beta )`$ provide lower bounds that strictly increase with increasing $`m`$.

We have done analogous summation calculations also for many other values of the coupling constants $`\lambda `$ and $`\beta `$. Of course, the performance of the Padé summations depend very much on the size of the coupling constants. For smaller values of $`\lambda `$ and the corresponding $`\beta `$, convergence is better than in Table 2, whereas for $`\lambda =1`$ and $`\beta =40`$ only the first digit of the summation results stabilize. For larger values of $`\lambda `$ and $`\beta `$, Padé summation produces only relatively crude upper and lower bounds. However, to emphasize that the typical qualitative features of the summation results in Table 2 – the same rate of divergence of two perturbation series and the occurrence of strictly decreasing upper bounds $`P_m^m`$ and strictly increasing lower bounds $`P_m^{m+1}`$ – were consistently observed in all cases considered.

Thus, Wynn’s epsilon algorithm is apparently unable to detect any substantial difference between the perturbation series (29) for $`\mathrm{\Delta }_0(\beta )`$, whose Stieltjes nature was established rigorously by Simon (Theorem IV.2.1 of ), and the perturbation series (28) for $`\mathrm{\Delta }E_0(\lambda ^2)`$, whose Stieltjes nature we conjecture.

## 4 Padé Predictions

As shown by countless articles from all branches of physics, Padé approximants have become the standard tool to overcome problems with slowly convergent or divergent power series. However, Padé approximants have other useful features that are not as well known yet. For example, Padé approximants can be used to make predictions for higher-order series coefficients that were not used for the construction of the approximant.

On a heuristic level the prediction capability of Padé approximants, which was apparently first noted and used by Gilewicz , can be explained quite easily. Let us assume that the partial sums (14) of the power series for some function $`f(z)`$ are to be converted to Padé approximants. Then, the *accuracy-through-order* relationship (12) implies that a Padé approximant $`P_m^l(z)`$ to $`f(z)`$ can be expressed as the partial sum $`f_{l+m}(z)`$ from which it was constructed plus a term $`z^{l+m+1}_m^l(z)`$, which was generated by the transformation of the partial sum to the rational approximant:

$$P_m^l(z)=\underset{\nu =0}{\overset{l+m}{}}\gamma _\nu z^\nu +z^{l+m+1}_m^l(z)=f_{l+m}(z)+z^{l+m+1}_m^l(z).$$
(32)
Similarly, the power series (10) can be expressed as follows:

$$f(z)=\underset{\nu =0}{\overset{l+m}{}}\gamma _\nu z^\nu +z^{l+m+1}\underset{\nu =0}{\overset{\mathrm{}}{}}\gamma _{l+m+\nu +1}z^\nu =f_{l+m}(z)+z^{l+m+1}_{l+m+1}(z).$$
(33)
Let us now assume that the indices $`l`$ and $`m`$ so large that the Padé approximant $`P_m^l(z)`$ provides a sufficiently accurate approximation to $`f(z)`$. Then, the Padé transformation term $`_m^l(z)`$ must also provide a sufficiently accurate approximation to the truncation error $`_{l+m+1}(z)`$ of the power series. In general, we have no reason to assume that

$$_m^l(z)=_{l+m+1}(z)$$
(34)
might hold *exactly* for finite values of $`l`$ and $`m`$. Consequently, Taylor expansions of $`_m^l(z)`$ and $`_{l+m+1}(z)`$, respectively, will in general produce different results. Nevertheless, the *leading* coefficients of the Taylor expansion for $`_m^l(z)`$ should in such a case provide sufficiently accurate approximations to the corresponding coefficients of the Taylor series for $`_{l+m+1}(z)`$.

It is important to note that this prediction capability of Padé approximants does not depend on the convergence of the power series expansions for $`_m^l(z)`$ and $`_{l+m+1}(z)`$, respectively, which was used implicitly in our heuristic reasoning given above. Padé approximants are able to make predictions about series coefficients even if the power series (10) for $`f`$ as well as the power series expansions for $`_l^m`$ and $`_{l+m+1}(z)`$ are only asymptotic as $`z0`$. This fact explains why the prediction capability of Padé approximants can be so very useful in the case of violently divergent perturbation expansions (see  and references therein).

Theoretically, very little is known about the prediction of series coefficients that were not used for the construction of the Padé approximant. A notable exception are again Stieltjes series for which inequalities are known.

Let us assume that the partial sums (22) of the moment expansion for some Stieltjes function $`F(z)`$ are to be converted to Padé approximants. Such a Padé approximant $`P_m^l(z)`$ possesses the following power series expansion:

$$P_m^l(z)=\underset{\nu =0}{\overset{\mathrm{}}{}}(1)^\nu \mu _\nu ^{[l/m]}z^\nu .$$
(35)
In Theorem 5.2.7 on p. 220 of  it was shown that for all $`n_0`$ and for all $`lm1`$ the coefficients $`\mu _n^{[l/m]}`$ in (35) are bounded in magnitude by the Stieltjes moments $`\mu _n`$ in (20) according to

$$0\mu _n^{[l/m]}\mu _n.$$
(36)
This inequality can be used to analyze the Stieltjes nature of a moment expansion of the type of (20). With the help of computer algebra systems like Maple or Mathematica it is possible to construct Padé approximants $`P_m^l(z)`$ in an unspecified symbolic variable $`z`$, and this can even be done free of rounding errors if the coefficients of the series to be transformed are exact integers like the coefficients $`b_n`$ in (2) or exact rational numbers like the coefficients $`_n`$ in (5). In the next step, a the leading part of a power series expansion of the Padé approximant must be constructed. The resulting series coefficients can then be compared with the corresponding coefficients of the moment expansion.

Again, this poses no principal problems for computer algebra systems like Maple and Mathematica. However, the accuracy-through-order relationship (12) implies that inequality (36), which is to be checked, is by default satisfied for all indices $`nl+m`$, and only for $`nl+m+1`$ we obtain useful information about the Stieltjes nature of the moment expansion, from which the Padé approximant was constructed. Thus, if $`l+m`$ becomes large, the brute-force approach based on computer algebra systems becomes very demanding both with respect to computer time and memory because it requires both the symbolic construction of complicated Padé approximants and also symbolic differentiations of very high orders.

These computational problems can be simplified considerably with the help of a recently derived recursive scheme (Section 3 of ) that permits a direct calculation of the transformation term $`_m^l(z)`$ in (32) if the corresponding Padé approximant $`P_m^l(z)`$ can be computed with the help of Wynn’s epsilon algorithm according to (15); that is, for Padé approximants of the type $`P_k^{k+n}(z)`$ with $`k,n_0`$.

It follows from the accuracy-through-order relationship (12) in combination with (15) that $`\epsilon _{2k}^{(n)}`$ can be expressed as follows if the partial sums (14) of the power series for some function $`f(z)`$ are used as input data for Wynn’s epsilon algorithm (13):

$$\epsilon _{2k}^{(n)}=f_{n+2k}(z)+z^{n+2k+1}\phi _{2k}^{(n)}(z).$$
(37)
The quantities $`\phi _{2k}^{(n)}(z)`$ can be computed with the help of the recursive scheme in Eq. (3.15) of , which uses the coefficients $`\gamma _n`$ of the power series (10) for $`f(z)`$ as input data, as follows:

$`\phi _0^{(n)}(z)`$ $`=`$ $`0,n_0,`$ (38a)
$`\phi _2^{(n)}(z)`$ $`=`$ $`{\displaystyle \frac{\left[\gamma _{n+2}\right]^2}{\gamma _{n+1}\gamma _{n+2}z}},n_0,`$ (38b)
$`\phi _{2k+2}^{(n)}(z)`$ $`=`$ $`\phi _{2k}^{(n+2)}(z)+{\displaystyle \frac{\alpha _{2k+2}^{(n)}(z)}{\beta _{2k+2}^{(n)}(z)}},k,n_0,`$ (38c)
$`\alpha _{2k+2}^{(n)}(z)`$ $`=`$ $`{\displaystyle \frac{\gamma _{n+2k+2}+\delta \phi _{2k}^{(n+1)}(z)}{\gamma _{n+2k+1}+\delta \phi _{2k}^{(n)}(z)}}{\displaystyle \frac{\gamma _{n+2k+2}+\delta \phi _{2k}^{(n+1)}(z)}{\gamma _{n+2k+1}+z\phi _{2k}^{(n+1)}(z)\phi _{2k2}^{(n+2)}(z)}},`$ (38d)
$`\beta _{2k+2}^{(n)}(z)`$ $`=`$ $`{\displaystyle \frac{1}{\gamma _{n+2k+2}+\delta \phi _{2k}^{(n+1)}(z)}}{\displaystyle \frac{z}{\gamma _{n+2k+1}+\delta \phi _{2k}^{(n)}(z)}}`$ (38e)
$`+{\displaystyle \frac{z}{\gamma _{n+2k+1}+z\phi _{2k}^{(n+1)}(z)\phi _{2k2}^{(n+2)}(z)}},`$
where

$$\delta \phi _{2k}^{(n)}(z)=z\phi _{2k}^{(n+1)}(z)\phi _{2k}^{(n)}(z).$$
(39)
The rational function $`\phi _{2k}^{(n)}(z)`$ can be expressed as a power series in $`z`$ according to

$$\phi _{2k}^{(n)}(z)=g_0^{(n,2k)}+g_1^{(n,2k)}z+g_2^{(n,2k)}z^2+\mathrm{}+g_\nu ^{(n,2k)}z^\nu +\mathrm{}.$$
(40)
The coefficients $`g_\nu ^{(n,2k)}`$ of this series expansion can be used to predict the coefficients $`\gamma _{n+2k+\nu +1}`$ with $`\nu 0`$ of the power series expansion (10) for $`f(z)`$ that were not used for the construction of either $`\epsilon _{2k}^{(n)}`$ or $`\phi _{2k}^{(n)}(z)`$.

Thus, we can compute the rational function $`\phi _{2k}^{(n)}`$ with the help of the recursive scheme (38) in the case of the perturbation series (28) for $`\mathrm{\Delta }E_0(\lambda ^2)`$ and (29) for $`\mathrm{\Delta }_0(\beta )`$. In the case of the $`𝒫𝒯`$-symmetric perturbation series we obtain the expansion

$$\phi _{2k}^{(n)}(\lambda ^2)=b_0^{(n,2k)}+b_1^{(n,2k)}\lambda ^2+b_2^{(n,2k)}\lambda ^4+\mathrm{}+b_\nu ^{(n,2k)}\lambda ^{2\nu }+\mathrm{},$$
(41)
and in the case of the perturbation series for the quartic anharmonic oscillator we obtain

$$\phi _{2k}^{(n)}(\beta )=_0^{(n,2k)}+_1^{(n,2k)}\beta +_2^{(n,2k)}\beta ^2+\mathrm{}+_\nu ^{(n,2k)}\beta ^\nu +\mathrm{}.$$
(42)
The coefficients $`b_\nu ^{(n,2k)}`$ and $`_\nu ^{(n,2k)}`$ with $`\nu 0`$ can be used to obtain predictions for the coefficients $`b_{n+2k+\nu +2}`$ and $`_{n+2k+\nu +2}`$, respectively, that were not used for the construction of the Padé approximants $`\epsilon _{2k}^{(n)}=P_k^{k+n}`$ or the transformation terms $`\phi _{2k}^{(n)}`$.

In Table 3 we compute the rational function $`\phi _{66}^{(0)}(\lambda ^2)`$ corresponding to the diagonal Padé approximant $`\epsilon _{66}^{(n)}=P_{33}^{33}(\lambda ^2)`$ with the help of the recursive scheme (38) from the coefficients $`b_\nu `$ with $`1\nu 67`$. The resulting expansion coefficients $`b_\nu ^{(0,66)}`$ defined in (41) with $`\nu 0`$ provide predictions to the coefficients $`b_{\nu +68}`$ of the perturbation series (28). All calculations for Table 3 were done free of rounding errors using the exact rational arithmetics of Maple. Only in the final step the coefficients were converted to floating point numbers for the sake of readability.

The results in Table 3 show that the expansion coefficients $`b_\nu ^{(0,66)}`$ with $`\nu 0`$ provide already remarkably accurate predictions to the corresponding coefficients $`b_{\nu +68}`$. Moreover, the Stieltjes inequality (36) is satisfied in all cases.

In Table 4 we do the same calculations as in Table 3, but this time for the ground-state energy shift $`\mathrm{\Delta }_0(\beta )`$ of the quartic anharmonic oscillator. Thus, the expansion of the rational function $`\phi _{66}^{(0)}(\beta )`$, which is computed from the coefficients $`_\nu `$ with $`1\nu 67`$ of the perturbation series (29), provides predictions $`_\nu ^{(0,66)}`$ to the coefficients $`_{\nu +68}`$.

A comparison of Tables 3 and 4 shows that their qualitative features are identical. In this context it is quite remarkable that although the coefficients $`b_n`$ grow significantly more rapidly in magnitude than the coefficients $`_n`$, which follows from the large-order estimates (3) and (6), Padé prediction nevertheless yields results of virtually identical quality.

The Padé prediction of unknown series coefficients based on the recursive scheme (38) is certainly computationally simpler than the straightforward symbolic computation and expansion of Padé approximants. The rational function $`\phi _{2k}^{(n)}(z)`$ has a simpler structure than $`\epsilon _{2k}^{(n)}=P_k^{k+n}(z)`$, and the first $`n+2k`$ symbolic differentiations can be avoided. Nevertheless, the recursive symbolic computation of the rational function $`\phi _{2k}^{(n)}(z)`$ from the coefficients $`\gamma _0`$, $`\gamma _1`$, …, $`\gamma _{n+2k}`$ of the power series (10) can become quite demanding, in particular if $`k`$ becomes large.

The problems connected with the computation of $`\phi _{2k}^{(n)}(z)`$ can largely be avoided if one only tries to compute a prediction

$$g_{2k}^{(n)}=g_0^{(n,2k)}$$
(43)
for the *first* term $`\gamma _{n+2k+1}`$ not used for the construction of either $`\epsilon _{2k}^{(n)}`$ or $`\phi _{2k}^{(n)}(z)`$. For that purpose, we have only to set $`z=0`$ in (38). This yields the following recursive scheme (Eq. (3.17) of ):

$`g_0^{(n)}`$ $`=`$ $`0,n_0,`$ (44a)
$`g_2^{(n)}`$ $`=`$ $`{\displaystyle \frac{\left[\gamma _{n+2}\right]^2}{\gamma _{n+1}}},n_0,`$ (44b)
$`g_{2k+2}^{(n)}`$ $`=`$ $`g_{2k}^{(n+2)}+{\displaystyle \frac{\left[\gamma _{n+2k+2}g_{2k}^{(n+1)}\right]^2}{\gamma _{n+2k+1}g_{2k}^{(n)}}}{\displaystyle \frac{\left[\gamma _{n+2k+2}g_{2k}^{(n+1)}\right]^2}{\gamma _{n+2k+1}g_{2k2}^{(n+2)}}},`$ (44c)
$`k,n_0.`$
The main advantage of this recursive scheme over the recursive scheme (38), from which it was derived, is that it only involves numbers but no symbolic expressions.

In Table 5 we present selected results for the Padé predictions of the *first* coefficients of the perturbation expansion (28) for $`\mathrm{\Delta }E_0(\lambda ^2)`$ and (29) for $`\mathrm{\Delta }_0(\beta )`$, which were not used in the Padé approximants $`\epsilon _{2[[n/2]]}^{(n2[[n/2]])}=P_{[[n/2]]}^{n[[n/2]]}`$ for $`2n191`$. The first predictions $`b_0^{(n2[[n/2]],2[[n/2]])}`$ to $`b_{n+2}`$ and $`_0^{(n2[[n/2]],2[[n/2]])}`$ to $`_{n+2}`$ were computed with the help of the recursive scheme (44). For the sake of readability, we present in Table 5 only the relative errors

$$R_n=\frac{b_0^{(n2[[n/2]],2[[n/2]])}b_{n+2}}{b_{n+2}}$$
(45)
and

$$_n=\frac{_0^{(n2[[n/2]],2[[n/2]])}_{n+2}}{_{n+2}}.$$
(46)
If the input data of the recursive scheme (44) are exact rational numbers as the coefficients $`b_n`$ or $`_n`$, then the predictions can be computed free of rounding errors. However, it turned out that that the predictions computed in this way were *huge* rational numbers which slowed down computation considerably. Therefore, we used the floating-point arithmetics of Maple with an accuracy of 600 decimal digits for the computation of results presented in Table 5.

The results in Table 5 show that the first coefficients not used for the construction of the Padé approximants $`\epsilon _{2[[n/2]]}^{(n2[[n/2]])}=P_{[[n/2]]}^{n[[n/2]]}`$ can be predicted with remarkable accuracy by the recursive scheme (44) if $`n`$ is sufficiently large. Moreover, the agreement of the relative errors $`R_n`$ and $`_n`$ is stunning. However, for our purposes most important is the observation that all relative errors in Table 5 are negative which is in agreement with the Stieltjes inequality (36).

Thus, the prediction results of this section show that there is no significant difference between the prediction of coefficients $`b_n`$ of the perturbation series (28), whose Stieltjes nature we want to establish, and the prediction of the coefficients $`_n`$ of the analogous quartic anharmonic oscillator perturbation series (29), whose Stieltjes nature was established rigorously by Simon (Theorem IV.2.1 of ). Accordingly, these results provide further numerical evidence that the perturbation series (28) for the ground-state energy shift $`\mathrm{\Delta }E_0(\lambda ^2)`$ is indeed a Stieltjes series.

## 5 Summary and Conclusions

As shown by countless monographs or articles, the mathematical theory of the conventional Hermitian Hamiltonians of quantum mechanics is well established and has reached a high degree of sophistication. Moreover, the divergence and the summation of the perturbation expansions resulting from these Hamiltonians is also comparatively well understood, in particular in the case of the Padé summation of Stieltjes series. In contrast, the rigorous mathematical theory of non-Hermitian, $`𝒫𝒯`$-symmetric Hamiltonians is virtually nonexistent. It is not known whether divergent perturbation expansions resulting from such Hamiltonians can be summed to yield a uniquely determined result. Consequently, the best we can do for the moment is to perform numerical studies from which we can try to draw general qualitative conclusions.

The main intention of this paper is to provide numerical evidence that the perturbation series (28) for the ground-state energy shift $`\mathrm{\Delta }E_0(\lambda ^2)`$ is a Stieltjes series, because this would guarantee that certain subsequences of the Padé table constructed from the partial sums of this perturbation series converge to a uniquely determined summation result, as discussed in Section 2.

If the Padé approximants are computed with the help of Wynn’s recursive algorithm (13) according to (15) and (17) – as it is done in this paper – and if the series to be transformed is the moment expansion of a Stieltjes function, then the Padé approximants must satisfy inequality (27). As shown in Table 2, the Padé summation results for the perturbation series (28) for the ground-state energy shift $`\mathrm{\Delta }E_0(\lambda ^2)`$ as well as for the analogous perturbation series (29) for the ground-state energy shift $`\mathrm{\Delta }_0(\beta )`$ of the quartic anharmonic oscillator, which is known to be a Stieltjes series, satisfy this inequality. Moreover, the divergence of the two perturbation expansions (28) and (29) as well as as the convergence of their Padé summation results is virtually indistinguishable if the two coupling constants $`\lambda `$ and $`\beta `$ satisfy (7).

If a Padé approximant $`P_m^l(z)`$ constructed from the partial sums (22) of the moment expansion for a Stieltjes function $`F(z)`$ is expanded in power series around $`z=0`$ according to (35), then this series must be strictly alternating for $`z>0`$ just like the Stieltjes series for $`F(z)`$. Moreover, the coefficients $`\mu _n^{[l/m]}`$ of this expansion are bounded in magnitude by the Stieltjes moments $`\mu _n`$ according to the inequality (36).

Thus, via the inequality (36) it can be checked whether a moment expansion of the type of (20) is a Stieltjes series. However, the accuracy-through-order relationship (12) implies that this inequality is by default satisfied for all indices $`nl+m`$, and only for $`nl+m+1`$ do we obtain useful information. In particular for large values of $`l`$ and $`m`$, the symbolic construction of Padé approximants and their subsequent expansion may become quite demanding both with respect to time and memory.

These computational problems can to some extent be overcome by expressing the Padé approximants $`P_k^{k+n}(z)`$, that can be computed with the help of Wynn’s epsilon algorithm, according to (37) by the partial sum from which it was constructed plus the transformation term $`z^{n+2k+1}\phi _{2k}^{(n)}(z)`$. The quantities $`\phi _{2k}^{(n)}(z)`$ can be computed recursively with the help of (38), and their computation is less demanding than the recursive computation of the corresponding Padé approximants $`P_k^{k+n}(z)`$. Moreover, the Taylor expansion of $`\phi _{2k}^{(n)}(z)`$ according to (40) yields the desired coefficients $`\mu _{n+2k+\nu +1}^{[k+n/k]}`$ with $`\nu 0`$.

In Tables 3 and 4, we proceed as described above and confirm the validity of the Stieltjes inequality (36) in the case of the first 15 coefficients of the quantities $`\phi _{66}^{(0)}(\lambda ^2)`$ and $`\phi _{66}^{(0)}(\beta )`$ corresponding to the Padé approximants $`P_{33}^{33}(\lambda ^2)`$ and $`P_{33}^{33}(\beta )`$, respectively. In this context, it is remarkable that the quality of the predictions is virtually indistinguishable although the coefficients $`b_n`$ of the $`𝒫𝒯`$-symmetric perturbation series grow much more rapidly than the coefficients $`_n`$ for the quartic anharmonic oscillator.

For larger values of $`n`$ and $`k`$ the symbolic computation of the quantities $`\phi _{2k}^{(n)}(z)`$ via the recursive scheme (38) becomes quite demanding. In such a case, it is much simpler to compute only the prediction for the first series coefficient not used for the computation of $`\phi _{2k}^{(n)}(z)`$. This can be done with the help of the recursive scheme (44). In Table 5 we show that all first predictions, which we can compute from the coefficients $`b_\nu `$ and $`_\nu `$ with $`1\nu 192`$, satisfy the Stieltjes inequality (36). Moreover, the quality of the prediction results is again virtually identical.

We are of course aware that numerical results cannot replace rigorous mathematical proofs. Nevertheless, we believe that our numerical experiments are both interesting and useful, and that they provide considerable evidence that the perturbation series (28) for the ground-state energy shift $`\mathrm{\Delta }E_0(\lambda ^2)`$ of the $`𝒫𝒯`$-symmetric Hamiltonian (1) is indeed a Stieltjes series, which would imply that the Padé summation of this divergent perturbation series converges.

## Acknowledgment

CMB is grateful to the U.S. Department of Energy for financial support and EJW thanks the Fonds der Chemischen Industrie for financial support.