text:
# I Introduction

## I Introduction

Collisions at relativistic heavy ion colliders like the Relativistic Heavy Ion Collider RHIC/Brookhaven and the Large Hadron Collider LHC/CERN (operating in its heavy ion mode) are mainly devoted to the search of the Quark Gluon Plasma. However, peripheral heavy ion collisions also open up a broad area of studies as advocated by Baur and collaborators . Examples are the possible discovery of an intermediate-mass Higgs boson  or beyond standard model physics  using peripheral ion collisions, which have been discussed at length in the literature. More promissing than these may be the study of hadronic physics, which will appear quite similarly to the two-photon hadronic physics at $`e^+e^{}`$ machines with the advantage of a huge photon luminosity peaked at small energies . Due to this large photon luminosity it will become possible to discover resonances that couple weakly to the photons .

Double-pomeron exchange will also occurs in peripheral heavy ion collisions and their contribution is similar to the two-photons one as discussed by Baur  and Klein . A detailed calculation performed by Müller and Schramm of Higgs boson production have shown that the diffractive contribution is much smaller than the electromagnetic one . We can easily understand this result remembering that the coupling between the Higgs boson and the pomerons is intermediated by quarks, and according to the pomeron model of Donnachie and Landshoff  when in the vertex pomeron-quark-quark any of the quark legs goes far “off-shell” the coupling with the pomeron decreases. Therefore, we do not need to worry about the pomeron-pomeron contribution in peripheral heavy ion collisions when heavy (or far “off-shell”) quarks are present. However, this is not what happens in the case of light resonances , where double diffraction were claimed to be as important as photon initiated processes. In particular, Engel et al.  have shown that at the LHC the diffractive production of hadrons may be a background for the photonic one.

In Ref.  it was remarked that the effect of removing “central collisions” should also be performed in the double-pomeron calculation, implying in a considerable reduction of the background calculated in Ref.. This claim is the same presented by Baur  and Cahn and Jackson  in the case of early calculations of peripheral heavy ion collisions. Roughly speaking we must enforce that the minimum impact parameter ($`b_{min}`$) should be larger than ($`R_1+R_2`$), where $`R_i`$ is the nuclear radius of the ion “$`i`$”, in order to have both ions coming out intact after the interaction.

In this work we will compute the production of resonances, pion-pairs and a hadron cluster with invariant mass $`M_X`$ through photon-photon and pomeron-pomeron fusion in peripheral heavy ion collisions at the energies of RHIC and LHC. We will take into account the effect of the impact parameter as discussed in the previous paragraph for photons as well as for pomerons. We also compare this approach to cut the central collisions with the use of an absorption factor in the Glauber approximation. The inclusion of pion-pairs production is important because they certainly will be studied at these colliders, and they also represent a background for glueball (and other hadrons) detection. The pomeron physics within the ion will be described by the Donnachie and Landshoff model . We will focus on the values of the cross sections that shall be measured in the already quoted ion colliders, and point out when pomeron-pomeron processes can be considered competitive or not with photon-photon collisions. The arrangement of our paper is the following: Section 2 contains a discussion of the photons and pomerons distributions in the nuclei. In Sect. 3 we introduce the cross section for the elementary processes. Finally, Sect. 4 contains the results and conclusions.

## II Photons and pomerons distribution functions

### A Photons in the nuclei

The photon distribution in the nucleus can be described using the equivalent-photon or Weizsäcker-Williams approximation in the impact parameter space. Denoting by $`F(x)dx`$ the number of photons carrying a fraction between $`x`$ and $`x+dx`$ of the total momentum of a nucleus of charge $`Ze`$, we can define the two-photon luminosity through

$`{\displaystyle \frac{dL}{d\tau }}={\displaystyle _\tau ^1}{\displaystyle \frac{dx}{x}}F(x)F(\tau /x),`$ (1)
where $`\tau =\widehat{s}/s`$, $`\widehat{s}`$ is the square of the center of mass (c.m.s.) system energy of the two photons and $`s`$ of the ion-ion system. The total cross section $`ZZZZ\gamma \gamma ZZX`$, where $`X`$ is the particle produced within the rapidity gap, is

$`\sigma (s)={\displaystyle 𝑑\tau \frac{dL}{d\tau }\widehat{\sigma }(\widehat{s})},`$ (2)
where $`\widehat{\sigma }(\widehat{s})`$ is the cross-section of the subprocess $`\gamma \gamma X`$.

There remains only to determine $`F(x)`$. In the literature there are several approaches for doing so, and we choose the conservative and more realistic photon distribution of Ref.. Cahn and Jackson , using a prescription proposed by Baur , obtained a photon distribution which is not factorizable. However, they were able to give a fit for the differential luminosity which is quite useful in practical calculations:

$$\frac{dL}{d\tau }=\left(\frac{Z^2\alpha }{\pi }\right)^2\frac{16}{3\tau }\xi (z),$$
(3)
where $`z=2MR\sqrt{\tau }`$, $`M`$ is the nucleus mass, $`R`$ its radius and $`\xi (z)`$ is given by

$$\xi (z)=\underset{i=1}{\overset{3}{}}A_ie^{b_iz},$$
(4)
which is a fit resulting from the numerical integration of the photon distribution, accurate to $`2\%`$ or better for $`0.05<z<5.0`$, and where $`A_1=1.909`$, $`A_2=12.35`$, $`A_3=46.28`$, $`b_1=2.566`$, $`b_2=4.948`$, and $`b_3=15.21`$. For $`z<0.05`$ we use the expression (see Ref. )

$$\frac{dL}{d\tau }=\left(\frac{Z^2\alpha }{\pi }\right)^2\frac{16}{3\tau }\left(\mathrm{ln}(\frac{1.234}{z})\right)^3.$$
(5)
The condition for realistic peripheral collisions ($`b_{min}>R_1+R_2`$) is present in the photon distributions showed above, and the applications of Sect. 4 are straightforward once we determine the cross sections for the elementary processes.

### B Pomerons in the nuclei

In the case where the intermediary particles exchanged in the nucleus-nucleus collisions are pomerons instead of photons, we can follow closely the work of Müller and Schramm  and make a generalization of the equivalent photon approximation method to this new situation. So the cross section for particle production via two pomerons exchange can be written as

$`\sigma _{AA}^{PP}={\displaystyle 𝑑x_1𝑑x_2f_P(x_1)f_P(x_2)\sigma _{PP}(s_{PP})},`$ (6)
where $`f_P(x)`$ is the distribution function that describe the probability for finding a pomeron in the nucleus with energy fraction $`x`$ and $`\sigma _{PP}(s_{PP})`$ is the subprocess cross section with energy squared $`s_{PP}`$. In the case of inclusive particle production we use the form given by Donnachie and Landshoff

$`f_P(x)={\displaystyle \frac{1}{4\pi ^2x}}{\displaystyle _{\mathrm{}}^{(xM)^2}}𝑑t|\beta _{AP}(t)|^2|D_P(t;s^{})|^2,`$ (7)
where $`D_P(t;s^{})`$ is the pomeron propagator

$`D_P(t;s)={\displaystyle \frac{(s/m^2)^{\alpha _P(t)1}}{\mathrm{sin}(\frac{1}{2}\pi \alpha _P(t))}}\mathrm{exp}\left({\displaystyle \frac{1}{2}}i\pi \alpha _P(t)\right),`$ (8)
with $`s`$ the total squared c.m. energy. The Regge trajectory obeyed by the pomeron is $`\alpha _P(t)=1+\epsilon +\alpha _P^{}t`$, where $`\epsilon =0.085`$, $`\alpha _P^{}=0.25`$ GeV-2 and $`t`$ is a small exchanged four-momentum square, $`t=k^2<<1`$, so the pomeron behaves like a spin-one boson. The term in the denominator of the pomeron propagator, $`[\mathrm{sin}(\frac{1}{2}\pi \alpha _P(t))]^1`$, is the signature factor that express the different properties of the pomeron under C and P conjugation. At very high c.m. energy this factor falls very rapidly with $`𝐤^2`$, whose exponential slope is given by $`\alpha _P^{}\mathrm{ln}(s/m^2)`$, $`m`$ is the proton mass, and it is possible to neglect this $`𝐤^2`$ dependence,

$`\mathrm{sin}{\displaystyle \frac{1}{2}}\pi (1+\epsilon \alpha _P^{}𝐤^2)\mathrm{cos}({\displaystyle \frac{1}{2}}\pi \epsilon )1.`$ (9)
If we define the pomeron range parameter $`r_0`$ as

$`r_0^2=\alpha _P^{}\mathrm{ln}(s/m^2),`$ (10)
the pomeron propagator can be written as

$`|D_P(t=𝐤^2;s)|=(s/m^2)^\epsilon e^{r_0^2𝐤^2}.`$ (11)
Since we are interested in the spatial distribution of the virtual quanta in the nuclear rest frame we are using $`t=𝐤^2`$.

The nucleus-pomeron coupling has the form

$`\beta _{AP}(t)=3A\beta _0F_A(t),`$ (12)
where $`\beta _0=1.8`$ GeV-1 is the pomeron-quark coupling, $`A`$ is the atomic number of the colliding nucleus, and $`F_A(t)`$ is the nuclear form factor for which is usually assumed a Gaussian expression (see, e.g., Drees et al. in )

$`F_A(t)=e^{t/2Q_0^2},`$ (13)
where $`Q_0=60`$ MeV.

Performing the $`t`$ integration of the distribution function in Eq.(7) we obtain

$`f_P(x)`$ $`=`$ $`{\displaystyle \frac{(3A\beta _0)^2}{(2\pi )^2x}}\left({\displaystyle \frac{s^{}}{m^2}}\right)^{2\epsilon }{\displaystyle _{\mathrm{}}^{(xM)^2}}𝑑te^{t/Q_0^2}`$ (14)
$`=`$ $`{\displaystyle \frac{(3A\beta _0Q_0)^2}{(2\pi )^2x}}\left({\displaystyle \frac{s^{}}{m^2}}\right)^{2\epsilon }\mathrm{exp}\left[\left({\displaystyle \frac{xM}{Q_0}}\right)^2\right].`$ (15)
The total cross section for a inclusive particle production is obtained with the above distribution and also with the expression for the subprocess $`PPX`$ as prescribed in Eq.(6). However, in Eq.(6) the cases where the two nuclei overlap are not excluded. To enforce the realistic condition of a peripheral collision it is necessary to perform the calculation taking into account the impact parameter dependence, $`b`$. It is straightforward to verify that in the collision of two identical nuclei the total cross section of Eq.(6) is modified to

$`{\displaystyle \frac{d^2\sigma _{AA}^{PPX}}{d^2b}}={\displaystyle \frac{Q^2}{2\pi }}e^{Q^2b^2/2}\sigma _{AA}^{PP},`$ (16)
where $`(Q^{})^2=(Q_0)^2+2r_0^2`$. The total cross section for inclusive processes is obtained after integration of Eq.(16) with the condition $`b_{min}>2R`$ in the case of identical ions.

For exclusive particle production the determination of the pomeron distribution function in the nuclei is slightly modified, because in this case it is necessary some specific assumption about the pomeron internal structure . Following Ref. the distribution function of pomerons is

$`f_P(x)={\displaystyle \frac{(3A\beta _0)^2}{(2\pi )^2x}}{\displaystyle _{\mathrm{}}^{(xM)^2}}𝑑t(tx^2M^2)F_A(t)^2|D(t)|^2,`$ (17)
and the cross section for a resonance production as a function of the impact parameter is

$`{\displaystyle \frac{d^2\sigma _{AA}^{PPR}}{d^2b}}`$ $`=`$ $`2\pi \left({\displaystyle \frac{3A\beta _0}{2\pi ^2}}\right)^4{\displaystyle \frac{dx_1}{x_1}\frac{dx_2}{x_2}Q_1^4Q_2^4\stackrel{~}{Q}^2e^{x_1^2M^2/Q_1^2}e^{x_2^2M^2/Q_2^2}}`$ (18)
$`\times `$ $`\left({\displaystyle \frac{x_1x_2s^2}{m^4}}\right)^{2\epsilon }\sigma _{AA}^{PPR}(x_1x_2s)b^2\stackrel{~}{Q}^2e^{b^2\stackrel{~}{Q}^2/2},`$ (19)
with $`\sigma _{AA}^{PPR}(x_1x_2s)`$ indicating the subprocess cross section (double pomeron fusion producing a resonance), and where

$`\stackrel{~}{Q}^2={\displaystyle \frac{1}{2}}(Q_1^2+Q_2^2),`$ (20)
with $`Q_i^2Q_0^2+2r_0^2`$ for idential ions. In the calculations we are going to perform we noticed that the aproximation $`Q_i^2Q_0^2`$ is quite reasonable, because for the energies that we shall consider the pomeron range parameter (Eq.(10)) is smaller than the width of the Gaussian form factor and consequently $`\stackrel{~}{Q}^2Q_0^2`$. Therefore, we obtain the final expression

$`{\displaystyle \frac{d^2\sigma _{AA}^{PPR}}{d^2b}}`$ $`=`$ $`2\pi \left({\displaystyle \frac{3A\beta _0Q_0^2}{2\pi ^2}}\right)^4{\displaystyle \frac{dx_1}{x_1}\frac{dx_2}{x_2}e^{x_1^2M^2/Q_0^2}e^{x_2^2M^2/Q_0^2}}`$ (21)
$`\times `$ $`\left({\displaystyle \frac{x_1x_2s^2}{m^4}}\right)^{2\epsilon }\sigma _{AA}^{PPR}(x_1x_2s)b^2Q_0^4e^{b^2Q_0^2/2}.`$ (22)
As discussed previously, to enforce the condition of peripheral collisions we integrate Eq.(22) with the condition $`b_{min}>2R`$.

Another way of to exclude events due to inelastic central collisions is through the introduction of an absortion factor computed in the Glauber aproximation . This factor modifies the cross section in the following form

$`{\displaystyle \frac{d\sigma _{AA}^{gl}}{d^2b}}`$ $`=`$ $`{\displaystyle \frac{d\sigma _{AA}^{PPR}}{d^2b}}\mathrm{exp}\left[A^2b\sigma _0{\displaystyle \frac{dQ^2}{(2\pi )^2}F_A^2(Q^2)e^{iQb}}\right]`$ (23)
$`=`$ $`{\displaystyle \frac{d\sigma _{AA}^{PPR}}{d^2b}}\mathrm{exp}\left[A^2b\sigma _0{\displaystyle \frac{Q_0^2}{4\pi }}e^{Q_0^2b^2/4}\right],`$ (24)
where $`\sigma _0`$ is the nucleon-nucleon total cross section, whose value for the different energy domains that we shall consider is obtained directly from the fit of Ref.

$`\sigma _0=Xs^ϵ+Y_1s^{\eta _1}+Y_2s^{\eta _2},`$ (25)
with $`X=18.256`$, $`Y_1=60.19`$, $`Y_2=33.43`$, $`ϵ=0.34`$, $`\eta _1=0.34`$, $`\eta _2=0.55`$, $`F_A(Q^2)=e^{Q^2/2Q_0^2}`$ and we exemplified Eq.(24) for the case of resonance production, i.e., $`\sigma _{AA}^{PPR}`$ is the total cross section for the resonance production to be discussed in the next section. The integration in Eq.(24) is over all impact parameter space and in the last section we discuss the differences between the two approaches showed above for removing central collisions.

## III Subprocesses initiated by photons and pomerons

### A Resonances

The main motivation to study resonance production in peripheral heavy ion collisions is that the high photon luminosity will allow us to observe resonances that couple very weakly to photons. The simplicity of this calculation also enable us to test the methods for removing central collisions, as well as to check up to which degree the double pomeron exchange is or not a background for the two photon physics.

To estimate the production of single spin-zero resonances, we note that these states can be formed by photon-photon fusion with a coupling strength that is measured by their photonic width

$`\widehat{\sigma }_{\gamma \gamma R}={\displaystyle \frac{8\pi ^2}{M_R\widehat{s}}}\mathrm{\Gamma }_{R\gamma \gamma }\delta \left(\tau {\displaystyle \frac{M_R^2}{\widehat{s}}}\right),`$ (26)
where $`M_R`$ is the ressonance mass and $`\mathrm{\Gamma }_{R\gamma \gamma }`$ its decay width in two photons. Using this expression into Eq.(2) we obtain the total cross section for the production of pseudo-scalar mesons.

To compute the cross section of the subprocess $`PPR`$ we can use the pomeron model of Donnachie and Landshoff . In this model it is assumed that the pomeron couples to the quarks like a isoscalar photon . This means that the cross sections of $`PPX`$ subprocesses can be obtained from suitable modifications on the cross-section for $`\gamma \gamma X`$. Another aspect to be considered is that the pomeron-quark-quark vertex is not point-like, and when either or both of the two quark legs in this vertex goes far off shell the coupling is known to decrease. So the quark-pomeron coupling $`\beta _0`$ must be replaced by

$`\stackrel{~}{\beta }_0(q^2)={\displaystyle \frac{\beta _0\mu _0^2}{\mu _0^2+Q^2}},`$ (27)
where $`\mu _0^2=1.2`$ GeV2 is a mass scale characteristic of the pomeron, in the case of resonance production $`Q=M_R/2`$ measures how far one of the quark legs is off mass shell and $`M_R`$ is the resonance mass. Therefore, the process $`PPR`$ is totally similar to the one initiated by photons unless from an appropriate change of factors. The cross section we are looking for is obtained changing the fine-structure constant $`\alpha `$ by $`9\stackrel{~}{\beta }/16\pi ^2`$, where $`\stackrel{~}{\beta }`$ is giving by Eq.(27) and $`9=3^2`$ is a color factor, leading to

$`\sigma _{PP}^R={\displaystyle \frac{9}{2}}{\displaystyle \frac{\stackrel{~}{\beta }^4}{\alpha ^2}}{\displaystyle \frac{\mathrm{\Gamma }(R\gamma \gamma )}{M_R}}\delta (x_1x_2sM_R^2).`$ (28)
Using this expression in Eq.(22) the total cross section is equal to

$`\sigma _{AA}^{PPR}`$ $`=`$ $`{\displaystyle \frac{9\pi }{8}}{\displaystyle \frac{(\stackrel{~}{\beta }Q_0)^4}{\alpha ^2}}\left({\displaystyle \frac{3A\beta _0Q_0}{2\pi }}\right)^4{\displaystyle \frac{\mathrm{\Gamma }(R\gamma \gamma )}{M_R}}\left({\displaystyle \frac{M_R^2s}{m^4}}\right)^{2ϵ}{\displaystyle \frac{Q_0^4}{M_R^2}}`$ (29)
$`\times `$ $`{\displaystyle \frac{dx}{x}\mathrm{exp}\left[\left(\frac{M_R^2M}{sQ_0x}\right)^2\frac{(xM)^2}{Q_0^2}\right]_{b_{min}}^{\mathrm{}}𝑑b2\pi b^3e^{Q_0^2b^2/2}},`$ (30)
where $`b_{min}=2R`$.

### B Pion pair production

The continuous production of pion pairs ($`\pi ^+\pi ^{}`$) is also an interesting signal to be observed in peripheral heavy ion collisions, mostly because they are a background for glueball and other resonances decays. Here we discuss the subprocess cross sections for two photons, $`ZZ\gamma \gamma ZZ\pi ^+\pi ^{}`$, and two pomerons exchange $`ZZPPZZ\pi ^+\pi ^{}`$.

The cross section for pion pair production by two photons can be calculated approximately by using a low energy theorem derived from partially-conserved-axial-vector-current hypothesis and current algebra and is equal to

$`\sigma (\gamma \gamma \pi ^+\pi ^{}){\displaystyle \frac{2\pi \alpha ^2}{s}}\left(1{\displaystyle \frac{4m_\pi ^2}{s}}\right)^{(1/2)}\left[{\displaystyle \frac{m_V^4}{\left(\frac{1}{2}s+m_V^2\right)\left(\frac{1}{4}s+m_V^2\right)}}\right]^2,`$ (32)
where $`m_\pi `$ is the pion mass an $`s`$ its squared energy, $`m_V`$ is a free parameter, whose value that provides the best fit to the experimental data is $`m_V1.4`$ GeV. This expression shows a nice agreement with the experimental data . For large values of $`s`$ it deviates from the Brodsky and Lepage formula . However, since most of the photon distribution is concentrated in the small $`x`$ region, i.e., the photons carry a small fraction of the momentum of the incoming ion, the difference is negligible.

Using Eqs. (32) and (2) we obtain

$`\sigma (s)={\displaystyle \frac{2\pi \alpha ^2}{s}}{\displaystyle _{\tau _{min}}^1}{\displaystyle \frac{d\tau }{\tau }}\left(1{\displaystyle \frac{4m_\pi ^2}{s\tau }}\right)^{(1/2)}\left[{\displaystyle \frac{m_V^4}{\left(\frac{1}{2}s\tau +m_V^2\right)\left(\frac{1}{4}s\tau +m_V^2\right)}}\right]^2{\displaystyle \frac{dL}{d\tau }}.`$ (33)
In the case of double pomeron exchange producing a pion pair we use once again the Donnachie and Landshoff model for the pomeron, obtaining the cross section for $`PP\pi ^+\pi ^{}`$ from the photonic one changing $`\alpha ^29\stackrel{~}{\beta }_0^4/16\pi ^2`$ in $`\sigma (\gamma \gamma \pi ^+\pi ^{})`$, and the resulting expression replaces $`\sigma _{AA}^{PPR}(x_1x_2s)`$ in Eq.(22). The total cross section appears after we perform the integration in the parameter space representation of the following equation

$`{\displaystyle \frac{d^2\sigma _{AA}^{PP\pi ^+\pi ^{}}}{db^2}}`$ $`=`$ $`\left({\displaystyle \frac{\pi ^2}{8}}\right){\displaystyle \frac{9}{4}}{\displaystyle \frac{(\stackrel{~}{\beta }_0Q_0)^4}{s}}\left({\displaystyle \frac{3A\beta _0Q_0}{2\pi ^2}}\right)^4{\displaystyle \frac{dx_1}{x_1^2}\frac{dx_2}{x_2^2}e^{(x_1M)^2/Q_0^2}e^{(x_2M)^2/Q_0^2}}`$ (34)
$`\times `$ $`\left({\displaystyle \frac{x_1x_2s^2}{m^4}}\right)^{2\epsilon }\left(1{\displaystyle \frac{4m_\pi ^2}{x_1x_2s}}\right)^{1/2}\left[{\displaystyle \frac{m_V^4}{\left(\frac{x_1x_2s}{2}+m_V^2\right)\left(\frac{x_1x_2s}{4}+m_V^2\right)}}\right]^2`$ (35)
$`\times `$ $`{\displaystyle _{b_{min}}^{\mathrm{}}}𝑑b2\pi Q_0^4b^3e^{b^2Q_0^2/2}.`$ (36)
### C Multiple particle production

The elementary cross section for multiple-particle production via two photons fusion can be described by the parametrization

$`\sigma _{\gamma \gamma hadrons}=C_1\left({\displaystyle \frac{s}{s_0}}\right)^ϵ+C_2\left({\displaystyle \frac{s}{s_0}}\right)^\eta ,`$ (37)
where $`C_1=173`$ nbarn, $`C_2=519`$ nbarn, $`s_0=1`$ GeV2, $`ϵ=0.079`$ and $`\eta =0.4678`$. The total cross section comes out from Eq.(2).

Within the Donnachie and Landshoff model it is straightforward to see that with the above parametrization the differential cross section to produce a cluster of particles with mass $`M_X`$ through double pomeron exchange is

$`{\displaystyle \frac{d\sigma }{dM_X}}`$ $`=`$ $`{\displaystyle \frac{(3A\beta _0\stackrel{~}{\beta }_0\mu _0)^4}{(2\pi )^4R_N^4(16\pi ^2\alpha ^2)}}{\displaystyle \frac{1}{2M_X}}{\displaystyle \frac{ds^{}}{s^{}}\left[C_1\left(\frac{s^{}}{s_0}\right)^ϵ+C_2\left(\frac{s^{}}{s_0}\right)^\eta \right]}`$ (38)
$`\times `$ $`\mathrm{exp}\left[\left({\displaystyle \frac{s^{}MR_N}{s}}\right)^2\left({\displaystyle \frac{M_X^2MR_N}{s^{}}}\right)^2\right]{\displaystyle _{b_{min}}^{\mathrm{}}}𝑑bb{\displaystyle \frac{e^{b^2/2R_N^2}}{R_N^2}},`$ (39)
to obtain this expression we used the pomeron distribution function in the nucleus for inclusive process (Eq.(7)).

To be possible a comparison with the work of Engel et al. , we also make use of the Ter-Martirosyan  model for diffractive multiparticle production. In this model the subprocess $`PPhadrons`$ is characterized by the cross section

$`\sigma _{PP}^{tot}(\mathrm{ln}(M_X^2/m^2),t_1,t_2)8\pi r(t_1)r(t_2),`$ (40)
which is a function of the triple-pomeron vertex $`r(t)`$, where t is the exchanged momentum. Using the value of $`r(0)`$ from Ref. , $`\sigma _{𝒫𝒫}^{tot}=8\pi r^2(0)140`$ $`\mu `$ barn. Note that we have clear differences between the approaches described above. Eq.(37) is a parametrization valid for a wide range of momenta, and with this one we naively apply the model of Ref.  to compute the total cross section for multiparticle production. On the other hand Eq.(40) is obtained in another specific model and it is not expected to be valid for the same range of energies as Eq.(37). This difference is going to be discussed in the last section.

Streng  applied the model of Ref.  for proton-proton collisions where the initial protons are scattered almost elastically, emerging with a very large fraction of the initial energy,

$`|x_1|,|x_2|c,c0.9.`$ (41)
The double pomeron exchange produce a particle cluster within a large rapidity gap and with a mass of the order

$`M_X^2s(1|x_1|)(1|x_2|),`$ (42)
where $`s`$ is the reaction energy squared. As the scattering is almost elastic, i.e., the emerging beam has approximately the same energy as the incident one, the following kinematical boundaries can be introduced

$`M_0M_X(1c)\sqrt{s},`$ (43)
$`{\displaystyle \frac{M_X^2}{(1c)}}s_1(1c)s,`$ (44)
where $`M_0=2`$ GeV and $`c=0.9`$. These limits have been translated for the case of heavy ions by Engel et al.  and we will proceed as them. If we consider Eq.(40), dress it with the pomeron distribution functions within the nuclei, and subtract the central collisions considering the absortion factor computed in the Glauber aproximation  we reproduce the results of Ref. .

## IV Results and Conclusions

Peripheral collisions at relativistic heavy ion colliders provide an arena for interesting studies of hadronic physics. Resonances coupling weakly to photons can be studied due to the large photon luminosity, the continuous production of pion pairs will be observed not only as a reaction of interest as well as a possible background for some resonance decays. A hadron cluster produced within a large rapidity gap will give information about photon-photon and double pomeron exchange. In this work we estimate the cross section for these processes. One of the main points is to verify if the double pomeron exchange is or not a background for the purely electromagnetic process. We discussed double pomeron exchange according to the Donnachie and Landshoff  model and calculated the cross sections in the impact parameter space. The condition for a realistic peripheral collision is imposed integrating the cross section with $`b_{min}>2R`$ in the case of two identical ions with radius $`R`$.

We considered the production of pseudoscalars resonances in the collision of 238U for energies available at RHIC ($`\sqrt{s}=200`$ GeV/nucleon), and collisions of 206Pb at energies available in LHC ($`\sqrt{s}=6.300`$ GeV/nucleon). Our results are shown in Table 1. Contrarily to the result of Ref.  the double pomeron exchange is not important when the cut in the impact parameter is introduced. For a realistic peripheral collision in the case of resonance production the pomeron-pomeron process is at least two orders of magnitude below the photon-photon one. Note in Table 1 that the rate of diffractive resonance production decreases with the increase of the meson mass. The main reason for this behavior lies in the fast decrease of the pomeron-quark coupling as shown in Eq.(27).

Note that the results of this table assume $`100\%`$ of efficiency in tagging the peripheral collision, even if we consider a small efficiency we recall that the cross section for light resonances imply in approximately billions of events/yr which easily survive the cuts for the background separation proposed by Nystrand and Klein (see the last paper of Ref. ). One of the most important cuts to separate inelastic nuclear reactions, associated with grazing collisions, is the small multiplicity of the final state, and this is exactly what we may expect in the final state of the particles discussed in Table 1.

The decays of $`\pi ^0`$, $`\eta `$, etc… will be dominated by two (or three) body decays in the central region of rapidity, and easily separated from the larger multiplicity common to inelastic collisions. It is interesting that in the case of $`\pi ^0`$ and $`\eta `$ production we may focus on the $`2\gamma `$ decay, and even if it is possible to separate the background from inelastic nuclear reactions, we still have the background of the photon-photon scattering through the QED box diagram producing the same final state. The box diagram will be dominated by light quarks, electron and muon, and for these we can use the asymptotic expression of $`\gamma \gamma `$ scattering ($`\sigma (s)20/s`$). Integrating this expression in a bin of energy (proportional to the resonance partial width into two-photons) centered at the mass of the resonance, we obtain a cross-section smaller than the resonant one with subsequent decay into two-photons. We do not considered the interference between the box and resonant diagram because on resonance the two processes are out of phase. It is opportune to mention that the decay products will fill the central region of rapidity, which is also one of the conditions proposed in Ref.  to isolate the peripheral collisions.

As discussed in Sect. 2 we have two different ways to enforce a realistic peripheral heavy ion collision. One is a geometrical cut in the impact parameter space where $`b_{min}>2R`$ is imposed, the other is through the introduction of the absorption factor in the Glauber approximation as given by Eq.(24). In Table 2 we compare the ratios between the total cross sections for diffractive resonance production computed with Eq.(24) and the one with the cut on the impact parameter (given by Eq.(LABEL:sigmaex)), in the collision of 238U for energies available at RHIC ($`\sqrt{s}=200`$ GeV/nucleon), and collisions of 206Pb at energies available in LHC ($`\sqrt{s}=6.300`$ GeV/nucleon).

The values of Table 2 show that the geometrical cut is less restrictive than the one given by the Glauber absorption factor. However, which one is more realistic also depends on the energy and on the ion that we are considering. In Table 3 we present the cross section for $`\pi ^0`$ production for different ions and at different energies. From Eq.(24) we notice that small variations in $`\sigma _0`$ (the nucleon-nucleon total cross section) are also promptly transmitted to the total cross section, and modify the ratios between the different methods to exclude inelastic collisions. Table 3 shows that the difference between the methods also become less important for light ions, but the most striking fact in this table is that for light ions double pomeron exchange starts becoming a background for photon-photon processes! According to Table 3 for 28Si the diffractive $`\pi ^0`$ production is a factor of 2 down the electromagnetic one (assuming the geometrical cut). This is not surprising because we know that for proton-proton the double pomeron exchange process should be larger than the electromagnetic one for producing a light resonance.

In Table 4 we show the pion pair cross section for different ions. The values were obtained using the geometrical cut, and even if with this procedure the diffractive cross section is a little bit overestimated for heavy ions we verify that photon-photon dominates. For light ions the diffractive process is already of the order of $`10\%`$ of the electromagnetic one.

The simulations discussed in the last paper of Ref.  have shown that the $`\gamma \gamma `$ interactions produce final states with small summed transverse momentum $`(|\overline{p}_T|)`$. Therefore, a cut of $`|\overline{p}_T|40\mathrm{\hspace{0.17em}\hspace{0.17em}100}MeV/c`$ can reduce considerably the background of non-peripheral collisions. In Table 4 we present the cross section for pion pair production through double photon interaction with $`|\overline{p}_T|100MeV/c`$. With this cut the cross section was reduced almost by a factor of $`4`$. The electromagnetic process with the restriction on $`p_T`$ is still larger than those of double pomeron exchange without this cut, and the introduction of this cut in the diffractive process produces a similar reduction.

The results for a hadron cluster production with invariant mass $`M_X`$ is depicted in Fig.1. In the figure it is shown the cross section for four different ions (Pb, Au, Ag, Ca) at energies that will be available at RHIC and LHC. The results were obtained integrating the cross sections with the condition $`b_{min}>2R`$. At LHC the photon-photon process will dominate the cross sections for heavy ions, whereas for light ions and small invariant mass they become of the same order. For heavy ions the diffractive process is indeed negligible. Note that our result for photons is similar to the one of Engel et al. , but the diffractive cross sections is slightly smaller than the one of Ref. . We credit this deviation to the differences in our approachs to calculate the subprocess cross section, mainly in the use of Eq.(37) with the changes prescribed by the Donnachie and Landshoff  instead of Eq.(40) given by the Ter-Martirosyan  model. They also use a value for $`\sigma _0`$ that is smaller than the one we considered here, which gives a smaller cut of the central collisions. We believe that the use of Eq.(37) and the model of Ref.  is more appropriate for the full range of momenta. Actually, diffractive models are plagued by uncertainties and the measurement of the double pomeron exchange in heavy ion colliders will provide useful information to distinguish between different models.

For multiple particle production we will not have the criteria of low multiplicity to help us to select the truly peripheral collisions, as well as it is far from clear if the cut in transverse momentum will be very effective to select the $`\gamma \gamma `$ events. However, we can separate the peripheral events on the basis of a clustering in the central region of rapidity, although an extensive and detailed simulation of the background processes will be necessary in order to set the precise interval of rapidity needed to cut the inelastic nuclear collisions.

As verified by Drees, Ellis and Zeppenfeld , Eq.(13) is a reasonable approximation for the form factor obtained from a Fermi or Woods-Saxon density distribution. However, their result shows that for heavy final states the photon-photon luminosity is slightly underestimated, and we can expect the same for the Pomeron one. A simple form factor expression consistent with the Fermi distribution has been recently obtained in Ref. , and its use would yield a few percent larger cross section in the case of a very heavy hadron cluster production.

In the case of peripheral heavy ion collisions at RHIC we surely cannot neglect the diffractive contribution, for light ions and a hadron cluster with low invariant mass it surely dominates photon-photon collisions. Notice that these results may change if we use the Glauber absorption factor to compute the cross section (depending on the energy, the ion and invariant mass), but the actual fact is that double pomeron exchange cannot be neglected at RHIC.

In conclusion, we estimated the production of resonances, pion pairs and a cluster of hadrons with invariant mass $`M_X`$ in peripheral heavy ion collisions at energies that will be available at RHIC and LHC. The condition for a realistic peripheral collision was studied with the use of a geometrical cut, where the minimum impact parameter was forced to be larger than twice the identical nuclei radius. The introduction of an absorption factor in the Glauber approximation to eliminate central collisions was also studied. We find out that the most restrictive method to account for inelastic collisions depended on the energy, the ion, as well as on the value of $`\sigma _0`$ (the nucleon-nucleon total cross section). The geometrical cut is not allways the most restrictive way to enforce peripheral collisions, an this is a topic that should be answered by the future experiments. In both cases we noticed that at energies of the LHC operating in the heavy ion mode and for very heavy ions the double pomeron exchange is not a background for the two photon process. The situation changes considerably for light ions and mostly for the energies available at RHIC, where double pomeron exchange cannot be neglected.

## Acknowledgments

One of us (C.G.R.) thanks Paulo S. R. da Silva for useful discussions, and C. A. Bertulani for a helpful remark. This research was supported in part by the Conselho Nacional de Desenvolvimento Cientifico e Tecnologico (CNPq) (AAN), Fundação de Amparo a Pesquisa do Estado de São Paulo (FAPESP) (CGR and AAN), and by Programa de Apoio a Núcleos de Excelência (PRONEX).


format: text

text:
# Exact and Approximate Dynamics of the Quantum Mechanical 𝑂(𝑁) Model

## I Introduction

Initial value problems in quantum field theory are of great interest in areas such as heavy ion collisions, dynamics of phase transitions, and early Universe physics. However, the solution of the corresponding functional Schrödinger equation is essentially impossible and one is forced to resort to approximate methods such as mean field approaches of the Hartree type or the large $`N`$ expansion. The application of variational techniques such as Hartree is limited in scope since the errors are uncontrolled. While $`1/N`$ methods promise better error control since they are based on a systematic expansion, at next-to-leading order these methods can become extremely complicated and expensive to implement. The motivation for our work in this paper is to implement the $`1/N`$ expansion at the first nontrivial order in a quantum mechanical example. Not only does this simplify the analysis but it also opens the possibility of comparing the approximate results with numerical simulations of the time dependent Schrödinger equation, a luxury not available in the field theoretic case. However, it should be kept in mind that quantum mechanics and quantum field theory are very different. For example, in the quantum mechanics applications discussed below, the $`O(1/N)`$ corrections do not correspond to inter-particle collisions (as they do in field theory) since we are restricting ourselves to one-particle quantum mechanics. Nevertheless, as discussed in more detail below, quantum mechanical examples provide excellent test-beds for key issues such as positivity violation and late-time accuracy of the approximations.

The $`O(N)`$ model has been extensively employed in time independent applications in statistical physics and quantum field theory  and several recent applications have studied time-dependent phenomena. The dynamics of the chiral phase transition following the expansion of a quark-gluon plasma produced during a relativistic heavy ion collision has been modeled by an $`O(4)`$ $`\sigma `$-model at leading order in $`1/N`$. The nonequilibrium dynamics of an $`O(N)`$-symmetric $`\lambda \varphi ^4`$ theory, again treated at leading order, has been investigated in detail . Even at leading order, the $`1/N`$ expansion captures the phase transition, but does not contain enough of the dynamics to allow for rethermalization, since direct scattering first occurs at next order. The $`O(N)`$ model has been used in inflationary models of the early Universe  with the scalar field often starting at the top of a hill in the potential and “rolling” down, giving rise to a quantum roll problem. It has also been applied to study primordial perturbations arising from defect models of structure formation .

The general method for obtaining the dynamical $`1/N`$ approximation via path integral techniques in quantum field theory was discussed earlier in Ref.  and applied later  to a quantum mechanical system of $`N+1`$ coupled oscillators (a one-dimensional version of scalar electrodynamics). Two different sets of approximate actions were considered, which differed by terms of order $`1/N^2`$, both of them being energy conserving. The first method in Ref.  was a perturbative expansion of the generating functional in powers of $`1/N`$. The second method was to first Legendre transform the action to order $`1/N`$, and then find the equations of motion. When these two methods diverged from each other, they also diverged from an exact solution for the case $`N=1`$. However, due to computational restrictions it was not possible to study numerically the accuracy of the approximation as a function of $`N`$. Remedying that deficiency is the main motivation of the present study, since for the quantum roll problem numerical solutions can be obtained for arbitrary $`N`$.

One of the subtle issues in expansions involving moment based truncation schemes such as $`1/N`$, which is present both in quantum field theory and in quantum mechanics, relates to the imposition of constraints arising from the positivity of the underlying probability density function or functional. The importance of these constraints is well known in areas such as turbulence and beam dynamics . In this paper, we show that possible violations of these constraints must be tamed in $`1/N`$ expansions if the approximation is at all expected to succeed at moderate values of $`N`$. This is possible by using certain resummation schemes which we will discuss elsewhere.

In this paper we show that the naive next-to-leading order $`1/N`$ expansion violates unitarity (or more generally, positivity) leading to an instability at least for $`N`$ less than some value $`N_T`$. We have numerical evidence for a sharp threshold at $`NN_T`$ beyond which we have not been able to detect an instability. This behavior appears to be related to the nature of the effective potential at next-to-leading order: At this order, the effective potential has the property of not being defined everywhere for values of $`N<N_c`$, where $`N_c`$ depends on the values of the parameters specifying the potential, and $`N_TN_c`$. However, for $`N>N_c`$, the effective potential is defined everywhere.

A comparison of the $`1/N`$ expansion and Hartree is of interest since both agree at infinite N. At finite $`N`$, the next-to-leading order large N and Hartree approximations differ and provide alternative routes to improving the leading order result which, for the quantum roll problem, consists of harmonic oscillations in $`r^2`$ where $`r`$ is the radial degree of freedom. At finite values of $`N`$, the inclusion of nonlinearities leads to amplitude modulation effects on top of the harmonic motion. The ability to capture this modulation is a good test for the next-to-leading order large N and Hartree approximations. Our numerical results provide evidence that neither of these methods are satisfactory at late times (relative to the oscillation time), though they work reasonably well at short to intermediate times.

Our results suggest that it is important to find ways to improve the naive $`1/N`$ expansion at next-to-leading order. Work using resummation schemes is in progress and short discussions of relevant issues are included in this paper.

The paper is organized as follows. In Section II we present the $`O(N)`$ model as it pertains to quantum mechanics and in Section III we derive equations of motion for the large-$`N`$ approximation to order $`1/N`$. We derive the corresponding equations for the time dependent Hartree approximation (TDHA) in Section IV. In section V we show how the same TDHA equations can be obtained from an equal-time Green’s function approach which is computationally more attractive. The energies for the various approximations are calculated in Section VI. Section VII describes the two initial conditions which preserve the $`O(N)`$ symmetry, namely a quantum roll, and the time evolution of an offset Gaussian centered at an $`O(N)`$ symmetric point. In Section VIII we determine the effective potential to both order $`1/N`$ and for the Hartree approximation. Numerical results and comparisons with the approximations are discussed in Section IX and our conclusions are discussed in Section X.

## II The $`O(N)`$ model

The Lagrangian for the $`O(N)`$ model in quantum mechanics is given by:

$$L(x,\dot{x})=\frac{1}{2}\underset{i=1}{\overset{N}{}}\dot{x}_i^2V(r),$$
(1)
where $`V(x)`$ is a potential of the form

$$V(r)=\frac{g}{8N}\left(r^2r_0^2\right)^2,r^2=\underset{i=1}{\overset{N}{}}x_i^2.$$
(2)
The time-dependent Schrödinger equation for this problem is given by:

$$i\frac{\psi (x,t)}{t}=\left\{\frac{1}{2}\underset{i=1}{\overset{N}{}}\frac{^2}{x_i^2}+V(r)\right\}\psi (x,t).$$
(3)
For arbitrary initial conditions, given present computational constraints, these equations can be numerically integrated only for small $`N4`$. The initial conditions for the quantum roll problem allow a numerical solution for all $`N`$, and in this case we can attempt to study fully the behavior of the large-$`N`$ expansion. (For the shifted Gaussian initial conditions, however, this is not possible, and we used numerical solutions obtained for $`N=1`$ and $`2`$ to benchmark the large-$`N`$ approximations and the TDHA solutions at short times.)

The symmetry of the quantum roll problem is such that only the radial part of the wave function is of interest. Assuming a solution of the form

$$\psi (r,t)=r^{(1N)/2}\varphi (r,t),$$
(4)
the time dependent Schrödinger equation for $`\varphi (r,t)`$ reduces to:

$$i\frac{\varphi (r,t)}{t}=\left\{\frac{1}{2}\frac{^2}{r^2}+U(r)\right\}\varphi (r,t)$$
(5)
with an effective one dimensional potential $`U(r)`$ given by

$$U(r)=\frac{(N1)(N3)}{8r^2}+\frac{g}{8N}\left(r^2r_0^2\right)^2.$$
(6)
It is further useful to make the rescaling:

$$r^2=Ny^2,r_0^2=Ny_0^2.$$
(7)
The potential (6) then becomes:

$$u(y,N)=\frac{U(y)}{N}=\frac{(N1)(N3)}{8N^2y^2}+\frac{g}{8}(y^2y_0^2)^2,$$
(8)
corresponding to the new Schrödinger equation,

$$i\frac{\varphi (y,\stackrel{~}{t})}{\stackrel{~}{t}}=\left\{\frac{1}{2N^2}\frac{^2}{y^2}+u(y,N)\right\}\varphi (y,\stackrel{~}{t})$$
(9)
where $`\stackrel{~}{t}=Nt`$.

The method of choice to investigate the long-time behavior of the exact solution is the split-operator method, which has been presented in detail in Ref. . The wave function is expanded as a Fourier series in the radial component, and the solution is obtained as the repeated application of a time-evolution operator in symmetrically split form. As a result, the use of a Fast-Fourier Transform algorithm is required. For the purpose of the present implementation, 256 radial grid points, a value of 20 for the radial grid boundary, and a time step size of 0.01, provides a conservation of the wave function unitarity to better than 9 significant figures. The accuracy of the method has been established by comparing results with a second method, where we first solve for the eigenvalues and eigenfunctions, and then use the expansion

$$\varphi (r,t)=\underset{n}{}C_ne^{iE_nt}\varphi _n(r),$$
(10)
where $`C_n`$ was determined from the initial conditions. This method was restricted to moderate values of $`N`$. Results from the two methods agreed in the cases where they were used together.

## III The Large-$`N`$ approximation

The large $`N`$ approximation has been worked out for the $`O(N)`$ model in $`1+3`$ dimensions in Ref. . The Lagrangian (1) with the potential function (2) is obtained from that paper by specializing to $`0+1`$ dimensions, and replacing $`\varphi _a(t)x_i(t)`$, $`vr_0`$, and $`\lambda g`$.

To implement the large $`N`$ expansion, it is useful to rewrite the Lagrangian in terms of the composite field $`\chi `$ by adding a constraint term to (1), given by:

$$\frac{N}{2g}\left[\chi \frac{g}{2N}(r^2r_0^2)\right]^2,$$
(11)
which yields an equivalent Lagrangian,

$$L^{}(x,\dot{x},\chi )=\underset{i}{}\frac{1}{2}\left(\dot{x}_i^2\chi x_i^2\right)+\frac{r_0^2}{2}\chi +\frac{N}{2g}\chi ^2.$$
(12)
The generating function $`Z[j,J]`$ is given by the path integral over the classical fields $`x_i(t)`$:

$`Z[j,J]=e^{iW[j,J]}`$ $`=`$ $`{\displaystyle d\chi \underset{i}{}\mathrm{d}x_i\mathrm{exp}\left\{iS[x,\chi ;j,J]\right\}},`$
$`S[x,\chi ;j,J]`$ $`=`$ $`{\displaystyle _𝒞}dt\left\{L^{}+{\displaystyle \underset{i}{}}j_ix_i+J\chi \right\}.`$
The effective action, to order $`1/N`$, is obtained by integrating the path integral for the generating functional for the Lagrangian (12), over the $`x_i`$ variables, and approximating the integral over $`\chi `$ by the method of steepest descent (keeping terms up to order $`1/N`$). A Legendre transform of the resulting generating functional then yields the effective action, which we find to be:

$`\mathrm{\Gamma }[q,\chi ]=`$ (15)
$`{\displaystyle _𝒞}\mathrm{d}t\{{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}[\dot{q}_i^2(t)\chi (t)q_i^2(t)]+{\displaystyle \frac{i}{2}}{\displaystyle \underset{i}{}}\mathrm{ln}[G_{ii}^1(t,t)]`$
$`+{\displaystyle \frac{r_0^2}{2}}\chi (t)+{\displaystyle \frac{N}{2g}}\chi ^2(t)+{\displaystyle \frac{i}{2}}\mathrm{ln}[D^1(t,t)]\},`$
where the integral is over the close time path $`𝒞`$, discussed in Ref.  and $`q(t)=x_i(t)`$. Here $`G_{ij}^1(t,t^{})`$ and $`D^1(t,t^{})`$ are the lowest order in $`1/N`$ inverse propagators for $`x_i`$ and $`\chi `$, given by

$`G_{ij}^1(t,t^{})`$ $`=`$ $`\left\{{\displaystyle \frac{\mathrm{d}^2}{\mathrm{d}t^2}}+\chi (t)\right\}\delta _𝒞(t,t^{})\delta _{ij}G^1(t,t^{})\delta _{ij},`$
$`D^1(t,t^{})`$ $`=`$ $`{\displaystyle \frac{N}{g}}\delta _𝒞(t,t^{})\mathrm{\Pi }(t,t^{}),`$
where

$`\mathrm{\Pi }(t,t^{})`$ $`=`$ $`{\displaystyle \frac{i}{2}}{\displaystyle \underset{i,j}{}}G_{ij}(t,t^{})G_{ji}(t^{},t)`$ (17)
$`+{\displaystyle \underset{i,j}{}}q_i(t)G_{ij}(t,t^{})q_j(t^{}).`$
Here $`\delta _𝒞(t,t^{})`$ is the closed time path delta function.

The equations of motion for the classical fields $`q_i(t)`$, to order $`1/N`$, are

$`\left\{{\displaystyle \frac{\mathrm{d}^2}{\mathrm{d}t^2}}+\chi (t)\right\}q_i(t)`$ (18)
$`+i{\displaystyle \underset{j}{}}{\displaystyle _𝒞}dt^{}G_{ij}(t,t^{})D(t,t^{})q_j(t^{})=0,`$ (19)
with the gap equation for $`\chi (t)`$ given by

$$\chi (t)=\frac{g}{2N}r_0^2+\frac{g}{2N}\underset{i}{}\left[q_i^2(t)+\frac{1}{i}𝒢_{ii}^{(2)}(t,t)\right].$$
(20)
The next-to-leading order $`x_i`$ propagator $`𝒢_{ij}^{(2)}(t,t^{})`$ and self energy $`\mathrm{\Sigma }_{ij}(t,t^{})`$ to order $`1/N`$ turn out to be

$`𝒢_{ij}^{(2)}(t,t^{})=G_{ij}(t,t^{})`$ (21)
$`{\displaystyle \underset{k,l}{}}{\displaystyle _𝒞}dt_1{\displaystyle _𝒞}dt_2G_{ik}(t,t_1)\mathrm{\Sigma }_{kl}(t_1,t_2)G_{lj}(t_2,t^{}),`$ (22)
$`\mathrm{\Sigma }_{kl}(t,t^{})=iG_{kl}(t,t^{})D(t,t^{})q_k(t)D(t,t^{})q_l(t^{}).`$ (23)
These equations agree with (2.18–2.22) of Ref. . We mention here that the actual equation for $`𝒢`$ which follows from the effective action differs from Eq. (21) in that the final $`G`$ in the integral equation is replaced by the full $`𝒢`$. This leads to a partial resummation of the $`1/N`$ corrections which which guarantees positivity of $`x^2(t)`$ (this restricted result does not imply that the full positivity problem for the density matrix has been solved). However, it does not improve the long-time accuracy of the results .

In order to solve for $`D(t,t^{})`$, we first write

$$\frac{N}{g}D(t,t^{})=\delta _𝒞(t,t^{})+\frac{N}{g}\mathrm{\Delta }D(t,t^{}).$$
(24)
Then $`\mathrm{\Delta }D(t,t^{})`$ satisfies the integral equation,

$$\frac{N}{g}\mathrm{\Delta }D(t,t^{})=\frac{g}{N}\mathrm{\Pi }(t,t^{})_𝒞dt^{\prime \prime }\mathrm{\Pi }(t,t^{\prime \prime })\mathrm{\Delta }D(t^{\prime \prime },t^{}),$$
(25)
in agreement with (2.13–2.16) of Ref. .

We are now in a position to solve these coupled equations for the motion of $`q_i(t)`$ and $`\chi (t)`$ for given initial conditions. For the initial conditions discussed in Sec. VII, we find:

$$G_{ij}(t,t^{})/i=\delta _{ij}f(t)f^{}(t^{}),$$
(26)
where $`f(t)`$ and $`f^{}(t)`$ satisfy the homogeneous equation,

$$\left\{\frac{\mathrm{d}^2}{\mathrm{d}t^2}+\chi (t)\right\}\left(\begin{array}{c}f(t)\\ f^{}(t)\end{array}\right)=0,$$
(27)
with initial conditions:

$$f(0)=\sqrt{G},\dot{f}(0)=1/(2\sqrt{G}).$$
(28)
However, $`\mathrm{\Delta }D(t,t^{})`$ cannot be factored into products of functions like $`G_{ij}(t,t^{})`$.

We solve (19) and (20) simultaneously with (21) and (25), using the Chebyshev expansion technique of Appendices A and B of Ref. .

## IV The time dependent Hartree approximation

It is useful to compare our results for the large-$`N`$ approximation to the time dependent Hartree approximation (TDHA) suitably formulated for the $`O(N)`$ problem. The static Hartree approximation is based on the idea of varying the parameters of a Gaussian wave function so as to minimize the energy (the generalization to the time-dependent case is given below). For the $`O(N)`$ problem this amounts to placing an $`N`$-dimensional Gaussian some (radial) distance away from the origin and then carrying out the minimization procedure. In contrast, the leading-order large $`N`$ wave function is a Gaussian which is locked at the origin. At infinite $`N`$ the TDHA becomes exact and equivalent to the leading order large-$`N`$ approximation, a well-known result. At finite $`N`$, the TDHA and the next-to-leading order large $`N`$ approximation can be thought of as two competing schemes to improve on the leading-order result.

There are several ways of implementing the Hartree approximation: The most common is by using the time-dependent variational principle of Dirac . This has the advantage of giving a classical Hamiltonian description for the dynamics of the variational parameters, which can be hidden in other formulations.

The idea behind this approach is that the variation of

$`\mathrm{\Gamma }[\psi ,\psi ^{}]`$ $`=`$ $`{\displaystyle dt\psi (t)|i\frac{}{t}H|\psi (t)}`$ (29)
is stationary for the exact solution of the Schrödinger equation. We consider Gaussian trial wave functions of the form:

$`\psi (x,t)=𝒩\mathrm{exp}[ip_i(t)z_i(t)`$ (30)
$`z_i(t)({\displaystyle \frac{G_{ij}^1(t)}{4}}\mathrm{\Pi }_{ij}(t))z_j(t)],`$ (31)
where $`𝒩`$ is the normalization constant, and we have set $`z_i(t)=x_iq_i(t)`$. Here $`q_i(t)`$, $`p_i(t)`$, $`G_{ij}(t)`$ and $`\mathrm{\Pi }_{ij}(t)`$ are time-dependent variational parameters, to be determined by minimizing the Dirac action. We note that $`\mathrm{\Pi }_{ij}(t)`$, which is used only in this section, is conjugate to $`G_{ij}(t)`$ and is not to be confused with the self energy $`\mathrm{\Pi }(t,t^{})`$ defined in Eq. (17).

The $`n`$-point functions can be calculated from the generating functional using the formula,

$$z_iz_j\mathrm{}z_n=\frac{^nZ[j]}{j_ij_j\mathrm{}j_n}|_{j=0},$$
(32)
where

$`Z[j]`$ $`=`$ $`𝒩^2{\displaystyle \underset{s}{}\mathrm{d}x_s}`$ (34)
$`\times \mathrm{exp}\left[{\displaystyle \frac{1}{2}}z_i(t)G_{ij}^1z_j(t)+j_iz_i(t)\right]`$
$`=`$ $`\mathrm{exp}\left[{\displaystyle \frac{j_iG_{ij}j_j}{2}}\right].`$ (35)
The expectation value of the time derivative is given by

$`i{\displaystyle \frac{}{t}}=p_i\dot{q}_iG_{ij}\dot{\mathrm{\Pi }}_{ij}`$ (36)
and the expectation value of the kinetic energy is

$$\frac{1}{2}\frac{^2}{x_i^2}=\frac{p_ip_i}{2}+\frac{1}{8}G_{ii}^1+\mathrm{\hspace{0.17em}2}\mathrm{\Pi }_{ij}G_{jk}\mathrm{\Pi }_{ki}.$$
(37)
For the expectation value of $`V`$ we first expand the potential in a Taylor series about $`z_i=0`$,

$`V(q,z)=V(q)+V_i(q)z_i+{\displaystyle \frac{1}{2}}V_{ij}(q)z_iz_j+\mathrm{}`$

where

$`V_i(q)`$ $`=`$ $`{\displaystyle \frac{g}{2N}}q_i\left(q_sq_sr_0^2\right),`$ (38)
$`V_{ij}(q)`$ $`=`$ $`{\displaystyle \frac{g}{2N}}\left[\delta _{ij}\left(q_sq_sr_0^2\right)+2q_iq_j\right],`$ (39)
$`V_{ijk}(q)`$ $`=`$ $`{\displaystyle \frac{g}{N}}\left(\delta _{ij}q_k+\delta _{ik}q_j+\delta _{jk}q_i\right),`$ (40)
$`V_{ijkl}(q)`$ $`=`$ $`{\displaystyle \frac{g}{N}}\left(\delta _{ij}\delta _{kl}+\delta _{il}\delta _{jk}+\delta _{ik}\delta _{jl}\right).`$ (41)
Thus

$`V(q,z)={\displaystyle \frac{g}{8N}}[(q_jq_jr_0^2)^2+(z_iz_i)^2+4(z_iz_i)(z_jq_j)`$
$`+4(z_iq_i)^2+2(z_iz_i)(q_jq_jr_0^2)+4(z_iq_i)(q_jq_jr_0^2)]`$
Taking the expectation value, we obtain:

$`V`$ $`=`$ $`{\displaystyle \frac{g}{8N}}[(q_jq_jr_0^2)^2+2G_{ii}(q_jq_jr_0^2)`$ (43)
$`+4G_{ij}q_iq_j+G_{ii}G_{jj}+2G_{ij}G_{ji}]`$
The Hartree equations of motion are Hamilton’s equations for the variational parameters:

$`\dot{q}_i`$ $`=`$ $`p_i,`$ (44)
$`\dot{p}_i`$ $`=`$ $`V_i{\displaystyle \frac{1}{2}}V_{ijk}G_{jk},`$ (45)
$`\dot{G}_{ij}`$ $`=`$ $`2\left(G_{ik}\mathrm{\Pi }_{kj}+G_{jk}\mathrm{\Pi }_{ki}\right),`$ (46)
$`\dot{\mathrm{\Pi }}_{ij}`$ $`=`$ $`{\displaystyle \frac{1}{8}}G_{ik}^1G_{kj}^12\mathrm{\Pi }_{ik}\mathrm{\Pi }_{kj}{\displaystyle \frac{1}{2}}V_{ij}{\displaystyle \frac{1}{4}}V_{ijkl}G_{kl}.`$ (47)
Solutions of this set of equations determine the time-dependent Hartree approximation to the true solution of the Schrödinger equation.

## V Method of equal-time Green’s functions

Solutions of the TDHA equation (44) require computing the matrix inverse of $`G_{ij}`$. This can be difficult to carry out in practice for large $`N`$. Fortunately, the method of equal-time Green’s functions provides a way to avoid this technical difficulty . We begin by considering the time-evolution of the one point functions:

$`\dot{q}_i`$ $`=`$ $`p_i,`$ (48)
$`\dot{p}_i`$ $`=`$ $`V_i{\displaystyle \frac{1}{2}}V_{ijk}G_{jk},`$ (49)
$`=`$ $`{\displaystyle \frac{g}{2N}}\left\{q_i(q_kq_k+G_{kk}r_0^2)+q_k(G_{ik}+G_{ki})\right\}`$ (50)
where $`V_i`$ and $`V_{ijk}`$ are given by (38), as well as the evolution of the two-point functions:

$`G_{ij}(t)`$ $`=`$ $`z_iz_j,K_{ij}=\dot{z}_i\dot{z}_j,`$ (51)
$`F_{ij}`$ $`=`$ $`{\displaystyle \frac{1}{2}}[z_i\dot{z}_j+\dot{z}_jz_i].`$ (52)
Here we have again set $`z_i(t)=x_iq_i(t)`$. All of the expectation values are taken with respect to the Gaussian trial wave function, Eq. (30). To obtain the equations of motion for the two-point functions, we use the exact equation of motion and the factorization resulting from the Gaussian approximation. This yields

$`\dot{G}_{ij}`$ $`=`$ $`F_{ij}+F_{ji},`$ (53)
$`\dot{F}_{ij}`$ $`=`$ $`K_{ij}z_i{\displaystyle \frac{V}{x_j}},`$ (54)
$`\dot{K}_{ij}`$ $`=`$ $`\left[{\displaystyle \frac{V}{x_i}}\dot{z}_j\dot{z}_i{\displaystyle \frac{V}{x_j}}\right].`$ (55)
Here we have used the Lagrange equations of motion

$$\ddot{x}_i+\frac{V}{x_i}=0,$$
(56)
where

$`{\displaystyle \frac{V}{x_i}}`$ $`=`$ $`V_{ij}z_j+{\displaystyle \frac{1}{6}}V_{ijkl}z_jz_kz_l`$ (58)
$`+\text{terms with even powers of }z_i,`$
and the fact that for our Gaussian wave packet, $`z_i=0.`$ The canonical commutation relations give:

$$z_i\dot{z}_j\dot{z}_jz_i=x_i\dot{x}_j\dot{x}_jx_i=i\delta _{ij},$$
(59)
and we have

$`z_iz_jz_kz_l=G_{ij}G_{kl}+G_{il}G_{jk}+G_{ik}G_{jl},`$ (60)
$`z_iz_jz_k\dot{z}_l=G_{ij}F_{kl}+F_{il}G_{jk}+G_{ik}F_{jl}`$ (61)
$`+i\left(G_{ij}\delta _{kl}+G_{jk}\delta _{il}+G_{ik}\delta _{jl}\right)/2,`$ (62)
$`\dot{z}_iz_jz_kz_l=F_{ji}G_{kl}+F_{li}G_{jk}+F_{ki}G_{jl}`$ (63)
$`i\left(G_{kl}\delta _{ji}+G_{jk}\delta _{li}+G_{jl}\delta _{ki}\right)/2.`$ (64)
Finally, from Eqs. (54) and (55) we get:

$`\dot{F}_{ij}`$ $`=`$ $`K_{ij}V_{jk}G_{ik}`$ (66)
$`V_{jklm}\left(G_{ik}G_{lm}+G_{im}G_{kl}+G_{il}G_{km}\right)/6,`$
$`\dot{K}_{ij}`$ $`=`$ $`V_{ik}F_{kj}V_{jk}F_{ki}`$ (69)
$`V_{iklm}\left(G_{kl}F_{mj}+G_{lm}F_{kj}+G_{km}F_{lj}\right)/6`$
$`V_{jklm}\left(F_{ki}G_{lm}+F_{mi}G_{kl}+F_{li}G_{km}\right)/6.`$
For Gaussian initial conditions, the equal-time Green’s function method is assured to give the same result as the Hartree method, if $`F_{ij}(0)`$ and $`K_{ij}(0)`$ satisfy the requirements:

$`F_{ij}(0)=0,K_{ij}(0)=G_{ij}^1(0)/4.`$

Choosing $`K_{ij}`$ independently of $`G_{ij}`$, corresponds to a mixed initial density matrix, rather than a pure state. If we further choose $`G_{ij}(0)`$ to be diagonal and equal to the same number $`G_0`$,

$`G_{ij}(0)=\delta _{ij}G_0,`$

then $`K_{ij}(0)`$ is given by:

$`K_{ij}(0)=\delta _{ij}/(4G_0).`$

For the initial conditions pertinent to the quantum roll, $`q_i(t)=0`$ for all $`t`$. Then $`G_{ij}(t)`$, $`F_{ij}(t)`$, and $`K_{ij}(t)`$ are all proportional to the unit matrix, and have no off-diagonal terms. For the offset initial condition, we choose large-$`N`$ symmetric initial conditions so that $`q_i(0)=q_0`$, and $`p_i(0)=0`$, and $`G_{ij}(0)=G_0\delta _{ij}`$. In that case, all the $`q`$’s and $`p`$’s are identical,

$`q_i(t)=q(t),p_i(t)=p(t),`$

and the matrices $`G_{ij}(t)`$, $`F_{ij}(t)`$ and $`K_{ij}(t)`$ become off-diagonal in a simple way so that all the diagonal elements are equal and all the off-diagonal elements are equal. That is, we can write

$`G_{ij}(t)`$ $`=`$ $`G(t)\delta _{ij}+\overline{G}(t)(1\delta _{ij}),`$
$`F_{ij}(t)`$ $`=`$ $`F(t)\delta _{ij}+\overline{F}(t)(1\delta _{ij}),`$
$`K_{ij}(t)`$ $`=`$ $`K(t)\delta _{ij}+\overline{K}(t)(1\delta _{ij}).`$
For this case, Eqs. (48), (53), (66), and (69) simplify to the following set of coupled equations:

$`\dot{q}`$ $`=`$ $`p,`$ (70)
$`\dot{p}`$ $`=`$ $`{\displaystyle \frac{g}{2N}}q\left\{Nq^2r_0^2+(N+2)G+2(N1)\overline{G}\right\},`$ (71)
$`\dot{G}`$ $`=`$ $`2F,\dot{\overline{G}}=2\dot{\overline{F}},`$ (72)
$`\dot{F}`$ $`=`$ $`K{\displaystyle \frac{g}{2N}}\{G[(N+2)(q^2+G)r_0^2]`$ (74)
$`+2(N1)\overline{G}(q^2+\overline{G})\},`$
$`\dot{\overline{F}}`$ $`=`$ $`\overline{K}{\displaystyle \frac{g}{2N}}\{\overline{G}[(3N2)q^2+(N+4)G`$ (76)
$`+2(N2)\overline{G}r_0^2]+2Gq^2\},`$
$`\dot{K}`$ $`=`$ $`{\displaystyle \frac{g}{N}}\{F[(N+2)(q^2+G)r_0^2]+`$ (78)
$`+\overline{F}\mathrm{\hspace{0.17em}2}(N1)(q^2+\overline{G})\},`$
$`\dot{\overline{K}}`$ $`=`$ $`{\displaystyle \frac{g}{N}}\{\overline{F}[(3N2)q^2+(N+2)G`$ (80)
$`+2(N2)\overline{G}r_0^2]+F\mathrm{\hspace{0.17em}2}(q^2+\overline{G})\}.`$
If we let $`r_0^2=Ny_0^2`$, and then take the limit $`N\mathrm{}`$, we recover the leading order in the large-$`N`$ result, as discussed in Ref. .

The equal time Green’s function method is easier to implement numerically than the Hamiltonian system described by Eq. (44), since no matrix inversion is involved. However, if one wants to find the wave function or the energy, instead of just obtaining the Green’s functions, matrix inversion is once again required.

## VI Energy

It is important to note that even though the Hartree and large $`N`$ approximations are truncations of the true dynamics, they are nevertheless energy conserving. In the large-$`N`$ approximation, to order $`1/N`$, the expectation value of the Hamiltonian is given by

$`E`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}\left[\dot{x_i}^2(t)+\widehat{\chi }(t)x_i^2(t)\right]`$ (82)
$`{\displaystyle \frac{r_0^2}{2}}\widehat{\chi }(t){\displaystyle \frac{N}{2g}}\widehat{\chi }^2(t).`$
We write these expectation values in terms of the CTP Green’s functions. By definition, the disconnected two-point Green’s functions are introduced as

$`D_{\mathrm{dis}}(t,t^{})`$ $`=`$ $`i\mathrm{T}_𝒞[\widehat{\chi }(t)\widehat{\chi }(t^{})]`$ (83)
$`=`$ $`i\chi (t)\chi (t^{})+𝒟(t,t^{}),`$ (84)
$`G_{ij,\mathrm{dis}}(t,t^{})`$ $`=`$ $`i\mathrm{T}_𝒞[x_i(t)x_j(t^{})]`$ (85)
$`=`$ $`iq_i(t)q_j(t^{})+𝒢_{ij}(t,t^{}),`$ (86)
where $`𝒟`$ and $`𝒢_{ab}`$ denote the connected two-point Green’s functions

$`𝒟(t,t^{})`$ $`=`$ $`\left[{\displaystyle \frac{\delta ^2W[J,j]}{\delta J(t)\delta J(t^{})}}\right]_{J,j=0},`$ (87)
$`𝒢_{ij}(t,t^{})`$ $`=`$ $`\left[{\displaystyle \frac{\delta ^2W[J,j]}{\delta j_i(t)\delta j_j(t^{})}}\right]_{J,j=0}.`$ (88)
To obtain the energy from (82), we require the expectation values

$`\widehat{\chi }(t)`$ $`=`$ $`\chi (t),`$
$`\widehat{\chi }^2(t)`$ $`=`$ $`\chi ^2(t)+D(t,t)/i,`$
$`x_i^2(t)`$ $`=`$ $`q_i^2(t)+𝒢_{ii}(t,t)/i,`$
$`\dot{x}_i^2(t)`$ $`=`$ $`\dot{q}_i^2(t)+{\displaystyle \frac{^2𝒢_{ii}(t,t^{})/i}{tt^{}}}|_{t=t^{}},`$
$`\widehat{\chi }(t)x_i^2(t)`$ $`=`$ $`\chi (t)\left[q_i^2(t)+𝒢_{ii}(t,t)/i\right]K_{ii}(t,t,t),`$
where $`K_{ij}(t_1,t_2,t_3)`$ is the 3-point Green’s function defined as

$`K_{ij}(t_1,t_2,t_3)={\displaystyle _𝒞}dtG_{ik}(t_1,t)G_{kj}(t_2,t)D(t,t_3).`$
The energy for the next-to-leading order large-N approximation is then given by:

$`E`$ $`=`$ $`{\displaystyle \frac{r_0^2}{2}}\chi (t){\displaystyle \frac{N}{2g}}\left\{\chi ^2(t)+\mathrm{\Delta }D(t,t)/i\right\}`$ (91)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}\left\{\dot{q}_i^2(t)+{\displaystyle \frac{^2𝒢_{ii}(t,t^{})/i}{tt^{}}}|_{t=t^{}}\right\}`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}\left\{\chi (t)\left[q_i^2(t)+𝒢_{ii}(t,t)/i\right]K_{ii}(t,t,t)\right\}.`$
Using the equations of motion, we can show directly that (91) is conserved.

It is easy to evaluate the energy at $`t=0`$ for the quantum roll problem using the initial conditions (28). The result for the leading and next-to-leading order large-N approximation is:

$`{\displaystyle \frac{E}{N}}=ϵ_0+{\displaystyle \frac{1}{N}}ϵ_1+\mathrm{}`$

where,

$`ϵ_0`$ $`=`$ $`{\displaystyle \frac{1}{8G}}+{\displaystyle \frac{1}{8}}gy_0^4{\displaystyle \frac{1}{4}}gGy_0^2+{\displaystyle \frac{1}{8}}gG^2,`$ (92)
$`ϵ_1`$ $`=`$ $`{\displaystyle \frac{1}{4}}gG^2.`$ (93)
Our initial wave function was chosen to be Gaussian, so that the parameters of the Hartree approximation agree exactly with the energy and parameters of the exact wave function at $`t=0`$. However, the leading order in the large-$`N`$ approximation for the same value of $`G`$ will disagree with the exact energy by $`ϵ_1`$. This discrepency dissapears when we include the $`1/N`$ corrections.

Since the Hartree approximation leads to a canonical Hamiltonian dynamical system, the corresponding energy in that approximation is also a constant of the motion. It is given by:

$`E`$ $`=`$ $`{\displaystyle \frac{1}{2}}p_i^2+{\displaystyle \frac{1}{8}}G_{ii}^1+2\mathrm{\Pi }_{ij}G_{jk}\mathrm{\Pi }_{ki}`$ (96)
$`+V(q)+{\displaystyle \frac{1}{2}}V_{ij}G_{ij}`$
$`+{\displaystyle \frac{1}{4!}}V_{ijkl}\left(G_{ij}G_{kl}+G_{il}G_{jk}+G_{ik}G_{jl}\right).`$
We used this expression to check the accuracy of our numerical solutions. As with the next-to-leading order $`1/N`$ expression, for the quantum roll initial condition, Eq. (96) agrees with the exact result.

## VII Initial Conditions

### A Quantum Roll

We wish to study initial conditions which are consistent with $`O(N)`$ symmetry. This implies immediately that all the $`x_i(t)`$ have to be identical, with $`x_i(0)=0`$, and $`G_{ij}(t)`$ must be diagonal. The quantum roll problem is defined by a Gaussian initial wave function that is centered on the origin:

$$\psi _0(r)=\frac{1}{(2\pi G)^{N/4}}\mathrm{exp}\left\{\frac{r^2}{4G}\right\}.$$
(97)
In this section, $`GG(0)`$.

One of the difficulties in studying the systematics of the $`1/N`$ expansion is the fact that, at next to leading order, every different value of $`N`$ (with all other parameters held constant) defines a different initial value problem. In this sense one cannot naively compare individual solutions, exact or approximate, at different values of $`N`$. In effect one has to tune the parameters of the problem at each $`N`$ in order to maintain certain invariance properties which allow different $`N`$ evolutions to be compared to each other. This parameter tuning process is described below.

Since the infinite $`N`$ limit has very precise properties, several technical issues arise when one wants to approach this limit starting at $`N=1`$ in a uniform manner. To study the large $`N`$ limit it is convenient to make a rescaling to the $`y`$ variables, given in (7). At very large $`N`$, the potential energy $`u(y,N)`$ is (8):

$`u(y,N)`$ $`=`$ $`{\displaystyle \frac{(N1)(N3)}{8N^2y^2}}+{\displaystyle \frac{g}{8}}(y^2y_0^2)^2,`$ (98)
$``$ $`{\displaystyle \frac{1}{8y^2}}+{\displaystyle \frac{g}{8}}(y^2y_0^2)^2.`$ (99)
In this limit, $`u(N,y)`$ has a minimum which is independent of $`N`$, and the large $`N`$ limit consists of harmonic oscillations about this minimum \[the reason for this is that the large $`N`$ limit also corresponds to an effectively large mass limit in the Schrödinger equation (9)\].

One way to uniformly study the motion of a wave packet as a function of $`N`$ is to choose initial conditions so that there is a uniform overlap of the initial wave function with the set of eigenfunctions of the Hamiltonian in the $`N\mathrm{}`$ limit. We can obtain this constant overlap if we allow the coupling constant $`g`$ to be a slowly varying function of $`N`$. This can be done in several ways that differ by terms of order $`1/N^2`$. The method presented below leads to uniform results even at $`N=1`$ as we change the parameters with $`N`$. Our method is to keep the distance between the centers of the initial wave function and the position of the minimum of the potential a constant as $`N`$ is varied.

Using (97), we define $`\stackrel{~}{r}`$ by:

$$\stackrel{~}{r}=r=\sqrt{2G}\left[\frac{\mathrm{\Gamma }((N+1)/2)}{\mathrm{\Gamma }(N/2)}\right],$$
(100)
and $`G_{\mathrm{}}`$ by the variance:

$`{\displaystyle \frac{G_{\mathrm{}}}{2}}=r^2r^2=G\left\{N2\left[{\displaystyle \frac{\mathrm{\Gamma }((N+1)/2)}{\mathrm{\Gamma }(N/2)}}\right]^2\right\}.`$

Solving the above equations for $`G`$, we have

$$G(N)=\frac{G_{\mathrm{}}}{2N4\left[\mathrm{\Gamma }((N+1)/2)/\mathrm{\Gamma }(N/2)\right]^2}.$$
(101)
Substitution of this expression into (100) yields

$$\stackrel{~}{r}(N)=\sqrt{\frac{G_{\mathrm{}}}{N[\mathrm{\Gamma }(N/2)/\mathrm{\Gamma }((N+1)/2)]^22}}.$$
(102)
In the limit when $`N`$ goes to infinity, we have

$`G(N)G_{\mathrm{}},\stackrel{~}{r}(N)\stackrel{~}{r}_{\mathrm{}}=\sqrt{(N1)G_{\mathrm{}}},`$

which defines $`\stackrel{~}{r}_{\mathrm{}}`$, and agrees with the asymptotic form of the rescaled version of the initial wave function (97):

$`\varphi _0(r)`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi G)^{N/4}}}\mathrm{exp}\left\{{\displaystyle \frac{r^2}{4G}}+{\displaystyle \frac{N1}{2}}\mathrm{ln}r\right\},`$
$``$ $`{\displaystyle \frac{1}{(2\pi G)^{N/4}}}\mathrm{exp}\left\{{\displaystyle \frac{(r\stackrel{~}{r}_{\mathrm{}})^2}{2G_{\mathrm{}}}}+O(1/\sqrt{N})\right\}.`$
In order to ensure that the initial wave function has a finite overlap with the energy eigenfunctions of the Schrödinger equation at large $`N`$, we will keep the value of $`G_{\mathrm{}}`$ (and not $`G`$) fixed in our simulations.

Another quantity that should be kept constant is the basic oscillation frequency. In order to do this, we first find the Gaussian oscillations about the minimum of the one dimensional potential, defined by Eq. (6). We expand $`U(r)`$ as

$$U(r)=U(\overline{r})+\frac{1}{2}\overline{m}^2(r\overline{r})^2+\mathrm{},$$
(103)
where $`\overline{r}`$ is given by the solution of the equation,

$$\frac{(N1)(N3)}{4\overline{r}^4}=\frac{g}{2N}(\overline{r}^2r_0^2),$$
(104)
and $`\overline{m}^2`$ by

$`\overline{m}^2`$ $`=`$ $`{\displaystyle \frac{3(N1)(N3)}{4\overline{r}^4}}+{\displaystyle \frac{g}{2N}}(3\overline{r}^2r_0^2),`$ (105)
$`=`$ $`{\displaystyle \frac{g}{N}}(3\overline{r}^22r_0^2).`$ (106)
The frequency of oscillation is determined by $`\overline{m}`$, and this is the quantity to be kept fixed as $`N`$ is changed.

The last technical issue is to keep the distance between the center of the initial wave function, $`\stackrel{~}{r}`$, and the minimum of the potential, $`\overline{r}`$, a constant as we vary $`N`$. That is, we keep

$`\delta r=\overline{r}\stackrel{~}{r},`$

constant for all $`N`$.

With this strategy of keeping $`G_{\mathrm{}}`$, $`\overline{m}`$ and $`\delta r`$ fixed, we can now determine how the coupling constant must vary with $`N`$. We first define $`m^2`$ to be the second derivative of $`V(r)`$ evaluated at $`r=\overline{r}`$,

$$m^2=\frac{\mathrm{d}^2V(r)}{\mathrm{d}r^2}|_{r=\overline{r}}=\frac{g}{2N}(3\overline{r}^2r_0^2).$$
(107)
Then, (105) becomes

$$m^2(N)=\overline{m}^2\frac{3(N1)(N3)}{4\overline{r}^4}.$$
(108)
Solving (107) for $`r_0`$, substituting into (104) and solving for $`g`$ gives:

$$g(N)=\frac{N}{\overline{r}^2}\left\{\overline{m}^2\frac{(N1)(N3)}{\overline{r}^4}\right\}.$$
(109)
with $`\overline{r}=\stackrel{~}{r}+\delta r`$. The value of $`r_0^2`$ is then determined by (104):

$$r_0^2(N)=\overline{r}^2\frac{2N}{g(N)}\frac{(N1)(N3)}{4\overline{r}^4}.$$
(110)
Thus, for fixed values of $`G_{\mathrm{}}`$, $`\overline{m}`$ and $`\delta r`$, Eqs. (101), (102), (109), and (110) determine values for $`G(N)`$, $`\stackrel{~}{r}(N)`$, $`g(N)`$, and $`r_0(N)`$ all values of $`N`$.

In the limit, $`N\mathrm{}`$, we find that

$`g(\mathrm{})`$ $`=`$ $`{\displaystyle \frac{1}{G_{\mathrm{}}}}\left(\overline{m}^2{\displaystyle \frac{1}{G_{\mathrm{}}^2}}\right)`$ (111)
$`m^2(\mathrm{})`$ $`=`$ $`\overline{m}^2{\displaystyle \frac{3}{4G_{\mathrm{}}^2}}`$ (112)
To summarize, in order to establish appropriate initial conditions for the quantum roll problem, we have kept the variance $`G_{\mathrm{}}`$ constant instead of $`G`$, and have allowed the parameters describing the potential function $`g`$ and $`r_0`$ to change with $`N`$ in order to compare solutions that have close to the same oscillation frequencies. In our numerical runs, we chose the values

$$G_{\mathrm{}}=1;\overline{m}^2=2;\delta r=2;g(\mathrm{})=1$$
(113)
Fig. 1 displays the variation of the potential parameters with $`N`$.

### B Shifted Gaussian Initial Conditions

The second $`O(N)`$ invariant initial condition we investigated had a wave function localized in a wave packet near the center of the valley of the classical potential at $`r=r_0`$. For $`N=1`$ this would be the standard double-well tunnelling problem; for higher values of $`N`$, tunneling is avoided by going around the barrier. Therefore this initial condition is qualitatively different from the roll problem and provides a different arena for testing approximations. However, since this initial condition violates the $`O(N)`$ symmetry of the potential, numerical solution is at present possible only for very small $`N`$.

We take the initial wave function to be a shifted Gaussian of the form:

$$\psi _0(x)=\frac{1}{(2\pi G)^{N/4}}\mathrm{exp}\left\{\underset{i}{}\frac{(x_ir_0/\sqrt{N})^2}{4G}\right\}.$$
(114)
The energy $`E`$ of this state can be determined from Eq. (96) by the substitutions,

$`G_{ij}\delta _{ij}G;q_ir_0/\sqrt{N};p_i0;\mathrm{\Pi }_{ij}0,`$

from which we find:

$$E=\frac{N}{8G}+\frac{g}{8N}\left\{N(N+2)G^22NGr_0^2+r_0^4\right\}.$$
(115)
On the other hand, the height of the classical potential barrier is given by:

$$E_b=\frac{g}{8N}r_0^4.$$
(116)
For $`N=1`$, the necessary requirement for tunneling is that $`E<E_b`$.

In the general case (arbitrary $`N`$), we have:

$$M^2=\frac{^2V(r)}{r^2}|_{r_0}=\frac{g}{N}r_0^2\mathrm{or}r_0^2=\frac{N}{g}M^2.$$
(117)
If the initial state is close to the ground state of a harmonic potential that approximates the potential at the bottom of the well, then the width $`G`$ of the wave function is, approximately,

$$G=\frac{1}{2\sqrt{M^2}},$$
(118)
which can be combined with Eq. (115) to give the desired energy of the initial state in terms of the values of $`N`$ and $`g`$.

We are interested in initial conditions where the energy per oscillator does not increase as a function of $`N`$. To implement this we fix $`M^2=1`$, which corresponds to $`G=1/2`$ for the initial width. The barrier height is then given by

$$E_b=\frac{N}{8g},$$
(119)
and the total energy by

$$E=\frac{N+1}{4}+\frac{N+2}{32}g$$
(120)
We explored three cases: $`E=0.5E_b`$, $`E=E_b`$, and $`E=2E_b`$. For each of these cases, Eqs. (119) and (120) determine $`g`$ for each $`N`$. In all cases we took $`x_i(0)=r_0/\sqrt{N}`$, $`\dot{x}_i=0`$, $`G_{ij}(0)=G\delta _{ij}`$, and $`\dot{G}_{ij}(0)=0`$. As a consequence, all of the oscillators $`x_i(t)`$ move identically.

## VIII Effective potential

It is well-known that the static effective potential is not always a useful guide to the true dynamics of the system (See, e.g., Ref. ). Nevertheless, one may seek to gain qualitative insight into some aspects of quantum dynamics this way, though care is certainly indicated (See, e.g., Ref.  for the Gaussian effective potential). Indeed, there appears to be an interesting connection with the properties of the effective potential at next-to-leading order and with the corresponding dynamical evolution (discussed in the next section).

The effective potential in the large $`N`$ approximation has been previously obtained by Root  to order $`1/N`$; however, we recalculate it here using our equations. When $`x_i`$ and $`\chi `$ are independent of time, we can ignore the closed time path ordering and use Fourier transforms, passing the poles by using the Feynman contour. Then, from the action given in (15), we find

$`V_{\text{eff}}^{[1]}(r,\chi )={\displaystyle \frac{N\chi }{g}}\left(\mu ^2{\displaystyle \frac{\chi }{2}}\right)+{\displaystyle \frac{1}{2}}\chi r^2`$ (122)
$`+{\displaystyle \frac{N}{2}}{\displaystyle \frac{dk}{2\pi i}\mathrm{ln}[\stackrel{~}{G}^1(k)]}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{dk}{2\pi i}\mathrm{ln}[\stackrel{~}{D}^1(k)]},`$
where $`\chi `$ satisfies the requirement

$`{\displaystyle \frac{}{\chi }}V_{\text{eff}}(r,\chi )=0.`$ (123)
In this section to make contact with Root , we have $`\mu ^2=gr_0^2/(2N)<0`$. In order to examine the large $`N`$ limit, we again rescale (7) to the $`y`$ variables. Then for the Green’s functions, we find:

$`\stackrel{~}{G}^1(k)`$ $`=`$ $`(k^2\chi ),`$
$`\stackrel{~}{D}^1(k)`$ $`=`$ $`{\displaystyle \frac{N}{g}}Ny^2\stackrel{~}{G}(k)+{\displaystyle \frac{iN}{2}}{\displaystyle \frac{dp}{2\pi }\stackrel{~}{G}(p)\stackrel{~}{G}(kp)}`$
$`=`$ $`{\displaystyle \frac{N}{g}}\left\{1g{\displaystyle \frac{y^2}{k^2\chi }}{\displaystyle \frac{g}{2\sqrt{\chi }}}{\displaystyle \frac{1}{k^24\chi }}\right\}`$
$`=`$ $`{\displaystyle \frac{N}{g}}{\displaystyle \frac{(k^2m_+^2)(k^2m_{}^2)}{(k^24\chi )(k^2\chi )}},`$
where $`m_\pm ^2=b\pm \sqrt{b^2c}`$, with

$`b`$ $`=`$ $`{\displaystyle \frac{5}{2}}\chi +{\displaystyle \frac{g}{2}}\left(y^2+{\displaystyle \frac{1}{2\sqrt{\chi }}}\right)`$
$`c`$ $`=`$ $`4\chi ^2+g\left(4y^2\chi +{\displaystyle \frac{1}{2}}\sqrt{\chi }\right).`$
For the Feynman contour, we have

$$\frac{dk}{2\pi i}\mathrm{ln}(k^2\chi )=\sqrt{\chi }+\text{constant terms}.$$
(124)
Thus the effective potential (122) becomes

$`{\displaystyle \frac{V_{\text{eff}}^{[1]}(y,\chi )}{N}}`$ $`=`$ $`{\displaystyle \frac{\chi }{2}}\left(y^2y_0^2\right){\displaystyle \frac{\chi ^2}{2g}}+{\displaystyle \frac{\sqrt{\chi }}{2}}`$ (126)
$`+{\displaystyle \frac{1}{2N}}\left(m_++m_{}3\sqrt{\chi }\right).`$
The gap equation which determines $`\chi `$ follows from (123)

$`\chi `$ $`=`$ $`{\displaystyle \frac{g}{2}}\left(y^2y_0^2\right)+{\displaystyle \frac{g(N3)}{4N\sqrt{\chi }}}+{\displaystyle \frac{g}{2N}}{\displaystyle \frac{(m_++m_{})}{\chi }}.`$ (127)
To leading order in the large $`N`$ expansion, Eqs. (126, 127) reduce to the parametric set

$`{\displaystyle \frac{V_{\text{eff}}^{[0]}(\chi )}{N}}`$ $`=`$ $`{\displaystyle \frac{\chi ^2}{2g}}+{\displaystyle \frac{\sqrt{\chi }}{4}},`$ (128)
$`y^2(\chi )`$ $`=`$ $`y_0^2+{\displaystyle \frac{2}{g}}\chi {\displaystyle \frac{1}{2\sqrt{\chi }}}.`$ (129)
Equations (126) and (127) agree with Root, however he used the leading order expression for $`\chi `$ in (129), rather than the full $`\chi `$ of (127).

There exist two real solutions of Eq. (127) for $`\chi `$ with $`y`$ greater than some minimum value $`y_{\text{min}}`$. The next-to-leading order large $`N`$ effective potential, from Eq. (126) is therefore double valued for $`y>y_{\text{min}}`$, and does not exist for smaller values of $`y`$. The physical solution branch corresponds to the one that matches on to the leading order result; the other branch is an unphysical solution. Since it follows from a Legendre transformation, the effective potential (at any order in $`1/N`$) has to be a convex function. The nonexistence of the effective potential at $`y<y_{\text{min}}`$ implies that no quantum state can be associated with the next-to-leading order large $`N`$ approximation in this range.

In Fig. 2, we plot the physical branch of the effective potential as a function of $`y`$, for values of $`N`$ from $`1`$ to $`100`$, for the case $`g=1`$ and $`y_0=2`$. For comparison, we also show in this figure the leading order potential function from Eq. (129), which is single-valued and finite for all $`y`$. (In contrast to the next-to-leading order case we can always associate a Gaussian wave function with the leading order approximation.)

In the case of the Hartree approximation, one can define an “effective potential” as the expectation value of the Hamiltonian using the variational wave function (30) for static configurations . Setting $`p_i(t)=0`$ and $`\mathrm{\Pi }_{ij}(t)=0`$, and putting $`_iq_i^2=r^2`$ and $`G_{ij}=\delta _{ij}G`$ in Eqs. (37) and (43), we find:

$`{\displaystyle \frac{V_{\text{eff}}^{[H]}(y,G)}{N}}`$ $`=`$ $`{\displaystyle \frac{1}{8G}}+{\displaystyle \frac{g}{8}}\left(y^2y_0^2\right)^2`$ (131)
$`+g{\displaystyle \frac{N+2}{4N}}\left(y^2+{\displaystyle \frac{1}{2}}G\right)G{\displaystyle \frac{g}{4}}y_0^2G,`$
The value of $`G`$ is fixed by the requirement that

$$\frac{V_{\text{eff}}^{[H]}(y,G)}{G}=0,$$

which gives the gap equation for the Hartree approximation:

$$\chi =\frac{g}{2}(y^2y_0^2)+\frac{g}{N}\left(y^2+\frac{1}{2\chi }\right)+\frac{g}{4\sqrt{\chi }},$$
(132)
where we have set $`G=1/2\sqrt{\chi }`$. Parametric equations for the Hartree effective potential are then given by:

$`{\displaystyle \frac{V_{\text{eff}}^{[H]}(\chi )}{N}}`$ $`=`$ $`{\displaystyle \frac{1}{2g}}\left({\displaystyle \frac{N}{N+2}}\right)^2\chi ^2{\displaystyle \frac{N}{(N+2)^2}}y_0^2\chi `$ (135)
$`+{\displaystyle \frac{N+4}{4(N+2)}}\sqrt{\chi }+{\displaystyle \frac{g}{2(N+2)^2}}y_0^4`$
$`+{\displaystyle \frac{g}{4(N+2)}}{\displaystyle \frac{y_0^2}{\sqrt{\chi }}}{\displaystyle \frac{g}{16N}}{\displaystyle \frac{1}{\chi }}`$
$`y^2(\chi )`$ $`=`$ $`{\displaystyle \frac{N}{N+2}}\left(y_0^2+{\displaystyle \frac{2}{g}}\chi {\displaystyle \frac{1}{2\sqrt{\chi }}}\right){\displaystyle \frac{1}{(N+2)\chi }}.`$ (136)
Note that in the limit $`N\mathrm{}`$, Eq. (135) reduces to Eq. (129), the leading order large $`N`$ result. Note also that the Hartree effective potential is not derived from a Legendre transform and hence is not subject to a convexity constraint.

In Fig. 3, we plot the Hartree effective potential from Eq. (135) as a function of $`y`$, for our chosen parameters $`g_{\mathrm{}}=1`$ and $`y_0=2`$, for different values of $`N`$. In contrast to the smooth behavior exhibited by the large $`N`$ effective potentials, the Hartree effective potential shows a “first-order transition” in the placement of the minimum of the potential as a function of $`N`$.

The minimum of the effective potential corresponds to a determination of the ground state energy. In Fig. 4, we show values of the minimum energies of the large $`N`$ and Hartree effective potentials as a function of $`N`$. For $`N2`$, the Hartree minimum is generally greater than that for the next-to-leading order large-$`N`$ approximation. Since the Hartree approximation is a variational ansatz it gives an upper bound to the minimum energy. The fact that the next-to-leading order large $`N`$ results are lower than this bound is encouraging, although no guarantee of absolute accuracy.

It is interesting to ask the question how the point $`y_{\text{min}}`$, below which the next-to-leading order large $`N`$ effective potential does not exist, changes as a function of $`N`$. We know that at “infinite $`N`$,” (leading order), $`y_{\text{min}}=0`$ but it is important to know how this limit is reached. For instance, is there a finite value of $`N`$ beyond which $`y_{\text{min}}=0`$? In Fig. 5, we plot $`y_{\text{min}}`$ as a function of $`N`$, for the Hartree and next-to-leading order large $`N`$ approximations. As already stated, the Hartree approximation displays a first order phase transition between the broken and unbroken symmetry solutions at $`N=6.2`$, whereas for the next-to-leading order large $`N`$ approximation a different type of behavior is found: for $`N18.6`$, $`y_{\text{min}}`$ is finite, but for $`N18.6`$, it hits the origin. Thus for $`N18.6`$, we can associate a quantum state (though not known explicitly) with the next-to-leading order approximation.

The critical value of $`N`$ is fixed by the value of $`\chi `$ at the gap equation at the inflection point. If we write the gap equation (127) as

$$f(\chi ,y^2,N)=f_0(\chi ,y^2)+\frac{1}{N}f_1(\chi ,y^2)=0,$$
(137)
where

$`f_0(\chi ,y^2)`$ $`=`$ $`{\displaystyle \frac{g}{2}}(y^2y_0^2)+{\displaystyle \frac{g}{4}}{\displaystyle \frac{1}{\sqrt{\chi }}}\chi ,`$
$`f_1(\chi ,y^2)`$ $`=`$ $`{\displaystyle \frac{3g}{4}}{\displaystyle \frac{1}{\sqrt{\chi }}}+{\displaystyle \frac{g}{2}}{\displaystyle \frac{(m_++m_{})}{\chi }}.`$
Then the critical point is determined by

$$N_\text{c}=\frac{f_1(\overline{\chi },0)}{f_0(\overline{\chi },0)}=\frac{f_1(\overline{\chi },0)/\overline{\chi }}{f_0(\overline{\chi },0)/\overline{\chi }},$$
(138)
where $`\overline{\chi }`$ is given by the solution of this system of equations. For the parameters of Eq. (113), we find numerically that $`N_\text{c}=18.60`$ in excellent agreement with the results shown in Fig. 5.

## IX Numerical Results

### A Quantum Roll

We begin with a discussion of our results for the quantum roll problem. We first examine the short time behavior, $`0<t<3`$, to see if the next-to-leading order large $`N`$ approximation gives an improvement over the leading order solutions. In Fig. 6, we plot the values of $`r^2/N`$ from the numerical solution, the leading and next-to-leading order large-$`N`$ approximations, and the Hartree approximation, for $`N=20`$. The next-to-leading order large $`N`$ approximation is clearly better than the leading order solution, and also better than the Hartree results. Similar behavior is seen for other values of $`N`$ (we also ran $`N=50,80,`$ and 100).

The long time behavior of these approximations is typically of much more interest. We examined behavior over times $`0<t<100`$ to see how long the approximations remained viable. Fig. 7 displays $`r^2/N`$ for the numerical solution, the leading and next-to-leading order large-$`N`$ approximations, and the Hartree approximation for $`N=20`$ and 100. The next-to-leading order large-$`N`$ approximation for $`N=20`$ blows up at $`t84`$. This instability is connected to a violation of unitarity in the particular implementation of this approximation and will be discussed in greater detail below. In general, at these moderate values of $`N`$, the approximations track the numerical solutions reasonably well though they do get out of phase as time progresses. As $`N`$ is increased, the phase errors are considerably reduced as is apparent in the results for $`N=100`$.

The energy of the next-to-leading order large $`N`$ and Hartree approximations is the same as the exact one, but the energy of the leading order large $`N`$ approximation differs from it by terms of order $`1/N`$. (This is because we need to keep the initial values of the parameters the same.) To make a comparison between the approximations this difference has to be compensated for; we do this by rescaling time by a constant multiplicative factor so as to match the last oscillation maxima. This effect is of order $`1/N`$. For $`N=100`$, we find that for $`0<t<100`$, the next-to-leading order $`1/N`$ approximation is always more accurate than the leading order; however, when comparing the next-to-leading order with the Hartree, although less accurate for $`t<50`$, the Hartree approximation starts becoming more accurate at $`t50`$ (however, the errors are very similar in magnitude, of the order of a few per cent).

We now return to a discussion of the blow ups first encountered in the $`N=20`$ case discussed above. The failure of truncation schemes (of which $`1/N`$ is an example) to maintain a rigid connection with the existence of a probability distribution function is a well-known problem in nonequilibrium statistical mechanics. It is often the case that, when this connection is lost (failure of reality or positivity conditions), instability soon follows. In our case, the violation of unitarity is manifest in that positive expressions such as $`r^2/N`$ can turn out to be negative. Note that since both the Hartree and leading order $`1/N`$ approximations are variational in nature, they can never violate unitarity. In Fig. 8, we show the blow up or failure time for the next-to-leading order large $`N`$ approximation as a function of $`N`$; the failure time being defined as the time at which $`r^2/N`$ becomes negative. It is interesting to note that for values of $`N`$ near the critical value $`N_\text{c}=18.60`$ where the effective potential extends down to the origin, the failure time starts to increase rapidly (consistent with our interpretation of the connection of the static effective potential to an associated quantum state). For values of $`N`$ greater than 21, we could not find failure for $`t<150`$. Thus, it may well be that above a certain value of $`N`$, the unitarity violation is pushed out to times of no practical significance, or may even disappear altogether.

We discussed previously how to choose initial conditions so as to reach the large $`N`$ limit in a controlled manner. Starting with these initial conditions, we show in Fig. 9 results from numerical solution of the Schrödinger equation for the time evolution of $`r^2/N`$. In the strict large $`N`$ limit one expects pure harmonic oscillations about the minimum of the infinite-$`N`$ effective potential. For any finite value of $`N`$, however, at sufficiently large times the pure harmonic motion is overcome by nonlinear effects, and interesting behavior is found, as shown in Fig. 10, indicating the presence of a very non-Gaussian wave function. The ability to capture this long time behavior is an important test of $`1/N`$ methods.

In order to test whether Hartree or the next-to-leading order approximation incorporate nonlinearities correctly so as to capture the late time modulation behavior, we ran a comparison against the numerical results for $`N=21`$, the results being displayed in Fig. (11). It is clear that both approximations do not give satisfactory results. This provides additional motivation for the development of alternative $`1/N`$ expansions which would incorporate selective resummations in order to reduce the coefficient of the error term at late times.

### B Shifted Gaussian Initial Conditions

We now discuss the time evolution of a quantum state having an initial wave function given by Eq. (114). For this problem, because of the lack of symmetry, exact solutions were only obtained for $`N2`$. For $`N=1`$, depending on whether the energy is above or below the barrier height, one observes either slow tunneling with rapid oscillations in one well, or slower oscillation in the complete range. At these low values of $`N`$, the large-$`N`$ expansion breaks down quickly, as in the quantum roll, but even here at $`N=1`$, the $`1/N`$ corrections improve the short time accuracy of $`q(t)`$.

A more relevant comparison is to consider larger values of $`N`$ at which the approximations have a better chance of capturing the exact behavior. In the next figure, we compare the Hartree with the leading and next-to-leading order large-$`N`$ approximation for $`q(t)`$, at $`N=50`$. Fig. 12 displays the results for a run with $`E>E_b`$ using the equal time Green’s function approximation (see Sec. V) method for obtaining the Hartree results. Unlike the situation for the roll initial condition where the Hartree and next-to-leading order large $`N`$ results are not dramatically different, here the qualitative behavior is quite dissimilar (whereas the Hartree and leading order results are in fact very close).

## X Conclusions

Testing the $`1/N`$ approximation in quantum mechanics has already enabled us to arrive at some useful conclusions. In order to interpret our results, it is important to keep in mind that $`1/N`$ approximations are a form of resummed perturbation theory and are therefore only valid at weak coupling. Thus for couplings of order unity, it is unrealistic to expect the approximation to give good results for small values of $`N`$. Our results have shown that at sufficiently large $`N`$ the next-to-leading order approximation is a clear improvement over the leading order approximation, however, at late times this approximation (as well as Hartree) fails to capture the nonlinear effects that lead to nontrivial amplitude modulation of the radial oscillations in the quantum roll problem.

We have noted the presence of a finite-time breakdown in the evolution given by the next-to-leading order approximation. This result is related to the fact that the large $`N`$ expansion for the expectation values does not necessarily correspond to a positive semi-definite density matrix when truncated at any finite order in $`1/N`$. (At lowest order, the $`1/N`$ approximation is equivalent to a Gaussian variational ansatz for the density matrix and does not have this problem.) This last aspect is already clear even in static situations such as the lack of a real effective potential for all $`N`$ in the next-to-leading order approximation. This type of finite-time breakdown induced by unitarity/positivity violation has also been noted in simulations of quantum systems where the coupled equal times Green’s function approach was truncated at fourth order  or where high order cumulant expansion methods were used .

Two aspects of this breakdown deserve further mention: First, the time at which breakdown occurs appears to be strongly connected with the behavior of the effective potential. For values of $`N`$ not very much bigger than the critical value $`N_c`$ (beyond which the effective potential exits over the entire range of $`y`$), the breakdown time increases extremely steeply and may even be pushed to times long enough to be no longer an obstacle to practical calculations ( this still needs to be demonstrated). Second, it is important to point out that avoiding the breakdown via a partial resummation does not automatically guarantee better late time accuracy (or convergence) since such a scheme is also only next-to-leading order accurate. However, it will help in the sense that one may carry out simulations at smaller values of $`N`$, thus making it easier to compare against the late-time numerical solutions of the corresponding Schrödinger equation.

One possible way of correcting the problem of a manifestly positive operator such as $`r^2`$ becoming negative is to solve for the full Green’s function $`𝒢_{ij}(t,t^{})`$:

$`𝒢_{ij}(t,t^{})=G_{ij}(t,t^{})`$ (139)
$`{\displaystyle \underset{k,l}{}}{\displaystyle _𝒞}dt_1{\displaystyle _𝒞}dt_2G_{ik}(t,t_1)\mathrm{\Sigma }_{kl}(t_1,t_2)𝒢_{lj}(t_2,t^{}),`$ (140)
rather than the next-to-leading order one as in (21). This equation is the exact equation one obtains by varying the effective action and it contains terms of all orders in $`1/N`$ (thus, strictly speaking, one is no longer truncating at some fixed order).

However, just making this correction does not increase the time period during which the approximation is accurate. In order to extend the accuracy of the $`1/N`$ approximation to late times, it appears necessary to use a more robust approximation based on the Schwinger-Dyson equations. Several approximations of this sort are possible, which may both cure the positivity problem as well as lead to accurate results at late times. These will be discussed separately .

## XI Acknowledgements

The authors acknowledge helpful conversations with Luis Bettencourt, Yuval Kluger, and Emil Mottola. S.H. acknowledges stimulating discussions with Larry Yaffe. The work of B.M. and J.F.D. at UNH is supported in part by the U.S. Department of Energy under grant DE-FG02-88ER40410. B.M. and J.F.D thank the Theory Group (T-8) at LANL and the Institute for Nuclear Theory at the University of Washington for hospitality during the course of this work. F.C. would like to thank the Physics Department at Yale University for their hospitality where some of this research was carried out.