Michal Slowinski

Franck Thibault

Yan Tan

Jin Wang

Agnès Pernet-Liu

Shui-Ming Hu

Samir Kassi

Alain Campargue

Magdalena Konefal

Hubert Jóźwiak

An-Wen Liu

Magdalena Konefa

Konrad Patkowski

Piotr Żuchowski

Roman Ciury Lo

Daniel Lisak

Piotr Wcis Lo

H-2-He collisions: Ab initio theory meets cavity-enhanced spectra

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Molecular hydrogen differs in many respects from most other diatomic molecules, revealing its nonclassical nature due to a large rotational constant (B = 60.853 cm -1 ) [1]. The most straightforward example can be traced back to the early twentieth century, when the heat capacity of H 2 at low temperatures was discovered to dramatically diverge from the expected classical behavior [2]. The lack of inelastic channels in H 2 scattering, up to relatively large temperatures of a few hundred kelvin [3], results in many other atypical features observed in H 2 collision studies. In the particular case of the optically excited H 2 molecule, the lack of inelasticity causes the rate of optical excitation damping to be much smaller than the rate of velocity-changing collisions. This property results in strong collisional Dicke narrowing of the H 2 rovibrational lines [4,5]. Furthermore, the opening of the first inelastic channel leads to a strong dependence of the generalized spectroscopic cross-section on the collision energy [3], which enhances the speed dependence of the line-shape parameters. Pronounced collisional effects make the molecular hydrogen well suited to study not only the potential energy surfaces (PESs) and quantum-scattering calculation methodology, but also the models describing the collision-perturbed velocity distribution of the optical coherence [3,[6][7][8][9], which constitute a long-term development at the frontier of quantum optics, collision theory and statistical physics [START_REF] Rautian | Kinetic problems of non-linear spectroscopy[END_REF][START_REF] Mccourt | Nonequilibrium Phenomena in Polyatomic Gases[END_REF].

Experimental studies on the interactions between two hydrogen molecules and between a hydrogen molecule * mslowinski@umk.pl † piotr.wcislo@fizyka.umk.pl and a noble gas atom have a several-decades-long history [6,7,[START_REF] Mckellar | [END_REF][13][14][15]. Recently, the rotationally inelastic scattering of the HD molecule colliding with D 2 [16,17] and with He [18] at 1 K was observed in coexpanded supersonic beams and theoretically handled by Croft et al. [19].

The influence of the molecular hydrogen collisions, with perturbers such as H 2 isotopologues or noble gas atoms, on the shapes of the rotational and rovibrational lines in H 2 were already established in a wide range of temperatures from 20 K to 1200 K [20][21][22][23][24][25][26]. In particular, the effects of the H 2 -He collisions on the widths and shifts of H 2 rovibrational lines were subjected to intense experimental and theoretical studies [3,25,[27][28][29][30][31][32]. However, the collision-perturbed line profiles calculated from first principles never reached agreement with experimental data at the subpercent level.

In this work, we report a full ab initio description of the collision-perturbed shapes of rovibrational lines for the simplest benchmark system, i.e., Heperturbed H 2 . We demonstrate, for the first time, agreement between measured and ab initio computed collision-perturbed shapes of molecular lines at the subpercent level: the root-mean-square difference (calculated within ±FWHM) between experimental and theoretical profiles is smaller than one-hundredth of the profile amplitude. We merged all of the theoretical tools that are necessary to reach this agreement; i.e., we used the state-of-the-art statistical model of the collision-perturbed shape of molecular lines [33][34][35][36], we derived all parameters of this model from quantum scattering calculations [30,[37][38][39], and we used the highly accurate ab initio PES [30,40], referred to as BSP3. We demonstrate that the fully ab initio calculations can provide reference spectra for atmospheric studies of Solar System planets and exoplanets at the required metrological (subpercent) level of accuracy. Therefore, this approach is well-suited for populating line-by-line spec-troscopic databases such as HITRAN [41], GEISA [42], and ExoMol [43], including the recent advanced database structures [44] involving beyond-Voigt line-shape parameters and their wide ranges of temperature representations. We show that this approach can be used to experimentally test the quality of quantum-chemical calculations of the PESs [30,38,40,45,46].

II. EXPERIMENTS

To validate our ab initio quantum scattering calculations, we performed measurements of the He-perturbed 3-0 S(1) H 2 line centered at 12265.59 cm -1 (815.33 nm) and the 2-0 Q(1) H 2 line centered at 8075.31 cm -1 (1238.39 nm) using the experimental setups in Hefei [47] and Grenoble [48][49][50][51], respectively. Both setups as seen in Fig. 1, are based on continuous-wave laser sources and high-finesse ring-down (RD) cavities. The two RD cavities are 1.0 and 1.4 m long, with finesse of 63000 and 131000, which results in effective path lengths of 20 km and 59 km for Hefei and Grenoble setups. The acousto-optic modulators (AOMs) are used as optical switches to generate the RD events. The frequency of the Hefei laser source (Ti:sappphire) is set by a temperature-stabilized (∆T ≈ 10 mK at 302 K) ultralow expansion (ULE) etalon reaching 10 MHz accuracy (see Refs. [52,53]), whereas for the Grenoble setup the external cavity diode laser (ECDL) is referenced to the optical frequency comb (OFC) reaching 300 kHz absolute accuracy (see Refs. [54][55][56]). The length of the Grenoble cavity is tunable with a piezoelectric actuator, which allows the spectra to be recorded on a much denser frequency grid. The temperatures of the cavities were controlled to within 1 K and 0.15 K for the Hefei and Grenoble setups, respectively. For the 3-0 S(1) line, the H 2 /He mixing ratio spans between 33% and 10%, with 1% accuracy, while it has a constant value of 4.9(1)% for the 2-0 Q(1) line. The spectra were recorded at four pressures ranging from 0.36 to 1.35 atm with an accuracy of 0.5% and signal-to-noise ratio (SNR) up to 370 for the 3-0 S(1) line, as well as at nine pressures ranging from 0.07 to 1.05 atm with an accuracy of 0.15% and SNR up to 2800 for the 2-0 Q(1) line.

Our line-shape calculations [34,35,[START_REF] Ciury Lo | AIP Conf. Proc[END_REF][START_REF] Lisak | [END_REF], originating from first principles, are based on the generalized Hess method (GHM) [37,59,60]. To better reproduce the physics of the velocity-changing collisions, we replace the commonly used hard-collision model [61,62], which is used in GHM, with the model based on the Boltzmann operator for rigid-sphere approximation of the potential [9]. The latter model based on the isotropic part of the PES properly reproduces the dependence of the colliding partners on the mass ratio and describes the relation between velocities before and after collision [9], which is not the case for the simple hard-collision model. Moreover, we introduce the full speed dependence [63] of the collisional broadening and shift [60,64,65] into GHM.

For the H 2 -He system, the PES is three-dimensional; i.e., it depends on the distance between the center of mass of the hydrogen molecule and the helium atom, R, the intramolecular distance, r, and the angle between the inter-and intramolecular axes, θ. To carry out the close-coupling quantum scattering calculations [START_REF] Hutson | Molscat computer code, version 14[END_REF], the PES is projected on Legendre polynomials and averaged over the rovibrational wave functions of the hydrogen molecule [3,32]. The generalized spectroscopic cross-sections, σ q λ , for the electric quadrupole lines considered here (q = 2), as functions of the initial kinetic energy, E KIN , are determined from S-matrix elements [32,60,64]. For the zero rank of the velocity tensor, λ = 0, σ q λ reduces to the pressure broadening and shift cross-sections (PBXS and PSXS, respectively) [3,65,[START_REF] Ben-Reuven | [END_REF]68], while for λ = 1 it provides the complex Dicke cross-section [31,32,59,60,64] (RDXS and IDXS for its real and imaginary parts, respectively). The resulting cross-sections are shown in Fig. 2, panels (a)-(d).

We use the generalized spectroscopic cross-sections to calculate the corresponding line-shape parameters at a given temperature T, i.e., the speed-dependent pressure broadening, Γ(v), and shift, ∆(v) [63,69,70], as well as the complex Dicke parameter, ν opt [30,39,44,71,72]:

Γ(v)+i∆(v) = n 2πc 2ṽ 2 p v √ π exp -v 2 /ṽ 2 p × × ∞ 0
x 2 e -x 2 sinh (2vx/ṽ p ) σ q 0 (xṽ p ) dx, (1a)

ν opt = n 2πc m p m + m p × × ∞ 0 v r f m (v r ) 2 3 v 2 r v2 r σ q 1 (v r ) -σ q 0 (v r ) dv r , (1b)
where v, v p and v r are the speed of the active molecule, the speed of the perturber with the most probable value ṽp = 2k B T /m p , and their relative speed with mean value vr , respectively, and x = v r /ṽ p . m, m p and µ are the masses of the active molecule and perturber and their reduced mass, respectively. f m (v) is the Maxwell-Boltzmann speed distribution, n is the number density of the perturber, k B is the Boltzmann constant and c is the speed of light. The ab initio speed-dependent broadenings and shifts are shown in Fig. 2, panels (e) and (f). The calculated line-shape parameters are gathered in Table I. We denote speed-averaged broadening and shift as Γ 0 and ∆ 0 , respectively: refer to the dashed horizontal lines in Fig. 2 (e) and (f). Following Ref. [35], we introduce two parameters, Γ SD = vm We calculate the He-perturbed shapes of the H 2 lines by averaging the velocity distribution of the optical coherence [33,78,79]

I(ω) = 1 π e d 3 v f m ( v) h(ω, v). ( 2
)
f m ( v)
1 = -i(ω -ω 0 -k • v)h(ω, v) -Ŝf h(ω, v), (3)
where ω and ω 0 are the angular frequencies of the electromagnetic radiation and the unperturbed frequency of the molecular transition, respectively, k is the wave vector and Ŝf is the total collision operator describing relaxation and dephasing of the optical coherence as well as The absorption axis has been normalized to the area of the spectral line, A. The red lines under the line profiles represent the differences between experiment and theory. rRMSE values on the plots describe the root-mean-square errors relative to the profile amplitude. Panels (b) and (d): The solid red, dotted blue and dashed green lines are the zoomed residuals for the BSP3 [30,40], SK [38], and mMR [45,46] potentials, respectively. The gray vertical lines indicate the theoretical unperturbed line position [76,77]. The pressure shift is atypically large for H2 and despite small pressure range can be directly seen in this figure; see the deviation of the line maximum from the gray vertical line. its flow between different velocity classes. This operator can be written as [34,35,[START_REF] Ciury Lo | AIP Conf. Proc[END_REF][START_REF] Lisak | [END_REF]

Ŝf = -Γ(v) -i∆(v) + ν opt M f BB , (4)
where ν opt M f BB is the velocity-changing operator for the rigid-sphere model of collisions [80,81]. By using the complex Dicke parameter derived from GHM, ν opt , we introduce the correlations between internal and translational degrees of freedom to the velocity-changing operator. The ultimate collision operator used to compare ab initio calculations with experiment is a sum of two contributions [9]: the H 2 -He operator discussed here, Eq. ( 4), and the H 2 -H 2 operator taken from Ref. [82].

IV. EXPERIMENTAL VALIDATION OF THE QUANTUM SCATTERING CALCULATIONS AND TEST OF THE AB INITIO POTENTIALS

The direct comparison between theory and experiment is shown in Fig. 3, panels (a) and (c). Despite the fact that none of the line-shape parameters were fitted, we achieved subpercent agreement between the simulated and raw experimental profiles; i.e., the root-mean-square error of the ab initio model relative to profile amplitude (rRMSE) calculated within ±FWHM is smaller than 1%. The pressure-average values of rRMSE are 0.99%and 0.33% for the 3-0 S(1) and 2-0 Q(1) lines, respec-tively. We do not know the reason why the discrepancy is larger for the 3-0 S(1) line. This is the first time that the fully ab initio line-shape calculations, includ-ing not only the collisional perturbations of both inter-nal and translational motions of molecules but also the correlations between them, are validated for highly accurate experimental spectra. In addition to the line-shape parameters, whose values are directly taken from our ab initio calculations, several other parameters (not related to the collisional effects) were fitted: the unperturbed line position frequency, ν 0 = ω 0 /2π, the line intensity, S, and the background slope and baseline. The fitted values of ν 0 and S are consistent with the theoretical values [51,76,77] within the declared 1σ uncertainties.

The advantage of the method employed in this work is that it allows the absolute value of the collisional cross-sections to be accurately tested experimentally, which is not the case for most of the other experimental techniques [16][17][18]. As an example, we consider the pressure broadening cross-section (PBXS in Fig. 2 panel (a)), which gives the Γ 0 parameter after thermal averaging: refer to the dashed horizontal lines in Fig. 2 (e). At high pressures (at which the collisional broadening dominates), the relative error of Γ 0 is approximately equal to the residual rRMSE. For instance, in the case of the 2-0 Q(1) line, rRMSE = 0.18% at the highest pressure corresponds to subpercent agreement between our ab initio and experimental values of Γ 0 .

In addition to demonstrating that the full theoretical methodology can reproduce the collision-perturbed experimental profiles at an unprecedented level of accuracy, we show that this approach can serve as a link between high-precision spectroscopy and quantum chemical calculations, thus enabling experimental tests of PESs in regions that are inaccessible with other techniques. In Figure 3, panels (b) and (d), we show the differences between the experimental and theoretical profiles for three different PESs: the Schaefer and Köhler (SK) PES [38], modified Muchnik and Russek (mMR) PES [45,46], and the PES used in our main analysis, namely the extended Bakr, Smith and Patkowski (BSP3) [30,40]. Contrary to the main analysis here we do not fit the line inten-sity but fix it to the theoretical value [51]. We do so because for the older less accurate PESs [38,45,46] the fitted line intensities, to some extent, artificially com-pensate inaccurate values of some of the line-shape pa-rameters. This comparison clearly demonstrates that the best agreement with the experiment is achieved for the BSP3 PES. This result is consistent with our expecta-tions, since the most recent BSP3 surface was calculated using the coupled-cluster method with single, double and perturbative triple excitations (CCSD(T)), augmented by full configuration interaction corrections for which the relative uncertainty in the minimum is estimated as 0.4% (± 0.045 cm -1 ). Two other potentials studied here, SK and mMR, were obtained using much smaller Gaussian basis sets and lower levels of electronic structure theory. Moreover, the BSP3 PES involves a consistent treatment of a large range of intramolecular H-H distances and ensures a proper asymptotic behavior for large H 2 -He separations [30].

V. CONCLUSION

We demonstrated for the first time that the calculations from first principles are able to reproduce the experimental collision-perturbed molecular spectra at subpercent level accuracy. We used quantum scattering calculations to construct a collision operator that describes not only the dephasing and damping of the optical coherence but also flows of an optical coherence between different velocity classes. We solved a Boltzmann-like transport-relaxation equation to determine the velocity distribution of optical coherence and, hence, the collisionperturbed shapes of molecular lines. The calculations were validated on the highly accurate experimental spectra of the H 2 S(1) 3-0 and Q(1) 2-0 electric quadrupole transitions perturbed by helium. In addition to validating the theoretical methodology of accurate determination of the collision-induced shapes of molecular lines, we demonstrated that this approach can be used to test the quantum-chemical calculations of the potential energy surfaces. This work demonstrates that theoretical calculations are capable of delivering line-shape parameters with the accuracy necessary to interpret experimental spectra for atmospheric, planetary, and astrophysical applications. With this approach, we show the possibility of calculating and validating beyond-Voigt line-shape parameters in order to populate line-by-line spectroscopic databases such as HITRAN [41], GEISA [42], or Exo-Mol [43] with requested accuracies in a wide temperature range for recent advanced database structures [44]. Such approach is important for the lines relevant for atmospheric and astrophysical applications, especially for those which have not been accessed yet with laboratory experiments.

FIG. 1 .

 1 FIG. 1. Cavity ring-down spectrometers in Hefei (a) and Grenoble (b) laboratories. AOM -acousto-optic modulator, ECDL -external cavity diode laser, EOM -electro-optic modulator, FC -fiber coupler, OFC -optical frequency comb, OI -optical isolator, PD -photodiode, PFC -polarization fiber coupler, PZT -piezoelectric transducer, Ti:Sapph -Ti:sapphire laser.

FIG. 2 .

 2 FIG. 2. (a-d) Generalized spectroscopic cross-sections for the 3-0 S(1) and 2-0 Q(1) lines in helium-perturbed H2: see the respective black and red lines, respectively. PBXS, PSXS, RDXS and IDXS are pressure broadening, pressure shift, and real and imaginary parts of complex Dicke cross-sections, respectively. The two curves in panel (c) are overlapping. (e, f) Ab initio speed dependence of the line broadening, Γ, and shift, ∆ (the same color notation), plotted for 1 amagat of He at the experimental temperatures. The black and red dashed lines are the speed-averaged collisional broadening, Γ0, and shift, ∆0, respectively. The gray line is the Maxwell Boltzmann distribution at 296 K in arbitrary units; the distribution of relative kinetic energy in panels (a)(d), and the distribution of active molecule speed in panels (e) and (f).

2 d

 2 dv Γ(v) v=ṽ and ∆ SD = vm 2 d dv ∆(v) v=ṽ , that quantify the strength of the speed dependence of the broadening and shift, respectively, where v m is the most probable active molecule speed. The Γ SD /Γ 0 and ∆ SD /∆ 0 ratios, larger than 0.4, reveal extraordinarily pronounced speed dependence for the H 2 case; for typical atmospheric systems, these ratios are at the level of 0.1[73][74][75]. TABLEI. Calculated values of the line-shape parameters for the 3-0 S(1) and 2-0 Q(1) H2 lines. The calculations are performed for T = 296.6 K for the 3-0 S(1) line and for T = 294.2 K for the 2-0 Q(1) line. All parameters are given in units of 10 -3 cm -1 for n = 1 amagat.

  is the Maxwell-Boltzmann distribution of the active molecule velocity, v, and f m ( v) h(ω, v) is a scalar function proportional to the velocity distribution of the optical coherence. The h(ω, v) function fulfills the transport-relaxation equation[33,79]

2 )FIG. 3 .

 23 FIG.3. Panels (a) and (c): Comparison of the raw experimental spectra (black points) and the simulated line profiles (red lines). The absorption axis has been normalized to the area of the spectral line, A. The red lines under the line profiles represent the differences between experiment and theory. rRMSE values on the plots describe the root-mean-square errors relative to the profile amplitude. Panels (b) and (d): The solid red, dotted blue and dashed green lines are the zoomed residuals for the BSP3[30,40], SK[38], and mMR[45,46] potentials, respectively. The gray vertical lines indicate the theoretical unperturbed line position[76,77]. The pressure shift is atypically large for H2 and despite small pressure range can be directly seen in this figure; see the deviation of the line maximum from the gray vertical line.

ACKNOWLEDGMENTS M.S. contribution is supported by the Polish National Science Centre Project No. 2017/26/D/ST2/00371 and by the A next-generation worldwide quantum sensor network with optical atomic clocks project carried out within the TEAM IV Programme of the Foundation for Polish Science cofinanced by the European Union under the European Regional Development Fund. P.W. and H.J. contributions are supported by the Polish National Science Centre Project Nos. 2015/19/D/ST2/02195 and 2018/31/B/ST2/00720. S.-M.H. acknowledges the support from National Natural Science Foundation of China (21688102). M.K. and D.L. contributions are supported by the Polish National Science Centre Projekt No. 2015/18/E/ST2/00585. K.P. is supported by the U.S. National Science Foundation CA-REER award CHE-1351978. P. Ż. contribution is supported by the Polish National Science Centre Project No. 2015/19/B/ST4/02707. The project is cofinanced by the Polish National Agency for Academic Exchange under the PHC Polonium program (dec.

PPN/X/PS/318/2018). The research is part of the program of the National Laboratory FAMO in Toruń, Poland. It is also supported by the COST Action CM1405 MOLIM.


Jean Zinn-Justin

Ulrich D Jentschura

Imaginary Cubic Perturbation: Numerical and Analytic Study

Keywords: numbers: 11.10.Jj, 11.15.Bt, 11.25.Db, 12.38.Cy, 03.65.Db

The analytic properties of the ground state resonance energy E(g) of the cubic potential are investigated as a function of the complex coupling parameter g. We explicitly show that it is possible to analytically continue E(g) by means of a resummed strong coupling expansion, to the second sheet of the Riemann surface, and we observe a merging of resonance and antiresonance eigenvalues at a critical point along the line arg(g) = 5π/4. In addition, we investigate the convergence of the resummed weak-coupling expansion in the strong coupling regime, by means of various modifications of order-dependent mappings (ODM), that take special properties of the cubic potential into account. The various ODM are adapted to different regimes of the coupling constant. We also determine a large number of terms of the strong coupling expansion by resumming the weak-coupling expansion using the ODM, demonstrating the interpolation between the two regimes made possible by this summation method.

I. INTRODUCTION

The Hamiltonian

H = - 1 2
∂ 2 ∂x 2 + 1 2 x 2 + i √ g 6 x 3 , (1)
with g real positive, has been a subject of a number of studies since it has been conjectured by Bessis and Zinn-Justin (1992) to have a real spectrum. Arguably, Eq. (1) provides the simplest example of a PT symmetric Hamiltonian, that is, of a Hamiltonian invariant under a simultaneous complex conjugation (i.e., time reversal) and parity transformation x → -x (see Ref. [START_REF] Bender | Real spectra in non-Hermitian Hamiltonian having PT symmetry[END_REF]). In particular, the reality of its spectrum was proven first in the strong coupling limit [START_REF] Dorey | Spectral equivalences, Bethe ansatz equations, and reality properties in PT -symmetric quantum mechanics[END_REF], then more generally in Ref. [START_REF] Shin | On the reality of eigenvalues for a class of PT -symmetric oscillators[END_REF] (together with the positivity of the eigenvalues).

Here, we continue our investigations [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF] regarding its numerical and analytic properties, including the analytic structure of its resonance energy eigenvalues as a function of the complex coupling parameter g, specially in the strong coupling regime and near the first level crossing singularity.

We denote by E(g) an energy eigenvalue of the Hamiltonian [START_REF] Bender | Real spectra in non-Hermitian Hamiltonian having PT symmetry[END_REF]. One verifies that its perturbative expansion contains only integer powers of g:

E(g) ∼ g→0 L=0 E L g L , (2)
with real coefficients E L . The series is divergent for all g, but is Borel summable [START_REF] Caliceti | Perturbation theory of odd anharmonic oscillators[END_REF][START_REF] Caliceti | Distributional Borel summability of odd anharmonic oscillators[END_REF], and a steepest descent calculation of the path integral representation of the corresponding quantum partition function [START_REF] Lipatov | Divergence of the perturbation-theory series and the quasi-classical theory[END_REF][START_REF] Brézin | Perturbation theory at large order. I. The ϕ N interaction[END_REF][START_REF] Zinn-Justin | Perturbation series at large orders in quantum mechanics and field theories: application to the problem of resummation[END_REF][START_REF] Jentschura | Generalized nonanalytic expansions, PT -symmetry and large order formulas for odd anharmonic oscillators[END_REF] yields a large order behaviour for the nth energy eigenvalue of the form

E L ∼ L→∞ (-1) L+1 6 π 3/2 288 n n! Γ L + n + 1 2 A L+n+1/2 , (3)
* Electronic address: jean.zinn-justin@cea.fr where A = 24 5 is the instanton action. More recently, Padé summability was also rigorously established in Ref. [START_REF] Grecchi | Padé summability of the cubic oscillator[END_REF]. The imaginary part of the energy levels on the cut for g = -|g| + i 0 is positive. This proof confirms numerical investigations based on the summation of the perturbative expansion by Padé approximants [START_REF] Bender | Numerical evidence that the perturbation expansion for a non-Hermitian PTsymmetric Hamiltonian is Stieltjes[END_REF].

More information on analytic and numerical properties was provided by an ODM summation [START_REF] Seznec | Summation of divergent series by order dependent mappings: Application to the anharmonic oscillator and critical exponents in field theory[END_REF][START_REF] Zinn-Justin | Summation of divergent series: Orderdependent mapping[END_REF] of the perturbative series for the ground state energy (see Ref. [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF]). It was based on the mapping g → λ,

g = ρ λ (1 -λ) 5/2 , ( 4
)
where ρ is an order-dependent parameter.

In particular, the leading term of the strong coupling expansion could be determined (and was found in agreement with results coming direct solutions of Schrödinger equation [START_REF] Fernandez | Strong-coupling expansions for the PT -symmetric oscillators V (x) = a(ix) + b(ix) 2 + c(ix) 3[END_REF][START_REF] Jentschura | Calculation of the Characteristic Functions of Anharmonic Oscillators[END_REF]) and the positivity of the imaginary part verified numerically for all g < 0. It was shown that along the negative real axis, the imaginary part above the cut is a simple positive decreasing function interpolating smoothly between a strong coupling power-law behaviour and a non-analytic exponential tunneling factor obtained from semi-classical instanton calculation for weak coupling. Finally, it was argued that the ground state energy is analytic up to arg(g) = 5π/4 in the second Riemann sheet, a property that we discuss in this article in more detail.

By a simple coordinate translation, one can relate the Hamiltonian (1) to another PT symmetric Hamiltonian,

H [qc] = - 1 2
∂ 2 ∂x 2 + i 1 6 x 3 + 1 2 χ x , (5)
where χ = g -4/5 . This correspondence suggests to sum the series by another ODM of the form

g = ρ λ (1 -λ) 5/4 , (6)
which is shown here to be much more efficient in the strong coupling regime, while the mapping (4) (and a variation thereof discussed below) is more efficient for weak coupling. The method [START_REF] Caliceti | Distributional Borel summability of odd anharmonic oscillators[END_REF] is used here in order to determine a number of coefficients of the large g expansion (small χ expansion). Finally, studying the limit of weak χ, we also construct a corresponding continued fraction (a subclass of Padé approximants), which enables us to obtain more insight into the analytic structure. Combining all methods, we have a precise control on the analytic continuation of the ground state energy in the first Riemann sheet of the uniformization variable χ = g -4/5 , except near the cut in the χ-plane. Based on our observation, we conjecture that the ground state energy of the Hamiltonian ( 5) is a real analytic function, with a cut on the real negative χ axis, which corresponds to a complex phase of arg(g) = 5π/4 for the original coupling parameter g and thus lies on the second sheet of the Riemann surface for the original coupling parameter g. This domain is not covered by Ref. [START_REF] Shin | On the reality of eigenvalues for a class of PT -symmetric oscillators[END_REF]. The imaginary part of the ground state energy on the cut for χ = -|χ|+ i 0 is positive. The small χ expansion of the resonance energies converges in a disk. From the point of view of the initial Hamiltonian [START_REF] Bender | Real spectra in non-Hermitian Hamiltonian having PT symmetry[END_REF], the singularity closest to the origin has the interpretation of a level crossing.

This can be seen as follows: the cut for negative χ corresponds to a cut of the energy level on the second sheet [START_REF] Bender | Anharmonic oscillator[END_REF][START_REF] Delabaere | Spectral analysis of the complex cubic oscillator[END_REF][START_REF] Kleinert | Convergence behavior of variational perturbation expansion: A method for locating Bender-Wu singularities[END_REF] of the Riemann surface ("Bender-Wu cut"). It starts from the limit of the circle of convergence of the strong coupling expansion, i.e., from a point on the boundary of the circle of convergence of the χ expansion. Two resonance energies are attached to points infinitesimally displaced above and below the cut. The functions describing the real parts of the two resonance energies are equal on the cut. They cross as we cross the cut. By contrast, the functions describing their imaginary parts are equal and opposite on the cut and vary smoothly in its immediate vicinity.

At the start of the cut, the imaginary part of the resonance and antiresonance energy is zero, but the real parts still cross, and this configuration therefore corresponds to a level crossing, given by the confluence of two energy levels on the second sheet of the Riemann surface characterizing the resonance energies as a function of g. For the Hamiltonian (5), it corresponds, in addition, to a spontaneous breaking of the PT symmetry (namely, the spectrum composed of the resonance and antiresonance energies with their equal and opposite imaginary parts corresponds to a manifestly non-PT symmetric time evolution). Starting from the region where the two eigenvalues merge and are real, we can say that after the singularity, they become complex conjugate.

To confirm this picture, that is, analyticity in a cutplane in the χ variable, sign of the imaginary part of the ground state energy on the cut and level merging singularity we investigate together the ground state and the first excited state, forming symmetric combinations: sum and difference squared of the corresponding energies. We verify that these functions are regular at the level merging singularity, showing that the level merging involves the ground state and the first excited state. This gives us a totally consistent picture and, moreover, allows us to determine with increased precision a number of terms of the small χ expansion. In turn, from the small χ expansion a continued fraction expansion is derived, whose coefficients are all positive and seem to converge to a value consistent with a square root singularity, within errors.

As a final remark, we notice that the Hamiltonian (1) can be considered as a toy model for the i φ 3 quantum field theory in two and three space dimensions, whose renormalization group functions describe the Lee-Yang edge singularity of Ising-like statistical models [START_REF] Fisher | Yang-Lee Edge Singularity and φ 3 Field Theory[END_REF][START_REF] Bender | Scalar Quantum Field Theory with Cubic Interaction[END_REF].

This paper is organized as follows. We first discuss the scaling of the coupling constant and the related strong coupling Hamiltonians in Sec. II. General properties of ODM summation methods are recalled in Sec. III. Optimized methods for the problem at hand are introduced in Sec. IV. Numerical results for g positive and negative are obtained and discussed in Sec. V. Sec. VI is devoted to the strong coupling expansion and level merging. Analytic properties of the ground state resonance in the second Riemann sheet and continued fraction expansion are discussed in Sec. VII. Conclusions are relegated to Sec. VIII.

II. STRONG COUPLING HAMILTONIANS

The time-independent Schrödinger equation

- 1 2
∂ 2 ∂x 2 + 1 2 x 2 + i 6 √ g x 3 ψ(x) = E(g) ψ(x) (7)
corresponds to the Hamiltonian [START_REF] Bender | Real spectra in non-Hermitian Hamiltonian having PT symmetry[END_REF]. With proper boundary conditions at infinity, this equation determines the eigenvalues E(g) of the operator. After rescaling the variable x as x → g -1/5 x, the eigenvalue equation can be rewritten as follows (the superscript [sc] stands for the strong-coupling limit)

H [sc] (ξ) ψ(x) = E [sc] (ξ) ψ(x) , (8)
where

ξ = g -2/5 , (9)
and the strong coupling Hamiltonian is

H [sc] (ξ) = - 1 2 ∂ 2 ∂x 2 + i 6 x 3 + 1 2 ξ x 2 . ( 10
)
The energies of the original and the strong coupling Hamiltonian are related by

E(g) = g 1/5 E [sc] (g -2/5 ) = g 1/5 E [sc] (ξ) . (11)
In particular, we have in leading order

E(g) ∼ g→∞ g 1/5 E [sc] (0) , (12)
and the strong coupling expansion

E [sc] (ξ) = ∞ L=0 E [sc] L ξ L . (13)
In the uniformizing variable √ ξ = g -1/5 , the singularities of E(g) in the complex plane correspond to level crossings, whereas the cut of the original resonance energy, defined as a function of g, is moved to ξ = ∞.

We then shift the coordinate x → x+i ξ. The equation becomes

- 1 2
∂ 2 ∂x 2 + i 6 x 3 + i 2 ξ 2 x - 1 3 ξ 3 ψ(x) = E [sc] (ξ) ψ(x). (14)
We conclude that

E [sc] (ξ) = - 1 3 ξ 3 + E [qc] (ξ 2 ) , (15)
where E [qc] (χ ≡ ξ 2 ) is a function of the linear coupling

χ = ξ 2 = g -4/5 , (16)
and thus an eigenvalue of the PT symmetric Hamiltonian (5) corresponding to the Schrödinger equation

H [qc] (χ) ψ(x) = E [qc] (χ) ψ(x) , (17)
where

H [qc] (χ) = - 1 2
∂ 2 ∂x 2 + i 1 6 x 3 + 1 2 χ x . (18)
A generalization of the considerations leading to Eq. ( 14) to the normalization of the cubic potential used in Ref. [START_REF] Jentschura | Structure, Time Propagation and Dissipative Terms for Resonances[END_REF] is given in Appendix C. This generalization is somewhat non-trivial because the precise form of the strong coupling Hamiltonian depends on the normalization of the initial perturbation. When translated to the conventions of Ref. [START_REF] Jentschura | Structure, Time Propagation and Dissipative Terms for Resonances[END_REF], the coefficient of relative order ξ 3 becomes 1/108 for all states of the spectrum; this affords an explanation for the observation made originally in Table 2 of Ref. [START_REF] Jentschura | Structure, Time Propagation and Dissipative Terms for Resonances[END_REF]. The expansion of the ground state resonance energy E [qc] (χ) in powers of χ reads

E [qc] (χ) = ∞ L=0 E [qc] L χ L . ( 19
)
Unlike the perturbative expansion in g, it is convergent in a disk. Therefore, E(g) has a convergent large g expansion of the form

E(g) = g 1/5 - 1 3 ξ 3 + E [qc] (χ) = - 1 3g + g 1/5 ∞ L=0 E [qc] L g -4L/5 . (20)
These simple transformations show that E [sc] (ξ) has only one term odd in ξ, which moreover can be read off from Eq. [START_REF] Fisher | Yang-Lee Edge Singularity and φ 3 Field Theory[END_REF]. As a consequence, for g large, the relevant expansion variable is g -4/5 rather than g -2/5 , an observation that influences the construction of optimized ODM summation methods for our problem. In the uniformizing variable g -1/5 , it is sufficient to determine the eigenvalues for

- π 4 ≤ arg g -1/5 ≤ π 4 , (21)
in order to know them in the whole first Riemann sheet of the χ variable (-π ≤ arg(χ) ≤ π). Conversely, knowledge of the resonance eigenvalues for -π ≤ arg(χ) ≤ π implies a determination for -5π/4 ≤ arg(g) ≤ 5π/4 which exceeds the first Riemann sheet. Indeed, the Hamiltonian ( 17) is PT symmetric for all real values of the "quadratic coupling" χ but corresponds to the Hamiltonian (1) only in the case χ ≥ 0. It is thus interesting to study its spectrum also in the region χ < 0, which corresponds to the analytic continuation g = |g| exp(5 i π/4). Equation ( 17) is also related to residues of poles of solutions of a linear system associated to the first Painlevé equation [START_REF] Masoero | Poles of Intégrale Tritronquéée and Anharmonic Oscillators[END_REF].

III. GENERAL SUMMATION METHODS

A. Idea of the method

We here briefly recall the ODM summation method. Details can be found elsewhere [START_REF] Seznec | Summation of divergent series by order dependent mappings: Application to the anharmonic oscillator and critical exponents in field theory[END_REF][START_REF] Zinn-Justin | Summation of divergent series: Orderdependent mapping[END_REF]. The method is based on some a priori knowledge, or educated guess, of the analytic properties of the function that is expanded. It applies both to convergent and divergent series, although it is mainly useful in the latter situation.

In what follows, we consider a function analytic in a sector and mappings g → λ of the form

g = ρ ζ(λ) , (22)
where ζ(λ) is a real analytic function increasing on 0

≤ λ < 1, such that ζ(λ) = λ + O(λ 2 ) and, for λ → 1, ζ(λ) ∝ (1 -λ) -α with α > 1. A possible choice is g = ρ λ (1 -λ) α . ( 23
)
The parameter α has to be chosen in accordance with the analytic properties of the function E (a rational number in our examples) and ρ is an adjustable parameter, whose interpretation is an order-dependent, "artificial, local radius of convergence" of the truncated, divergent perturbative expansion (which becomes smaller as the order of the transformation is increased). An important property that singles out relevant mappings here is the following: one chooses mappings such that, for g → ∞ and thus λ → 1, the quantity g -1/α has a regular expansion in powers of 1 -λ. A more general discussion of the method can be found in Ref. [START_REF] Seznec | Summation of divergent series by order dependent mappings: Application to the anharmonic oscillator and critical exponents in field theory[END_REF][START_REF] Zinn-Justin | Summation of divergent series: Orderdependent mapping[END_REF].

After the mapping ( 22), E is given by a Taylor series in λ of the form

E g(λ) = ∞ L=0 P L (ρ) λ L , (24)
where the coefficients P L (ρ) are polynomials of degree L in ρ. Since the result is formally independent of the parameter ρ, the parameter can be chosen freely. At ρ fixed, the series in λ is still divergent, but it has been verified on a number of examples (all Borel summable), and proven in certain cases [START_REF] Guida | Improved Convergence Proof of the Delta Expansion and Order Dependent Mappings[END_REF] that, by adjusting ρ order by order, one can devise a convergent algorithm. The Kth approximant E (K) (g) is constructed in the following way: one truncates the expansion at order K and chooses ρ as to cancel the last term. Since P K (ρ) has K roots (real or complex), one chooses, in general, for ρ the largest possible root (in modulus) ρ = ρ K with P K (ρ K ) = 0, for which P ′ K (ρ) is small. This leads to a sequence of approximants E(g) ≈ E (K) (g), where

E (K) (g) = K L=0 P L (ρ K ) λ L (g, ρ K ) ; P K (ρ K ) = 0 . (25)
The term P K+1 (ρ K ) then gives an order of magnitude of the error.

In the case of convergent series, it is expected that ρ K has a non-vanishing limit for K → ∞. By contrast, for divergent series it is expected that ρ K goes to zero for large K as

ρ K = O E -1/K K . ( 26
)
The intuitive idea here is that ρ K corresponds to a 'local' radius of convergence. Let us conclude this section with a few remarks. Above, we have used the selection criterion P K (ρ K ) = 0 for the order-dependent ρ K which are used to calculate the polynomials P L (ρ K ), for L = 0, . . . , K -1. Alternatively, one can choose the largest roots ρ K of the polynomials P ′ K (ρ K ) = 0 for which P K is small. Other mixed criteria involving a combination of P K and P ′ K can also be used. Indeed, one empirically observes that the approximant is not very sensitive to the precise value of ρ K , within errors. In the ODM method, the determination of the sequence of the ρ K 's is the most time-consuming computational step. Indeed, once the ρ K are known, for each value of g, the calculation reduces to inverting the mapping [START_REF] Jentschura | Structure, Time Propagation and Dissipative Terms for Resonances[END_REF] and simply summing the Taylor series in λ to the relevant order.

B. Convergence analysis

We here give a heuristic analysis [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF][START_REF] Seznec | Summation of divergent series by order dependent mappings: Application to the anharmonic oscillator and critical exponents in field theory[END_REF] of the convergence of the ODM method that shows how the convergence can be optimized. This will justify the choice of the class of zeros of the polynomials P K or P ′ K and provide a quantitative analysis of the corresponding convergence.

In order to focus the analysis, we assume that we are dealing with functions that have the properties of the eigenvalues of the Hamiltonian (1). We consider only real functions analytic in a cut-plane with a cut along the real negative axis and a Cauchy representation of the form

E(g) = E(0) + g π 0- -∞ dg ′ ∆(g ′ ) g ′ (g ′ -g) , (27)
where the subtraction ensures the convergence of the integral for g → -∞.

For the example we discuss here, one can derive by a steepest descent calculation [START_REF] Lipatov | Divergence of the perturbation-theory series and the quasi-classical theory[END_REF][START_REF] Brézin | Perturbation theory at large order. I. The ϕ N interaction[END_REF][START_REF] Zinn-Justin | Perturbation series at large orders in quantum mechanics and field theories: application to the problem of resummation[END_REF] an asymptotic behaviour of the form

∆(g) ∝ g→0- 1 (-g) b+1 exp A g , A > 0 . ( 28
)
The function E(g) can be expanded in powers of g,

E(g) = ∞ L=0 E L g L , E L = 1 π 0- -∞ dg g L+1 ∆(g) , (29) the Cauchy representation of E L being valid for L > 0.
The asymptotic behaviour (28) then implies a large order behaviour of the form (3),

E L ∝ L→∞ (-A) -L Γ(L + b + 1) ∝ k→∞ (-A) -L L b L! . (30)
A remark is in order. When g becomes negative in Eq. ( 7), the expression i √ g becomes real, and the corresponding resonance and antiresonance eigenvalues are attached to the upper and lower side of the integration contour which extends from g ′ = -∞ to g ′ = 0. The discontinuity ∆(g ′ ) is the difference of the imaginary parts of the resonance and antiresonance eigenvalues. For the potential we study here, it has been proved that ∆(g ′ ) is positive [START_REF] Grecchi | Padé summability of the cubic oscillator[END_REF], and this is consistent with Eq. ( 28). In expression (27), for g > 0 the integrand is positive throughout the integration contour (the denominator contains a product of two negative quantities). The eigenvalues E(g) of the PT symmetric Hamiltonian (1) thus are positive since E(0) = 1/2 is positive. Similarly, the expression (29) shows that the sign of E L is (-1) L+1 , a result consistent with the estimate given in Eq. (30).

We now introduce the order-dependent mapping [START_REF] Jentschura | Structure, Time Propagation and Dissipative Terms for Resonances[END_REF],

g = ρ ζ(λ) , (31)
with the understanding that the point λ = 1 corresponds to the point g → ∞. Because we can alternatively formulate the dispersion relation in terms of λ, the Cauchy representation then can be written as where Γ initially is the image of the cut on the real negative axis (see Fig. 1 and2 for examples). The contour Γ can then be deformed if the function has analyticity properties beyond the first Riemann sheet.

E g(λ) = 1 π Γ dλ ′ ∆ g(λ ′ ) λ ′ -λ , (32)
We expand

E g(λ) = ∞ K=0 P K (ρ) [λ(g)] K (33)
with

P K (ρ) = 1 π Γ dλ λ K+1 ∆ g(λ) . ( 34
)
For K → ∞, the factor λ -K favours small values of |λ|, but for too small values of λ, the exponential decay of ∆(g(λ)) takes over. Thus, P K (ρ K ) can be evaluated by the steepest descent method. The ansatz

ρ K ∼ K→∞ R K , R > 0 , (35)
implies that in the limit K → ∞, the saddle point values of λ are independent of K and that g(λ) → 0. The former statement can easily been seen by taking the logarithm of the integrand in Eq. ( 34) and observing that K factors out in the equation defining the saddle point. Thus, ∆(g) can be replaced by its asymptotic form (28) for g → 0 -, except for g(λ) of order one (or larger) and thus λ close to the initial endpoint λ = 1 for K → ∞.

The contribution of the endpoint depends on the analytic properties of ∆(g) but is bounded by a quantity of order exp(C

1 K 1-1/α ) [see Eq. (46)].
In what follows we set

R A = µ , ( 36
)
since µ is the only parameter (and it is independent of the normalization of g). The behaviour of P K (ρ K ) is then given by the sum of leading saddle contributions. Each saddle point contribution is of order

{P k (ρ K )} saddle point = O(e Kσ ) , σ ≡ 1 µ ζ(λ)
-ln(λ) ,

(37) where λ is determined by the saddle point equation

σ ′ (λ) = - 1 µ ζ ′ (λ) ζ 2 (λ) - 1 λ = 0 . ( 38
)
The analysis is simple only if λ at the saddle point is real negative. By contrast, if λ is complex, the analytic properties of ∆(g) and the possibility of deforming the contour Γ become relevant. For µ large, the equation has a unique solution (λ ∼ -1/µ), which is real negative. Then we note that at this saddle point, as a function of µ,

∂σ ∂µ = - 1 µ 2 ζ(λ) > 0 , (39)
as long as the saddle point value of λ is negative. This suggests decreasing µ as much as possible to improve the convergence. The exponential rate σ corresponding to the saddle point vanishes for

0 = Re σ(λ) = 1 µ ζ(λ) -ln(-λ) , (40)
and this defines a special value µ = µ c of the parameter µ.

In [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF][START_REF] Seznec | Summation of divergent series by order dependent mappings: Application to the anharmonic oscillator and critical exponents in field theory[END_REF][START_REF] Zinn-Justin | Summation of divergent series: Orderdependent mapping[END_REF], we have discussed in some detail mappings of the form

g = ρ λ (1 -λ) α , (41)
for various values of α. Most notably, α = 5/2 has been applied to the present problem (for a motivation regarding this mapping, see Appendix A).

The system of Eqs. ( 38) and ( 40) can be solved for the variables µ = µ c and λ = λ c . We obtain mapping (a):

g = ρ λ (1 -λ) 5/4 , (42a)
⇒ µ (a) c = 3.811 522 . . . , λ (a) c = -0.259 901 . . . ,

mapping (b): g = ρ λ (1 -λ/2) (1 -λ) 5/2 , (42b) ⇒ µ (b) c = 4.445 762 . . . , λ (b) c = -0.216 262 . . . , mapping (c): g = ρ λ (1 -λ) 5/2 , (42c)
⇒ µ (c) c = 4.895 690 . . . , λ (c) c = -0.189 645 . . . . Then, we expect the error ε K at order K to be of the same order as the leading saddle point contribution, that is, of order P K (ρ K ) λ(g, ρ K ) K or in view of Eq. (37),

ε K = O(e Kσ λ(g, ρ K ) K ) . ( 43
)
For K → ∞, at g fixed, λ goes to 1 as

1 -λ(g, ρ K ) ∼ K→∞ C 2 (Kg) -1/α , (44)
where for µ = µ c and methods (42a), (42b) (42c), respectively

C (a) 2 = (A µ (a) c ) 4/5 ≈ 10.23 , (45a)
C (b) 2 = A µ (b) c 2 2/5 ≈ 3.40 , (45b)
C (c) 2 = A µ (c) c 2 2/5 ≈ 3.53 . (45c)
Then,

λ K (g, ρ K ) ∼ K→∞ exp -C 2 K 1-1/α g -1/α . (46)
Then, two cases are possible: Case (i): if the contribution to the integral corresponding to g(λ) of order unity decreases exponentially with K (which implies the possibility of deforming the contour Γ), we can also choose µ < µ c and the ODM method converges for all g = 0. Moreover, this in turn implies that the function E(g) is an entire function of g -1/α . This clearly is a non-generic situation, but several examples have been met in the form of simple integrals.

Case (ii): in a generic situation where the mapping removes the singularity at infinity but other singularities are present, as explained in Refs. [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF][START_REF] Seznec | Summation of divergent series by order dependent mappings: Application to the anharmonic oscillator and critical exponents in field theory[END_REF][START_REF] Zinn-Justin | Summation of divergent series: Orderdependent mapping[END_REF] we then expect the optimal mapping to correspond to µ = µ c and P K (ρ K ) (or P K+1 (ρ K ) if we choose for ρ K a zero of P K ) to be, asymptotically for K → ∞, itself of order

P K (ρ K ) = O exp C 3 K 1-1/α ( 47
)
where α = 5/4 for method (42a) and α = 5/2 for methods (42b) and (42c), respectively. The domain of convergence then depends on the sign of the constant C 3 . For C 3 > 0, the domain of convergence is

|g| < C 2 C 3 α cos arg(g) α α . (48)
For α = 5/2, this domain extends beyond the first Riemann sheet and requires analyticity of the function E(g) in the corresponding domain. For C 3 < 0, the domain of convergence is the union of the sector | arg(g)| < πα/2 and the domain

|g| > C 2 C 3 α -cos arg(g) α α . ( 49
)
Again for α = 5/2, this domain extends beyond the first Riemann sheet. We will see that this is the situation realized in the three methods (42a), (42b) and (42c).

One might now wonder how this analysis is related to the zeros of P K or P ′ K . Choosing the relevant zeros is a natural way of finding the region of minimal values. In fact, the relevant zeros lie in the region where two leading saddle points cancel or the leading saddle point(s) cancel the contribution coming from the remaining part of the Cauchy integral.

Indeed, in the discussed examples of case (i), one could show that the zeros correspond to a cancellation between saddle points. In the case (ii), empirically one observes that the zeros are consistent with the asymptotic behaviour (35) but, in general, with power-law corrections:

ρ K = R K (1 + δ K ) , (50)
where δ K → 0 for K → ∞. A comparison with Eq. (37) then leads to the following correction factor to the saddle point contribution

exp - K δ K µ c ζ(λ c ) = exp Kδ K ln(-λ c ) . (51)
Consistency with the preceding analysis, cancellation of the endpoint contribution in the integral and convergence analysis, then requires

δ K = O(K -1/α ) ( 52
)
and this is also what is generally found.

IV. OPTIMIZED SUMMATION METHODS

A. ODM summation with α = 5/4

We now consider the spectrum of the Hamiltonian (1), which can be obtained by solving the time-independent Schrödinger equation

- 1 2
∂ 2 ∂x 2 + 1 2 x 2 + i 6 √ g x 3 ψ(x) = Eψ(x) (53)
with appropriate boundary conditions. The spectrum has a perturbative expansion in integer powers of g with real (rational) coefficients alternating in sign, of the form (2), which can be derived up to large orders by standard techniques. Instanton calculus, based on a steepest descent evaluation of the corresponding path integral, allows to calculate the classical action A relevant for the large order behaviour (3). One finds

A = 24 5 . (54)
According to Eq. ( 15), the strong coupling expansion of an eigenenergy E(g) of the cubic contains a coefficient g 1/5 (-1/3) g -6/5 = -1 3 g -1 . We can subtract this term and introduce the function

F (g) = 1 3 + g E(g) , (55)
which has a regular small g expansion and a large g expansion of the form

F (g) = g 6/5 E [qc] (g -4/5 ) = g 6/5 ∞ K=0 E [qc] K g -4K/5 , (56)
where the E

[qc]

K coefficients have been defined in Eq. ( 19). Our aim is to investigate to which extent the strongcoupling expansion is described by an ODM summation of weak-coupling perturbation theory. To this end, according to the discussion made in Appendix A of Ref. [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF], we now introduce the mapping (42a) which reads

g = ρ λ (1 -λ) 5/4 , (57)
and set

F g(λ) = (1 -λ) -3/2 φ(λ, ρ). (58)
The function φ has a regular expansion both at λ = 0 (g = 0) and λ = 1 (g → ∞). The analysis of Sec. III then yields µ Although the summation methods apply to all eigenvalues, from now on we consider only the ground state energy. The first terms of its perturbative expansion are

E(g) = 1 2 + 11 288 g - 930 (288) 2 g 2 + O(g 3 ) . ( 59
)
The convergence analysis of Ref. [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF], partially reproduced in Sec. III, clearly indicates that in this problem, one expects a convergence of the ODM method only in e -CK ζ with ζ < 1, and this is confirmed here. We find strong evidence for a finite distance singularity. It is consistent with this property that we can fit the zeros ρ K of P ′ K (ρ) using the functional form

ρ K = Aµ (a) c K 1 - 12.94 (K + 3) 4/5 + 11.97 . ( 60
)
The zeros of P k can be fitted by an analogous formula but with slightly different coefficients for the correction term. It has proven convenient to use the latter expression in the ODM summation in the place of the exact values.

There is strong numerical evidence for a singularity at g ≈ 0.687 exp(5 i π/4). (see Sec. VII B). One thus expects the ODM approximants of order K to converge to the energy eigenvalues as (see Eq. (45a)) constant × exp -13.8 -10.23 Re (g -4/5 ) K 1/5 , (61) consistent with a domain of convergence Re g -4/5 = Re χ > -1.351 .

(62)

Indeed, this expression reproduces well the region where the convergence factor becomes close to unity. The domain of convergence contains a segment on the real negative χ axis not contained in the convergence domain obtained by the ODM method (42c) of Ref. [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF], which is

Re g -2/5 = Re √ χ > -0.0826 . ( 63
)
The expression (51) combined with a fit of the relevant zeros of the polynomials P ′ K yields a coefficient of about 17.4 instead of 13.8 in Eq. ( 61). Finally, a rough numerical analysis of the convergence toward E

[qc] 0 rather yields exp(-18 K 1/5 ). This implies that either the asymptotic regime is reached very slowly, and this is the most likely explanation, or our heuristic convergence analysis in this particular example is oversimplified. This problem requires a further analysis.

Compared with the ODM method of Ref. [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF], the convergence is expected here to be slower at small g and, eventually, at extremely large K. At order 150, for g real and positive, the transition occurs for |g| slightly larger than 1. Of course, the method with α = 5/2 also converges in a larger sector, up to |arg(g)| = 5π/4 instead of |arg(g)| = 5π/8 for α = 5/4.

Finally, because the convergence is smooth (unlike what happens with α = 5/2 methods) the precision of the results can still be improved (see Tables I, II) by using repeatedly the algorithm that to any sequence S n substitutes

S n → S n S n+2 -S 2 n+1 S n + S n+2 -2S n+1 . (64)
Empirically, we find that it is most advantageous to apply this method to the odd and even orders separately. This is a variant of Aitken's ∆ 2 process as it is known in the literature. Then at order 150, the method (42a) remains the most efficient one for large coupling |g| ≥ 1.

B. ODM summation with α = 5/2

In Ref. [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF], the ODM method with the mapping (42c),

g = λ (1 -λ) 5/2 , ( 65
)
was applied directly to the energy eigenvalue E(g), which has a large-order expansion of the form

E(g) = g 1/5 ∞ K=0 E [sc] K g -2k/5 . (66)
This expansion effectively becomes a large-order expansion in powers of g -4/5 if one subtracts the term of order -1 3 g -6/5 . The latter is the only term not accounted for by the powers of g -4/5 .

The numerical results in Sec. V justify the introduction of the function F (g), which is defined in Eq. ( 56), and the mapping (42a) that precisely takes the special form of the large g expansion into account. However, while the method (42a) is more efficient for g large, the method (42c) yields more precise results for K very large and converges in a larger sector of the complex plane.

It is thus natural to investigate an ODM method, still with α = 5/2, that incorporates in a better way the properties of the large g expansion than the method given by Eq. (42c). To this end, we again introduce the function (55) and modify the form of the mapping to transfer the property of even powers of g -2/5 into even powers of 1-λ. A possible choice for a modified mapping is the mapping (42b),

g = ρ λ(1 -λ/2) (1 -λ) 5/2 (67)
⇒ 2g ρ -2/5 = (1 -λ) 1 -(1 -λ) 2 -2/5 . Also, (1 -λ) -3 = 2g ρ 6/5 1 -(1 -λ) 2 -6/5 . ( 68
)
With this mapping, the function

φ(λ) = (1 -λ) 3 F g(λ) (69)
is indeed only a function of (1-λ) 2 . However, in contrast to the previous mapping, the symmetry g -2/5 → -g -2/5 is only implemented approximatively, at a finite order, because the series in powers of λ is truncated. Since one cannot expect geometric convergence of the ODM method defined by Eq. (42b), the relevant value of µ is µ For the ground state energy, we empirically find that the location of the zeros of P ′ K can be fitted by

ρ K = A µ c K 1 - 5.0 (K + 2) 4/5 + 4.6 . ( 70
)
One also expects a term decreasing only like K -2/5 but its coefficient is apparently too small to be detected compared to the larger K -4/5 contribution. The same remark actually applies to the method (42c). We have not compared the convergence of methods (42b) to method (42c) systematically, but methods (42b) is definitively better for g → ∞ than with the initial α = 5/2 mapping of Eq. (42c) (see Tables I andII). For example, at order 55, with the new mapping (42b), one finds the value E

[qc] 0 = 0.372 545 [START_REF] Zinn-Justin | Perturbation series at large orders in quantum mechanics and field theories: application to the problem of resummation[END_REF] for the leading strong-coupling coefficient defined according to Eq. ( 19), corresponding to a relative error of the order of 10 -7 , that is, two order of magnitudes better as compared to the mapping (42c). At order 150, one finally obtains E

[qc] 0 = 0.372 545 78(6), a value that differs by 4 × 10 -9 from the strong coupling limit in (92).

Values of the ground state energy E(g) for various real, positive and negative values of g, and order 55 and 150, are reported in Tables I-IV.

The general conclusion is: at the order we calculate, as expected, the results for the modified method (42b) are more precise than for the more straightforward method (42c) if g is large, but comparable or worse for g small. However, for |g| ≥ 1, results given by the method (42a) with α = 5/4 are always the best. This is consistent with the fact that higher-order corrections of the form g -4K/5 are well represented by the corresponding expansion for E [qc] (g). Note that the correction of relative order g -6/5 for E [sc] (g) provides for a certain irregularity in the expansion [START_REF] Seznec | Summation of divergent series by order dependent mappings: Application to the anharmonic oscillator and critical exponents in field theory[END_REF]. As a conclusion, the advantages of the modified version (42b) seem to be somewhat limited, and we have not investigated the method more thoroughly, but we have included the results in the tables to stress the general consistency of the different implementations of the ODM summation methods.

V. NUMERICAL RESULTS: GROUND STATE ENERGY

We first display a few typical results obtained for g finite, positive and negative, by the various methods we discuss in this article: the ODM method with α = 5/4, the methods of Sec. IV B and Ref. [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF] and, finally, the continued fraction of Sec. VII B below.

A. Real positive axis

For g finite and real positive, we have chosen the same values as in Ref. [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF] and added results for the case g = 0.5. Most results are displayed in Table I for the results at order 55, to allow for a direct comparison with [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF] and give an idea about the rate of convergence, and in Table II at order 150. We use everywhere the convention that the numerical uncertainty applies to the digit in parentheses (this may also affect the preceding digit by a shift ±1).

For g = 0.5, at order 150, with the ODM method given I, but for a transformation order of 150. For weak and moderate coupling, the method (42c) now takes over. For strong coupling, the convergence of method (42a) with α = 5/4 clearly is superior to that of the other ODM methods. For completeness, we give also the results obtained from the simple continued fraction method (99). For the method defined by (42a), the numerical results are post-accelerated using the algorithm (64). Clearly, for small values of g, that is, g < 1, at order 150 the method of [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF] gives the most precise results. The entries for E(1.0) exceed the maximum number of columns available in Table II; for results obtained at transformation order 55, see Table I. Finally, to compare with the Padé summation of [START_REF] Bender | Numerical evidence that the perturbation expansion for a non-Hermitian PTsymmetric Hamiltonian is Stieltjes[END_REF], we give results for g = 288/49, which corresponds to λ 2 = 1/7 in [START_REF] Bender | Numerical evidence that the perturbation expansion for a non-Hermitian PTsymmetric Hamiltonian is Stieltjes[END_REF]. In Ref. [START_REF] Bender | Numerical evidence that the perturbation expansion for a non-Hermitian PTsymmetric Hamiltonian is Stieltjes[END_REF], at order 192, the reported result is E = 0.61273 81063 88986. In Ref. [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF], using the method (42c), at order 55, one finds

g 0.
E( 288 49 ) = 0.61273 81063 89(1) . ( 78
)
With α = 5/4, using method (42a) and convergence acceleration, at order 55, one obtains E( 288 49 ) = 0.61273 81063 88984 124 [START_REF] Lipatov | Divergence of the perturbation-theory series and the quasi-classical theory[END_REF] .

(

) 79
This improves over the ODM method (42c) in Ref. [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF] by six orders of magnitude. Then, at order 150, the result becomes

E( 288
by a numerical solution of the Schrödinger equation. The general conclusion is: for values g ≤ 1, at order 150 the method with the mapping (42c) of Ref. [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF] is the most precise, and the advantage increases with increasing order and decreasing parameter g. By contrast, for g > 1, at least up to order 150, the method with α = 5/4 takes over. Note that from now on and, in particular in the tables, we only quote the results obtained by the method (42a) after extrapolation.

Finally, the continued fraction of Sec. VII B, constructed from the strong coupling expansion (but which incorporates additional information about level merging), provides a rather good representation of the function in a wide domain.

For completeness, in the spirit of reference [START_REF] Bender | Numerical evidence that the perturbation expansion for a non-Hermitian PTsymmetric Hamiltonian is Stieltjes[END_REF], we give some indication about the summation of the perturbative series by a continued fraction expansion. We define (for the ground state)

h 0 (g) = E(g)/E(0), (83)
and introduce the relation

h p-1 (g) = 1 + κ p g h p (g) , h p (0) = 1 , (84)
which allows calculating h p from h p-1 and determining recursively the coefficients κ p . The truncated continued fraction (obtained by replacing, at some order p, h p (g) by 1) generates, alternatively, [n + 1/n] and [n/n] Padé approximants. For a Stieltjes function, all coefficients κ p are positive. Moreover, for a divergent series with a large order behaviour of the form (3) one expects the coefficients κ p to grow asymptotically linearly with p and the error for the continued fraction truncated at order p to behave like exp(-C p/g). This is indeed what is observed. The coefficient of the linear term is compatible with 5/48. More precisely, a good fit is κ p = (10 p + 3 × (-1) p )/96 + O(1/ √ p) for p → ∞. Also, the coefficient of p/g extrapolates to C = 2 48/5 = 6.19 . . . with a good precision. Compared to the various ODM summations, the convergence is limited to the first Riemann sheet, is poorer than the ODM method with α = 5/2 for g small and is much poorer than the method with α = 5/4 for g large. For example, at order 150, for g = 0.5 the error is about 1.5 × 10 -42 , for g = 1 about 6 × 10 -31 , for g = 5 about 10 -14 and about 2 × 10 -7 for g = 21.6.

B. Negative real axis

Results for four values of g < 0 are displayed in Tables III and IV. From the numerical evidence, we conclude that ODM methods with α = 5/2 converge also for g < 0 but, from the theoretical analysis, the convergence is expected to be poorer than for g > 0, in agreement with the data. The ODM method (42b) with α = 5/4 converges also on the real negative axis for |g| large enough, that is, |g| > 0.53 approximately.

From the analysis of Ref. [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF], we know that the imaginary part itself is an analytic function with singularities on the real negative axis only at g = 0 and at infinity. Moreover, it is a simple positive decreasing function. Compared with the method with α = 5/2, the positivity of the imaginary part for g = -|g| + i0 can again be verified, but with higher precision for g large. For |g| → ∞ and Im (g) = 0 + , from the leading term in Eq. ( 92) suitably rotated into the complex plane, one now infers

E(g)/|g| 1/5 = g→-∞ 0.3013958756586835717823(7) + 0.2189769214314493762936(0) i , (85)
a result that uses the perturbative series up order 150. For g = -0.5, the method (42b) with α = 5/4 is not to converge and, indeed, numerical indications are that this is the case. The result is somewhat equivalent to a direct summation of the initial asymptotic series, with an error of the order of the imaginary part TABLE III: On the real negative axis, at order 150, the real part Re E(g) is approximated by several variants of ODM methods. The straightforward ODM method (42c) provides the best results for weak coupling, whereas method (42a) with α = 5/4, post-accelerated, is superior in the moderate and strong coupling domain. Note that the simple continued fraction method (99) also provides for a well converged answer.

-g 0.5 on the cut. However, compared to the direct summation of the series, in some range the successive results first oscillate around the exact value with the order K with a rather slowing increasing amplitude rather than blowing up. Otherwise, we notice a situation comparable to the real positive axis. For |g| < 1, the ODM method of [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF] yields the most precise results, while for |g| ≥ 1, the method (42a) with α = 5/4 and convergence acceleration takes over. However, the continued fraction of Sec. VII B gives the most precise results in the intermediate range and is equivalent for |g| large.

VI. STRONG COUPLING EXPANSION AND LEVEL CROSSING

We here demonstrate that considerable information on subleading corrections to the strong coupling asymptotics can be obtained by investigating the ODM resummed weak-coupling expansion of the energy levels.

In order to explore more thoroughly the strong coupling expansion and the analytic properties of the ground state energy of the cubic Hamiltonian, it is very useful to also consider the linear coupling Hamiltonian [START_REF] Bender | Anharmonic oscillator[END_REF],

H [qc] = - 1 2 ∂ 2 ∂x 2 + i 1 6 x 3 + 1 2 χ x , (86)
because a perturbative expansion in the variable

χ = g -4/5 , (87)
of the energy level of the Hamiltonian H [qc] , defined in Eq. ( 17), is equivalent to a large g expansion of the corresponding energy of the Hamiltonian (1). We recall here that the two Hamiltonians are only equivalent for χ > 0, the negative χ axis corresponding to arg(g) = 5π/4. For real g, the ground state energy E(g) of the Hamiltonian (1) is a real analytic function. From the numerical evidence provided by the ODM summation in Ref. [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF], we conjecture that, that, in the variable g -1/5 , the eigenvalues can be continued up to | arg g -1/5 | = π/4. This implies that in the variable χ, the eigenvalues E [qc] (χ) are real analytic with singularities only on the negative real axis. The series in powers of χ are convergent in a disk. In the case of the ground state, at the singularity nearest to the origin (numerical results indicate χ = χ c = -1.3510 . . .). Note that, by contrast with a Hermitian Hamiltonian, PT symmetric Hamiltonians do not experience eigenvalue repulsion, which explains why it is possible for eigenvalues to merge for χ real.

A. Strong coupling expansion

We now concentrate on the behaviour of E(g) for g → ∞, or, equivalently, on the small χ expansion given in Eq. ( 19) of E [qc] (χ):

E [qc] (χ) = K E [qc] K χ K . (88)
First, we have calculated the values of the coefficients of the large coupling expansion as determined by the ODM method (42a) with α = 5/4, with the improvement by the algorithm (64). The successive coefficients of the small χ expansion are related to φ(λ, ρ) and its derivatives taken at λ = 1 (see Appendix A for the first terms). At leading order, one finds

φ(1, ρ) = ρ 1/5 E [qc] 0 . (89)
For the leading term, we obtain

E

[qc] 0 = 0.37254 57904 52207 09825 0601(1).

Note that even at order 150, with the ODM method (42c) for the coefficient E

[qc] 0 one still finds an uncertainty of 4.0 × 10 -6 . With the method (42b), one obtains E

[qc] 0 = 0.372 545 78 [START_REF] Zinn-Justin | Perturbation series at large orders in quantum mechanics and field theories: application to the problem of resummation[END_REF].

From the ODM method (42a) with α = 5/4 applied to the perturbative expansion of the function (56), at order 150, we have then determined, using the algorithm described in Appendix B, with decreasing relative precision (about 10 -6 for the last term), all terms up to order K = 20. The analysis of the behaviour of the coefficients with increasing order strongly suggests the existence of a square root singularity located at χ = χ c = -1.351 ± 0.002. A square root singularity is consistent with the existence of a level merging. To confirm this analytic structure and locate the singularity more precisely, we have assumed that the point χ c corresponds to a merging between the ground state energy and the first excited state energy.

B. Level merging and strong coupling expansion

In addition to the ground state, as described, we have also generated the perturbative series for the first excited state. We do not report here the details of the numerical study of the energy of the first excited state. The behaviour of the series with respect to the ODM summation methods is very similar. The coefficients of the strong coupling expansion can be determined with comparable precision. The convergence properties of the ODM method (42c) indicate that, again, the eigenvalue is analytic in a cut-plane and that the first singularity corresponds to the level merging with the ground state energy.

We have then formed the two symmetric combinations, the half sum S 01 and the half difference squared ∆ 01 (GS and ES stand for ground state and first excited state, respectively):

∆ 01 = 1 4 E [qc] GS -E [qc] ES 2 , S 01 = 1 2 E [qc] GS + E [qc] ES . (90)
If the ansatz is correct, the two functions are not singular at χ c and ∆ 01 must vanish linearly at χ c . This is indeed what is observed (see Figs. 3 and4). The result confirms that at χ c the eigenvalues corresponding to the ground state and the first excited state, merge and for χ < χ c become complex conjugate.

More precisely, the direct summation of the series for ∆ 01 with the ODM method of Sec. IV A shows that ∆ 01 vanishes linearly with χ at the point

χ c = -1.3510415966(3) , (91)
a result fully consistent with the direct study of the ground state. At this point, E [qc] (χ c ) = 0.41330579447 [START_REF] Shin | On the reality of eigenvalues for a class of PT -symmetric oscillators[END_REF]. As a by-product, one also obtains for the ground state, the value of E [qc] (χ = -2 4/5 ) = 0.38985020(5) -0.364427(9)i, which confirms the result coming from the continued fraction in Table VI.

By the same method as for the ground state, we have determined the strong coupling expansions of S 01 and ∆ 01 . From these strong coupling expansions, one can recalculate the strong coupling expansion of the ground state energy. The results are completely consistent with those from the direct expansion. However, the precision is improved for the higher order coefficients as expected since the singularity at χ c is now explicitly generated.

As a necessary ingredient for a more precise determination of the coefficients of the continued fraction of Sec. VII B for the ground state, we have calculated the strong-coupling expansions of the difference ∆ 01 and the sum S 01 up to to order 28. We have then inferred the strong coupling expansion of the ground state at the same order although we here give only the first 20 terms. These are useful for reference purposes and read as follows: (92)

E [qc] (χ) = 0.
Note that the errors are strongly correlated. This expansion is also consistent with the 20 first coefficients reported in Ref. [START_REF] Fernandez | Strong-coupling expansions for the PT -symmetric oscillators V (x) = a(ix) + b(ix) 2 + c(ix) 3[END_REF], which have been determined with a 10 -10 -10 -11 relative precision by a numerical solution of the eigenvalue equation. It is also consistent with results obtained for the few first terms from a numerical solution of the Schrödinger equation [START_REF] Jentschura | Calculation of the Characteristic Functions of Anharmonic Oscillators[END_REF], which have relative errors of order 10 -9 .

VII. CONTINUED FRACTION AND ANALYTIC PROPERTIES A. General remarks

To discuss more thoroughly the analytic properties of the ground state energy, it is convenient to consider the Hamiltonian [START_REF] Bender | Anharmonic oscillator[END_REF] and parameterize the energy in terms of the coupling χ.

Inverting the relations [START_REF] Grecchi | Padé summability of the cubic oscillator[END_REF] and [START_REF] Kleinert | Convergence behavior of variational perturbation expansion: A method for locating Bender-Wu singularities[END_REF], one finds a form relevant for the large χ behaviour of E [qc] (χ):

E [qc] (χ) = 1 3 χ 3/2 + χ 1/4 E(χ -5/4 ) . (93)
The perturbative expansion of E(g) in g translates into a large χ expansion, because large χ corresponds to weak coupling g,

χ 1/4 E(χ -5/4 ) = 1 2 χ 1/4 + 11 288 χ -1 - 155 13824 χ -9/4 + • • • • (94)
In particular, the imaginary part for χ = -|χ| + i 0 has the expansion

Im E [qc] (χ) = - 1 3 |χ| 3/2 + 1 2 √ 2 |χ| 1/4 + 155 13824 √ 2 |χ| -9/4 + • • • • (95)
For |χ| ≫ 1, the imaginary part is clearly negative. The question that arises is its sign on the cut for |χ| of order unity.

Note that the large χ behaviour (93) implies that the Cauchy representation for the perturbative expansion coefficient E

[qc]

K can be written as

E [qc] K = 1 π χc -∞ Im E [qc] (χ) χ K+1 dχ , (96)
converges only for K > 1. The conjecture Im E [qc] (χ) ≤ 0, which will be substantiated in Sec. VII B, is thus consistent with the signs of coefficients of the expansion (92). Finally, some additional information can be obtained directly from the ODM method with α = 5/2 (see Table VI) and also by summing the small χ expansion (92) as we will show.

B. Continued fraction

To continue the expansion (92) outside the circle of convergence, we introduce the continued fraction expansion of

Ẽ[qc] (χ) = E [qc] (χ) -E [qc] (0) E [qc] ′ (0)χ (97) motivated by the prejudice that Ẽ[qc] is a Stieltjes func- tion since Im Ẽ[qc] (χ) = Im E [qc] (χ) E [qc] ′ (0) χ (98) is positive if Im E [qc] 0
is negative and Ẽ[qc] behaves like χ 1/2 for |χ| → ∞. We then construct the continued fraction expansion for the function Ẽ[qc] by introducing the recursion relation

f p-1 (χ) = 1 + a p χ f p (χ) , f p (0) = 1 , (99)
with the initial condition f 0 (χ) = Ẽ[qc] (χ). It allows us to calculate the coefficients a p recursively. The truncated continued fractions, obtained by replacing at some order f p (χ) by 1, generate, alternatively, [n + 1/n] and [n/n] Padé approximants, depending on whether p is odd or even.

For a Stieltjes function, all coefficients a p are positive. Indeed, using the coefficients (92), which are obtained from the series at order 150, we find a p > 0 (see Table V) for p ≤ 27 and the coefficients seem to converge slowly, within errors, toward the limit expected for a square root singularity at χ = χ c : 1/(4χ c ) = -0.18504241514 . . .. For example, the two last coefficients determined with sufficient precision are a 24 = 0.1850464 [START_REF] Brézin | Perturbation theory at large order. I. The ϕ N interaction[END_REF] and a 25 = 0.185042 [START_REF] Brézin | Perturbation theory at large order. I. The ϕ N interaction[END_REF].

To test the idea further, we have substituted in the infinite continued fraction

f pmax (χ) → 1 2 + 1 2 1 + 4γχ , (100)
p max varying from 9 to 27. This amounts to taking a p constant, a p = γ ≈ 0.1850 for p > p max . Remarkably enough, for χ → ∞, the corresponding approximation for E [qc] (χ) behaves like χ 3/2 with a coefficient close to the exact value 1/3. We have thus adjusted γ more precisely to get 1/3 exactly. This yields a very stable sequence up to p = 25, converging slowly toward the expected limit 0.18504 . . .: C. Spectrum for χ < 0

The convergence of the various methods is rather poor for χ < 0, specially in the neighbourhood of the cut. Table VI displays a few results for E [qc] (χ). According to Table VI, the most precise results are in general obtained from the modified continued fraction, with a precision that for |g| > 1 is comparable with the ODM method defined by Eq. (42a). Note that for the first value, which corresponds to |g| = 0.5, the ODM method with α = 5/4 does not converge, as expected, and does not even provide a reliable estimate. The method of Sec. IV B is not useful either. Finally, in the neighbourhood of the singularity χ = χ c a calculation based on the determination of ∆ 01 and S 01 gives the most precise results (see Figs. 3 and4).

The imaginary part after the cut for χ = -|χ| + i 0 is necessarily negative when calculated from the approximated continued fraction. Therefore, the general consistency of the results coming from the continued fraction and the other ODM methods provides an additional confirmation of our conjecture about the sign of the discontinuity.

VIII. CONCLUSIONS

Let us briefly summarize the findings of the current numerical study of properties of the cubic anharmonic oscillator. For g ≥ 0, the Hamiltonian (1) is PT symmetric and has a real positive spectrum. The energy eigenvalues have divergent, Borel summable expansions in powers of the coupling constant g. The imaginary parts of the eigenvalues on the cut on the real negative axis are positive. Starting from the PT symmetric case and developing a strong coupling expansion in the variable √ ξ = g -1/5 , one sees that eigenvalues are univalued functions in some neighbourhood of the origin in the ξ plane, [qc] (χ), obtained by various ODM summation methods, on the real negative axis, at transformation order 150. Method (42a) and the continued fraction approximation (99) provide the best numerical convergence. with a simple pole at z = √ ξ = g -1/5 = 0 as a unique singularity.

In Sec. III, we have discussed various possible summation methods for the determination of energy eigenvalues of the cubic potential. Notably, we explore the numerical properties of three such methods here, given in Eqs. (42a), (42b) and (42c). The latter was used in our previous paper [START_REF] Zinn-Justin | Order-dependent mappings: strong coupling behaviour from weak coupling expansions in non-Hermitian theories[END_REF]. Summing the perturbative expansion of the ground state energy E(g) by various implementations of the ODM method, we uncover additional properties. From the apparent convergence of the ODM method for α = 5/2 [see Eq. (42c)], we conclude that in the variable z the function E(g(z)) has no other singularity in the sector |arg z| < π/4. In the variable χ = g -4/5 = z 4 , this translates into the property that the ground state energy E [qc] (χ) of the Hamiltonian (5) also is a real analytic function in a cut plane with a cut along the real negative axis.

The ODM summation method (42a) converges very well in the strong coupling regime and thus allows us to determine a number of terms of the large g, thus small χ expansion. The precision can be further improved by the convergence acceleration algorithm (64) and, finally, by a combined calculation of the ground state and first excited energies. The first 20 terms are given in Eq. (92). A direct analysis of the strong coupling expansion indicates the existence of a square root singularity located at χ = χ c ≈ -1.351 [g ≈ 0.687 exp(5iπ/4) in the second sheet]. These properties are confirmed by the direct calculation of the difference (E GS -E ES )

2 , which vanishes linearly at the point χ c = -1.3510415966(3), as given in Eq. ( 91). The strong coupling expansion, in turn, can be summed by an expansion in a continued fraction (Sec. VII B). The calculated coefficients of the continued fraction display an asymptotic behavior consistent with a square root singularity. The extrapolated continued fraction yields results in remarkable agreement with the more direct ODM calculations.

As a function of z, E [qc] (χ(z)) is invariant under the rotation z → z exp(iπ/2). It is thus entirely defined by its values for -π/4 ≤ arg(z) ≤ +π/4, that is, -5π/4 ≤ arg(g) ≤ +5π/4. As a consequence, as a function of z, E has cuts only on the lines arg z = nπ/4 where level merging can occur. Note that, from the point of view of Eq. ( 17), for χ < 0 the Hamiltonian is still PT symmetric and, thus, the singularity at χ = χ c corresponds to a kind of spontaneous PT symmetry breaking.

All coefficients of the associated expansion in a continued fraction are found to be positive as long as we can estimate them with enough precision, that is, up to order 27 (corresponding to a strong coupling expansion up to order 28), and seems to converge toward a positive value consistent with the square root singularity of the first level merging. This gives a very strong argument that the eigenvalue of the strong coupling Hamiltonian Im E [qc] (-|χ| + i 0) is negative on the cut. This conjecture is also consistent with the large χ behaviour, the Cauchy representation, the coefficient of the square root singularity. However, this property cannot be shared by all eigenvalues and may even be unique to the ground state.

To go beyond this study, one would have to apply these ODM summation methods more directly to the spectral equation, in order to eliminate all level merging singularities (for the pure cubic potential see Ref. [START_REF] Voros | Airy function -exact WKB results for potentials of odd degree[END_REF]). Finally, we would like to emphasize the empirical evidence gathered during the current study, for the robustness of the ODM summation methods (42a), (42b) and (42c), which applied to different functions with varying implementation always give consistent results.

FIG. 1 :

 1 FIG.1:(Color online.) Image in the complex λ plane of the negative real g axis under the mapping g = ρ λ/(1 -λ) 5/2 , for ρ = 1/2. A point infinitesimally displaced above the negative real axis is mapped onto a point in the upper complex λ plane.

= 3 .

 3 811 522 . . . and thus R = A µ (a) c = 18.295 306 . . . in the notation of Sec. III B [see Eq. (42a)].

c = 4 .

 4 445 762 607 . . . (see Sec. III B). It follows that C (b) 2 = (R/2) 2/5 = (A µ (b) c /2) 2/5 = 3.401 535 . . . [see Eq. (45b)].

  49 ) = 0.61273 81063 88984 12476 20895 52(6) . (80) With the method (42c), at the same order 150, one obtains E( 288 49 ) = 0.61273 81063 88984 12476(3) .(81)We have also verified that E( 288 49 ) = 0.61273 81063 88984 12476 . . .

13 + 5 . 14 - 3 . 15 + 2 . 16 - 1 . 17 + 1 .

 135143152161171 37254 57904 52207 09825 06011(5) + 0.36753 58055 44193 60353 04(6) χ + 0.14378 77004 15066 51583 39(0) χ 2 -0.02658 61056 27059 38713 52(9) χ 3 + 0.00988 71650 79200 88729 05(5) χ 4 -0.00461 00192 93623 15160 2(3) χ 5 + 0.00240 93426 35048 47521 1(7) χ 6 -0.00134 88515 29319 85498(8) χ 7 + 0.00079 06119 76816 97837(2) χ 8 -0.00047 88478 41414 5725(4) χ 9 + 0.00029 72375 58426 7145(5) χ 10 -0.00018 80657 95326 713(9) χ 11 + 0.00012 08255 49560 587(6) χ 12 -7.86045 58627 946(5) × 10 -5 χ 16744 64642 199(1) × 10 -5 χ 42729 47828 030(3) × 10 -5 χ 29050 20869 87(5) × 10 -5 χ 54091 59219 76(3) × 10 -5 χ 04266 03452 042(4) × 10 -5 χ 18 -0.70913 87535 56(4) × 10 -5 χ 19 + 0.48450 81660 66(0) × 10 -5 χ 20 + O(χ 21 ) .

γ

  17 = 0.1850687 , γ 18 = 0.1850261 , γ 19 = 0.1850603 , γ 20 = 0.1850239 , γ 21 = 0.1850562 , γ 22 = 0.1850310 , γ 23 = 0.1850562 , γ 24 = 0.1850367 , γ 25 = 0.185048(4).(101) corresponding to a singularity located at χ = χ c ≈ -1.3504, consistent with a direct analysis of the behaviour of the coefficients in the expansion (92). As a consequence, in all Tables and Figures we have reported results obtained from this approximated continued fraction.

χ -2 4 / 5 = 1 4 / 5 =ES ) 2 / 4

 4514524 FIG. ∆01 = (E [qc]GS -E [qc] ES ) 2 /4 as a function of χ for -1.85 ≤ χ ≤ -0.85 [the value at the origin is ∆01(χ = 0) = 0.2263073 . . .].

FIG. 4 :

 4 FIG. 4: S01 = (E [qc] GS + E [qc] ES )/2 as a function of χ for -1.85 ≤ χ ≤ -0.85 [the value at the origin amount to S01(χ = 0) = 0.8482634 . . .].

TABLE I :

 I The ODM mappings given in Eq. (42a), (42b) and (42c) are compared with respect to their convergence properties on the real positive axis at a transformation order of 55. Method (42a) with α = 5/4 and extrapolation clearly provides the best answers for large and moderate coupling. The modified method (42b) with α = 5/2 interpolates between method (42a) and (42c); indeed, the application of the straightforward ODM method (42c) still provides for the best convergence at weak coupling.

	g	0.5	1
	Method of Eq. (42a) 0.51689 17642 53171 97821 11588 95(6)	0.53078 17593 04176 67113 556(1)
	Method of Eq. (42b) 0.51689 17642 53171 97821 11588(9)	0.53078 17593 04176 67113(4)
	Method of Eq. (42c) 0.51689 17642 53171 97821 11588 9(6)	0.53078 17593 04176 671135(6)
	g	5.0	21.6
	Method of Eq. (42a) 0.60168 39332 05191 96159(1)	0.73340 99204 85427 96(4)
	Method of Eq. (42b) 0.60168 39332 05(2)	0.73340 9920(5)
	Method of Eq. (42c) 0.60168 39332 05(2)	0.73340 992(1)

TABLE II :

 II Same as Table

TABLE IV :

 IV The imaginary part Im E(g) at transformation order 150 shows the same characteristics as the real part discussed in TableIII.

	-g	0.5	1.0
	Method of Eq. (42a)	0.00026 6662(1)	0.01551 79258 2059(4)
	Method of Eq. (42b)	0.00026 66618 8(2)	0.01551 79258 2(0)
	Method of Eq. (42c)	0.00026 66618 82408(1)	0.01551 79258 20(6)
	Continued fraction	0.00026 66618 824(6)	0.01551 79258 20594 2572(2)
	-g	5.0	21.6
	Method of Eq. (42a)	0.18385 80861 86171 172893(3)	0.35140 17775 93691 93624 451(6)
	Method of Eq. (42b)	0.18385 80861(9)	0.35140 1777(6)
	Method of Eq. (42c)	0.18385 8086(0)	0.35140 18(0)
	Continued fraction	0.18385 80861 86171 17289 33(1) 0.35140 17775 93691 93624 45(1)

TABLE V :

 V Coefficients ap of the continued fraction for K = 150. The singularity is consistent with a limit 0.1850424 . . ..

	p = 1	p = 2	p = 3	p = 4
	0.39122 09320 72635 98993 0.18489 83296 22869 56168 0.18699 38616 55333 76095 0.18768 40668 91881 09149
	p = 5	p = 6	p = 7	p = 8
	0.18519 10686 95077 9010(9) 0.18477 49761 41944 774(1) 0.18470 22248 17836 84(5) 0.18507 38286 03338 2(9)
	p = 9	p = 10	p = 11	p = 12
	0.18513 42403 94894(8)	0.18510 36391 9575(7)	0.18499 98937 209(1)	0.18501 58945 34(4)
	p = 13	p = 14	p = 15	p = 16
	0.18504 17233 0(4)	0.18506 43232 9(6)	0.0.18504 19101(5)	0.18503 66878(1)

TABLE VI :

 VI Values for E

Acknowledgements

J.Z.-J. gratefully acknowledges CERN's hospitality, where a major part of this work was completed. U.D.J. acknowledges support by a Grant from the Missouri Research Board and by the National Science Foundation (Grant PHY-8555454).

Appendix A: ODM and strong coupling expansion

We assume that we know a perturbative expansion for E(g) in powers of the coupling constant g, and that E(g) also has a strong-coupling asymptotic expansion of the form,

a property shared by the example [START_REF] Bender | Anharmonic oscillator[END_REF] we discuss here. It is also shared by the quartic anharmonic oscillator [START_REF] Seznec | Summation of divergent series by order dependent mappings: Application to the anharmonic oscillator and critical exponents in field theory[END_REF] and all x N perturbations to the quantum harmonic oscillator. We then consider the conformal mapping (41),

This transformation maps the real positive g axis onto the finite λ interval [0, 1]. For g → ∞, λ → 1 and g has an expansion at λ = 1 of the form

with Λ 0 = ρ -1/α . The function

then has a Taylor series expansion at λ = 0,

as well as at λ = 1,

with ϕ 0 = ǫ 0 ρ β , where ǫ 0 is the coefficient defined in (A1). This last property explains, to a large extent, the good convergence of the method even for g → ∞.

Appendix B: Large g expansion: a few terms

In order to determine the successive terms of the large g, small χ expansion we do not have to differentiate numerically the values of E(g) calculated at g large but finite. In the ODM method, the various terms have explicit analytic forms and we give here the first terms. First, we set

Then we write the relation (A2) as

The function ψ(λ) can then be expanded successively in powers of s and z:

where the ǫ n coefficients are given in Eq. (A1). The first terms yield

)

Appendix C: Numerical evidence for strong coupling

In Ref. [START_REF] Jentschura | Structure, Time Propagation and Dissipative Terms for Resonances[END_REF], we started from the Hamiltonian

which entails both a change in the normalization of the coupling term and also a change in the complex phase of the coupling. The transformation x → g -1/10 x brings the Hamiltonian into the form

We now define g = u -5/2 and write

A shift of the variable x → x -1 6 u then results in

g -4/5 x + 1 180 g -6/5 .

(C4)

If one now writes the strong-coupling expansion of the N th energy level in the form of Eq. ( 17) of Ref. [START_REF] Jentschura | Structure, Time Propagation and Dissipative Terms for Resonances[END_REF],

then the terms [see Eq. (C4)]

generate a convergent series in even powers of g -4/5 , and the term 1 180 g -6/5 (C7) in Eq. (C4) is the only term of odd power in g -2/5 in Eq. (C5). From the point of view of the strong coupling expansion, this term L N,3 = 1/108 represents a shift, independent of the quantum number N , of all levels of the cubic potential. This affords an explanation for the observations made in Table 2 of Ref. [START_REF] Jentschura | Structure, Time Propagation and Dissipative Terms for Resonances[END_REF], where on the basis of numerical evidence, it was conjectured that the terms of order K = 1 according Eq. (C5) vanish, and that the numerical value for the coefficient of the term of order K = 3 is 1/108 uniformly for all levels.