Have you read the basic paper: "Anharmonic Oscillator. II. A Study of Perturbation Theory in Large Order," by Carl M. Bender and Tai Tsun Wu, Phys. Rev. D 7, 1620 (1973)?

That paper shows you how to obtain the series for the quartic case. We were able to find several hundred terms for other QM perurbation series by their method.

See the appendix in: J.Reeve, M. Stone "Late terms in the asymptotic expansion for the energy levels of a periodic potential" Phys.Rev.D 18 (1978) 4746

I have not done it for the cubic, but it is probably not too hard.

In your example, you can do this analytically as the unperturbed potential in the harmonic oscillator, for which there are analytical solutions for the eigenenergies and eigenfunctions.

General formulae

Use perturbation theory and a sensible choice of your unperturbed basis.

The energy $E_n$ will be written as: $$E_n = E_n^{(0)} + E_n^{(1)} + E_n^{(2)} + dots ,$$ where the LHS is the true value and the LHS terms are $n^{mathrm{th}}$ corrections.

Same thing for the wavefunction $psi_n$: $$psi_n = psi_n^{(0)} + psi_n^{(1)} + psi_n^{(2)} + dots .$$

The first order correction for the energy $E_n^{(1)}$ is given by: $$E_n^{(1)} = int psi_n^{(0)*} hat{H}^{(1)} psi_n^{(0)}, $$ the second order is: $$E_n^{(2)} = int psi_n^{(0)*} hat{H}^{(1)} psi_n^{(1)}, $$ and so on.

The first order correction for the wavefunction is: $$ psi_n^{(1)} = sum_{ineq n} psi_n^{(0)} frac{int psi_i^{(0)*} hat{H}^{(1)} psi_n^{(0)}}{E_n^{(0)} - E_i^{(0)}} .$$

$H^{(1)}$ is the perturbation to the Hamiltonian.

You can combine the previous two formulae to re-write the second-order energy correction as: $$ E_n^{(2)} = sum_{ineq n} frac{| int psi_i^{(0)*} hat{H}^{(1)} psi_n^{(0)}|^2}{E_n^{(0)} - E_i^{(0)}} .$$

Your example

I am going to use $hat{H}_1 = gx^3$ for the perturbed potential. The unperturbed potential is the harmonic potential so will use the analytical eigenergies and eigenfunctions for $E_n^0$ and $psi^0_n$.

In your example, let's take $x^2/2$ to the the "basis" potential, which is none other than the harmonic oscillator which known and analytical wavefunctions (so that you can compute the integrals easily). The perturbed contribution $H^{(1)}$ is $gx^3$.

So the total energy $E$ will be the energy from the harmonic oscillator plus the correction due to the perturbed potential.

For the ground state then ($n=0$):

• Zero-th order: $$ E^{(0)} = frac{1}{2}, $$ which is the harmonic (unperturbed potential) contribution.

• First order: $$ E^{(1)} = 0, $$ because $gx^3$ is odd.

• Second order: I did this with Mathematica using the last equation from the previous section. It converges after 3 terms in the sum and it gives: $$ E^{(2)} = -1.375 g^2 = -frac{11}{8}g^2.$$

Then you go on.

--- Addition: ---

Higher-order

I then copied all the 3rd, 4th and higher terms from wikipedia into Mathematica and got:

• Third order: $$ E^{(3)} = 0.$$

• Fourth order (converged in 8 terms): $$ E^{(4)} = -14.5313 g^4 = -frac{465}{32}g^4.$$

1. Here we will for fun try to reproduce the first few terms in a perturbative series for the ground state energy $E_0$ of the 1D TISE $$begin{align} Hpsi_0~=~&E_0psi_0, cr H~=~&frac{p^2}{2}+frac{omega^2}{2}q^2+V_{rm int}(q), cr V_{rm int}(q)~=~&gq^3, cr g~=~&frac{lambda}{6},end{align} tag{A}$$ using an Euclidean path integral in 0+1D $$begin{align} e^{W_c[J]/hbar}~=~&Z[J]cr ~=~&int!{cal D}q ~expleft{ frac{1}{hbar}int_{[0,T]}!dt~(-L_E +J q) right} cr ~=~&expleft{-frac{1}{hbar}int_{[0,T]}!dt~V_{rm int}left(hbar frac{delta}{delta J}right) right} Z_2[J],end{align} tag{B}$$ cf. Refs. 1-3. The Euclidean Lagrangian is $$begin{align} L_E~=~&frac{1}{2}dot{q}^2+frac{omega^2}{2}q^2+V_{rm int}(q), cr q(T)~=~&q(0),end{align} tag{C}$$ with periodic boundary condition.

2. The free quadratic part is the harmonic oscillator (HO) $$ begin{align}Z_2[J]&~=~cr Z_2[J!=!0]&expleft{frac{1}{2hbar}iint_{[0,T]^2}!dt~dt^{prime} J(t) Delta(t,t^{prime})J(t^{prime}) right}.end{align} tag{D}$$ The partition function for the HO can be calculated either via path integrals or via its definition in statistical physics:
   $$begin{align} Z_2[J!=!0]~=~& sum_{ninmathbb{N}_0}e^{-(n+1/2)omega T}cr ~=~&left(2sinhfrac{omega T}{2}right)^{-1}.end{align} tag{E}$$ The free propagator is $$ begin{align} Delta(t,t^{prime})~=~&frac{1}{2omega}e^{-omega |t-t^{prime}|}, cr left(-frac{d^2}{dt^2}+omega^2right)Delta(t,t^{prime}) ~=~&delta(t!-!t^{prime}).end{align}tag{F}$$

3. The main idea is to use the fact that the ground state energy can be inferred from the connected vacuum bubbles
   $$ W_c[J!=!0]~sim~ -E_0 Tquadtext{for}quad Ttoinfty.tag{G}$$ Here we are using the linked-cluster theorem. The Feynman rule for the cubic vertex is $$ -frac{lambda}{hbar}int_{[0,T]}!dt. tag{H}$$

4. The unperturbed zero-mode energy for the HO is the well-known $$E_0(g!=!0)~=~frac{hbaromega}{2},tag{I} $$ cf. eq. (E). Set $hbar=1=omega$ to compare with OP's eq. (3).

5. There are two 2-loop vacuum bubbles: the dumbbell diagram $O!!-!!O$ and the sunset diagram $theta$ with symmetry factor $S=8$ and $S=12$, respectively, cf. Fig. 1.

   

   $uparrow$ Fig. 1. (From Ref. 4.) The two 2-loop vacuum bubbles: the dumbbell diagram $O!!-!!O$ and the sunset diagram $theta$.

   They make up the next-to-leading (NL) order contribution $$ -frac{11hbar^2 g^2}{8omega^4} tag{J} $$ to the ground state energy $E_0$, cf. OP's eq. (3) & Ref. 5.

   Proof of eq. (J): The dumbbell Feynman diagram $O!!-!!O$ is$^1$ $$ begin{align} frac{(-lambda/hbar)^2}{8} left(frac{hbar}{2omega}right)^{3}& iint_{[0,T]^2}!dt~dt^{prime}e^{-omega |t-t^{prime}|}cr ~=~& frac{hbarlambda^2}{32omega^4}T + {cal O}(T^0)cr ~=~& frac{9hbar g^2}{8omega^4}T + {cal O}(T^0).end{align}tag{K}$$ The sunset Feynman diagram $theta$ is$^1$ $$ begin{align} frac{(-lambda/hbar)^2}{12} left(frac{hbar}{2omega}right)^{3}& iint_{[0,T]^2}!dt~dt^{prime}e^{-3omega |t-t^{prime}|}cr ~=~& frac{hbarlambda^2}{144 omega^4}T + {cal O}(T^0)cr ~=~& frac{hbar g^2}{4omega^4}T + {cal O}(T^0).end{align}tag{L}$$ $Box$

6. Five 3-loop vacuum bubbles make up the next-to-next-to-leading (NNL) order, cf. Fig. 2.

   

   $uparrow$ Fig. 2. (From Ref. 4.) The five 3-loop vacuum bubbles.

7. It is in principle possible to calculate to any order by drawing Feynman diagrams. An $n$-loop-integral is analytically doable by breaking the integration region $[0,T]^n$ into $n$-simplexes.

References:

1. M. Marino, Lectures on non-perturbative effects in large $N$ gauge theories, matrix models and strings, arXiv:1206.6272; section 3.1.

2. R. Rattazzi, The Path Integral approach to Quantum Mechanics, Lecture Notes for Quantum Mechanics IV, 2009; subsection 2.3.6.

3. R. MacKenzie. Path Integral Methods and Applications, arXiv:quant-ph/0004090, section 6.

4. M. Srednicki, QFT, 2007; figures 9.1 + 9.2. A prepublication draft PDF file is available here.

5. I. Gahramanov & K. Tezgin, A resurgence analysis for cubic and quartic anharmonic potentials, arXiv:1608.08119, eqs. (3.1) + (3.2).

--

$^1$ It is straightforward to check that $$ iint_{[0,T]^2}!dt~dt^{prime}e^{-omega |t-t^{prime}|} ~=~frac{2T}{omega} + {cal O}(T^0). tag{M}$$

1. In this answer, we will sketch the derivation of a large-$n$ expression$^1$ for the $a_n$-coefficients in the perturbative series $$E_0(g)~=~sum_{ninmathbb{N_0}} a_n g^n tag{A}$$ for the ground state energy by considering a single instanton bounce.

2. The Euclidean Lagrangian is $$begin{align}L_E~=~&frac{1}{2}dot{q}^2+V(q), cr q(T)~=~&q(0),end{align} tag{B}$$ with periodic boundary condition. Let us be a bit more general and assume that the potential is of the form $$ begin{align} V(q)~=~&V_2(q)+V_{rm int}(q),cr V_2(q)~=~&frac{1}{2}omega^2q^2, cr V_{rm int}(q)~=~&g~{rm sgn}(q) {cal V}(|q|)cr ~=~&-V_{rm int}(-q),cr 0~leq~& {cal V}(|q|)~=~{cal O}(|q|^3).end{align} tag{C}$$ This means the potential is unstable for $qtopminfty$, and therefore has an instanton bounce for $q>0$ (or $q<0$) if $g<0$ (or $g>0$), respectively, cf. e.g. Refs. 1-3 and this related Phys.SE post.

3. The Euclidean partition function $Z$ is defined by analytical continuation of the position variable $q$ in the complex plane. Remarkably we only need to assume that this is possible; the details are largely irrelevant. Since initially all the parameters in the Euclidean partition function $Z$ are manifestly real, we expect the analytical continuation to satisfy the Schwarz reflection principle $$2i{rm Im} Z(g)~=~Z(g+i0^+)-Z(g-i0^+).tag{D} $$

4. Let us focus on the case $g<0$ so that the instanton $q_{rm cl}geq 0$ is in the positive $q$-direction. The classical energy is conserved: $$ begin{align}frac{1}{2}dot{q}_{rm cl}^2-V(q_{rm cl})~=~&E_{rm cl} crcr Updownarrow ~& crcr pmdot{q}_{rm cl}~=~&v(q_{rm cl})cr ~=~&sqrt{2(V(q_{rm cl})+E_{rm cl})}.end{align}tag{E}$$ The classical action contribution for the instanton bounce is $$begin{align} frac{S_{rm cl}}{2} ~stackrel{(B)}{=}~& int_0^{T/2}! dtleft(frac{1}{2}dot{q}_{rm cl}^2+V(q_{rm cl})right)cr ~stackrel{(E)}{=}~& int_{q_-}^{q_+} ! frac{dq}{v(q)}left(2V(q)+E_{rm cl}right) cr ~stackrel{(E)}{=}~& int_{q_-}^{q_+} ! dqfrac{2V(q)+E_{rm cl}}{sqrt{2(V(q)+E_{rm cl})}}cr ~stackrel{E_{rm cl}lesssim 0}{simeq}&~int_0^{q_+} ! dqsqrt{2V(q)}.end{align}tag{F}$$

5. A simple scaling argument shows that the coupling constant $g$ effectively plays the role of $hbar$ in the WKB/stationary phase approximation. The leading asymptotic WKB approximation for $gto 0$ of the Euclidean partition function is $$begin{align} {rm Re}Z(T) ~stackrel{rm WKB}{sim}&~ frac{1}{2sinhfrac{omega T}{2}}cr~approx~& e^{-omega T/2}end{align}tag{G}$$ and $$begin{align} {rm Im}Z(T) ~stackrel{rm WKB}{sim}&~crcr-frac{T}{2}sqrt{frac{1}{2pi hbar}frac{partial E_{rm cl}}{partial T}}& expleft(-frac{S_{rm cl}-TE_{rm cl}}{hbar}right)~ll ~1.end{align}tag{39.95}$$ Eq. (39.95) from Ref. 1 follows essentially via the Gelfand-Yaglom formula for functional determinants, cf. e.g. this related Phys.SE post.

6. The (imaginary part of the) ground state energy is $$ begin{align}-{rm Im}E_0 ~=~&lim_{Tto infty}frac{{rm Im} W_c(T)}{T}cr ~=~&lim_{Tto infty}frac{hbar~{rm Im} ln Z(T)}{T}cr ~approx~&lim_{Tto infty}frac{hbar~{rm Im}Z(T)}{T~{rm Re}Z(T)}.end{align}tag{H}$$ The sought-for $a_n$-coefficients can then in principle be extracted from $$begin{align} a_n ~=~&frac{1}{2pi i}oint_{mathbb{R}} ! mathrm{d}g frac{E_0(g)}{g^{n+1}}cr ~=~&frac{1}{2pi i}int_{mathbb{R}} ! mathrm{d}g frac{E_0(g-i0^+)-E_0(g+i0^+)}{g^{n+1}}cr ~=~& -frac{1}{pi}int_{mathbb{R}} ! mathrm{d}g frac{{rm Im} E_0(g)}{g^{n+1}}.end{align}tag{I}$$ This is the main answer to OP's question.

7. Let us provide some more details. We are interested in the limit $Tto infty$ where $q_-to 0$ and $$begin{align} 0 ~<~& -E_{rm cl}cr ~=~& V(q_+)cr ~=~& V(q_-) cr ~simeq~& V_2(q_-)cr ~=~& frac{1}{2}omega^2 q_-^2to 0.end{align}tag{J} $$ We calculate $$begin{align} frac{T}{2} ~=~& int_{q_-}^{q_+} ! frac{dq}{v(q)}cr ~stackrel{(E)}{=}~& int_{q_-}^{q_+} ! frac{dq}{sqrt{2(V(q)+E_{rm cl})}}cr ~stackrel{(L)+(M)}{simeq}&~ frac{I+I_2}{omega},end{align}tag{K}$$ where $$begin{align} I ~=~& int_0^{q_+}!dq~left(frac{omega}{sqrt{2V(q)}} -frac{omega}{sqrt{2V_2(q)}}right)cr ~=~& int_0^{q_+}!dq~left(frac{omega}{sqrt{2V(q)}} -frac{1}{q}right) end{align}tag{L}$$ and where $$begin{align} I_2 ~=~& omegaint_{q_-}^{q_+} ! frac{dq}{sqrt{2(V_2(q)+E_{rm cl})}}cr ~=~& int_{q_-}^{q_+} ! frac{dq}{sqrt{q^2 -q_-^2}} cr ~=~& left[ lnleft(sqrt{q^2 -q_-^2}+q right)right]_{q_-}^{q_+}cr ~simeq~&lnfrac{2q_+}{q_-}.end{align}tag{M}$$ This implies $$begin{align} q_-~stackrel{(M)}{simeq}~&2q_+e^{-I_2} crcr Downarrow ~& crcr -E_{rm cl} ~simeq~& frac{1}{2}omega^2 q_-^2cr ~stackrel{(K)}{simeq}~& 2omega^2q_+^2e^{2I-omega T} crcr Downarrow ~& crcr frac{partial E_{rm cl}}{partial T} ~simeq~&-omega E_{rm cl}.end{align}tag{N}$$
   This leads to the main formula $$begin{align} {rm Im}E_0~stackrel{(39.95)+(N)}{sim}&~crcr -frac{q_+omega^{3/2}}{2}sqrt{frac{hbar}{pi}} & e^Iexpleft(-frac{1}{hbar} S_{rm cl}right).end{align} tag{39.103}$$

8. Example: odd quartic interaction potential: $$begin{align} {cal V}(|q|)~=~& frac{1}{4}|q|^4 crcr Downarrow ~& crcr q_+~=~&omegasqrt{frac{2}{|g|}},end{align}tag{O}$$ $$begin{align} S_{rm cl}~stackrel{(F)}{=}~& 2int_0^{q_+}!dq~sqrt{2V(q)} cr ~=~& int_0^{q^2_+}!d(q^2)~sqrt{omega^2 - frac{|g|}{2}q^2} cr ~=~& frac{2omega^2}{|g|}int_0^1!du~sqrt{1-u}cr ~=~& frac{4omega^2}{3|g|},end{align}tag{P}$$ $$begin{align} I ~stackrel{(L)}{=}~& int_0^{q_+}!dq~left(frac{omega}{sqrt{2V(q)}} -frac{1}{q}right)cr ~=~ & int_0^{q_+}!frac{dq}{q}left(frac{1}{sqrt{1 -frac{|g|}{2omega^2}q^2}} -1right)cr ~=~& left[-lnleft(sqrt{1 -frac{|g|}{2omega^2}q^2} +1right)right]_0^{q_+} cr ~=~& ln 2,end{align}tag{Q}$$ $$begin{align} {rm Im}E_0~stackrel{(39.103)}{sim}&~crcr -sqrt{frac{2hbar}{pi |g|}}omega^{5/2}&expleft(-frac{4omega^2}{3|g|hbar}right).end{align}tag{R}$$

9. Example: odd cubic interaction potential: $$begin{align} {cal V}(|q|)~=~& |q|^3 crcr Downarrow ~& crcr q_+~=~&frac{omega^2}{2|g|}, end{align}tag{S}$$ $$begin{align} S_{rm cl}~stackrel{(F)}{=}~& 2int_0^{q_+}!dq~sqrt{2V(q)} cr ~=~&2 int_0^{q_+}!dq~qsqrt{omega^2 -2|g|q} cr ~=~& frac{1}{2g^2} int_0^{omega^2}!du~(omega^2 - u)sqrt{u}cr ~=~& frac{2omega^5}{15g^2}, end{align}tag{T}$$ $$begin{align} I ~stackrel{(L)}{=}~& int_0^{q_+}!dq~left(frac{omega}{sqrt{2V(q)}} -frac{1}{q}right)cr ~=~& int_0^{q_+}!frac{dq}{q}left(frac{1}{sqrt{1 -frac{2|g|}{omega^2}q}} -1right) cr ~=~& left[-2lnleft(sqrt{1 -frac{2|g|}{omega^2}q} +1right)right]_0^{q_+}cr ~=~& ln 4,end{align}tag{U}$$ $$begin{align} {rm Im}E_0~stackrel{(39.103)}{sim}&~crcr -frac{omega^{7/2}}{|g|}sqrt{frac{hbar}{pi}}&expleft(-frac{2omega^2}{15|g|^2hbar}right).end{align}tag{V}$$

References:

1. J. Zinn-Justin, QFT & Critical Phenomena, 2002; chapter 39.

2. M. Marino, Instantons and large $N$; chapter 2.

3. M. Marino, Lectures on non-perturbative effects in large $N$ gauge theories, matrix models and strings, arXiv:1206.6272; section 3.1.

4. I. Gahramanov & K. Tezgin, A resurgence analysis for cubic and quartic anharmonic potentials, arXiv:1608.08119; section 3.

--

$^1$ It is in principle possible to derive exact expressions for finite $n$ by considering multi-instantons and their interactions, cf. e.g. Ref. 4.