Radial Anharmonic Oscillator: Perturbation Theory, New Semiclassical Expansion, Approximating Eigenfunctions. I. Generalities, Cubic Anharmonicity Case

TL;DR

This paper develops a comprehensive perturbation and semiclassical expansion framework for the radial anharmonic oscillator, introducing an approximant that accurately models eigenfunctions across all coupling regimes and dimensions.

Contribution

It introduces a novel interpolating approximant combining multiple expansions, achieving high-precision eigenfunction and energy calculations for the cubic anharmonic oscillator.

Findings

Approximant deviates less than 10^{-4} from exact eigenfunctions.

Achieves 7-8 significant digit accuracy in variational energies.

Provides a unified approach valid for all coupling constants and dimensions.

Abstract

For the general $D$-dimensional radial anharmonic oscillator with potential $V(r)= \frac{1}{g^2}\,\hat{V}(gr)$ the Perturbation Theory (PT) in powers of coupling constant $g$ (weak coupling regime) and in inverse, fractional powers of $g$ (strong coupling regime) is developed constructively in $r$-space and in $(gr)$ space, respectively. The Riccati-Bloch (RB) equation and Generalized Bloch (GB) equation are introduced as ones which govern dynamics in coordinate $r$-space and in $(gr)$-space, respectively, exploring the logarithmic derivative of wavefunction $y$. It is shown that PT in powers of $g$ developed in RB equation leads to Taylor expansion of $y$ at small $r$ while being developed in GB equation leads to a new form of semiclassical expansion at large $(g r)$: it coincides with loop expansion in path integral formalism. In complementary way PT for large $g$ developed in RB equation leads to an expansion of $y$ at large $r$ and developed in GB equation leads to an expansion at small $(g r)$. Interpolating all four expansions for $y$ leads to a compact function (called the {\it Approximant}), which should uniformly approximate the exact eigenfunction at $r \in [0, \infty)$ for any coupling constant $g \geq 0$ and dimension $D > 0$. Free parameters of the Approximant are fixed by taking it as a trial function in variational calculus. As a concrete application the low-lying states of the cubic anharmonic oscillator $V=r^2+gr^3$ are considered. It is shown that the relative deviation of the Approximant from the exact ground state eigenfunction is $\lesssim 10^{-4}$ for $r \in [0, \infty)$ for coupling constant $g \geq 0$ and dimension $D=1,2,\ldots$. In turn, the variational energies of the low-lying states are obtained with unprecedented accuracy 7-8 s.d. for $g \geq 0$ and $D=1,2,\ldots$.

Tables

Table 1. Table 1 : Ground state energy for the cubic potential ( r 2 + g ​ r 3 ) superscript 𝑟 2 𝑔 superscript 𝑟 3 (r^{2}\,+\,g\,r^{3}) for D = 1 , 2 , 3 , 6 𝐷 1 2 3 6 D=1,2,3,6 and g = 0.1 , 1 , 10 𝑔 0.1 1 10 g=0.1,1,10 , labeled by quantum numbers (0,0) for D > 1 𝐷 1 D>1 . Variational energy E 0 ( 1 ) superscript subscript 𝐸 0 1 E_{0}^{(1)} , the first correction E 2 subscript 𝐸 2 E_{2} (rounded to three s.d.) found with use of Ψ 0 , 0 ( t ) superscript subscript Ψ 0 0 𝑡 \Psi_{0,0}^{(t)} , see text, the corrected energy E 0 ( 2 ) = E 0 ( 1 ) + E 2 superscript subscript 𝐸 0 2 superscript subscript 𝐸 0 1 subscript 𝐸 2 E_{0}^{(2)}=E_{0}^{(1)}+E_{2} shown. E 0 ( 2 ) superscript subscript 𝐸 0 2 E_{0}^{(2)} coincides with Lagrange Mesh results (see text) in nine displayed decimal digits, hence all printed digits are exact.

$g$	$D=1$			$D=2$
	$E_0^{\left(1\right)}$	$-E_2$	$E_0^{\left(2\right)}$	$E_0^{\left(1\right)}$	$-E_2$	$E_0^{\left(2\right)}$
0.1	1.053120300	$5.39\times {10}^{-7}$	1.053119761	2.124027648	$4.40\times {10}^{-7}$	2.124027208
1.0	1.387428891	$4.00\times {10}^{-8}$	1.387428851	2.877490906	$3.76\times {10}^{-8}$	2.877490868
10.0	2.729533139	$6.56\times {10}^{-7}$	2.729532483	5.794213459	$5.58\times {10}^{-7}$	5.794212901
$g$	$D=3$			$D=6$
	$E_0^{\left(1\right)}$	$-E_2$	$E_0^{\left(2\right)}$	$E_0^{\left(1\right)}$	$-E_2$	$E_0^{\left(2\right)}$
0.1	3.208922743	$4.00\times {10}^{-7}$	3.208922343	6.528432540	$2.02\times {10}^{-7}$	6.528432338
1.0	4.442965260	$3.15\times {10}^{-8}$	4.442965229	9.465319951	$1.85\times {10}^{-8}$	9.465319933
10.0	9.094985589	$4.23\times {10}^{-7}$	9.094985166	19.981458504	$1.96\times {10}^{-7}$	19.981458308

Table 2. Table 2 : The first excited state energy for the cubic potential ( r 2 + g ​ r 3 ) superscript 𝑟 2 𝑔 superscript 𝑟 3 (r^{2}\,+\,g\,r^{3}) for different D 𝐷 D and g 𝑔 g . They labeled by quantum numbers (0,1) for D > 1 𝐷 1 D>1 , as for D = 1 𝐷 1 D=1 it corresponds to 1st excited state. Variational energy E 0 ( 1 ) superscript subscript 𝐸 0 1 E_{0}^{(1)} , the first correction E 2 subscript 𝐸 2 E_{2} found with use of Ψ 1 , 0 ( t ) superscript subscript Ψ 1 0 𝑡 \Psi_{1,0}^{(t)} for D = 1 𝐷 1 D=1 and Ψ 0 , 1 ( t ) superscript subscript Ψ 0 1 𝑡 \Psi_{0,1}^{(t)} for D > 1 𝐷 1 D>1 , see text, the corrected energy E 0 ( 2 ) = E 0 ( 1 ) + E 2 superscript subscript 𝐸 0 2 superscript subscript 𝐸 0 1 subscript 𝐸 2 E_{0}^{(2)}=E_{0}^{(1)}+E_{2} shown, the correction E 2 subscript 𝐸 2 E_{2} rounded to three s.d. All printed digits for E 0 ( 2 ) superscript subscript 𝐸 0 2 E_{0}^{(2)} are exact.

$g$	$D=1$			$D=2$
	$E_0^{\left(1\right)}$	$-E_2$	$E_0^{\left(2\right)}$	$E_0^{\left(1\right)}$	$-E_2$	$E_0^{\left(2\right)}$
0.1	3.208922765	$4.21\times {10}^{-7}$	3.208922343	4.305557665	$3.55\times {10}^{-7}$	4.305557309
1.0	4.442965265	$3.59\times {10}^{-8}$	4.442965229	6.068723537	$2.92\times {10}^{-8}$	6.068723507
10.0	9.094985630	$4.64\times {10}^{-7}$	9.094985166	12.579594377	$3.48\times {10}^{-7}$	12.579594029
$g$	$D=3$			$D=6$
	$E_0^{\left(1\right)}$	$-E_2$	$E_0^{\left(2\right)}$	$E_0^{\left(1\right)}$	$-E_2$	$E_0^{\left(2\right)}$
0.1	5.412425220	$2.86\times {10}^{-7}$	5.412424933	8.784695351	$1.21\times {10}^{-7}$	8.784695230
1.0	7.745092165	$2.41\times {10}^{-8}$	7.745092141	13.018486318	$1.49\times {10}^{-8}$	13.018486303
10.0	16.215748127	$2.66\times {10}^{-7}$	16.215747861	27.841430199	$1.37\times {10}^{-7}$	27.841430061

Table 3. Table 3 : The second excited state energy for the cubic potential r 2 + g ​ r 3 superscript 𝑟 2 𝑔 superscript 𝑟 3 r^{2}\,+\,g\,r^{3} for different D 𝐷 D and g 𝑔 g , they labeled by quantum numbers (0,2) for D > 1 𝐷 1 D>1 , as for D = 1 𝐷 1 D=1 it corresponds to 2nd excited state n r = 2 subscript 𝑛 𝑟 2 n_{r}=2 . Variational energy E 0 ( 1 ) superscript subscript 𝐸 0 1 E_{0}^{(1)} found with use of Ψ 2 , 0 ( t ) superscript subscript Ψ 2 0 𝑡 \Psi_{2,0}^{(t)} for D = 1 𝐷 1 D=1 and Ψ 0 , 2 ( t ) superscript subscript Ψ 0 2 𝑡 \Psi_{0,2}^{(t)} for D > 1 𝐷 1 D>1 , see text, the first correction E 2 subscript 𝐸 2 E_{2} and the corrected energy E 0 ( 2 ) = E 0 ( 1 ) + E 2 superscript subscript 𝐸 0 2 superscript subscript 𝐸 0 1 subscript 𝐸 2 E_{0}^{(2)}=E_{0}^{(1)}+E_{2} shown. Correction E 2 subscript 𝐸 2 E_{2} rounded to 3 s.d. All 9 printed decimal digits in E 0 ( 2 ) superscript subscript 𝐸 0 2 E_{0}^{(2)} are exact.

$g$	$D=1$			$D=2$
	$E_0^{\left(1\right)}$			$E_0^{\left(1\right)}$	$-E_2$	$E_0^{\left(2\right)}$
0.1	5.436849553			6.528432582	$2.43\times {10}^{-7}$	6.52843233834
1.0	7.879141644			9.465319955	$2.21\times {10}^{-8}$	9.46531993256
10.0	16.641305904			19.981458531	$2.23\times {10}^{-7}$	19.98145830814
$g$	$D=3$			$D=6$
	$E_0^{\left(1\right)}$	$-E_2$	$E_0^{\left(2\right)}$	$E_0^{\left(1\right)}$	$-E_2$	$E_0^{\left(2\right)}$
0.1	7.652743974	$1.87\times {10}^{-7}$	7.652743787	11.069434802	$6.73\times {10}^{-8}$	11.069434735
1.0	11.224406591	$1.87\times {10}^{-8}$	11.224406573	16.699837135	$1.22\times {10}^{-8}$	16.699837123
10.0	23.860743313	$1.78\times {10}^{-7}$	23.860743135	36.070426676	$1.01\times {10}^{-7}$	36.070426576

Table 4. Table 4 : Variational energy E 0 ( 1 ) superscript subscript 𝐸 0 1 E_{0}^{(1)} of the third excited state with quantum numbers (1,0) - the first radial excitation for the potential ( r 2 + g ​ r 3 ) superscript 𝑟 2 𝑔 superscript 𝑟 3 (r^{2}+gr^{3}) for D = 2 , 3 , 6 𝐷 2 3 6 D=2,3,6 and at g = 0.1 , 1 , 10 𝑔 0.1 1 10 g=0.1,1,10 and its node r 0 ( 0 ) superscript subscript 𝑟 0 0 r_{0}^{(0)} , found with use of Ψ 1 , 0 ( t ) superscript subscript Ψ 1 0 𝑡 \Psi_{1,0}^{(t)} , see text. Correction E 2 subscript 𝐸 2 E_{2} (not shown) contributes systematically to the 7th d.d.

$g$	$D=2$		$D=3$		$D=6$
	$E_0^{\left(1\right)}$	$r_0^{\left(0\right)}$	$E_0^{\left(1\right)}$	$r_0^{\left(0\right)}$	$E_0^{\left(1\right)}$	$r_0^{\left(0\right)}$
0.1	6.570942086	0.953377788	7.709696613	1.162457356	11.15814973	1.626236134
1.0	9.690374810	0.780305457	11.517370500	0.941956538	17.128462944	1.289494458
10.0	20.681623429	0.532055331	24.758598615	0.638726047	37.346045552	0.865045854

Table 5. Table 5 : Different partial sums ( III.18 ) for the leading term ε ~ 0 subscript ~ 𝜀 0 \tilde{\varepsilon}_{0} of the strong coupling expansion ( V.19 ) for the one-dimensional cubic anharmonic oscillator in PT with trial function ( V.21 ) as zero approximation. Exact ε ~ 0 subscript ~ 𝜀 0 \tilde{\varepsilon}_{0} obtained via Lagrange mesh method with 50 mesh points.

${\tilde{\epsilon }}_{0,0}$	${\sum }_{k=0}^1{\tilde{\epsilon }}_{0,k}$	${\sum }_{k=0}^2{\tilde{\epsilon }}_{0,k}$	${\sum }_{k=0}^3{\tilde{\epsilon }}_{0,k}$	${\sum }_{k=0}^4{\tilde{\epsilon }}_{0,k}$	${\sum }_{k=0}^5{\tilde{\epsilon }}_{0,k}$	${\sum }_{k=0}^6{\tilde{\epsilon }}_{0,k}$	Exact ${\tilde{\epsilon }}_0$
0	1.053006976	1.021174929	1.022989568	1.022956899	1.022946414	1.022947763	1.022947875

Table 6. Table 6 : Ground state energy ε ~ 0 subscript ~ 𝜀 0 \tilde{\varepsilon}_{0} (see ( V.19 )) for the potential W = r 3 𝑊 superscript 𝑟 3 W=r^{3} for D = 1 , 2 , 3 , 6 𝐷 1 2 3 6 D=1,2,3,6 found in PT with the Approximant Ψ 0 , 0 ( t ) superscript subscript Ψ 0 0 𝑡 \Psi_{0,0}^{(t)} as the entry: ε ~ 0 ( 1 ) superscript subscript ~ 𝜀 0 1 \tilde{\varepsilon}_{0}^{(1)} is the variational energy, ε ~ 2 subscript ~ 𝜀 2 \tilde{\varepsilon}_{2} is the second PT correction and ε ~ 0 ( 2 ) = ε ~ 0 ( 1 ) + ε ~ 2 superscript subscript ~ 𝜀 0 2 superscript subscript ~ 𝜀 0 1 subscript ~ 𝜀 2 \tilde{\varepsilon}_{0}^{(2)}=\tilde{\varepsilon}_{0}^{(1)}+\tilde{\varepsilon}_{2} is corrected variational energy. Eight decimal digits in ε ~ 0 ( 2 ) superscript subscript ~ 𝜀 0 2 \tilde{\varepsilon}_{0}^{(2)} confirmed in Lagrange mesh calculation.

$D=1$			$D=2$
${\tilde{\epsilon }}_0^{\left(1\right)}$	$-{\tilde{\epsilon }}_2$	${\tilde{\epsilon }}_0^{\left(2\right)}$	${\tilde{\epsilon }}_0^{\left(1\right)}$	$-{\tilde{\epsilon }}_2$	${\tilde{\epsilon }}_0^{\left(2\right)}$
1.022948250	$3.75\times {10}^{-7}$	1.022947875	2.187461809	$3.09\times {10}^{-7}$	2.187461499
$D=3$			$D=6$
${\tilde{\epsilon }}_0^{\left(1\right)}$	$-{\tilde{\epsilon }}_2$	${\tilde{\epsilon }}_0^{\left(2\right)}$	${\tilde{\epsilon }}_0^{\left(1\right)}$	$-{\tilde{\epsilon }}_2$	${\tilde{\epsilon }}_0^{\left(2\right)}$
3.450562918	$2.29\times {10}^{-7}$	3.450562689	7.647118254	$1.01\times {10}^{-7}$	7.647118153

Table 7. Table 7 : Subdominant (next-to-leading) term ε ~ 1 subscript ~ 𝜀 1 \tilde{\varepsilon}_{1} in the strong coupling expansion ( V.19 ) of the ground state energy for the cubic anharmonic radial potential for D = 1 , 2 , 3 , 6 𝐷 1 2 3 6 D=1,2,3,6 .

$D=1$			$D=2$
${\tilde{\epsilon }}_1^{\left(1\right)}$	${\tilde{\epsilon }}_{1,1}$	${\tilde{\epsilon }}_1^{\left(2\right)}$	${\tilde{\epsilon }}_1^{\left(1\right)}$	${\tilde{\epsilon }}_{1,1}$	${\tilde{\epsilon }}_1^{\left(2\right)}$
0.410598524	$6.78\times {10}^{-7}$	0.410599202	0.766573847	$5.24\times {10}^{-7}$	0.766574371
$D=3$			$D=6$
${\tilde{\epsilon }}_1^{\left(1\right)}$	${\tilde{\epsilon }}_{1,1}$	${\tilde{\epsilon }}_1^{\left(2\right)}$	${\tilde{\epsilon }}_1^{\left(1\right)}$	${\tilde{\epsilon }}_{1,1}$	${\tilde{\epsilon }}_1^{\left(2\right)}$
1.092125224	$2.67\times {10}^{-7}$	1.092125491	1.967599668	$1.42\times {10}^{-7}$	1.967599810

Table 8. Table 8 : Parameters of fit ( V.23 ) for different D 𝐷 D .

Parameter	$D=1$	$D=2$	$D=3$	$D=6$
$a$	3.281	3.922	4.823	5.994
$b$	1.023	1.094	1.150	1.275

Equations

$$\hat{H}=-\frac{{\hslash }^2}{2M}\sum _{k=1}^D{\partial }_{x_k}^2+V(r),V(r)=\frac{1}{g^2}\hat{V}(gr)\equiv \frac{1}{g^2}\sum _{k=2}^m{a}_k{g}^k{r}^k,$$

$$L=\frac{1}{2}\dot{\vec{\phi }}-\frac{1}{g^2}\sum _{k=2}^m{a}_k{g}^k(\vec{\phi }{}^2{)}^{k/2},$$

$$\Psi =e^{-\frac{1}{\hslash }\Phi },$$

$$\hslash y^{\prime }-y^2=2M(E-V),y={\Phi }^{\prime },$$

$$y\ =\ y_{0}\ +\ \hbar y_{1}\ +\ \hbar^{2}y_{2}\ +\ \ldots\$$

$$V=V_0+\lambda V_1+{\lambda }^2{V}_2+\dots ,$$

$$\left[-\frac{{\hslash }^2}{2M}{\nabla }_D^2+V(r)\right]\psi =E\psi ,$$

$${\nabla }_D^2={\partial }_r^2+\frac{D-1}{r}{\partial }_r+\frac{{\Delta }_{S^{D-1}}}{r^2},{\partial }_r\equiv \frac{\partial }{\partial r},$$

$${\psi }_{n_r,\ell ,\{m_{\ell }\}}(r,\Omega )={\Psi }_{n_r,\ell }(r){\xi }_{\ell ,\{m_{\ell }\}}(\Omega ),$$

$$\mathcal{N}(D,\ell )=\frac{(2\ell +D-2)(\ell +D-3)!}{\ell !(D-2)!},$$

$$-{\Delta }_{S^{D-1}}{\xi }_{\ell ,\{m_{\ell }\}}={\hat{\mathbf{L}}}^2{\xi }_{\ell ,\{m_{\ell }\}}=\ell (\ell +D-2){\xi }_{\ell ,\{m_{\ell }\}}.$$

$$\left[-\frac{{\hslash }^2}{2M}\left({\partial }_r^2+\frac{D-1}{r}{\partial }_r-\frac{\ell (\ell +D-2)}{r^2}\right)+V(r)\right]{\Psi }_{n_r,\ell }(r)=E_{n_r,\ell }{\Psi }_{n_r,\ell }(r),{\partial }_r\equiv \frac{d}{dr}.$$

$${\hat{h}}_r=-\frac{{\hslash }^2}{2M}\left({\partial }_r^2+\frac{D-1}{r}{\partial }_r\right)+V(r),$$

$${\hat{h}}_r\Psi (r)=E\Psi (r),{\int }_0^{\infty }{\Psi }^2{r}^{D-1}dr<\infty ,\Psi (0)=1,\Psi (\infty )=0.$$

$$\hbar\,\partial_{r}y(r)\ -\ y(r)\,\left(y(r)-\frac{\hbar(D-1)}{r}\right)\ =\ 2M\,\bigg{[}E\ -\ V(r)\bigg{]}\ ,$$

$$y(r)=-\hslash {\partial }_r\left(\log \Psi (r)\right)=-\hslash \frac{{\partial }_r\Psi (r)}{\Psi (r)}={\partial }_r\Phi (r).$$

$$\displaystyle V(r\,;a_{2},a_{3},...,a_{m})\$$

$$\hslash {\partial }_{r}y(r)-y(r)\left(y(r)-\frac{\hslash (D-1)}{r}\right)=2M\left[E-\frac{1}{g^2}\hat{V}(gr)\right],$$

$$v\ =\ \bigg{(}\frac{2M}{\hbar^{2}}\bigg{)}^{\frac{1}{4}}\,r\ ,$$

$$y=(2M{\hslash }^2{)}^{\frac{1}{4}}\mathcal{Y}\left(v\right),E=\frac{\hslash }{(2M{)}^{\frac{1}{2}}}\epsilon .$$

$${\partial }_v\mathcal{Y}(v)-\mathcal{Y}(v)\left(\mathcal{Y}(v)-\frac{D-1}{v}\right)=\epsilon \left(\lambda \right)-\frac{1}{{\lambda }^2}\hat{V}\left(\lambda v\right),{\partial }_v\equiv \frac{d}{dv},$$

$$\lambda ={\left(\frac{{\hslash }^2}{2M}\right)}^{\frac{1}{4}}g,$$

$$\partial_{v}\mathcal{Y}(v)\Bigr{|}_{v=0}\ =\ \frac{\varepsilon}{D}\quad,\quad\mathcal{Y}(+\infty)\ =\ +\infty\ .$$

$$\mathcal{Y}(v)=\frac{\epsilon }{D}v+\frac{\left({\epsilon }^2-a_2{D}^2\right)}{D^2(D+2)}{v}^3-\frac{a_3\lambda }{D+3}{v}^4-\frac{a_4{D}^3(D+2){\lambda }^2-2\epsilon ({\epsilon }^2-a_2{D}^2)}{D^3(D+2)(D+4)}{v}^5+\dots ,$$

$$\mathcal{Y}(v)\ =\ a_{m}^{1/2}\lambda^{(m-2)/2}v^{m/2}\ +\$$

$$\frac{a_{m-1}{\lambda }^{(m-4)/2}}{2a_m^{1/2}}{v}^{(m-2)/2}+\frac{(4a_m{a}_{m-2}-a_{m-1}^2){\lambda }^{(m-6)/2}}{8a_m^{3/2}}{v}^{(m-4)/2}+\dots .$$

$$\mathcal{Y}(v)=\frac{\epsilon }{D}v,\epsilon =a_2^{1/2}D.$$

$$\frac{1}{\hslash }\Phi =\frac{\epsilon }{2D}{v}^2+\frac{\left({\epsilon }^2-a_2{D}^2\right)}{4D^2(D+2)}{v}^4-\frac{a_3\lambda }{5(D+3)}{v}^5-\frac{a_4{D}^3(D+2){\lambda }^2-2\epsilon ({\epsilon }^2-a_2{D}^2)}{6D^3(D+2)(D+4)}{v}^6+\dots ,$$

$$\frac{1}{\hslash }\Phi =\frac{2a_m^{1/2}}{m+2}{\lambda }^{(m-2)/2}{v}^{(m+2)/2}+\frac{a_{m-1}{\lambda }^{(m-4)/2}}{m{a}_m^{1/2}}{v}^{m/2}+\frac{(4a_m{a}_{m-2}-a_{m-1}^2){\lambda }^{(m-6)/2}}{4(m-2)a_m^{3/2}}{v}^{(m-2)/2}+\dots ,$$

$$u=gr,$$

Full text

Introduction

The Hamiltonian of the $D$-dimensional radial anharmonic oscillator is given by

[TABLE]

where $D=1$, $2,$ $...$, $r=\left(\sum_{k=1}^{D}x_{k}^{2}\right)^{1/2}$ is the radius and $M$ is the mass of the particle, $\partial_{x_{k}}\equiv\frac{\partial}{\partial{x_{k}}}$. Choosing the units in such a way that $\hbar=1$ and $M=1$ the coupling constant $g\geq 0$ is of dimension $[cm]^{-1}$, $a_{k},k=2,\ldots m$ are real dimensionless parameters; it is assumed that $a_{2},a_{m}>0$, hence, $V(r)$ has minimum at $r=0$. At $m>3$ the potential can have additional minima. For the sake of simplicity, we assume that none of them is degenerate with one at $r=0$ and the potential is positive, $V(r)>0$ at $r>0$, thus, the minimum at $r=0$ is global. Evidently, the Hamiltonian (1) has the infinite discrete spectra. Needless to say, its study is of fundamental interest in different branches of physics, it plays an important role in atomic-molecular physics, plasma, nuclear and solid state physics. In particular, at $D=3$ with a proper choice of parameters $\{a_{k}\}$ and relaxing the requirement that at $r=0$ there is a global minimum, (1) can be used to describe the low-lying ro-vibrational states of a diatomic molecule. Also, the Hamiltonian of the anharmonic oscillator represents a scalar quantum field theory in $(0+1)$ dimensions with an associated Lagrangian of the form

[TABLE]

where $\vec{\varphi}=(\varphi_{1},...,\varphi_{D})$ is a $D$-dimensional vector field. Many properties that appear in (0+1) dimensions like asymptotic perturbation theory can also be present in a realistic field theory Bender and Wu (1969). It is well known that quartic field theory at $m=4$ is renormalizable in $(3+1)$ dimensions, see e.g. Lipatov (1976), while cubic ($m=3$) and sextic ($m=6$) ones are renormalizable in six and three dimensions, respectively.

There were published many papers devoted to exploration of the spectra of the Hamiltonian (1) for some concrete potentials at $m=4,6$. It is worth mentioning several methods, which have been used in the past: Variational method Turbiner (2005, 2010), Rayleigh-Ritz method Taşeli and Eid (1996), Perturbation Theory (PT) in both the weak and strong coupling regimes Bender and Wu (1969, 1973); Weniger (1996); Graffi et al. (1970); Silverstone et al. (1985); Eletsky and Popov (1980); Fernández and Guardiola (1997), Padé Approximants Loeffel et al. (1969), Hill Determinant method Witwit (1993); Biswas et al. (1971); Tater and Turbiner (1993), WKB method Lu (1967); Hioe et al. (1978), Self-Similar Approximations Yukalova and Yukalov (1993), the method of Characteristic Functions Ari and Demiralp (1985) among many others.

The present paper aims primarily to study the ground state of the Hamiltonian (1) using as tools PT in both the weak and strong coupling regimes, the asymptotic series expansions of the eigenfunction at small and large distances and a new version of the semiclassical approximation. Our ultimate goal is to construct (for any $r\in[0,\infty)$) a locally-accurate approximation of the ground state wave function in arbitrary (integer) dimension $D>0$. We call this approximation the Approximant and it is denoted by $\Psi_{0,0}^{(t)}$. Moreover, we will show how the Approximant $\Psi_{0,0}^{(t)}$ can be modified to construct Approximants for the wave functions of excited states. It will be demonstrated that any of the above Approximants being used as trial function in variational calculation leads to highly accurate variational energy. Furthermore, if the trial function is used as the zeroth order approximation in PT scheme, a fastly-convergent PT occurs for both energy and wave function. A similar program was implemented in Turbiner (2005, 2010) for the low-lying states of the one-dimensional quartic anharmonic oscillator leading to an extremely high accuracy for the energy: $\sim 10$ correct significant digits (s.d.) for any $g\geq 0$ and also a locally accurate wave function (not less than 6 decimal digits (d.d.) of relative accuracy for any $g\geq 0$ and $r\in[0,\infty)$). In fact, the approximation of the wave function given in Turbiner (2005, 2010) can be regarded as a solution of a non-solvable, one-dimensional quartic anharmonic oscillator. It will be shown in this paper that a similar situation occurs for the $D$-dimensional cubic radial anharmonic potential.

Let us take the ground state wave function $\Psi$ in the exponential representation

[TABLE]

and call $\Phi$ the phase of the wave function $\Psi$. It is well known in one-dimensional case that the exponential representation (3) being substituted into the Schrödinger equation leads to the nonlinear Riccati equation for the derivative of the phase Landau and Lifshitz (1977),

[TABLE]

where $E$ is the energy and $V$ is the potential. Similar equation appears for the $D$-dimensional radial case. The Riccati equation (sometimes also called the Bloch equation), by choosing the zero order approximation $y_{0}$ accordingly, can be solved iteratively keeping $E$ fixed, by considering the expansion in powers of Planck constant,

[TABLE]

looking for corrections $y_{n}$, $n=0,1,2,...$ . Such an iterative procedure leads to the well-known semiclassical WKB approximation, see e.g. Landau and Lifshitz (1977). It is well-known that this expansion has zero radius of convergence, since $\hbar$ stands in front of leading derivative in (4).

Another iterative approach is represented by the so-called Non-Linearization Procedure Turbiner (1984). Assuming a potential of the form

[TABLE]

we construct PT looking for both $y$ and $E$ in the form of Taylor series in powers of the (formal) parameter $\lambda$. Evidently, this approach can be straightforwardly generalized to radial potentials in arbitrary dimension $D$. It is well known that the Non-Linearization Procedure seems highly appropriate to study the bound states for the general potential in (1).

For the one-dimensional anharmonic oscillator, a novel iterative approach was recently discover based on transforming the Riccati equation (4) into GB Equation, see Escobar-Ruiz et al. (2016, 2017); Shuryak and Turbiner (2018). It turned out this approach is also appropriate for the construction of a type of semiclassical expansion for path integral in Euclidian time for the density matrix in coinciding points with respect to quantum and thermal fluctuations around classical path (flucton). Later it will be revealed a certain connection between the standard semiclassical WKB approximation, the Non-Linearization Procedure and the PT for the GB equation. This connection, which seems already highly non-trivial even for one-dimensional case Shuryak and Turbiner (2018), will be established for the $D$-dimensional radial anharmonic oscillator (1). Note that the interpretation of a new semiclassical expansion in terms of the loop expansion for path integral in Euclidian time is not yet clear to the present authors and it waits to be uncovered. We discuss and exploit the above connection in order to design a prescription for construction of the Approximant $\Psi_{0,0}^{(t)}$ - the function which provides uniform approximation of the ground state function in $r,g,D$.

As illustration we present detailed results for the cubic anharmonic oscillator, seemingly for the first time. Based on the Non-Linearization Procedure we study the structure of the perturbative corrections when PT in powers of $g$ is developed. Using the Approximant as a trial function, we perform numerical variational calculations to compute the energy (and estimate its accuracy) for some low-lying states in different $D$. The deviation of the Approximant from the exact ground state function is also estimated in framework of Non-Linearization Procedure. In addition, we calculate the first two dominant terms in the strong coupling expansion of the ground state energy using the Approximant.

The paper is organized in the following way. In the Section I starting from the Riccati equation for logarithmic derivative of the wave function the RB and the GR equations are introduced. Section II is dedicated to study the domain of weak coupling regime based on RB and GR equations. The strong coupling regime is explored in Section III, in Section IV the uniform approximation of the phase and ground state eigenfunction is built. The cubic anharmonic oscillator is described in details in Section V.

I Generalities

The Schrödinger equation for a spherical-symmetric potential $V(r)$ in a $D$-dimensional space has the form

[TABLE]

where $D=1,2,\ldots$ and $\nabla^{2}_{D}=\sum_{k=1}^{D}\partial^{2}_{x_{k}}$ is the $D$-dimensional Laplacian. Following the spherical symmetry of the potential, we introduce in (I.1) the $D$-dimensional hyperspherical coordinates $\{r,\Omega\}$, where $r$ is the hyperradius and $\Omega$ is solid angle parametrized by $(D-1)$ Euler angles. In these coordinates the Laplacian $\nabla^{2}_{D}$ takes the form

[TABLE]

where $\Delta_{S^{D-1}}$ is the Laplacian on the $S^{D-1}$ sphere. Indeed, $\Delta_{S^{D-1}}\,=\,-\hat{\mathbf{L}}^{2}$, where $\hat{\mathbf{L}}$ is the $D$-dimensional angular momentum operator. Two remarks in row: $(i)$ for any spherically symmetric potential the $D$-dimensional angular momentum $\hat{\textbf{L}}$ is conserved, $(ii)$ in hyperspherical coordinates we can separate the hyperradial coordinate $r$ from the angular coordinates $\Omega$. Hence, in $D>1$ the wave function $\psi$ can be labeled by radial quantum number $n_{r}$, angular quantum momentum $\ell$ and $(D-2)$ magnetic quantum numbers, it is represented as the product of two functions

[TABLE]

where $n_{r}=0,1,2,\ldots$ , $\ell=0,1,2,...$ . The set of $(D-2)$ magnetic quantum numbers $\{m_{\ell}\}$, where each of them takes values from $-\ell$ to $+\ell$, in total contains $\mathcal{N}(D,\ell)$ different configurations, where

[TABLE]

and $\mathcal{N}(D,0)\ =\ 1$. The angular part of the wave function (for $D>1$) $\xi_{\ell,\{m_{\ell}\}}(\Omega)$ corresponds to a $D$-dimensional spherical harmonic Müller (1966); Efthimiou and Frye (2014), which satisfies the eigenvalue equation

[TABLE]

Hence, for given $\ell$ this equation is satisfied by $\mathcal{N}(D,\ell)$ different orthogonal spherical harmonics, which correspond to the same eigenvalue $\ell\,(\ell+D-2)$. Its degeneracy is equal to $\mathcal{N}(D,\ell)$. Note that for $D=3$ the degeneracy is given by the familiar formula $\mathcal{N}(3,\ell)=2\ell+1$. Since the theory of spherical harmonics is well established, we should focus solely on determining the radial part of the wavefunction $\Psi_{n_{r},\ell}(r)$.

After substituting (I.2), (I.3) and (I.5) to (I.1) we eventually arrive at the radial Schrödinger equation, which determine $\Psi_{n_{r},\ell}(r)$,

[TABLE]

For any radial potential, the energy $E_{n_{r},\ell}$ of excited state is always degenerate with respect to the quantum numbers $\{m_{\ell}\}$. Specifically, for given $n_{r}$, $\ell$ and $D>1$ we have $\mathcal{N}(D,\ell)$ different wave functions with the same energy. For this reason since now on we omit the label $\{m_{\ell}\}$ in the energy and the wave function.

The ground state is non-degenerate (and nodeless), it has to zero angular momentum $\ell=0$ thus being $S$-state and zero radial quantum number $n_{r}=0$. The angular part corresponds to the zero harmonics (which is a constant) with zero eigenvalue in (I.5). From (I.6) it can be seen that the ground wave function (as well as all $S$ state eigenfunctions) depends only on radial coordinate $r$ that now on we denote by $\Psi(r)$ without labels. It corresponds to the lowest energy eigenfunction of the radial Schrödinger operator

[TABLE]

for which the corresponding eigenvalue equation with boundary conditions read

[TABLE]

Here $E$ denotes the ground state energy, assuming that $\hbar,m$ are dimensionless, it has dimension $[cm]^{-2}$ as well as the potential. From now on $D$ can be considered as a continuous parameter. Interestingly, some non-trivial analytical properties of the ground state energy as a function of $D$ appear at both non-physical dimensions $D\leq 0$ and in the physical ones, $D>0$, see e.g. Dolgov and Popov (1979) and Doren and Herschbach (1986).

Taking the ground state wave function $\Psi$ in the exponential representation (3) and substituting it into the eigenvalue problem (I.8) with radial operator (I.7) we arrive at the Riccati equation

[TABLE]

where function $y(r)$ is the logarithmic derivative of $\Psi(r)$ or, equivalently, the derivative of the phase $\Phi(r)$,

[TABLE]

We focus on case when the potential $V(r)$ has the form of finite-degree polynomial in $r$,

[TABLE]

see (1). Hence, the Riccati equation reads

[TABLE]

c.f. (4). This is the basic equation in the present work.

I.1 Riccati-Bloch Equation

The Riccati equation (I.12) can be transformed into one without explicit dependence on $\hbar$ and $M$ by introducing new variable

[TABLE]

and making the replacements

[TABLE]

Then, the equation (I.12) becomes

[TABLE]

where

[TABLE]

plays the role of the effective coupling constant replacing $g$. It is evident that $\lambda v=gr$. The boundary conditions those we impose in equation (I.15) are

[TABLE]

Evidently, if $\hbar=1$ and $M=1/2$, the equation (I.15) coincides with (I.12) and $\lambda=g$. We call the equation (I.15) the Riccati-Bloch (RB) Equation.

The RB equation (I.15) allows us to construct the asymptotic expansions of $\mathcal{Y}(v)$, and ultimately of $y(r)$, at small and large $v$. For fixed $\lambda\neq 0$ and small $v$ we have the Taylor series

[TABLE]

while for large $v$,

[TABLE]

[TABLE]

It is easy to check that for $\lambda=0$, hence, for the harmonic oscillator case $V=a_{2}\,r^{2}$, the equation (I.15) has the exact solution

[TABLE]

Three remarks in a row: $(i)$ the expansion (I.18) at $\lambda=0$ is terminated and consists of the first term alone, leading to $\mathcal{Y}(v)=(\varepsilon/D)\,v$, in agreement with (I.20), $(ii)$ in general, any coefficient in (I.18) depends on energy explicitly, $(iii)$ the fourth and subsequent (unwritten) terms in (I.19) can acquire explicit dependence on $\varepsilon$. It is worth mentioning that for fixed $M$, large $v$ behavior can be obtained in two situations: large $r$ and fixed $\hbar$, or at fixed $r$ and small $\hbar$, see (I.13). Hence, the expansion (I.19) describes the semiclassical limit $\hbar\rightarrow 0$.

Integrating (I.18) and (I.19) in $v$, both expansions are converted into expansions for the phase $\Phi$, namely,

[TABLE]

and

[TABLE]

for small and large $v$, respectively. Later these expansions will be used to make interpolation between large and small $v$ to construct a (uniform) approximation for the phase.

I.2 Generalized Bloch Equation

There exists an alternative way to transform the original Riccati equation (I.12) into equation without the explicit dependence on $\hbar$ and $M$. This is achieved by using a change of variable

[TABLE]

and introducing a new unknown function

[TABLE]

where $\mathcal{Z}(u)$ satisfies the non-linear differential equation

[TABLE]

c.f.(I.15), with boundary conditions

[TABLE]

Here $\varepsilon$ and $\lambda$ are the same as in (I.14) and (I.16). They play the role of energy and effective coupling constant in both $v$-space (I.13) and $u$-space (I.23) dynamics.

At $D=1$ the equation (I.25) was called in Escobar-Ruiz et al. (2016, 2017); Shuryak and Turbiner (2018) the (one-dimensional) GB Equation. Here, the equation (I.25) is a natural extension of one-dimensional GB equation to $D$-dimensional case, for a $D$-dimensional radial potential. For this reason we continue to call it the Generalized (radial) Bloch Equation. Naturally, (I.25) can also be used to construct the asymptotic expansions similar to (I.18) and (I.19) for small $u$,

[TABLE]

and large $u$,

[TABLE]

respectively.

II The Weak Coupling Regime

In this Section we describe how to solve both the RB equation and the GB equation using PT in powers of the effective coupling constant $\lambda$. In particular, for the RB equation the perturbative solution can be obtained applying the so-called Non-Linearization Procedure Turbiner (1984). Due to the importance of this procedure, we begin Section by giving a brief description of this method taken as example the $D$-dimensional radial potential.

II.1 The Non-Linearization Procedure

Let us take the RB equation (I.15)

[TABLE]

with potential $V(v;\lambda)$ which admits a Taylor expansion

[TABLE]

Here $\lambda$ is formal parameter and the coefficient functions $V_{n}(v)$, $n=0,1,...,$ are real functions in $v$. Let us now assume that the unperturbed equation at $\lambda=0$,

[TABLE]

can be solved explicitly. It can always be achieved via the inverse problem: we take some function $\mathcal{Y}_{0}(v)$ and then calculate r.h.s.: the potential $V_{0}(v)$ and $\varepsilon_{0}$. Evidently, once we know $\mathcal{Y}_{0}(v)$, the wave function can be found,

[TABLE]

It leads to a constraint on the choice of $\mathcal{Y}_{0}(s)$: the function $\Psi_{0}(v)$ should be normalizable. Now, we can develop PT in powers of $\lambda$,

[TABLE]

and

[TABLE]

Substituting (II.5) and (II.6) into the RB equation (II.1) it is easy to see that the $n$th correction $\mathcal{Y}_{n}(v)$ satisfies the first order linear differential equation

[TABLE]

where

[TABLE]

Evidently, the solution for $\mathcal{Y}_{n}(v)$ is given by

[TABLE]

Since we are interested in finding bound states we have to impose the condition of the absence of current of particles for both $v\rightarrow 0$ and $v\rightarrow\infty$ as boundary condition Turbiner (1984),

[TABLE]

It can be easily checked that if $v$ tends to zero the condition (II.10) is satisfied automatically while at $v=\infty$ the correction $\varepsilon_{n}$ should be chosen accordingly

[TABLE]

where

[TABLE]

to satisfy (II.10). For $D=1$, this perturbative approach was called the Non-Linearization Procedure Turbiner (1984). Here we have given the straightforward extension of it to radial potentials in arbitrary $D$. We will continue to call it the Non-Linearization Procedure. In contrast with the Rayleigh-Schrödinger PT, the knowledge of entire spectrum of the unperturbed problem is not required to find constructively perturbative corrections in (II.5) and (II.6). It is sufficient to know the unperturbed ground state wave function of the unperturbed problem to which we are looking for corrections. This approach gives the closed analytic expression for both corrections $\varepsilon_{n}$ and $\mathcal{Y}_{n}(v)$ in form of nested integrals. Therefore, this procedure is the efficient method to calculate several orders in PT, see e.g. Turbiner (1984) and del Valle and Nader (2018). Interestingly, in this framework the convergence of the perturbation series (II.5) and (II.6) is guaranteed if the first correction $\mathcal{Y}_{1}(v)$ is bounded

[TABLE]

for any $D$, for discussion see Turbiner (1984).

We have presented a brief review of the Non-Linearization Procedure applied to the ground state. This approach can be modified to study the excited states by admitting a number of simple poles with residues equal to one in $\mathcal{Y}_{n}(v)$. The pole positions are found in PT in $\lambda$. Explicit formulas can be found in Turbiner (1984) for one-dimensional case. The generalization for the $D$-dimensional radial potential case is straightforward.

II.2 The Weak Coupling Expansion from the Riccati-Bloch Equation

Now, let us consider the $D$-dimensional radial anharmonic oscillator potential (I.11) in weal coupling regime assuming that the effective coupling constant $\lambda=\left(\frac{\hbar^{2}}{2M}\right)^{\frac{1}{4}}\,g$ is small. In order to apply the Non-Linearization Procedure of Subsection II.1 to construct the weak coupling expansion we choose in (II.2)

[TABLE]

hence, our potential is a finite degree polynomial in both $\lambda$ and $v$. The corresponding unperturbed equation (II.3) describes the $D-$dimensional spherical harmonic oscillator: an exactly solvable problem in any $D>0$, see e.g. Lynch and Mavromatis (1990). The unperturbed ground state wave function is given by

[TABLE]

In general, the explicit calculation of the perturbative corrections $\mathcal{Y}_{n}(v)$ and $\varepsilon_{n}$ can be carried out analytically by using (II.11), final expressions involve explicitly the incomplete gamma function. Let us find the asymptotic expansions of $n$th correction $\mathcal{Y}_{n}(v)$ in two limits: $v\rightarrow 0$ and $v\rightarrow\infty$. For small $v$ the correction $\mathcal{Y}_{n}(v)$ is given by the Taylor expansion

[TABLE]

where $b_{k}^{(n)}$ are some coefficients. It can be easily shown that the first coefficient $b_{0}^{(n)}$ is related to the energy correction $\varepsilon_{n}$,

[TABLE]

while the next coefficient always vanishes,

[TABLE]

due to the absence of the linear term in the potential (I.11). As for large $v$ at $n>0$ we have a Laurent series expansion,

[TABLE]

where $c_{k}^{(n)}$ are some coefficients. Interestingly, the leading coefficient $c_{0}^{(n)}$ does not depend on $D$ while the next-to-leading one always vanishes,

[TABLE]

At $D=0$ all corrections $\varepsilon_{n}$ vanish,

[TABLE]

It can be demonstrated by using the expression (II.12) for $\varepsilon_{n}$, and the asymptotic expansions (II.16) and (II.19). When $D\rightarrow 0$ the numerator of (II.12) is bounded,

[TABLE]

while the denominator has the asymptotic series expansion

[TABLE]

Consequently, when $D\rightarrow 0$ the correction $\varepsilon_{n}$ vanishes linearly, the coefficient $b_{0}^{(n)}$ remains finite, see (II.17). If $\varepsilon_{n}$ vanishes for all $n$, then their formal sum (II.5) results in $\varepsilon=0$ and, ultimately, $E=0$. In general, for $D\neq 0$ it can be shown that expansion for $\varepsilon$ is asymptotic: $\varepsilon_{n}$ grows factorially as $n\rightarrow\infty$, see e.g. Bender and Wu (1969, 1973); Banks and Bender (1972); Silverstone et al. (1985). Therefore, series (II.5) and (II.6) are divergent in $\lambda$. For some particular cases of (1), it has been proved that perturbation series for the energy can be summed to the exact one by using different techniques: calculating the Borel sum Graffi et al. (1970), taking Padé approximants Loeffel et al. (1969) and via renormalization Killingbeck (1981).

An interesting situation occurs when all odd monomial terms in $r$ in potential (I.11) are absent, i.e. the potential is (formally) an even function, $V(r)=V(-r)$. In this case, all odd corrections $\mathcal{Y}_{2n+1}(v)$ and $\varepsilon_{2n+1}$ vanish. The even correction $\mathcal{Y}_{2n}(v)$ has the form of a polynomial of finite degree

[TABLE]

where

[TABLE]

This implies that $\varepsilon_{2n}$ and $\mathcal{Y}_{2n}(v)$ can be calculated by linear algebra means. In turn, the energy correction $\varepsilon_{2n}$ is a finite-degree polynomial in $D$,

[TABLE]

where $d_{k}^{(2n)}$ are real coefficients. However, it is enough to have a single odd monomial term $a_{2q+1}v^{2q+1}$ in the potential (I.11) (naturally, $q$ is integer) to break the features (II.23) - (II.25). Hence, all coefficients in front of the singular terms in the expansion (II.19), as well as all higher order terms in the expansion (II.16) at $k>n$, are proportional to $a_{2q+1}$.

II.3 The Weak Coupling Expansion from the Generalized Bloch Equation

Taking in the GB equation (I.25) the effective coupling constant $\lambda=\left(\frac{\hbar^{2}}{2M}\right)^{\frac{1}{4}}\,g$ as a formal parameter we can develop PT in powers of $\lambda$. The expansion of the eigenvalue $\varepsilon$ in powers of $\lambda$ remains the same as for the RB equation, see (II.5),

[TABLE]

while $\mathcal{Z}(u)$ is of the form of the Taylor series

[TABLE]

We assume that the perturbative energy corrections $\varepsilon_{n}$ are already known: they can be found using the Non-Linearization Procedure, the standard Rayleigh-Schrödinger PT or any other suitable method. It is immediate to see that correction $\mathcal{Z}_{n}(u)$ is calculated by algebraic means and depends on corrections of smaller order,

[TABLE]

Note that the boundary condition $\mathcal{Z}(\infty)\,=\,+\infty$ implies the positive sign in the expression for $\mathcal{Z}_{0}(u)$ should be chosen.

Making the analysis of $\mathcal{Z}_{n}(u)$ correction (II.27) one can see that the boundary condition at $u=0$, see (I.26), can not be fulfilled in the case of arbitrary odd anharmonic radial potential (I.11) once the expansion (II.26) for $\mathcal{Z}(u)$ is used. In general, any $\mathcal{Z}_{n}(u)$ (and its derivative) at $n>1$ diverges at small $u$. For instance, for the cubic potential where $a_{3}\neq 0$ the $n$th correction $\mathcal{Z}_{n}(u)$ with $n>1$ behaves like

[TABLE]

when $u\rightarrow 0$. In turn, for the quintic potential with $a_{3}=0$ but $a_{5}>0$ the $n$th correction at $n>3$ behaves like

[TABLE]

when $u$ tends to zero. For the polynomial in $u$ potential of degree $(2k+1)$, if all other odd terms are absent $a_{3}=a_{5}=\ldots=a_{2k-1}=0$ but $a_{2k+1}>0$, the function $\mathcal{Z}_{n}(u)$ behaves like

[TABLE]

for small $u$ as long as $n>2k-1$. It is the clear indication that the radius of convergence of the expansion of $\mathcal{Z}(u)$ (see (II.26)) in $1/u$ is finite. However, if the anharmonic potential $V(r)$ is even, $V(r)=V(-r)$ the boundary condition at $u=0$ can be satisfied. It might be an indication that the radius of convergence in $1/u$ is infinite. It allows us to determine the correction $\varepsilon_{n}$ by imposing the boundary condition

[TABLE]

At this point, it is worth remarking that the expansion in powers of $\lambda$ (II.26) as well as (II.6) is divergent for any anharmonic oscillator due to the so-called Dyson instability, see Dyson (1952) and as for the discussion Turbiner (1984). This fact can be seen in the partial sums of (II.26), see Figs. 1 and 2.

Let us emphasize that the presence of a single odd degree monomial term in $u$ in the potential $V(u)$ breaks the feature (II.31): the expansion for $\mathcal{Z}(u)$ in powers of $\lambda$ only satisfies the single boundary condition at $u=\infty$. To explain why this happens, it is necessary to make a connection of the expansions (II.6) and (II.26).

II.4 Connection between $\mathcal{Y}\ \mbox{and}\ \mathcal{Z}$ expansions

Thus far we have constructed two different representations for the logarithmic derivative $y(r)$, see (I.10). From one side we have

[TABLE]

from the another side

[TABLE]

where the coupling constant $\lambda$ is defined in (I.16). The connection between them can be established if we use explicitly the expansion of $\mathcal{Y}_{n}(v)$ at large $v$, see (II.19). In this case $y(r)$ is given by

[TABLE]

or, equivalently

[TABLE]

Comparing the expansions (II.34) and (II.35) we can conclude that

[TABLE]

due to the uniqueness of the Taylor series. Therefore, $\mathcal{Z}_{n}(u)$ is a generating function (!) for the coefficients $c_{k}^{(n)}$, $k=0,1,...$ , see for graphical illustration Fig.3.

Equation (II.36) displays the reason why we cannot satisfy, in general, the boundary condition at $u=0$: $\mathcal{Z}_{n}(u)$ is constructed from the expansion (II.19) which is valid in the limit $r\rightarrow\infty$ ($u\rightarrow\infty$ for fixed $g$). In fact, some interesting properties of the corrections $\mathcal{Z}_{n}(u)$ occur when considering $r\rightarrow\infty$. For example, the first term in (II.33) has the expansion

[TABLE]

This expansion reproduces exactly the dominant asymptotic behavior of the function $y(r)$ up to the term $O\left(r^{(m-4)/2}\right)$, c.f. (I.19). Note that the next generating function $\mathcal{Z}_{1}(u)=0$ in agreement with (II.20). It is evident that the next generating functions in (I.19): $\mathcal{Z}_{2}(u)$, $\mathcal{Z}_{3}(u)$, …, also contribute to the expansion of $y(r)$. It will be discussed in detail for the cubic anharmonic oscillator, see below.

Let us conclude this Section by clarifying the perturbative approach for the GB equation and its connection with the semiclassical WKB approximation. By integrating in $r$, the expansion (II.33) can be converted into an expansion for the phase

[TABLE]

where $G_{n}(r)$ is equal to

[TABLE]

Note that keeping $g$ fixed, (II.38) can be regarded as a semiclassical expansion of the phase in powers of $\hbar^{\frac{1}{2}}$, see (I.16). The leading order term is

[TABLE]

for $V\geq 0$ it is the classical action at $E=0$ and it corresponds to the zero order of the standard one-dimensional WKB method (developed in the classically forbidden region). This is the only term in expansion (II.38) which contains no dependence on $D$. The second order term in (II.38) is

[TABLE]

Except for the integral - the third term - it looks like the first order correction to WKB at $E=0$ and $D=1$. Thus, it defines the determinant Escobar-Ruiz et al. (2017). The appearance of the extra term in the form of integral can be explained as follows. Let us take the zero order term in the standard WKB method (in the forbidden region) and expand it in powers of $\hbar^{1/2}$ using (I.14) and (I.16), thus

[TABLE]

Note that the first term is nothing but $G_{0}(r)$ while the second term is exactly the integral in (II.41). Consequently, one can see that (II.38) is the standard semiclassical WKB expansion in the classically-forbidden region, re-expanded in powers of $\hbar^{1/2}$. Needless to say, higher order generating functions $G_{3}(r),G_{4}(r),\ldots$ are related with higher order corrections in the WKB expansion in a similar way. They define two-, three-, etc loop contributions in the path integral formalism in Euclidian time for the density matrix at coinciding points.

Note that this is alternative consideration to the one-dimensional case presented in Escobar-Ruiz et al. (2016, 2017); Shuryak and Turbiner (2018). In these papers it was described how to obtain expansion (II.38) using a path integral formulation (in Euclidean time) considering a special classical solution at zero energy $E=0$ called flucton. In quite straightforward way this consideration can be extended to the $D$-dimensional radial case, where the flucton path appears in radial direction. Hence, the perturbative approach developed for the GB Equation can be called the true semiclassical approximation: both $\varepsilon(\lambda)$, and $\mathcal{Z}(u;\lambda)$ and $\Phi(r;\lambda)$ are expansions in powers $\hbar^{1/2}$.

III The Strong Coupling Expansion

In the previous Section II assuming an small effective constant $\lambda=\left(\frac{\hbar^{2}}{2M}\right)^{\frac{1}{4}}g$ (I.16) we studied the ground state energy $E$ and the phase of the wave function in the form of a Taylor expansion at $\lambda=0$. As a result the harmonic oscillator was taken as unperturbed problem. We now consider the opposite limit of large $\lambda$ taking as the unperturbed potential the leading monomial term $a_{m}g^{m-2}r^{m}$ in (I.11). This situation can be achieved if we consider the expansion in fractional powers $1/\lambda$ - the strong coupling expansion. It implies that all other terms in the potential (I.11), except for the leading one are taken as perturbation.

III.1 The Strong Coupling Expansion from the Riccati Equation

In order to develop the strong coupling expansion in the framework of the Non-Linearization Procedure, let us first take the radial Schrödinger equation (I.8) with the polynomial potential (I.11) and perform the Symanzik scaling transformation

[TABLE]

with idea to remove the coupling constant in front of the leading term in potential. As the result we arrive at the radial Schrödinger equation,

[TABLE]

with a potential $W(r)$ given by

[TABLE]

and rescaled energy $\tilde{E}$, where

[TABLE]

Note that now the parameter $\gamma^{-1}$ plays the role of the coupling constant for $W(r)$. By using the exponential representation (3), the equation (III.2) can be transformed into the Riccati equation (I.9) with potential $W(r)$ and energy $\tilde{E}$ , explicitly,

[TABLE]

where $y(r)$ is defined as in (I.10). With appropriate choice of variable

[TABLE]

with new unknown function and energy

[TABLE]

the equation (III.5) becomes the RB Equation

[TABLE]

c.f. (I.15) with different r.h.s., where

[TABLE]

see (III.4), plays a role of effective coupling constant, and $\tilde{\varepsilon}(\tilde{\lambda})$ plays a role of energy. Two remarks should be made: $(i)$ this equation does not contain explicit dependence on $\hbar$ and $M$; $(ii)$ the parameter $\tilde{\lambda}^{-1}$ is the effective inverse coupling constant, it appears instead of $\gamma^{-1}$. The latter implies that for $\tilde{\mathcal{Y}}(w)$ and $\tilde{\varepsilon}(\tilde{\lambda})$ it should eventually be developed the expansions

[TABLE]

and

[TABLE]

respectively, where $\tilde{\varepsilon}_{n}$, $n=0,1,...$ , are the coefficients of the strong coupling expansion. It is evident that similar expansion (III.11) occurs for the excited states. Contrary to the weak coupling expansion it has been demonstrated explicitly that, at least, for some particular cases of the potential (I.11) PT (III.11) has finite radius of convergence, see e.g. Turbiner and Ushveridze (1988) and references therein. It is related with the fact that the dominant term in the potential (I.11) is $r^{m}$ at $r\rightarrow\infty$ while other monomials are subdominant. Therefore, the phenomenon of the Dyson instability Dyson (1952) does not occur.

It has to be emphasized that the equation (III.8) has the same form as the Riccati equation (I.9), this implies that all the formulas derived in Section II.1 in the framework of the Non-Linearization Procedure are valid and can be used to calculate $\tilde{\mathcal{Y}}_{n}(w)$ and $\tilde{\varepsilon}_{n}$. Let us choose the potential (II.2) in the following form

[TABLE]

In this case, the unperturbed equation (II.3) is given by

[TABLE]

If the solution of this equation is known, the corresponding square-integrable unperturbed wave function $\Psi_{0}$ is obtained using (II.4). Then we use the formulas (II.11) and (II.12) to construct PT. As a result we obtain the coefficients $\tilde{\varepsilon}_{k}$ of the strong coupling expansion. However, since the exact solution of (III.13) for $m>2$ is unknown, we need to find $\tilde{\mathcal{Y}}_{0}(w)$ approximately. In the spirit of the present work, there are two possible ways to find it: $(i)$ one can choose a physically relevant trial function for the potential,

[TABLE]

and use the Non-Linearization Procedure assuming the convergence of PT; $(ii)$ we can construct a locally-accurate approximation of $\tilde{\mathcal{Y}}_{0}(w)$ for all $w\geq 0$. Let us describe both approaches.

III.1.1 Physically Relevant Trial Function

Probably, the simplest physically relevant trial function, see Turbiner (1984) for discussion, for the potential (III.14) was proposed for the first time in Turbiner (1979, 1981) handling the one-dimensional case. Undoubtedly, the similar trial function is also appropriate for the $D$-dimensional case of the radial potentials,

[TABLE]

where the second subscript marks the order of the correction to $\Psi_{0}$. By solving the inverse problem one can find the potential $W_{0,0}(w)$ for which $\Psi_{0,0}$ is the exact solution,

[TABLE]

Taking the difference between the potentials (III.14) and (III.16),

[TABLE]

one develop PT as in (II.5) and (II.6) with $W_{0,1}(w)$ as the perturbative potential for $W_{0,0}(w)$. In this manner, it allows us to calculate the zeroth order coefficient $\tilde{\varepsilon}_{0}$ in strong coupling expansion (III.11) iteratively, in the form of an expansion

[TABLE]

The first two coefficients can be found analytically,

[TABLE]

while several next ones can be obtained numerically. The function $\tilde{\mathcal{Y}}_{0}(w)$ is expanded in a similar way as in (III.18),

[TABLE]

where $\tilde{\mathcal{Y}}_{0,0}(w)$ is given by (III.15). Evidently, this expansion can be converted into the expansion of the phase,

[TABLE]

where $\Phi_{0,0}=a_{2p}^{1/2}\,\hbar\,{\dfrac{w^{p+1}}{p+1}}$. Using $\Psi_{0,0}$ as entry one can also calculate explicitly the first approximation to $\tilde{e}_{1}$ in the strong coupling expansion,

[TABLE]

The next terms, e.g. $\tilde{e}_{1,1}$, $\tilde{e}_{1,2}$, $...$, can be only written in the form of nested integrals, similar to those as in expression (II.12) and their computation is a numerical procedure. Of course, this procedure can be extended to calculate higher coefficients in the expansion (III.11).

III.1.2 Interpolation

One can construct another physically relevant trial function for the RB Equation (III.13). Let us find the asymptotic behavior of $y(r)$ at small $w$,

[TABLE]

[TABLE]

and large $w$

[TABLE]

[TABLE]

The first term in (III.25) chosen as zero approximation corresponds to (III.15). We have to note that large $w$ behavior can be obtained in two ways: taking either large $r$ at fixed $\hbar$, or small $\hbar$ at fixed $r$; therefore the expansion (III.25) should be valid in the semiclassical regime, see (III.6). A more sophisticated approximation for $y(r)$, and consequently for $\tilde{\mathcal{Y}}_{0,0}(w)$, is the result of the interpolation between expansions (III.24) and (III.25). The construction of such an interpolation is presented in the Section IV. As for the cubic anharmonic oscillator, we show how this interpolation leads to a locally accurate approximation.

III.2 The Strong Coupling Expansion from the Generalized Bloch Equation

There is an interesting connection between the strong coupling expansion developed for the GB equation (I.25) and the expansion of $y(r)$ at small $r$, see (I.18). For the case of polynomial potential of degree $m$ (I.11), let us take the equation (I.9), make change of variable (I.23), $u=gr$, then, following (III.1) and (III.2), make rescaling the energy

[TABLE]

and unknown function

[TABLE]

Finally, we arrive at

[TABLE]

c.f. (I.25), where the effective coupling constant $\tilde{\lambda}$ is given by (III.9),

[TABLE]

with

[TABLE]

see (III.4). This equation is the GB Equation for the strong coupling regime: after multiplication by $\tilde{\lambda}^{2}$ the l.h.s. is the same as in (I.25) - the GR equation for weak coupling regime - while the r.h.s. is different,

[TABLE]

Now we can develop PT in powers of inverse effective coupling constant $\tilde{\lambda}^{-1}$

[TABLE]

see (III.11), and

[TABLE]

hence, making expansions at $\tilde{\lambda}=\infty$. We assume the coefficients $\tilde{\varepsilon}_{n}$ are known a priori, in particular, from PT developed for the RB equation at strong coupling regime (III.8): (III.10) and (III.11). Inserting these two series (III.11) and (III.30) in (III.29) one can determine iteratively corrections $\tilde{\mathcal{Z}}_{n}(u),n=0,1,2\ldots$ by solving the same first order differential equation with different r.h.s.

[TABLE]

and

[TABLE]

Imposing the boundary condition $\tilde{\mathcal{Z}}(0)=\tilde{\varepsilon}/D$ to the original equation (III.29), we can find the boundary condition for the $n$th correction,

[TABLE]

The equations (III.31) - (III.32) are solved in elementary way, in particular,

[TABLE]

while the general $n$th correction is given by

[TABLE]

at $n\geq 2$. Note that it is not guaranteed the boundary condition $\tilde{\mathcal{Z}}(\infty)=+\infty$ is fulfilled for all these corrections. It implies that PT in inverse powers of $\tilde{\lambda}$ for $\tilde{\mathcal{Z}}(u)$ has a finite radius of convergence in $u$. Thus, the corrections $\tilde{\mathcal{Z}}_{n}(u)$ make sense in a bounded domain in $u$ only, see discussion below. From (III.34) one can see that any correction is a finite-degree polynomial in $u$,

[TABLE]

where $\alpha_{k}^{(n)}$ are real coefficients and the symbol $[\,\,]$ denotes the integer part. It is easy to see that

[TABLE]

following (III.33), while the next coefficient always vanishes,

[TABLE]

for any $n$. The coefficient $\alpha_{2}^{(n)}$ can be found in the form of finite sum

[TABLE]

while the next coefficient again vanishes,

[TABLE]

for any $n$, c.f. (III.37). From (III.35) one can see that $\tilde{\mathcal{Z}}_{n}(u)$ for $n<m$ and $n=m+1$ is always an odd polynomial in $u$:

[TABLE]

while for $n=m$ this property can be lost. In the case of an even anharmonic potential, $V(r)=V(-r)$, all odd corrections $\tilde{\mathcal{Z}}_{2n+1}(u)$ and $\tilde{\varepsilon}_{2n+1}$ vanish. In turn, all even corrections $\tilde{\mathcal{Z}}_{2n}(u)$ are odd polynomials in $u$ of degree $(2n+1)$ if $2n>m+1$.

The perturbative approach, which was used to calculate the strong coupling expansion of $\tilde{\mathcal{Z}}(u)$, is related to the asymptotic behavior of $\tilde{\mathcal{Z}}(u)$ at small $u$. Once we established the polynomial structure of the corrections $\tilde{\mathcal{Z}}_{n}(u)$, one can sum up some subseries of the perturbative expansion (III.30) keeping the degree of $u$ fixed. For example, the sum

[TABLE]

corresponds to summing up all terms $O(u)$ in (III.30). The next sum is carried out with all terms $O(u^{3})$. Using (III.38) one can show that it corresponds to

[TABLE]

The next sum contains coefficients in front of $u^{4}$ terms. Following (III.39), this sum has a single term,

[TABLE]

This summation procedure can be extended to higher order terms $u^{n},n>4$. After performing such a summations we obtain the coefficients in front of all terms in the Taylor expansion in powers of $u$ for the function $\tilde{\mathcal{Z}}(u)$, and, eventually, for $y(r)$. From another side one can construct the expansion in powers of $u$ directly from (III.28),

[TABLE]

[TABLE]

One can note the exact correspondence between the first three terms of the expansion (III.43) and the subseries previously defined in (III.40), (III.41) and (III.42). In this way we can see the connection between the strong coupling expansion, developed on the basis of the GB equation, and the small $u$ one: both expansions lead to the same representation of $\tilde{\mathcal{Z}}(u)$ and, eventually, of $y(r)$. Finally, this connection clarifies why the boundary condition $\tilde{\mathcal{Z}}(\infty)=+\infty$ can not always be fulfilled.

IV The Approximant

Now we formulate a prescription for the construction of a locally-accurate, uniform approximation of the ground state wave function. We denote this approximation by $\Psi_{0,0}^{(t)}$ and call it the Approximant. Here the subscripts indicate explicitly the quantum numbers of the ground state $(n_{r}=0\,,\,l=0)$. It is convenient to follow the exponential representation of the wavefunction, i.e. $\Psi_{0,0}^{(t)}=e^{-\frac{1}{\hbar}\Phi_{t}}$, since, the phase $\Phi$ of the exponential representation is not that sharp as the wave function $\Psi$, being sufficiently smooth. Hence, we focus on the construction of the trial phase $\Phi_{t}$.

We follow the prescription that the phase $\Phi_{t}$ interpolates between small and large distance $r$ expansions (I.21) and (I.22) and also the expansions of the strong coupling regime (III.24) and (III.25). One can see immediately that the simplest interpolation based on that prescription is of the following form

[TABLE]

cf. Turbiner (2005) for $D=1$ at $m=4$. Without loss of generality we put $\tilde{b}_{2}=1$. Here ${\hat{V}}(r;\{\tilde{a}\})$ and ${\hat{V}}(r;\{\tilde{b}\})$ are modified potentials of the form of the original potential (1), (I.11): instead of the external parameters $\{a\}$ some free parameters $\{\tilde{a}\}$ and $\{\tilde{b}\}$ are taken, respectively. Logarithmic terms can depend on free parameters $\{\tilde{c}\}$ as well. All those parameters will later be fixed by making the variational calculation with (IV.1) taken as the phase of the trial function, see a discussion below. At $g=0$ the phase (IV.1) becomes $a_{0}+a_{2}r^{2}$ that correspond to radial harmonic oscillator. Formula (IV.1) is the key result of this article.

The phase $\Phi_{t}$ has the outstanding property that with an appropriate choice of parameters $(\tilde{a}_{0},\tilde{a}_{1};\tilde{a}_{2},\tilde{b}_{2},\ldots\tilde{a}_{k},\tilde{b}_{k};\{\tilde{c}\})$ the generating function $G_{0}(r)$ is partially (or exactly) reproduced. In particular, in order to reproduce exactly the dominant term in expansion (I.22) the condition

[TABLE]

is imposed. To mimic possible appearance of logarithmic terms in expansion (I.22), in particular, in $G_{0}(r)$ and $G_{2}(r)$, some logarithmic terms to phase $\Phi_{t}$ should be added, see (II.41). These terms, mostly coming from $G_{0}(r)$ and $G_{2}(r)$, can be modified by replacing the concrete (numerical) coefficients by set $\{\tilde{c}\}$ of free parameters. For example, for a concrete two-term anharmonic oscillator potential

[TABLE]

the generating function $G_{0}(r)$

[TABLE]

contains logarithmic terms emerging in the hypergeometric function ${}_{2}F_{1}$ Bateman (1953), while the generating function $G_{2}(r)$ contains logarithmic terms explicitly

[TABLE]

see (II.39). In this case logarithmic terms added in (IV.1) should be a certain minimal modification of the logarithmic terms that appear in $G_{0}(r)$ and $G_{2}(r)$. For instance, a minimal modification of $g^{2}G_{2}(r)/(2M)^{1/2}$ is of the form

[TABLE]

where $\tilde{c}$ is a free parameter, if $\tilde{c}=1$ it coincides with (IV.5).

In general, some constraints on parameters have to be imposed to guarantee that $\Phi_{t}$ has the structure of (I.21). The remaining free parameters are fixed by giving them values of optimal parameters in a trial function used in variational calculus. In order to realize it, we use the Variational Method for the radial Hamiltonian (I.8) taking the Approximant $\Psi_{0,0}^{(t)}$ as a trial function. It is worth noting that the accuracy of variational calculations can be estimated using the connection between Variational Method and PT Turbiner (1980), it will be discussed in next Subsection.

It is worth mentioning that an approximation for the wave functions of excited states can be constructed using as a building block the phase for the ground state $\Phi_{t}$. In $D=1$, the $n$th wave function is labeled by a single quantum number $n$, $n=0,1,2,...$ . One can use an approximant of the form

[TABLE]

where $P_{n}$(r) is a polynomial of degree $n$ with real roots: its coefficients can be used as variational parameters. In $D>1$, due to the separation of variables in spherical coordinates for the Schrödinger equation (I.1), it is sufficient to approximate the radial part of the wave function. One can label the eigenstates by the pair of quantum numbers $(n_{r},\ell)$ omitting the magnetic quantum number. In this case we use an approximant of the form

[TABLE]

where $P^{(\ell)}_{n_{r}}$ is a polynomial of degree $n_{r}$ with positive roots with coeffcients as variational parameters. In both cases, (IV.7) and (IV.8), $\Phi_{t}$ has the same form of the phase, constructed for the ground state (IV.1), but with a different set of parameters $\{\tilde{a}\},\{\tilde{b}\},\{\tilde{c}\}$.

Another way to fix the coefficients in $P_{n}(r)$ is to impose orthogonality constraints: for chosen $n$ $P_{n}(r)$ (IV.7) should be orthogonal to all $k$ previously constructed functions, $k=0,1,...,(n-1)$. It fixes all parameters in (IV.7) while the remaining free parameters are treated as variational. For fixed angular momentum $\ell$ a similar procedure is used for $P^{(\ell)}_{n_{r}}(r^{2})$ and its parameters, the orthogonality constraint is imposed in such a way that the eigenstates with $k_{r}=0,1,\ldots(n_{r}-1)$ are orthogonal to the state with quantum number $n_{r}$. Again, the remaining free parameters are used as variational ones.

IV.1 Variational Method and PT

The realization of the Non-Linearization Procedure does not require the knowledge of the whole spectra of the unperturbed problem. This property gives a freedom to choose the unperturbed potential $V_{0}$ of equation (II.3) by solving the inverse problem. Exploiting this property we are able to find a relationship between PT and variational method Turbiner (1980) (for discussion Turbiner (1984)).

Let us consider square-integrable nodeless trial function $\Psi^{(t)}(r)$ at $r\in[0,\infty)$, thus it may correspond to a ground state, such that

[TABLE]

Notice that $\Psi^{(t)}$ is an eigenfunction of the operator (I.7) with potential

[TABLE]

We can always represent the original potential $V$ as the sum

[TABLE]

where formal parameter $\lambda$ is placed to one, $\lambda=1$. One can consider a formal PT in powers of $\lambda$ and eventually set $\lambda=1$. In this context, the variational energy for the radial operator $\hat{h}_{r}$ (I.7) is given by

[TABLE]

This result suggests to interpret the variational energy as the first two terms in a perturbation theory where $(V-V_{t})$ plays a role of perturbation potential. If the function $\Psi^{(t)}$ depends on certain parameters, we have a parameter-dependent variational energy, which can be minimized with respect to them. The variational principle guarantees that $E_{var}$ is above the ground state energy. The parameters that minimize $E_{var}$ are called optimal parameters. Calculating the next terms $E_{2}$, $E_{3}$, … in the expansion (IV.13) using formulas (II.11) and (II.12), one can estimate the accuracy of the variational calculation and perturbatively improve it, if it is convergent. Choosing different trial function $\Psi_{t}$, we define different approximations for the energy $E$.

Taking partial sums we define approximations of different orders. In particular, the variational energy

[TABLE]

and

[TABLE]

correspond to first and second order approximations, respectively. In general, the partial sum

[TABLE]

defines the $n$th approximation. If we calculate corrections for $y_{0}=\Phi_{t}^{\prime}$ in the framework of the Non-Linearization procedure, we can estimate the accuracy of $\Phi_{t}$ and ultimately the accuracy of the Approximant used as trial function. Moreover, if the correction $y_{1}$ is bounded we know that (IV.16) converges as $n\rightarrow\infty$ Turbiner (1984). Needless to say, this connection between the Variational Method and PT is also valid for excited states Turbiner (1980).

V Cubic Anharmonic Oscillator

The simplest radial anharmonic oscillator is characterized by a cubic anharmonicity,

[TABLE]

c.f. (1), (IV.3) at $m=3$. This potential is the subject of the forthcoming Section. It is worth noting before to proceeding to presentation that many features that the cubic anharmonic potential exhibits are also present in the general potential $V(r)\neq V(-r)$.

V.1 PT in the Weak Coupling Regime

For the cubic anharmonic oscillator the perturbative expansions of $\varepsilon$ and $\mathcal{Y}(v)$ (II.6), (II.5) remain the same functionally as for the general oscillator,

[TABLE]

with $\varepsilon_{0}=D$ and

[TABLE]

respectively, see (I.13), (I.14). If it is assumed $2M=\hbar=1$, then $\lambda=g$, $v=r$ and $\varepsilon=E$, $\mathcal{Y}=y$. Since potential (V.1) contains odd degree monomial $r^{3}$, the correction $\mathcal{Y}_{n}(v)$ is characterized by the infinite series (II.16) and (II.19) at small and large $v$, respectively. The first three terms in expansion of $\mathcal{Y}_{n}(v)$ at $n=1,2,3$ at both $v\rightarrow 0$ and $v\rightarrow\infty$ asymptotics are presented in Appendix A. The first corrections $\varepsilon_{1}$ and $\mathcal{Y}_{1}(v)$ can be calculated in closed analytic form using (II.11) and (II.12),

[TABLE]

where $\Gamma(a)$ and $\gamma(a,b)$ denote the (in)complete gamma functions, respectively, see Bateman (1953). Note that the higher energy corrections $\varepsilon_{n}$ can be computed only numerically. However, the behavior of $\varepsilon_{n}$ at large $D$ can be obtained in $1/D$-expansion analytically. This expansion can be constructed via the Non-Linearization Procedure, it is not be presented in this paper.

V.2 Generating Functions

We can determine the coefficients of $\mathcal{Y}_{n}(v)$, i.e. $c_{k}^{(n)}$ in (II.19), by algebraic means. In general, they are written in terms of the corrections $\varepsilon_{0}$, $\varepsilon_{1}$, $\ldots$, $\varepsilon_{n}$, see Appendix A for explicit formulas for $n=1,2,3$.

It was pointed out in Dolgov and Popov (1979) as well as in Dolgov and Popov (1978), see also Turbiner (1980, 1984), that the general expression for several first coefficients $c_{0}^{(n)}$, $c_{2}^{(n)},\ldots$ can be obtained by solving recurrence relations. For example, $c_{0}^{(n)}$ satisfies non-linear recurrence relation

[TABLE]

as for $c_{2}^{(n)}$ the recurrence relation is linear

[TABLE]

the solution of (V.4) requires the knowledge of the coefficients $c_{0}^{(n)}$. It is evident that in order to find $c_{3}^{(n)}$, the coefficients $c_{0}^{(n)}$, $c_{2}^{(n)}$ are needed, recurrence relation remains linear, etc. The simplest way to solve these recurrence relations is by using the generating functions. Surprisingly, we can calculate these generating functions straightforwardly via the algebraic, iterative procedure derived from the GB equation, see (II.27): we demonstrated that $k$th generating function is nothing but $\mathcal{Z}_{k}(u)$! For example, the first two generating functions are

[TABLE]

where $u=(gr)$. From these expressions, the explicit solutions of the equations (V.3) and (V.4) can be derived,

[TABLE]

Note that the coefficients $c_{0}^{(n)}$ are $D$-independent, $c_{2}^{(n)}$ depends on $D$ linearly. It can be demonstrated by induction that $c_{k}^{(n)}$ is polynomial in $D$ of degree $(k-1)$.

As mentioned before, an important property of the generating functions $\mathcal{Z}_{0}(gr)$, $\mathcal{Z}_{2}(gr)$, $\mathcal{Z}_{3}(gr)$, … is related to the asymptotic behavior of the function $y$ at large $v$. From (I.19) we have

[TABLE]

[TABLE]

see (I.13). Note that the first four terms of this expansion are independent of $\varepsilon$. The expansion of $(2M)^{1/2}\mathcal{Z}_{0}(gr)$ at large $r$ can be transformed into an expansion for large $v$, namely

[TABLE]

It reproduces exactly the expansion (V.8) up to $O(v^{-1/2})$. At large $r$, the function $(2M)^{1/2}\lambda^{2}\mathcal{Z}_{2}(gr)$ contributes to the expansion (V.8) in a similar way starting from $O(v^{-1})$, since

[TABLE]

Note that in the expansion of $(2M)^{1/2}\mathcal{Z}_{0}(gr)$ plus $(2M)^{1/2}\lambda^{2}\mathcal{Z}_{2}(gr)$ at large $v$ reproduces exactly the coefficient in front of $O(v^{-1})$ in expansion (V.8). The expansion of next function $(2M)^{1/2}\lambda^{3}\mathcal{Z}_{3}$ starts from $O(v^{-3/2})$. Together with the first two $\mathcal{Z}_{0,2}$ it reproduces the coefficient in front of $O(v^{-3/2})$ term exactly.

V.3 The Approximant and Variational Calculations

The first two generating functions, $G_{0}(r)$ and $G_{2}(r)$ in expansion (II.38) are given by

[TABLE]

The next two functions $G_{3}(r)$ and $G_{4}(r)$ are presented explicitly in Appendix A. The function $G_{0}(r)$ contains no logarithmic terms.

The generating functions $G_{0}(r)$ and $G_{2}(r)$ serve the building blocks for the construction of the Approximant. Following the general prescription (IV.1) we write the Approximant in the exponential representation $\Psi_{0,0}^{(t)}=e^{-\Phi_{t}}$ with

[TABLE]

where we set

[TABLE]

in order to reproduce exactly the leading term in asymptotic behavior of phase at $r\rightarrow\infty$. Here, the logarithmic terms in $\Phi_{t}$ represent a minimal modification of ones in $G_{2}(r)$, cf. (V.12). Additionally, we impose a constraint

[TABLE]

in order to guarantee vanishing $\Phi_{t}^{\prime}$ at $r=0$, thus, no linear term in expansion of (V.13). Eventually, the Approximant in its final form depends on three free parameters $\{\tilde{a}_{0},\tilde{a}_{2},\tilde{b}_{3}\}$ only, it reads

[TABLE]

where $\tilde{a}_{3,1}$ are defined by constraints (V.14), (V.15). This is the eventual expression for the approximate, 3-parametric ground state wave function, which is key result of Section on cubic radial anharmonic oscillator. As the first step this function is used as a trial function in variational calculation in order to fix 3 free parameters.

As for excited states at $D>1$ $(n_{r},\ell)$ the Approximant has the form

[TABLE]

where $P^{(\ell)}_{n_{r}}(r^{2})$ is the polynomial of degree $n_{r}$ with $n_{r}$ real positive roots and with leading term $r^{2\,n_{r}}$. $n_{r}$ coefficients of $P^{(\ell)}_{n_{r}}(r^{2})$ are found by imposing the $n_{r}$ orthogonality conditions:

[TABLE]

After imposing the orthogonality conditions and constraints (V.14), (V.15) on $a_{3,1}$, the function (V.17) depends eventually on 3 free parameters $\{\tilde{a}_{0},\tilde{a}_{2},\tilde{b}_{3}\}$, which are found by making minimization of the variational energy. For all studied states the optimal parameters $\{\tilde{a}_{0},\tilde{a}_{2},\tilde{b}_{3}\}$ demonstrate smooth behavior versus $g$ and $D$.

Variational energy calculation with trial function $\Psi_{0,0}^{(t)}$ (V.16) requires to perform a numerical integration of the integrals in numerator and denominator (IV.13), and also a numerical minimization. Computational code was written in FORTRAN 90 with use of the integration routine D01FCF from the NAG-LIB, which was built using the algorithm described in Genz and Malik (1980). The optimization routine to find the variational parameters was performed using the program MINUIT of CERN-LIB. Let us mention that the use of Approximant $\Psi_{n_{r},\ell}^{(t)}$ for excited states as trial functions in the form (V.17), c.f. (IV.7) or (IV.8), also require numerical integration and minimization.

Without loss of generality we set $\hbar=1$ and $M=1/2$, thus, putting $v=r$, $\varepsilon=E$, $\lambda=g$ and $\mathcal{Y}=y$, see (I.13), (I.14) and (I.16). The calculations of variational energy for four low-lying states with quantum numbers (0,0), (1,0), (0,1) and (0,2) for different values of $D$ and $g$ are presented in Tables 1 - 4. In general, variational parameters are smooth, slow changing functions in $g$ without strong $D$ dependence. As example, for the ground state (0,0) the plots of the parameters ${\tilde{a}}_{0,2,3}$ vs $g$ for $D=2,3,6$ are shown in Fig. 4.

For some states the variational energy $E_{var}=E_{0}^{(1)}$, the first correction $E_{2}$ to it as well as its corrected value $E_{0}^{(2)}=E_{var}+E_{2}$ are shown. All states are studied for different values of dimension $D$ and coupling constant $g$. For all cases the variational energy $E_{0}^{(1)}$ is obtained with absolute accuracy $10^{-7}-10^{-8}$ (7-8 s.d.). This accuracy is found by calculating the second correction $E_{2}$ (the first correction to the variational energy). Furthermore, making comparison of energies $E_{0}^{(2)}$ with ones coming from the Lagrange Mesh method (see below) one can see that the accuracy 8 d.d. is reached. This accuracy is confirmed independently by calculating the second correction $E_{3}$ to variational energy: this correction is $\lesssim 10^{-9}$ for any $D$ and $g$ we have studied. It indicates to very fast rate of convergence in Non-Linearization Procedure with trial function (V.16) taken as zero approximation.

To check the energies obtained in Non-Linearization Procedure we calculated them in the Lagrange Mesh Method Baye (2015). It is known that the Lagrange Mesh Method in the form proposed and developed by D. Baye with co-authors, see Baye (2015) and references therein, for various problems of atomic and molecular physics, provides easily 12 - 13 s.d. in energies. In present work we apply this method to study cubic radial anharmonic oscillator in different dimensions $D$. These studies are performed for one-dimensional case $D=1$ using the Hermite mesh, while for larger dimensions $D>1$ the Laguerre mesh is used. In both occasions 50 mesh points are used as maximum number of points. Details will be published elsewhere.

These numerical calculations show that all digits presented for $E_{0}^{(2)}$ - the variational energy with first (second order in PT) correction $E_{2}$ taken into account - in Tables 1 - 4 are exact(!): we obtain not less than nine decimal digits correctly in Non-Linearization Procedure. In general, the energies grow with increase of $D$ and/or $g$. Classifying states by radial quantum number $n_{r}$ and angular momentum $\ell$ as $(n_{r},\ell)$, we mention the hierarchy of eigenstates which holds for any fixed integer $D$ and $g$: (0,0), (0,1), (0,2), (1,0).

It has to be noted that to the best of the authors knowledge these results represent the first accurate variational calculations (and corrections to them found in the Non-Linearization Procedure) of the energy of the low-lying states for the $D$-dimensional cubic anharmonic oscillator.

The Non-Linearization Procedure allows us to estimate a deviation of the Approximant $\Psi_{n_{r},l}^{(t)}$ from the exact wave function $\Psi_{n_{r},l}$. As for the ground state the relative deviation is bounded and very small,

[TABLE]

in the whole range of $r\in[0,\infty)$ at any dimension $D$ we explored and at any coupling constant $g\geq 0$ studied. Therefore, the Approximant leads to a locally accurate approximation of the exact wave function $\Psi_{0,0}(r)$ once the optimal parameters for $\Psi_{0,0}^{(t)}$ are chosen. These optimal parameters are smooth slow-changing functions in $D$ and $g$. Similar situation appears for the excited states for different $D$ and $g$.

The Approximant $\Psi_{1,0}^{(t)}$ also provides an accurate estimate of the position of radial node of the exact wave function, see (V.17). The orthogonality condition imposed to $\Psi_{0,0}^{(t)}$ and $\Psi_{1,0}^{(t)}$ provides a simple analytic expression for the zero order approximation $r_{0}^{(0)}$ to $r_{0}$. Making comparison of $r_{0}^{(0)}$ with numerical estimates which come from the Lagrange Mesh Method with 50 mesh points, we can see the coincidence of $r_{0}$ and $r_{0}^{(0)}$ for at least 5 d.d. at integer dimension $D$ at any coupling constant $g\geq 0$. Results are presented in Table 4.

In all cases studied the first correction $y_{1}$ to the logarithmic derivative of the ground state function is not bounded, however, the ratio $|y_{1}/y_{0}|$ is bounded and small, thus, $y_{1}$ is a small function in comparison with $y_{0}$. In Figs. 4 - 7 it is presented $y_{0}$ and $y_{1}$ vs. $r$ for $g=1$ in physical dimensions $D=1,2,3$. Similar plots appear for $D=6$. Making analysis of these plots one see that in domain $1\gtrsim r\geq 0$, which gives dominant contribution to integrals defining the energy corrections, the $|y_{1}|$ is extremely small comparing to $|y_{0}|$ being close to zero. It explains why the energy correction $E_{2}$ is small being the order of $\sim 10^{-7}$, or $\sim 10^{-8}$. In similar way one can show numerically that the higher corrections $y_{2},y_{3},\ldots$ drop down to zero in the domain $1\gtrsim r\geq 0$ even faster indicating the convergence both $y_{n}$ and $E^{(n)}$ as $n\rightarrow\infty$.

V.4 The Strong Coupling Expansion

Here we present the results for the first two terms in the strong coupling expansion (III.11) for the ground state energy of cubic oscillator,

[TABLE]

assuming for simplicity $2M=\hbar=1$. It is worth mentioning that similar expansion holds for any excited state. Evidently, it has finite radius of convergence. This expansion corresponds to PT in powers of ${\hat{\lambda}}$ for potential

[TABLE]

defined in $w\in[0,\infty)$.

As the first step we focus on calculation of ${\tilde{\varepsilon}}_{0}$ in (V.19), which is, in fact, the ground state energy in the pure cubic radial potential $V=gr^{3}$ (ultra-strong coupling regime). In order to do it we begin with one of the simplest physically relevant trial function (III.15),

[TABLE]

for the ground state in potential (V.20) at ${\hat{\lambda}}=0$. Taking $\Psi_{0,0}^{(t)}$ as zero approximation we calculate variational energy analytically and develop the PT procedure, described in Section III.1, numerically. In this manner we obtain different PT corrections $\tilde{\varepsilon}_{0,k}\,,\ k=0,1,2,3,\ldots$ to the exact value of $\tilde{\varepsilon}_{0}$ and then form the partial sums (III.18). The first non-trivial partial sum $(\tilde{\varepsilon}_{0,k}+\tilde{\varepsilon}_{0,1})$ is equal to variational energy. Explicit results including partial sums up to sixth order are presented in Table 5 for the one-dimensional case, $D=1$. We have to note a slow convergence to the exact result with increase of the number of perturbative terms in expansion (III.18) taken into account. Even the partial sum $\sum_{k=0}^{6}\tilde{\varepsilon}_{0,k}$ with seven terms included allows us to reproduce 7 s.d. only of the highly-accurate result obtained in Lagrange mesh method with 50 mesh points. It is natural to accelerate convergence taking more advanced trial function as entry. As for the second term in the expansion (V.19) it can be estimated by calculating with (V.21) the normalized expectation value of $w^{2}$,

[TABLE]

This straightforward estimate with simplest trial function (V.21) provides ${\tilde{\varepsilon}}_{1}$ with accuracy $\sim 20\%$, see below. In similar way how it was done for ${\tilde{\varepsilon}}_{0}$ one can develop PT for ${\tilde{\varepsilon}}_{1}$ in (V.19). This PT is fast convergent: the first PT correction to (V.22) is equal to $-0.086$: it leads to ${\tilde{\varepsilon}}_{1}=0.409$ which differs from the exact value in $\sim 0.5\%$, see Table 7.

As alternative to the trial function (V.21) let us use the Approximant (V.13) in order to calculate the first two terms in the strong coupling expansion (V.19), see also (III.11), for the ground state. In Table 6 we present for different $D$ the leading coefficient of the strong coupling expansion $\tilde{\varepsilon}_{0}$ found variationally $\tilde{\varepsilon}_{0}^{(1)}$ and the second PT correction $\tilde{\varepsilon}_{2}$ calculated to it via the Non-Linearization Procedure. We introduce the partial sum $\tilde{\varepsilon}_{0}^{(2)}=\tilde{\varepsilon}_{0}^{(1)}+\tilde{\varepsilon}_{2}$. The final results are verified using the Lagrange mesh method with 50 mesh points, see discussion above. One can see that systematically the correction $\tilde{\varepsilon}_{2}$ is of order of $\sim 10^{-7}$. Hence, the first six decimal digits in variational energy are correct and it defines the accuracy of variational calculations of ${\tilde{\varepsilon}}_{0}$ with the Approximant (V.13) as the trial function. One can estimate the order of the third PT correction: all $D$ studied $\tilde{\varepsilon}_{3}\sim 10^{-2}\,\tilde{\varepsilon}_{2}$. It indicates to extremely high rate of convergence $\sim 10^{-2}$ of our PT. In Table 7 we present the first two approximations of the coefficient $\tilde{\varepsilon}_{1}$ in (V.19).

Summarizing the studies of $D$-dimensional radial cubic anharmonic oscillator for a given $D$, one can imagine that the numerical results for given eigenvalue can be described via a simple analytical formula. Inspiration comes from the fact that the eigenvalue vs coupling constant $g$ is very smooth, slow growing function. Basic idea is to interpolate the expansions in $g$ for the weak and strong coupling regimes, see (II.5) and (III.2). It has been already shown that this approach is appropriate and successful in description of various physical systems of different nature, e.g. del Valle and Nader (2018) and Olivares-Pilón and Turbiner (2018). In particular, for the ground state of the cubic radial anharmonic potential, the simplest interpolations is of the form

[TABLE]

where $a$ is a free parameter; we set $b=(\tilde{\varepsilon}_{0}/D)$, it is slow-changing parameter with $D$, see Table 6. This interpolation reproduces exactly both leading terms in (II.5) and (V.19). Parameter $a$ is fixed by fit of numerical data with minimal $\chi^{2}$. Optimal parameter $a$ and the parameter $b$ are presented in Table 8. For integer $D$ the above simple formula describes the ground state energy $E$ for any $g\in[0,\infty)$ with accuracy $\lesssim 2\%$. In a similar manner a straightforwardly modified formula works for energies of excited states.

For one-dimensional quartic oscillator $V=x^{2}+gx^{4}$ it was discovered long ago by Bender and Wu Bender and Wu (1969) and recently proved rigorously by Eremenko and Gabrielov Eremenko and Gabrielov (2009) that even (odd) parity eigenvalues form infinitely-sheeted Riemann surface in space of coupling constant $g$. This surface is characterized by infinitely many square-root branch points having a meaning of the points of level crossings at complex $g$. Each Riemann sheet contains infinitely many such square-root branch points manifesting that every two energy levels intersect. Seemingly, the same phenomenon holds for any one-dimensional anharmonic oscillator with two term potential $V=x^{2}+gx^{2m}$ defined in the whole line $x\in(-\infty,+\infty)$. We guess that for cubic anharmonic oscillator $V=r^{2}+gr^{3}$ at $r\in[0,\infty)$ for fixed integer dimension $D>1$ and each value of the angular momentum $\ell=0,1,\ldots$ the energies with quantum numbers $(n_{r},\ell),n_{r}=0,1,2,\ldots$ form infinitely-sheeted Riemann surface in variable ${\tilde{\lambda}}=g^{-4/5}$ with square-root branch points. It is evident that due to the Hermiticity of the original Hamiltonian (1) each branch point should have its complex-conjugate. It could be a definite challenge to calculate the distance ${\tilde{\lambda}}^{\star}$ from the origin to the nearest square-root branch point in ${\tilde{\lambda}}$ and find its $D$ dependence, since it defines the radius of convergence of strong coupling expansion for two intersecting energy levels $(n_{r},\ell)$ and $({\tilde{n}_{r}},\ell)$. The most important particular case is to find ${\tilde{\lambda}}^{\star}$ for the ground state $(0,0)$ and the first excited state $(1,0)$: it probably defines the radius of convergence of the strong coupling expansion (V.19). It is not clear what intersection of the ground state eigenvalue with one of the excited state is closest to the origin in ${\tilde{\lambda}}$. This challenge will be tackled elsewhere.

Conclusions

There are two views of anharmonic, polynomial, radial oscillator potential: (i) as in coordinate $r$-space it is a finite degree perturbation in coupling constant $g$ being polynomial in $r$ perturbation of the harmonic oscillator and (ii) as finite-degree polynomial in classical coordinate $(gr)$ with coupling constant dependence $1/g^{2}$ in front of the potential. These views suggest to make study the dynamics in these two spaces: $r$-space and in $(gr)$-space. It leads to two different non-linear differential equations for the logarithmic derivative of the wavefunction, both are of the first order with quadratic non-linearity: the RG equation and the GB equation with the same effective coupling constant $\lambda$ and energy $\varepsilon$. Perturbation theory in powers of $\lambda$ developed in RG equation leads effectively to a description of small $r$ domain, while developed in GR equation leads effectively to a description of large $(gr)$ domain. Latter generates a new version of semiclassical expansion in fractional powers of $\hbar$ in a very compact way. In complementary way the perturbation theory in inverse fractional powers of $\lambda$ developed in RG equation leads effectively to a description of large $r$ domain, while developed in GR equation leads effectively to a description of small $(gr)$ domain. Latter generates anti-semiclassical expansion in inverse fractional powers of $\hbar$.

Combining all four expansions together, we are able to construct for logarithmic derivative $y$ (4) the interpolating function $y_{0}$ (IV.1) for both the general case (I.11) and the particular case (V.13) of cubic oscillator (V.1). Eventually, we arrive to 3-parametric trial function for any eigenstate of cubic oscillator. It seems evident that, in general, $y_{0}$ should be very close to the exact $y$ - the first correction $y_{1}$ is small, $|y_{1}/y_{0}|\ll 1$ for any $r\in[0,\infty)$. We show that by taking $y_{0}$ as zero approximation the Non-Linearization procedure appears fast convergent. It is explicitly checked for cubic oscillator that the relative deviation of the Approximant (V.16) from exact ground state eigenfunction is less than $10^{-4}$, see (V.18), for any $r\in[0,\infty)$. Simultaneously, the absolute accuracy in energy reaches the extremely high value $\sim 10^{-7}$ at any dimension $D$ and coupling constant $g\in[0,\infty)$. Similar accuracy can be reached for any expectation value for which the dominant contribution to integrals comes from the domain $r\in[0,1]$. It manifests a solution of the problem of ground state of cubic radial anharmonic oscillator as well as its low-lying excited states! Based on this formalism similar accuracy is reached for quartic and sextic radial anharmonic oscillators that will be presented in subsequent paper del Valle and Turbiner .

All studies performed lead to smooth dependence on dimension $D$. We did not see any indication that physics dimension $D=3$ is special.

Acknowledgments

The authors thank J.C. López Vieyra and H. Olivares Pilón for their interest to the work and useful remarks, and especially for help with numerical and computer calculations. J.C. del V. is supported by a CONACyT PhD Grant No.570617 (Mexico). A.V.T. thanks the Stony Brook University and Simons Center for Geometry and Physics at Stony Brook, NY, and especially, E Shuryak for kind hospitality extended to him, where some of present work was carried out. This work is partially supported by CONACyT grant A1-S-17364 and DGAPA grant IN108815 (Mexico).

Appendix A First PT Corrections and Generating Functions $G_{3,4}$ for the Cubic Anharmonic Oscillator

The expansions in powers of $v$ for the first three corrections $\mathcal{Y}_{n}(v),n=1,2,3$ for the cubic anharmonic oscillator potential (V.1), see (II.16) and (II.19) are presented.

• •

At $v\rightarrow 0$

[TABLE]

In these formulas we have denoted $\epsilon_{n}=\frac{\varepsilon_{n}}{D}$.

• •

At $v\rightarrow\infty$

[TABLE]

A straightforward application of formulas (II.27) and (II.39) allows us to calculate the generating functions. In particular, $G_{3}$ and $G_{4}$ are given by

[TABLE]

where $w=(1+gr)^{1/2}$.