From divergent perturbation theory to an exponentially convergent self-consistent expansion

For many nonlinear physical systems, approximate solutions are pursued by conventional perturbation theory in powers of the nonlinear terms. Unfortunately, this often produces divergent asymptotic series, collectively dismissed by Abel as “an invention of the devil.” Although a lot of progress has been made on understanding the mathematics and physics behind this, new approaches are still called for. A related method, the self-consistent expansion (SCE), has been introduced by Schwartz and Edwards. Its basic idea is a rescaling of the zeroth-order system around which the solution is expanded, to achieve optimal results. While low-order SCEs have been remarkably successful in describing the dynamics of nonequilibrium many-body systems (e.g., the Kardar-Parisi-Zhang equation), its convergence properties have not been elucidated before. To address this issue we apply this technique to the canonical partition function of the classical harmonic oscillator with a quartic $g{x}^4$ anharmonicity, for which perturbation theory’s divergence is well-known. We obtain the $N$th order SCE for the partition function, which is rigorously found to converge exponentially fast in $N$, and uniformly in $g\ge 0$. We use our results to elucidate the relation between the SCE and the class of approaches based on the so-called “order-dependent mapping.” Moreover, we put the SCE to test against other methods that improve upon perturbation theory (Borel resummation, hyperasymptotics, Padé approximants, and the Lanczos $\tau$-method), and find that it compares favorably with all of them for small $g$ and dominates over them for large $g$. The SCE is shown to converge to the correct partition function for the double-well potential case, even when expanding around the local maximum. Our treatment is generalized to the case of many oscillators, as well as to any nonlinearity of the form $g|x{|}^q$ with $q\ge 0$ and complex $g$. These results allow us to treat the Airy function, and to see the fingerprints of Stokes lines in the SCE.

Figure 1

Performance of the SCE versus the parameter $\alpha$ and order $N$. (a) Relative error of the SCE as a function of $\alpha =M(N)/N$. $g$ is fixed at 1. Shown are measured errors for $N=20$ and 21. Also plotted is the analytical bound from Eq. (34), and the last, $n=N$ term in Eq. (19) (after summation over $l$, for $N=20$). The analytical bound is quite loose, while the $n=N$ term is a much better estimator of the error of the approximation, especially for $\alpha <{\alpha }^{*}$—this property is explored in Appendix pp1. All the graphs exhibit an optimum choice in the vicinity of $\alpha =1.3$, close to our estimated ${\alpha }^{*}\approx 1.317$ [Eq. (38)], marked by the vertical dashed line; note that this value was calculated in the limit of $N\gg 1$. The different behavior of $N=20$ and $N=21$ at optimum is explained in the text. (b) Convergence of the SCE at high order. Plotted are the relative errors of the SCE expansion of the partition function ${\mathcal{Z}}^{(N)}$ versus the order $N$ (markers). These are plotted for $M=\alpha N$ with $\alpha =1$, 2, and $\frac{4}{3}$, in order of merit (blue, purple, and yellow, respectively). The coupling $g$ is fixed at 1. The vertical axis reflects the number of significant digits successfully reproduced by the expansion. The data was fitted with a trend of the form ${10}^{-AN-B\sqrt{N}+C}$, corresponding to our analytical bound (continuous lines). The fits achieve a ${\chi }^2$ value of $\sim {10}^{-2}$ (${\chi }^2$ is defined in the text). Though even orders admit better accuracy, we evaluate at odd orders to mitigate the sensitivity of the fitted parameters to the exact choice of $\alpha$ near its optimum (cf. the left panel). The fitted values for the coefficient $A$ are 0.072, 0.200, and 0.283, for $\alpha =1$, 2, and $\frac{4}{3}$, respectively, representing the minimal amount of decimal places obtained per increment of $N$. The lower bounds for $A(\alpha )$ as predicted by Eq. (41) are 0.018, 0.176 and 0.243, in agreement with the numerical results. The minimal error shown is roughly $1.5\times {10}^{-92}$. Note that in order to obtain this many accurate digits, the intermediate calculations had to be performed with over 300 digit precision.

Figure 2

A comparison of the SCE with other asymptotic and numeric schemes, in order of merit (see Appendix pp2): Superasymptotics (i.e., truncation of PT at the least term) [1]; Borel resummation of the superasymptotic remainder [27 Sec. 7.2]; Padé approximants of PT [8]; 3rd level hyperasymptotics [5, 54]; SCE (with $\alpha =4/3$); and Lanczos’s Chebyshev $\tau$ approximations of $\mathcal{Z}(g)$ [55]. Since the Padé, SCE, and Lanczos methods converge and thus can be evaluated to arbitrary order, we choose to compare them with the asymptotic methods at the same order $N=N_0=\frac{1}{16g}$ as that used in the superasymptotic expansion. (a) The relative error of the SCE. The horizontal axis is the singulant $F=\frac{1}{16g}$ (cf. Subsection B.2), and at each point the SCE is evaluated to order $N_0=F$, which is the same truncation as the superasymptotic scheme. The line segments are a guide to the eye. (b) The errors produced by the other methods, relative to that produced by the SCE (i.e., the data was normalized by that shown in the top panel). The line segments are a guide to the eye. When evaluated at the same order, the SCE outperforms superasymptotics, Borel tail resummation, and the Padé approximants. Note the extremely similar errors produced by the SCE and hyperasymptotics, which is formally of order $2N_0$; this is discussed in the text. Lastly, the Chebyshev $\tau$ method appears more accurate than the SCE, but only because of the specific truncation at $1/16g$. In the large-$N$ limit, the SCE is more rapidly convergent, as demonstrated in Fig. 3. The dominance of the SCE in the strongly coupled regime is not manifest in this comparison, since for $g\gg 1$ the superasymptotic truncation occurs at $N_0=0$; however, for fixed and finite $N$, the SCE is vastly superior in this scenario, as depicted in Fig. 3.

Figure 3

Comparison of the SCE against numerical methods. (a) The relative error produced by SCE, Padé approximants [8] and the Chebyshev $\tau$ method [55] versus $g$ (see Appendix pp2). All methods are evaluated at $N=13$; the SCE is evaluated with $\alpha =\frac{4}{3}$. The SCE outperforms the Padé approximants across a wide range of values for $g$, and outperforms the $\tau$ method for $g\gtrsim 0.1$. This holds for other values of $N$ as well. As $g$ grows larger, the relative errors of the Padé and $\tau$ approximations diverge, though the approximations themselves are finite. The dashed red line shows the ultimate error of the SCE for $g\to \infty$, about $7.83\times {10}^{-6}$. This value was calculated by comparing the coefficients of $g^{-\frac{1}{4}}$ in the large-$g$ expansions of Eq. (4) [which is $\frac{1}{2}\mathrm{\Gamma }(\frac{1}{4})$] and Eq. (19) [obtained by taking $1-1/G\to 1$ and $\sqrt{\frac{2}{G}}\to {(gK)}^{-\frac{1}{4}}$]. The sharp dip in the error of the $\tau$ method near $g\approx 150$ is due to the error switching signs across that value, so at the crossover the relative error vanishes. This does not occur with the SCE since we evaluate at an odd order $N$ [cf. Fig. 1]. (b) The relative error produced by the same methods versus the approximation order $N$, up to order 200. All methods are evaluated at $g=0.01$; the SCE is evaluated with $\alpha =\frac{4}{3}$. The SCE is consistently more accurate than the Padé approximants, and it supersedes the $\tau$ method at roughly $N=180$, which depends on the value of $g$.

Figure 4

The convergence of the SCE for the case of a double-well potential. We examine the barely binding case of $g=0.01$, so that the differences between this case and the previous sections are accentuated. (a) The double-well counterpart of Fig. 1: Shown is the numerical accuracy of the SCE, evaluated for $\alpha =1$, 2, and $\frac{4}{3}$ (markers). Note that all three trends converge, though for the $\alpha =1$ case, the exponential convergence has barely begun to rein in the divergent stretched exponent at $N\approx 200$, showing that a proper choice of $\alpha$ is critical to the effectiveness of the approximation. All points are evaluated at odd order [cf. Fig. 1]. Against the data we fit a trend line of the form $10^{C-AN+B\sqrt{N}}$ (solid lines), and the fitted values of $A$ are 0.075, 0.197, and 0.282 for $\alpha =1$, 2, and $\frac{4}{3}$, respectively, essentially the same as those fitted in Fig. 1. (b) The Boltzmannian (normalized) distributions $e^{-\beta V(x)}/\mathcal{Z}$ of the double-well potential $V(x)=-\frac{1}{2}{x}^2+g{x}^4$, and the zeroth-order approximation $V(x)=+\frac{1}{2}G{x}^2$ around which the SCE is expanded. This plot corresponds to the point marked by a red star in the left panel, where $N=301$ and $\alpha =\frac{4}{3}$. Note how the SCE’s reference system bears little resemblance to the approximated one.

Figure 5

Trajectories of the three solutions of Eq. (90) for $y_{+}=G_{+}^{\frac{1}{4}}$ and varying $M$. (a) The trajectories for $z=1$. The arrowheads point in the direction the solution moves as $M$ is increased, up to 1000; the base of the arrow is the value of each root at $M=1$. Note that only one particular root is in the vicinity of $z^{\frac{1}{4}}$ (marked by a red star), tending to it asymptotically if $|z|$ is increased with $M$ held fixed. As $M$ increases the solutions tend towards the three cubic roots of $C_3(M)/3iM$, which lie on the positive imaginary axis and along the $\pm e^{\mp \frac{i\pi }{6}}$ directions, marked by the dashed lines. Note that the solution of interest (orange) tends to the principal branch of the cubic root, $\sqrt[3]{C_3(M)/3iM}$; this behavior experiences a crossover when transitioning across the Stokes line at $\arg z=\frac{2\pi }{3}$, after which the root which stems from $z^{\frac{1}{4}}$ tends to the imaginary axis, cf. the right panel. (b) The same trajectories for $z=16e^{\frac{2\pi i}{3}}$. Note that initially all roots lie on the line parallel to $e^{\frac{i\pi }{6}}$. The “primary” root coming from $z^{\frac{1}{4}}$ “collides” with the root coming from 0, after which they become a complex-conjugate pair of solutions for $r={(\frac{G}{z})}^{\frac{1}{4}}$, and continue to tend to the same asymptotes as before. For $\arg z$ any larger, the asymptotes of these two roots will be swapped.

Figure 6

Behavior of the SCE in the complex plane, for fixed $|z|$ and varying $\arg z$. (a) Trajectories of the factor $1-\sqrt{z}/\sqrt{G_{+}}$, as $M$ is increased up to $10^4$, for the listed values of $\arg z$ and $|z|=8$. The unit circle is plotted for reference. As the argument of $z$ is increased, the trajectories sweep outwards. For $|\arg z|<\frac{\pi }{3}$, the trajectories are contained inside the unit circle. For $\frac{\pi }{3}<|\arg z|<\frac{2\pi }{3}$, they exit it and eventually converge to $+1$ from outside. However, the closer $\arg z$ is to $\frac{2\pi }{3}$, the “sharper” the angle of protrusion through the circle, and the further outside the trajectory reaches before attaining its maximal magnitude. At $\arg z=\frac{2\pi }{3}$, the trajectory heads radially outwards on the negative real axis before turning sharply (red star) at $(-2)$, a value independent of $|z|$ which corresponds to the scattering point visible in Fig. 5. This sharp turn occurs at $M\approx 28.09$, as predicted by Eq. (95) for $|z|=8$. This point corresponds to the cusp in the right panel (red star). Lastly, for $|\arg z|>\frac{2\pi }{3}$, the origins of the trajectories jump abruptly to far away from the unit circle. Another arrowhead was placed on the points corresponding to $M=30$ for reference; for any larger $M$, the trajectories are continuous as a function of $\arg z$. (b) The relative error of the SCE versus the order $N$ for the same values of $z$ as in the left panel, with $M(N)=N$. For all arguments under $\frac{2\pi }{3}$, the error is roughly identical until order $\sim 10$, at which point the convergence rate is slowed down with increasing phase. As the argument is increased towards $\frac{2\pi }{3}$, a local bump in the error forms. For $\arg z=\frac{2\pi }{3}$, a cusp occurs at $N=29$ (red star). This is the first integer value of $N$ for which $M(N)$ is greater than the value at which the collision in Fig. 5 occurs, corresponding to the red star in the left panel. After this bump, the expansion returns to converge exponentially for all $\arg z$. The errors for $\arg z>\frac{2\pi }{3}$ exhibit the diverging nature of the Airy SCE in this Stokes wedge at low orders, but eventually the SCE starts to converge for sufficiently large $N$. Note the eye formed for $N<30$; for $N>30$, the error transitions smoothly across both Stokes lines, showing that the SCE smooths out the Stokes lines at large order. Lastly, the lower and upper dashed lines represent the uniform exponential convergence component in the $|\arg z|<\frac{\pi }{3}$ and $\frac{\pi }{3}<|\arg z|<\frac{2\pi }{3}$ Stokes sectors, respectively. The lower line was obtained by fitting an exponential profile to the SCE of $\mathrm{Ai}(0)$, which does not exhibit a stretched-exponential component (cf. Fig. 7). This line bounds the error for all $|\arg z|<\frac{\pi }{3}$, while for $\frac{\pi }{3}<|\arg z|<\frac{2\pi }{3}$, the error will eventually exceed it once the corresponding trajectory leaves the unit circle; for example, the yellow diamonds with $\arg z=\frac{4\pi }{9}$ cross this line at order $N\approx 200$. The upper dashed line was obtained by multiplying the lower line by $2^N$, as a bound on the factor ${(1-\sqrt{z}/\sqrt{G_{+}})}^N$ in the $\frac{\pi }{3}<|\arg z|<\frac{2\pi }{3}$ wedge. Note that to a good approximation, the cusp lies on this line (cf. the left panel).

Figure 7

Convergence of the SCE to Ai(0). Shown is the relative error of the expansion (dots) versus the order $N$, along with an exponential fit, $R^{(N)}=10^{B-AN}$ (solid line). We evaluate the series with $M=N$. The fitted value for $A$ is 0.492, showing that the SCE recovers another significant digit of Ai(0) every other order, greatly exceeding the convergence rate estimated in Table 1 for $q=3$.

Figure 8

Comparison of the Airy SCE with other asymptotic methods. (a) The relative error of the Airy SCE. The horizontal axis is the singulant $F=\frac{4}{3}{z}^{\frac{3}{2}}$, and at each point the SCE is evaluated to order $N_0=F$, which is the same truncation as the superasymptotic scheme. The solid lines are a guide to the eye. (b) The errors produced by superasymptotics, hyperasymptotics, Padé approximants and the $\tau$ method, relative to that produced by the SCE (i.e., the data was normalized by that shown in the top panel). The solid lines are a guide to the eye. Evaluated at the same order, the SCE consistently outperforms superasymptotics and the Padé approximants, and again produces an error very similar to hyperasymptotics (cf. Fig. 2 and subsequent discussion), although the hyperasymptotic series is formally of order $2N_0$. The Chebyshev $\tau$ method appears to give a few additional significant digits, but the SCE again converges faster in the large-$N$ limit—this is demonstrated in Fig. 9.

Figure 9

Comparison of the Airy SCE versus the Padé and Chebyshev $\tau$ approximations (a) The relative accuracy of the three methods as a function of $z$. All are evaluated at fixed order $N=7$. Also plotted is the estimated best achievable hyperasymptotic error, as predicted by Eq. (b25). By construction, the hyperasymptotic approximation order is $2N_0=2\times \frac{4}{3}{z}^{\frac{3}{2}}$, which gives $2N_0=7$ at $z\approx 1.9$. Much like the anharmonic case, SCE dominates the other methods for “stronger” couplings—this time, at smaller $z$. The noise seen near $z=10$ in the $\tau$ method is numerical noise due to evaluation at machine precision. The sharp dips in the SCE error are points where the error switches signs, i.e., points where ${\mathrm{Ai}}^{(N)}(z)$ and $\mathrm{Ai}(z)$ intersect. (b) The relative error of the SCE and the Lanczos Chebyshev $\tau$ method versus $N$, for $z=8$. As before [cf. Fig. 3], the SCE surpasses the $\tau$ method for sufficiently large $N$. Both methods were fitted with trend lines of the form $10^{C-AN-B{N}^{\frac{2}{3}}}$. The SCE yields the values $A=0.492$ and $B=1.446$, greatly exceeding our minimal bounds of uniform $A=0.305$ and $z$-dependent $B=1.28$. The $\tau$ method yields a negligible $A=0.007$, suggesting that it only contains the ${10}^{-B{N}^{\frac{2}{3}}}$ component, with $B=4.00$.

Figure 10

Phase contours of $k(f(z_n)-f_n)$. The orange stars are the saddles $z_n$, numbered from bottom to top. The red lines depict the steepest descent paths through $z_n$, where $k(f(z_n)-f_n)$ is real and positive. The colored contours map the absolute value of the argument of this quantity. Left panel: The steepest descent path through $z_2$, with ${\theta }_k=0$. Right panel: The paths through $z_{1,3}$ with ${\theta }_k=\pm \pi$. All contours point from left to right. Saddle 2 is adjacent to both 1 and 3, which are not mutually adjacent.

Figure 11

Plotted are the absolute values of the successive terms in the hyperasymptotic approximation of $\mathcal{Z}$, after multiplication by their symmetry factors as in (b23), for $g=\frac{1}{160}$. For this value, the hyperseries naturally halts after three scatterings. The zeroth and first levels are calculated exactly, the second level by numerical integration, and the last stage is approximated. The final relative error of the trans-series after it is summed is $5.7\times {10}^{-12}$.