text:
# Incorporating Chiral Symmetry and Heavy Quark Theory in Extrapolations of Octet Baryon Charge Radii

## I Introduction

Lattice QCD is so far the most successful method of calculating hadronic observables from the theory of QCD. However, computational limitations mean that hadronic observables are calculated on the lattice at quark masses larger than their physical values. Hence results from lattice simulations cannot be directly compared with experimental data. Although, with improvements in actions, algorithms and computer speed, future lattice calculations will be performed much closer to the physical regime, these improvements will proceed over many years. Therefore, to make sense of any lattice results produced to date, and to compare them with experiment, one must understand how to extrapolate lattice results, obtained at large quark masses, to the physical world.

One of the difficulties with calculating hadronic observables at heavy quark masses on the lattice is that chiral perturbation theory cannot be applied in this quark mass regime. However, chiral expansions of hadronic observables contain important non-analytic terms as a function of the quark mass, $`m_q`$ (or equivalently of $`m_\pi ^2`$, as $`m_qm_\pi ^2`$ in this range). It is vital that this non-analytic behaviour is included in any extrapolation to the physical regime .

The chiral expansion of the squared electric charge radius of a spin-1/2 octet baryon includes non-analytic behaviour in the form of logarithmic terms in $`m_\pi `$ (or $`m_q`$). To extrapolate the lattice results for electric charge radii we incorporate these logarithmic terms in our extrapolation formulae, while ensuring that the correct heavy quark behaviour is also maintained. A similar approach  has been successful in explaining why lattice simulations of pion and proton charge radii are similar in size, while experimental measurements reveal a significant difference. The dramatic differences in the chiral behaviour of the pion and proton charge radii account for the similarity of the lattice results at moderately heavy pion masses, while allowing good agreement with experiment at the physical pion mass. This illustrates the importance of including meson cloud effects in extrapolations of lattice results to the physical regime. In this paper, we improve the formalism of chiral extrapolations by incorporating both chiral symmetry and heavy quark effective theory.

Our lattice “data” for the electric charge radii is taken from the calculations of Ref. . The contributions to the charge radii from individual quark flavours are also given there. These results are the only available lattice calculations of the electric charge radii of spin-1/2 octet baryons. For each baryon, the electric charge radius was calculated at three different quark masses, corresponding to rather heavy pion masses (all exceeding 600 MeV). Here we extrapolate these results as functions of the squared pion mass, $`m_\pi ^2`$, to obtain predictions for the charge radii at the physical pion mass (139 MeV). Because the lattice calculations are quenched, we expect that there are errors in the lattice data which we have been unable to take into account. However, we expect that errors from the quenching approximation will be rather small. At the quark masses considered on the lattice, the dominant effect is a simple renormalization of the strong coupling constant, accounted for in setting the lattice spacing scale. When results from unquenched simulations become available, the formalism presented here may be readily applied.

## II Extrapolations

Dynamical breaking of chiral symmetry in the QCD lagrangian results in the formation of an octet of (pseudo-) Goldstone bosons. Goldstone boson loops give rise to significant non-analytic behaviour in hadronic observables, such as $`r_i^{2}`$ and magnetic moments, as a function of the quark mass, $`m_q`$. Using an expansion about the chiral SU(3) limit gives the following expression for the squared electric charge radius, $`r_i^{2}`$, of a spin-1/2 octet baryon (labelled by $`i`$)

$$r_i^{2}=\gamma _i+\underset{X=\pi ,K}{}\frac{6\alpha _i^{(X)}}{\left(4\pi f_\pi \right)^2}\mathrm{ln}\left(\frac{m_X}{\lambda }\right)+\mathrm{}.$$
(1)
Here $`f_\pi `$ is the pion decay constant (93 MeV) and $`\lambda `$ is the scale of the dimensional regularization. (The value of $`\gamma _i`$ is clearly correlated with the choice of $`\lambda `$.) Unlike SU(2)-flavour symmetry, SU(3) is significantly broken in the physical world, with the strange quark mass the same order of magnitude as $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ for low-energy phenomenology. The squared kaon mass exceeds the squared pion mass by over an order of magnitude. Given that the source of the meson cloud associated with a baryon is of a finite size, one might anticipate that the role of the kaon cloud will be suppressed away from the SU(3) chiral limit . The form factor describing the finite size of the kaon source will act to suppress kaon loop effects like $`m_X^4`$ at large $`m_X`$ (comparable to $`m_K`$). This was demonstrated within a particular chiral quark model for the nucleon magnetic moments in Ref. . Despite the model-dependence associated with the form factors, the lattice results themselves are very slowly varying functions of $`m_X`$ at values of the order $`m_K`$ or higher, thus supporting the general conclusion. Therefore kaon loop effects are expected to be small and slowly varying as a function of $`m_q`$. Hence we do not explicitly include kaon contributions in our extrapolation formulae. Conversely, since $`m_\pi `$ varies rapidly with $`m_q`$, the leading non-analytic behaviour in $`m_\pi `$ must be included explicitly in an extrapolation to the chiral regime.

To extrapolate the lattice calculations of the electric charge radii of the spin-1/2 octet baryons, we consider two distinct fitting procedures. Both these extrapolation schemes satisfy the constraints of chiral perturbation theory and heavy quark effective theory. The first extrapolation procedure we investigate is given simply by the formula

$$r_i^{2}=\frac{c_1+\chi _i\mathrm{ln}(m_\pi /\mathrm{\Lambda })}{1+c_2m_\pi ^2},$$
(2)
where $`r_i^{2}`$ are the lattice QCD results (at several values of $`m_\pi `$) extracted from Ref. , $`c_1`$ and $`c_2`$ are fit parameters chosen to best fit these results, $`\chi _i`$ (corresponding to the $`i^{\mathrm{th}}`$ baryon) is fixed (model-independently) by chiral perturbation theory and $`\mathrm{\Lambda }`$, which is directly correlated with $`c_1`$, is fixed at 1 GeV. This extrapolation procedure is not feasible for the neutral baryons because $`r_i^{2}0`$ as $`m_\pi `$ becomes large and thus sensitivity to the $`c_2`$ fit parameter is lost. In order to extrapolate the neutral baryon charge radii results we consider a second extrapolation procedure, focusing on individual quark “sector” (or quark flavour) contributions, as discussed below.

Clearly the extrapolation formula given in Eq. (2) builds in the correct chiral behaviour in the SU(2) limit, since in the limit $`m_\pi 0`$ it can be expanded as follows

$$r_i^{2}=c_1+\chi _i\mathrm{ln}(m_\pi )c_2m_\pi ^2+\mathrm{}$$
(3)
(Recall that the scale $`\mathrm{\Lambda }`$, in Eq. (2), has been set to 1 GeV. This choice is also implicit in Eqs. (6)-(10), below, where $`m_\pi `$ in the logarithm must be in GeV.) This agrees with the chiral SU(2) expansion of the squared electric charge radius (see Eq. (1)), provided we fix the coefficient $`\chi _i`$ to $`6\alpha _i^{(\pi )}/\left(4\pi f_\pi \right)^2`$ for the $`i^{\mathrm{th}}`$ baryon. The coefficients $`\alpha _i^{(\pi )}`$ and $`\chi _i`$ are given in Table I.

In the large $`m_q`$ limit, we expect that the quarks behave non-relativistically, and the squared charge radius falls off as $`m_q^2`$, as it does in non-relativistic quantum mechanics. In the region where $`m_q`$ is very large, $`m_qm_\pi `$, and hence we require that

$$r_i^{2}\frac{1}{m_\pi ^2},$$
(4)
as $`m_\pi `$ becomes extremely large. This is clearly satisfied by our first extrapolation formula, Eq. (2), since the logarithm is very slowly varying.

In the second extrapolation procedure we deal separately with the individual quark sector contributions to the baryon charge radii. For example, in the case of the nucleons, we extrapolate the up and down sector contributions separately. For the hyperons the strange and light sector results are extrapolated separately. This avoids the problem encountered with the neutral baryons which was mentioned previously, because now all the quantities being extrapolated are charged, even if the overall charge on the baryon is zero. This separation is valid because the squared electric charge radius can be decomposed as

$$r_i^{2}=\underset{q=u,d,s}{}e_qr_i^{(q)2},$$
(5)
where $`r_i^{(q)2}`$ is the contribution from the $`q^{\mathrm{th}}`$ quark sector and $`e_q`$ is the charge of this quark sector. Therefore, provided that the extrapolation formulae from each sector add so that the chiral and heavy quark limits of the sum are in agreement with Eqs. (1) and (4) respectively, this method contains the same physics as the first method, but simply makes use of the extra information contained in the individual quark sector results. Not only does this second extrapolation procedure solve the neutral baryon extrapolation difficulty, it also provides predictions for individual quark sector radii, which will be resolved at Jefferson Lab for the nucleon , and perhaps future accelerator facilities for hyperons.

Isolation of the individual quark sector contributions to the charge radii is relatively straightforward from the theoretical point of view. For example, to isolate the $`u`$-sector contribution to the charge radius of the proton one simply sets the $`d`$-quark charge to zero and calculates the proton charge radius as if only the $`u`$ quark carried charge. In the chiral expansion of the proton charge radius, the coefficient of the logarithm, $`\chi _p^{(\pi )}`$, originates from the pion loop, $`pn+\pi ^+`$ (see Fig. 1), and includes the charge of the pion cloud. Therefore, to extrapolate the $`u`$-sector contribution to the proton charge radius, the appropriate coefficient of the logarithm is $`\frac{2}{3}\chi _p^{(\pi )}`$, since the pion now carries charge $`+2/3`$. Thus we extrapolate the $`u`$-sector results of the proton according to

$$2e_ur_p^{(u)2}=\frac{c_1+\frac{2}{3}\chi _p^{(\pi )}\mathrm{ln}(m_\pi )}{1+c_2m_\pi ^2}\text{ ,}$$
(6)
where $`\chi _p^{(\pi )}`$ is the full chiral coefficient of the proton, given in Table I and $`r_p^{(u)2}`$ is the squared charge radius of a single $`u`$-quark of unit charge. Similarly, the $`d`$-sector results are extrapolated according to

$$e_dr_p^{(d)2}=\frac{c_1^{}+\frac{1}{3}\chi _p\mathrm{ln}(m_\pi )}{1+c_2^{}m_\pi ^2}\text{ ,}$$
(7)
where the factor of $`1/3`$ originates from the $`\overline{d}`$ contribution to the pion cloud. Clearly adding the left hand sides of Eqs. (6) and (7) yields the full expression for $`r_p^2`$. In the chiral limit the right hand sides add so that the correct chiral form for $`r_p^2`$, given in Eq. (3), is retained. The sum of Eqs. (6) and (7) also obeys the correct heavy quark behaviour, given in Eq. (4). Since the parameters $`c_2`$ and $`c_2^{}`$ are not necessarily the same, the individual quark sector extrapolation formulae (Eqs. (6) and (7)) cannot be added directly to give Eq. (2). Therefore, in general we do not expect the two extrapolation procedures to give exactly the same results. For the charged baryons this may be used to help quantify the systematic error of the approach. For the neutron, the analogous extrapolation functions are given by

$$e_ur_n^{(u)2}=\frac{c_1+\frac{2}{3}\chi _n^{(\pi )}\mathrm{ln}(m_\pi )}{1+c_2m_\pi ^2},$$
(8)
for the $`u`$-quark sector, and

$$2e_dr_n^{(d)2}=\frac{c_1^{}+\frac{1}{3}\chi _n^{(\pi )}\mathrm{ln}(m_\pi )}{1+c_2^{}m_\pi ^2},$$
(9)
for the $`d`$-quark sector results, where $`\chi _n^{(\pi )}`$ is given in Table I.

We now consider extrapolating the hyperon charge radii results using the second extrapolation procedure, i.e., extrapolating the strange and light quark sector results separately. In extrapolating the strange quark sector the charges of the light quarks are set to zero. Since the logarithmic term considered here originates from pionic corrections to the charge radius (and pions do not contain strange quarks), the coefficient of the logarithm in the strange sector extrapolation will be zero. Similarly, for the light sector extrapolation, strange quarks do not carry charge and hence the coefficient of the logarithm in this extrapolation will be the full coefficient, $`\chi _i^{(\pi )}`$. This results in the following extrapolation formulae for the hyperon quark sector contributions

$$e_lr_i^{(l)2}=\frac{c_1+\chi _i\mathrm{ln}(m_\pi )}{1+c_2m_\pi ^2}$$
(10)
and

$$e_sr_i^{(s)2}=\frac{c_1^{}}{1+c_2^{}m_\pi ^2}\text{ ,}$$
(11)
where $`i`$ runs over the hyperons only, and $`l`$ corresponds to the light-quark ($`u`$ and/or $`d`$) sector. Since the strange quark mass is held fixed in the light quark mass extrapolation, any variation in the strange quark sector is purely an environment effect from the surrounding light quarks. As such, the functional form for the strange quark sector is constrained by neither leading order chiral perturbation theory nor heavy quark effective theory. As we shall see, $`c_2^{}`$ is small and negative for each hyperon, which suggests that a simple linear ansatz for the strange quark sector extrapolation could also have been used.

## III Results

The lattice QCD simulations were performed on a $`24\times 12\times 12\times 24`$ periodic lattice using standard Wilson actions at $`\beta =5.9`$. Dirichlet boundary conditions were used for fermions in the time direction. Twenty-eight quenched gauge configurations were generated by the Cabibbo-Marinari  pseudo-heat-bath method. The conserved vector current was derived from the Wilson fermion action via the Noether procedure. The associated lattice Ward identity protects this vector current from renormalization. The radii were produced by fitting the electric form factor to dipole and monopole forms, allowing a charge radius to be extracted in each case. Since it is known from experiment that the dipole form is more suitable for parameterizing the electric form factor, we consider only the dipole results here. Statistical uncertainties in the lattice simulation results are calculated in a third-order, single elimination jackknife . Further details may be found in Ref. .

The extrapolations of lattice calculations for the charge radii of the spin-1/2 baryon octet are shown in Figs. 27. Extrapolations of baryon charge radii results, performed according to Eq. (2), are indicated by the solid lines, where the full circles ($``$) represent the baryon charge radii from lattice QCD and the extrapolated value at $`m_\pi =139`$ MeV. The individual quark sector extrapolations are shown by the dashed and dot-dashed lines, and the baryon charge radius predicted by this method is indicated by a full square ( ). Experimental measurements are indicated at the physical pion mass by an asterisk ($``$). Note that for the charged baryons, two extrapolation schemes and two corresponding predicted physical values are shown, whereas (for the reasons explained in section II) only one extrapolation procedure is shown for each neutral baryon.

In the case of the proton, the two extrapolated results agree very well with the experimental measurement. It can be seen that a traditional linear extrapolation in $`m_\pi ^2`$ would significantly underestimate the experimental result. Similarly the predicted charge radius for the neutron (produced by separate extrapolations of the $`u`$\- and $`d`$-sector results) agrees with the experimental data significantly better than the prediction from a conventional linear extrapolation in $`m_\pi ^2`$. Both predicted values for $`\mathrm{\Sigma }^{}`$ also agree very well with the two experimental measurements. These baryons are currently the only baryons of the spin-1/2 octet whose electric charge radius has been measured.

From simple quark model arguments , where the heavier strange quark has a smaller distribution than the light quarks, we expect the hierarchy of electric charge radii of charged octet baryons to be given as follows

$$\left|r_{\mathrm{\Sigma }^+}^2\right|\left|r_p^2\right|\left|r_\mathrm{\Sigma }^{}^2\right|\left|r_\mathrm{\Xi }^{}^2\right|$$
(12)
Clearly the results of the lattice extrapolations shown in Table I are in qualitative agreement with this expectation. Indeed, in the regime of the actual lattice data ($`m_\pi 600`$ MeV) the argument is even quantitative. However, as the chiral limit (and physical pion mass) is approached, the simple quark model description is no longer adequate and chiral physics gives rise to dramatic effects. For example, in the extrapolation of the $`d`$-quark sector of the proton (Fig. 2), chiral effects mean that the $`d`$-quark sector can actually make a positive contribution to the charge radius, via the $`\overline{d}`$ contribution in $`\pi ^+`$ – even though the total contribution is negative at the physical pion mass. This behaviour is not anticipated by the simple quark model.

For the neutral baryons the sign of the squared charge radius is important. In the neutron, the two $`d`$ quarks are most likely to be found in a spin 1 configuration, where they will undergo hyperfine repulsion. This leads to a small, negative charge radius squared. However, as one approaches the chiral limit, spontaneous chiral symmetry breaking, in particular the process $`np\pi ^{}`$, which carries $`d`$-quarks to larger radii and screens the $`u`$-quark contribution via the $`\overline{u}`$ in the $`\pi ^{}`$, leads to an enhancement of the negative charge radius. The remaining neutral baryons, $`\mathrm{\Sigma }^0`$, $`\mathrm{\Lambda }`$ and $`\mathrm{\Xi }^0`$, have a positive squared charge radius. This is because in each case the strange quark distribution is more localized than the light quark charge distribution (due to the larger mass of the strange quark). Therefore on average the light quark charge distribution occurs at a larger radius, resulting in a positive charge radius (since the light quark charge is positive in each case).

The lattice results used here were calculated with a strange quark mass of approximately 250 MeV . This is much larger than the physical strange quark mass of $`115\pm 8`$ MeV at a scale 2 GeV, taken from a careful analysis of QCD sum rules for $`\tau `$ decay . In an earlier study of lattice results for octet baryon magnetic moments  (where the results were extracted from the same lattice simulation ), it was found that the heavy strange quark had a significant effect on the predictions of the $`\mathrm{\Xi }`$ moments. Here we expect that the heavy strange quark should also have some effect on the $`\mathrm{\Xi }`$ charge radii. With a strange quark mass closer to the physical mass the strange quark contribution would be increased. This would result in a lower predicted charge radius for the $`\mathrm{\Xi }^0`$ and a larger (in magnitude) charge radius for the $`\mathrm{\Xi }^{}`$. In the absence of experimental measurements we will not attempt to correct for the effect of the strange quark mass here.

As we see from Table I, the extrapolated mean square charge radii obtained from both extrapolation procedures agree quite well for each charged octet baryon. For the proton, $`\mathrm{\Sigma }^{}`$ and $`\mathrm{\Xi }^{}`$ the reconstructed values completely cover the result from the original extrapolations of Eq. (2). In the case of the $`\mathrm{\Sigma }^+`$ the two values overlap only on the error bars. This is due to the small variation in the strange quark contribution (which is caused by an environment effect). When this environment effect is included in the baryon charge radius, the magnitude of the slope is increased, resulting in a larger charge radius after extrapolation.

In turning the dimensionless masses calculated on the lattice to physical units, the lattice spacing, $`a`$, was set in the traditional manner by fixing the nucleon mass, obtained by a naive linear extrapolation in $`m_\pi ^2`$, equal to the observed mass. Of course, such a linear extrapolation is known  to be inconsistent with chiral symmetry. We have checked that applying a more consistent chiral extrapolation would systematically lower values of $`r^2`$ obtained for the charged octet by of the order 15%. (The effect on neutral baryons is much smaller.) On the other hand, the data which we are forced to use is quenched data which omits some pion corrections. Although these are expected to be suppressed at the large values of $`m_\pi ^2`$ for which the data is available (c.f. Ref.), the associated systematic error would tend to increase the calculated values of $`r^2`$, perhaps by 5-10%. Rather than attempt to repair these deficiencies in the present data, we feel it would be more reasonable to simply accept that there is an additional systematic error of the order 15% associated with the extrapolated values shown in Table I. Adding this systematic in quadrature means that the values in Table I would become, for example, $`r_{}_{}{}^{}p^2=0.68\pm 0.14`$ fm<sup>2</sup>, $`r_{}_{}{}^{}\mathrm{\Sigma }_{}^{}^2=0.54\pm 0.09`$ fm<sup>2</sup>, $`r_{{}_{}{}^{}\mathrm{\Sigma }_{}^{+}}^2=0.77\pm 0.14`$ fm<sup>2</sup>. We look forward to repeating our analysis with unquenched data at lower quark mass, which is the best way to overcome these problems.

## IV Conclusion

In this paper we have investigated two methods of extrapolating lattice results for the electric charge radii of octet baryons to the physical regime. These procedures build in the correct leading non-analytic behaviour of the electric charge radii in the chiral limit, as well as the correct heavy quark behaviour. Both extrapolation procedures were performed for the charged octet baryons, with the predicted values agreeing very well. The extrapolation formulae seem to be very successful, as good agreement with experiment was obtained for the nucleons and the $`\mathrm{\Sigma }^{}`$. We await further experimental measurements of the baryon charge radii in order to test our predictions. In the future we hope to perform similar extrapolations of electric charge radius lattice results calculated with a more realistic strange quark mass and eventually with lighter, dynamical $`u`$ and $`d`$ quarks.

## Acknowledgement

This work was supported by the Australian Research Council and the University of Adelaide. DBL would like to thank Frank Lee and the Center for Nuclear Studies at George Washington University for their kind hospitality during which some of this research was carried out. We would also like to thank Pierre Guichon for helpful comments on the manuscript.


format: text

text:
# Quantum Phase Transitions and the Extended Coupled Cluster Method

## Abstract

We discuss the application of an extended version of the coupled cluster method to systems exhibiting a quantum phase transition. We use the lattice O(4) non-linear sigma model in (1+1)- and (3+1)-dimensions as an example. We show how simple predictions get modified, leading to the absence of a phase transition in (1+1) dimensions, and strong indications for a phase transition in (3+1) dimensions.

preprint: UMIST/Phys/TP/00-10

The nature of the ground-state wave function for (infinite) quantum systems is one of the underlying problems in many areas of physics, from atomic and condensed matter to high-energy physics. Special problems are encountered when one studies the transition between different ground states in cases where the balance of forces in a system with competing interactions is shifted. Since this is such a basic problem, a large number of calculational schemes exist to study the properties of the ground state. These schemes are either imaginary-time approaches, based on the fact that in imaginary time the state at lowest energy has the slowest rate of decay and can thus be filtered out, or are real-time schemes, usually based on the Rayleigh-Ritz variational principle. All of these methods have some weaknesses. In high-energy physics, for example, imaginary-time Monte-Carlo calculations can only be performed on finite lattices. As a real phase transition can only occur in an infinite system, one has to find ways to extrapolate to infinite lattices, and hope no finite-lattice artifacts arise. Even if this can be done, an even more fundamental problem in the Monte-Carlo integrations arises, due to the uncontrollable fluctuations introduced by fermionic degrees of freedom. This problem limits the validity of Monte-Carlo approaches in condensed matter as well as those in lattice-(gauge)-field theory. Real-time methods, on the other hand, are usually based on rather complex many-body techniques, and lack the inherent simplicity of the implementations of Monte-Carlo approaches. Also, even though most real-time methods start from the variational principle, this principle is often violated in an attempt to obtain a viable calculational scheme.

An important difference between real- and imaginary-time approaches is the way they deal with excitations – whereas in imaginary-time methods excited states need to be found by looking at exponentials with faster decay than the ground state, or via complicated Gramm-Schmidt orthogonalisation procedures, the calculation of excitations is often much more straight-forward for the real-time methods. Once we have constructed the ground state we can build excited states on top of this ground state, using sets of excitation operators. The calculational scheme used is then very similar for both the ground state and the excited states.

In this letter we shall concentrate on zero-temperature, also called quantum, phase transitions. The most important exact statement for such situations is the Coleman no-go theorem  which states that, in one space dimension and under rather generous conditions, there is no phase transition to a state with long-range order. As is well known from the exact solution of the $`XXZ`$ model in (1+1)D , this does not preclude a phase transition with an apparent gap-less mode, but without real long-range order. Notwithstanding such special cases, the theorem is a strong indictment against mean-field calculations, usually the lowest order in the real-time approaches, since these often predict phase transitions with long-range order even in one dimension. The reason is that mean-field theory ignores the strong interactions of the massless modes that accompany such a phase transition, which pushes up the mass of the would-be zero-mass modes.

One way to avoid this phenomenon is to concentrate on the phase with unbroken symmetry, and to construct a calculational scheme for that phase. Such approaches have been used in spin-lattice systems and in lattice field theories , but it has been found that they sometimes lead to termination points, probably indicating a phase transition, for the (1+1)-dimensional system as well as in higher dimensions. The excitation spectrum supports this behaviour, and suggests a phase transition to a state with long-range order.

The most complicated part of the real-time approaches is the choice of a truncation scheme. One must make a selection which has enough degrees of freedom not to restrict the physical behaviour of the system, nor to impose any assumptions of the nature of the phase transition on the system, and not too many to make the calculation infeasible. A good choice of wave function is key to any quantum many-body calculation, and more importantly, there must be a way to extend the calculations so as to test any assumptions being made.

Our approach is an example of such a formalism. The approach we use is, however, systematic and will not give rise to spurious phase transitions. The method is a recent offshoot of the coupled cluster method, called the extended coupled cluster method (ECCM) . It has several advantages over earlier versions. First of all, the mean-field approximation is a natural part of the hierarchy of approximations inherent to the method. It is thus suitable to deal with states at both sides of a phase transition, where one exists. We shall show in this letter that an (in principle infinite) set of improvements removes the spurious transition in one dimensional systems. For the moment we shall not include the mean-field term in our calculations, but concentrate on two-site correlations instead. While this restricts us to the symmetric phase, it makes it easier to perform numerical calculations, which can be quite substantial.

We shall illustrate our method with a simple field theory, the $`O(4)`$ non-linear sigma (NLS) model. We have chosen to discretise the Lagrangian for this model on a lattice in space, keeping time continuous. We then construct a Hamiltonian, which is the quantity being studied. The basic degrees of freedom are unitary two-by-two matrices, but these can also be represented by a four-dimensional unit vector $`(𝐧,n_4)`$, using a basis of Pauli matrices $`\tau `$,

$$U=n_4I+i\tau 𝐧.$$
(1)
The Hamiltonian is most succinctly expressed in the angle between nearest-neighbour vectors, and the on-site generalised angular momenta,

$$H=\frac{1}{2}\underset{𝐢}{}𝐈_𝐢^2+\lambda \underset{\mathrm{𝐢𝐣}}{}(1\mathrm{cos}\theta _{\mathrm{𝐢𝐣}}),$$
(2)
where the sum over $`\mathrm{𝐢𝐣}`$ runs over all nearest-neighbour pairs, and counts each pair (or lattice link) only once. The kinetic energy is proportional to the square of the (four-dimensional) angular momentum, $`𝐈_𝐢^2`$. This model shows great similarity with some Heisenberg-like models used in condensed-matter physics.

We expect a phase transition due to the competition between the kinetic and potential terms. For small $`\lambda `$, the weak-coupling limit, the first term is dominant and the system consists of a set of independent rotors. In the strong coupling limit the potential dominates and all vectors tend to align. As we increase $`\lambda `$ the excitation spectrum is also expected to change from a free rotor spectrum at small $`\lambda `$, to a spectrum containing three Goldstone modes, corresponding to the dynamical breaking of chiral symmetry, $`O(4)O(3),`$ at large $`\lambda `$.

This behaviour can most easily be seen in the mean-field approximation. Let us assume that on average all the vectors align along the 4-axis. It can be shown that in that case we obtain  a non-linear version of the Mathieu equation, which has a non-trivial solution (with broken symmetry) only for $`\lambda D>3`$, with a standard first-order phase transition at $`\lambda =3/D`$, where $`D`$ is the dimensionality of space. Below the phase transition point the wave function at each site is just the spherically symmetric “$`S`$-wave” of the free rotor, and above the phase transition we find a two-dimensional degeneracy, corresponding to positive and negative alignment with the 4-axis.

In the symmetric phase, one can ignore the one-body alignment, and concentrate on two-site correlations. We have investigated these processes using the normal coupled-cluster method , and have found that without the one-body term there is a termination point for one, two and three dimensions. This has been interpreted as a signature of a phase transition. In one dimension this cannot be true, since the Coleman theorem forbids such a phase transition. This makes one distrust the other results as well. We show that the extended CCM is able to give much more reliable results.

In the extended coupled cluster method the bra- and ket-states are parametrised independently, corresponding to the use of a bi-orthogonal basis. We define two operators $`\widehat{S}`$ and $`\widehat{S}^{\prime \prime }`$, from which we we build up a variational functional $`I`$ which is the ground state expectation value of the Hamiltonian in a special state,

$$I[S^{\prime \prime },S]=\mathrm{\Phi }_0|e^{\widehat{S}^{\prime \prime }}e^{\widehat{S}}He^{\widehat{S}}|\mathrm{\Phi }_0=\stackrel{~}{\mathrm{\Psi }}|H|\mathrm{\Psi },$$
(3)
where $`|\mathrm{\Phi }_0`$ is the model, or reference, state, that must be chosen separately. A key ingredient in the definition of the correlations is a set of generalised creation operators, that must be chosen with reference to the state $`|\mathrm{\Phi }_0`$, i.e., the hermitean conjugates of the creation operators annihilate $`|\mathrm{\Phi }_0`$. In the symmetric phase we shall use the free-rotor vacuum, the product of “$`S`$”-waves at each lattice site, as reference state. It can be shown that the use of single-site operators is a way to construct the mean-field approximation. When the mean-field is zero, the lowest non-trivial correlations are two-site operators,

$`\widehat{S}`$ $`=`$ $`{\displaystyle \underset{𝐢,𝐣}{}}\widehat{S}_{\mathrm{𝐢𝐣}},`$ (4)
where the sums run over the infinite lattice, and similarly for $`\widehat{S}^{\prime \prime }`$. The operators $`\widehat{S}_{\mathrm{𝐢𝐣}}`$ create a simultaneous excitation at the lattice sites $`𝐢`$ and $`𝐣`$. The precise definition of these operators is not crucial to the current discussion, and can be found in our previous work . By going to the coordinate representation, we find that specification of the operators $`\widehat{S}_{\mathrm{𝐢𝐣}}`$ corresponds to specification of a single function $`S(\mathrm{cos}\theta _{ij})`$. We choose to work in this representation. In the numerical calculation we have to truncate up to a certain length of the relative distance between $`𝐢`$ and $`𝐣`$. Both the bra- and ket-state correlations have to be truncated in the same manner to ensure the conservation of certain analytical properties of the variational principle underlying the method. The functions $`S(\mathrm{cos}\theta _{\mathrm{𝐢𝐣}})`$ can be expanded in Gegenbauer polynomials of the relative angle $`\theta _{\mathrm{𝐢𝐣}}`$ between the unit vectors at lattice sites $`𝐢`$ and $`𝐣`$. Because of the translational symmetry of the lattice the correlation functions depend only on the relative vector $`𝐢𝐣`$. For an optimal choice of $`S`$ and $`S^{\prime \prime }`$ the functional $`I`$ must be stationary with respect to independent variations of the bra- and ket-states. This leads to the non-linear ECCM equations

$$\frac{\delta I[\{S^{\prime \prime },S\}]}{\delta S_{𝐢𝐣}^{\prime \prime }}=0,\frac{\delta I[\{S^{\prime \prime },S\}]}{\delta S_{𝐢𝐣}}=0.$$
(5)
Since the ket state is defined as an exponential it contains automatically all possible independent combinations of correlations on top of the (lowest-order) sum of independent correlations. This is a natural choice, since the correlation operators are combined in linked objects, such as appear in Goldstone’s linked cluster theorem , and, indeed the coupled cluster method itself provides by far the easiest proof of this important theorem . Due to the similarity transform exponential $`e^{\widehat{S}}`$ on the left-hand side of the Hamiltonian in the functional Eq. (3), $`I`$ is a polynomial in the correlation functions. This has large advantages for high-order computer implementations to solve the Eqs. (5).

We shall contrast the results of the extended coupled cluster method to those of the normal coupled cluster method (NCCM), where the bra-state is parametrised linearly in $`S^{\prime \prime }`$,

$$\stackrel{~}{\mathrm{\Psi }}|=\mathrm{\Phi }_0|(1+\widehat{S}^{\prime \prime })e^{\widehat{S}}.$$
(6)
Therefore, in the NCCM, the overlap of the bra and ket states extends only as far as the correlation operator $`\widehat{S}^{\prime \prime }`$. Exponentiating $`\widehat{S}^{\prime \prime }`$ extends the overlap over the whole lattice and also breaks certain symmetries of the functional which allows one to cross a phase transition, instead of breaking down at the phase transition point, which often occurs for the NCCM.

Using the time-dependent variational principle we can find low-lying excitations of the system by considering harmonic fluctuations about the ground state, also called the RPA approximation. Since we have the functional in an analytical form we are able to calculate the excitation energies relatively straightforwardly.

We have applied the ECCM to the NLS model on a (1+1)- and a (3+1)-dimensional cubic lattice. We have studied various truncations, and found that for the expansion of $`S_{\mathrm{𝐢𝐣}}`$ in Gegenbauer polynomials five such functions appear to be enough. For the truncation of the range of the functions $`S`$ we have proceeded as far as practicable, by including all correlations that fit within a cube of size $`L`$. For the (1+1)-dimensional case we can almost use arbitrarily large values of $`L`$, say including a hundred lattice sites. For the $`(3+1)`$-dimensional result the best calculation was limited to those correlations that fit inside a cube of length 5.

For a linear chain the ground-state energy has converged completely with correlations which extend over up to 75 lattice sites. The susceptibility, the second-order derivative of the ground-state energy, changes smoothly with the coupling constant, indicating the absence of a transition region. This is confirmed by the excitation energy. The lowest excitation energy seems to have converged with 5 base functions and 75 correlations, although the global trend sets in at much lower orders. The mass gap remains finite. Clearly, these results are in agreement with the Coleman theorem and are a step beyond the mean-field, and NCCM, results.

For the cubic lattice the ground-state energy has more or less converged for 5 base functions and all 55 correlations that fit in a $`5\times 5\times 5`$ box. The susceptibility shows two regions: a linear, almost constant, susceptibility for the weak-coupling regime, and a steeper behaviour for the strong-coupling regime. In low-order calculations the behaviour between $`\lambda 0.6`$ and $`\lambda 0.7`$ changes drastically, while the trend on both sides outside this region is already predicted by low-order results. Even though the results do not allow us the make a definite statement about the nature of the transition, it seems weaker than a second-order phase transition. The excitation energies again confirm these results. Although the convergence is too slow to determine whether a massless mode appears, it dips at the phase transition and rises fast in the strong-coupling regime. Several low-lying excitation energies turn complex beyond the phase transition, indicating that excitations in this, distinctly different, phase are not properly described without the one-site correlations necessary to describe symmetry breaking. The absence of this phenomenon in (1+1)D is another indication that no phase transition occurs for this case.

The slow convergence in the transition region also indicates that many long-range correlations play an important role as is the case in a real phase transition. The absence of this behaviour in (1+1)D is due to the special properties of one-dimensional system that keep the correlation length finite. The ECCM has the property that both the ket and bra states are expressed in linked-cluster coefficients. Therefore, the length over which ECCM correlation operators extend is a good indication of the length of the physical correlations.

Since we have a dual expansion, in the extent of the correlations, and in each correlation in a set of basis functions, we must look at the truncation systematics for both of these. For the (1+1)-dimensional system, beyond a minimum length of the correlations (about 30 lattice units), the results only depend on the number of basis functions, which indicates a finite, physical correlation length. The convergence with the number of basis functions is quite fast, indicating that the correlations are smooth functions of the relative angle, which indicates that the large number of correlations are essential to describe the phase transition properly since single correlation operators contribute only very little. This also holds for the (3+1)-dimensional system, but we never fully reach convergence with the range of the correlations. Both expansions seems to interplay in a complicated way in the results. However, the ground-state energies generally convergence much faster than the excitation energies.

With the ECCM we find signatures for a phase transition of the $`O(4)`$ non-linear sigma model in (3+1)-dimensions, which confirm our naive expectations. In (1+1)-dimensions the phase transition does not occur, as expected from Coleman’s theorem. These results show clearly that ECCM is a proper many-body technique for the microscopic study of phase transitions and other critical phenomena in cases inaccessible with other methods, and we intend to use it in further studies.

We acknowledge support of a research grant (GR/L22331) from the Engineering and Physical Sciences Research Council (EPSRC) of Great Britain.


text:
# Extremal black holes and the limits of the third law

## Appendix: Asymptotic worldlines for incipient extremal black holes

To derive the late-time worldline for the centre of a body collapsing to an extremal black hole, we start with the usual RN geometry in the case $`Q^2=M^2`$,

$$\mathrm{d}s^2=\left(1\frac{M}{r}\right)^2\mathrm{d}t^2+\left(1\frac{M}{r}\right)^2\mathrm{d}r^2+r^2\mathrm{d}\mathrm{\Omega }^2,$$
(A.1)
where $`\mathrm{d}\mathrm{\Omega }^2`$ is the metric on the unit sphere. This line element, however, is singular on the horizon ($`r=M`$), so we want to transform to Kruskal coordinates. Normally, one does this by first constructing null coordinates $`u=tr_{}`$ and $`v=t+r_{}`$, then going to Kruskal coordinates $`𝒰`$ and $`𝒱`$ which are regular and finite on the horizon. (Here $`r_{}`$ is the usual tortoise coordinate.) However, the usual Kruskal transformation is $`𝒰=e^{\kappa u}`$ and $`𝒱=e^{\kappa v}`$. Because $`\kappa `$ vanishes in the limit of extremality, the Kruskal transformation results in constant values of $`𝒰`$ and $`𝒱`$, which are thus ill-defined in that limit.

This situation can be remedied, however, by a simple generalization of the Kruskal transformation. If one defines $`𝒰`$ and $`𝒱`$ through

$$\begin{array}{c}u=\psi (𝒰)\hfill \\ v=\psi (𝒱)\hfill \end{array}\},$$
(A.2)
where

$$\psi (\xi )=4M\left(\mathrm{ln}\xi \frac{M}{2\xi }\right),$$
(A.3)
then it is not difficult to show that the extremal RN metric (A.1) can be rewritten in terms of $`𝒰`$ and $`𝒱`$ as

$$\mathrm{d}s^2=\frac{(rM)^2}{r^2}\psi ^{}(𝒰)\psi ^{}(𝒱)\mathrm{d}𝒰\mathrm{d}𝒱+r^2\mathrm{d}\mathrm{\Omega }^2.$$
(A.4)
This metric is apparently degenerate on the horizon, but in fact near $`r=M`$ we have $`𝒰=\psi ^1\left(u\right)\psi ^1\left(\psi (rM)\right)=\left(rM\right)`$ and

$$\psi ^{}(𝒰)\frac{4M}{rM}+\frac{2M^2}{(rM)^2}\frac{2M^2}{(rM)^2}.$$
(A.5)
Furthermore, since $`𝒱`$ is everywhere nonzero and finite then $`\psi ^{}(𝒱)`$ is regular there. Now the form taken by the metric (A.4) is asymptotically

$$\mathrm{d}s^2\frac{2M^2}{r^2}\psi ^{}(𝒱)\mathrm{d}𝒰\mathrm{d}𝒱+r^2\mathrm{d}\mathrm{\Omega }^2,$$
(A.6)
because the $`(rM)^2`$ in the numerator of Eq. (A.4) is killed by the $`(rM)^2`$ in the denominator of Eq. (A.5). Consequently, $`𝒰`$ and $`𝒱`$ are good Kruskal-like coordinates.111Double-null coordinates that are regular on general degenerate Killing horizons have been constructed by Lake . In the present case, Lake’s coordinates would essentially correspond with keeping only the last term inside brackets in Eq. (A.3). We prefer to use instead the function $`\psi `$ given by Eq. (A.3) because it allows us to define coordinates $`𝒰`$ and $`𝒱`$ that cover not only a region near the horizon, but the whole exterior spacetime. Thus, our coordinates $`𝒰`$ and $`𝒱`$ can be regarded as a deformation of Lake’s double-null coordinates, with the difference becoming relevant far from the horizon.

The next step is to use this result to calculate the asymptotic worldline of the centre of the star. First, we match the coordinates $`𝒰`$ and $`𝒱`$ onto interior null coordinates, say, $`U`$ and $`V`$. In particular, if two nearby outgoing rays differ by $`\mathrm{d}U`$ inside the star, then they will differ by $`\mathrm{d}U=\beta (𝒰)\mathrm{d}𝒰`$, with $`\beta `$ a regular function, outside. Similarly, $`\mathrm{d}V=\zeta (v)\mathrm{d}v`$, where $`\zeta `$ is another regular function. Then if $`v=\overline{v}`$ represents the last ray that passes through the center of the star before the formation of the horizon, to first order we have $`\mathrm{d}V=\zeta (\overline{v})\mathrm{d}v`$, where $`\zeta (\overline{v})`$ is constant.

We can write near the horizon

$$\mathrm{d}U=\beta (0)\frac{\mathrm{d}𝒰}{\mathrm{d}u}\mathrm{d}u.$$
(A.7)
Since for the center of the star $`\mathrm{d}U=\mathrm{d}V=\zeta (\overline{v})\mathrm{d}v`$, this immediately integrates to

$$\zeta (\overline{v})(v\overline{v})=\beta (0)𝒰(u)=\beta (0)\psi ^1\left(u\right)2\beta (0)\frac{M^2}{u}.$$
(A.8)
The last approximation follows from Eq. (A.3) where $`\xi \psi ^1(2M^2/\xi )`$ near the horizon.

Thus the late-time worldline for the center of the star is, finally, represented by the equation

$$v\overline{v}\frac{A}{u},u+\mathrm{},$$
(A.9)
where $`A=2\beta (0)M^2/\zeta (\overline{v})`$ is a positive constant that depends on the dynamics of collapse. This is a hyperbolic trajectory, equivalent to that of a uniformly accelerated mirror in two-dimensional Minkowski space. From this point on, one can use results already in the literature, which establish that the particle production in such a scenario is nonthermal.


text:
# The Redshift and the Host Galaxy of GRB 980613: A Gamma-Ray Burst From a Merger-Induced Starburst? 1footnote 11footnote 1Partially based on the observations obtained at the W. M. Keck Observatory which is operated by the California Association for Research in Astronomy, a scientific partnership among California Institute of Technology, the University of California and the National Aeronautics and Space Administration.

## 1. Introduction

Studies of the cosmic gamma-ray bursts (GRBs) are currently one of the most active and exciting fields of astrophysics. A deluge of important clues in the long-standing puzzle of GRBs began with the discovery of long-lived X-ray afterglows by the BeppoSAX satellite (Costa et al. 1997). This was followed by the discovery of optical (van Paradijs et al. (1997)) and radio afterglows (Frail et al. 1997), and then the first redshift measurement which demonstrated the cosmological nature of GRBs (Metzger et al. 1997). To date, about a dozen redshift measurements of GRBs have been obtained. Studies of afterglows confirmed the synchrotron shock model (e.g.Wijers, Rees & Mészáros (1997)), and their physics now seems to be reasonably well understood. For a recent review of observational results, see, e.g., Kulkarni et al. (2000).

What is still not known is what are both the nature of the progenitors and the trigger that causes GRBs. The two currently popular models include mergers of neutron stars or other compact stellar remnants, leading to a creation of a black hole, and a collapsar or hypernova model in which an explosion of a massive star produces a black hole remnant. In both cases, spin energy of a debris torus or of the black hole itself is used to power the GRB. While other models are still viable, there is now a growing evidence in favor of the collapsar/hypernova model, at least for the long-duration bursts detected by BeppoSAX. The most direct evidence favoring this was the probable detection of supernovæ associated with GRB 980326 (Bloom et al. (1999)) and GRB 970228 (Reichart (1999); Galama et al. (2000)), as well as the still controversial case of GRB 980425 and SN 1998bw (Galama et al. (1998); Kulkarni et al. (1998)). All models involving massive stars and their remnants suggest that GRBs should be closely related to the massive star formation in galaxies.

Study of GRB host galaxies can provide valuable clues which can constrain the models, and redshifts are necessary in order to derive the physical parameters of the GRBs itself, primarily the energy scale, as well as the hosts: their luminosities, star formation rates, morphology, locations of GRBs within them, etc.

In this $`Letter`$ we present deep imaging and spectroscopic observations of the host galaxy of GRB 980613. GRB 980613 was detected and localized by BeppoSAX (Piro, Costa et al. (1998); Piro et al. (1998)). Following several unsuccessful attempts an optical transient (OT) was discovered by Hjorth et al. (1998). The detection of the host galaxy was reported by Djorgovski et al. (1998a). Its redshift determination was first reported by Djorgovski et al. (1999), and is described in more detail here, with additional data.

## 2. Observations and Data Reductions

In what follows we assume the Galactic foreground extinction of $`A_V=0.27`$ mag (Schlegel, Finkbeiner & Davis (1998)), and use the Galactic extinction curve from Cardelli, Clayton, & Mathis (1998) with $`R=A_V/E_{BV}=3.1`$.

Our initial imaging date were obtained on the W. M. Keck Observatory 10-m telescope (Keck II) by Dr. H. Ebeling, on 16.30 June 1998 UT, in the $`R`$ band, using the Low-Resolution Imaging Spectrometer (LRIS; Oke et al. (1995)); two 300 s exposures were obtained. We reduced these images in the standard manner and confirmed that no cosmic rays affected the immediate area of the transient. After the transient had faded we reobserved the field in $`R`$ band (29 November 1998 UT; 900s), $`I`$ band (24 March 1999 UT; 1000s), and $`K`$ band (7 February 1999 UT; 2040s). The final combined images in each filter is shown in fig. 1. The morphology of the system surrounding the GRB is complex. We label 5 distinct components in fig. 2. Table 1 provides their magnitudes and offsets relative to the position of the OT.

We referenced our $`R`$ and $`I`$ band magnitudes using star 2 (Halpern et al. 1998; Diercks et al. 1998). For each component we used a circular aperature and performed an aperature correction using the curve-of-growth of star 2. The total magnitude of system A+B+C is $`R=23.11\pm 0.05`$. We estimate the zeropoint systematic uncertainty due to an uncertain color correction at 0.15 mag. The 2-$`\sigma `$ upper limit to a point source detection is $`R=25.9`$ mag and $`I=24.9`$ mag. Our $`K`$ band imaging was taken under photometric conditions and we used SJ 9134 (Persson et al. 1998) to obtain a zeropoint for the night. The $`K`$ band aperature magnitudes of the 5 components are listed in Table 1.

We registered the images from 16 June 1998 and 29 November 1998 $`R`$ band images using 7 unsaturated stellar objects common to both images within 30 arcsec from the host galaxy. We used the optimum-filter technique in the IRAF package CENTER to position the astrometric tie objects. We determined the rotation, shift, and relative scale of the two images using GEOMAP and then registered the early-time image to the late-time image. The peak of the OT and the putative host were then estimated. Including the peak center errors and the registration uncertainty we find the OT was offset from the brightest $`R`$ band peak (A) by 0.52′′$`\pm `$ 0.13′′ East and 0.83′′$`\pm `$ 0.14′′S. This amounts to a radial angular projected offset of 0.98′′$`\pm `$ 0.14′′, and we interpret component A as the most likely host galaxy of the GRB. The $`I`$ and $`K`$ band images of the host were also registered to the late-time $`R`$ band image. Fig. 2 shows a color composite 15.3′′$`\times `$ 15.3′′region around the GRB and its host using the late-time $`R`$, $`I`$, and $`K`$ band images. The ellipse is a 3-$`\sigma `$ error contour about the position of the GRB.

Spectra of the host galaxy were obtained using LRIS on 17 December 1998, 14 January 1999, and 16 February 1999 UT, all in good observing conditions. We used a 1.5 arcsec wide long slit, always at a position angle close to the parallactic. On December 14, we used a 300 lines mm-1 grating, giving an effective instrumental resolution FWHM $`16`$Åand an approximate wavelength coverage 3950–8950 Å, and obtained 6 exposures totaling 8300 s. On January 14 and February 16, we used a 600 lines mm-1 grating, giving an effective instrumental resolution FWHM $`8`$ Åand an approximate wavelength coverage 5970–8540 and 5700–$`8270`$Å, respectively; 3 exposures of 1800 s were obtained on each night. The object was dithered on the spectrograph slit by several arcsec between the exposures. Exposures of an internal flat-field lamp and arc lamps were obtained at comparable telescope pointings immediately following the target observations. Exposures of standard stars from Oke & Gunn (1983) and Massey et al. (1998) were obtained and used to measure the instrument response curve. We estimate the net flux zero-point uncertainty, including the slit losses, to be about 10–20%.

Wavelength solutions were obtained from arc lamps in the standard manner, and then a second-order correction was determined from the wavelengths of isolated strong night sky lines, and applied to the wavelength solutions. This procedure largely eliminates systematic errors due to the instrument flexure, and is necessary in order to combine the data obtained during separate nights. The final wavelength calibrations have the r.m.s. $`0.20.5`$ Å, as determined from the scatter of the night sky line centers. All spectra were then rebinned to a common wavelength scale with a sampling of 2.5 Å using a Gaussian with a $`\sigma =2.5`$ Å as the interpolating/weighting function. This is effectively a very conservative smoothing of the spectrum, since it is smaller than the instrumental resolution. Individual spectra were extracted and combined using a statistical weighting based on the signal-to-noise ratio determined from the data themselves (rather than by the exposure time).

The final combined spectrum of the galaxy is shown in fig. 3. A strong \[O II\] 3727 line emission is present, and a weaker \[Ne III\] 3869 line is also seen. The weighted mean redshift is $`z=1.0969\pm 0.0002`$. The observed continuum is moderately blue, with the spectrophotometric colors $`(BV)0.45`$ and $`(VR)0.35`$ mag.

The corrected \[O II\] 3727 line flux is $`(5.25\pm 0.15)\times 10^{17}`$ erg cm-2 s-1 Hz-1, and its observed equivalent width is $`W_\lambda =125\pm 5`$ Å, i.e., $`60\pm 2.4`$ Å in the restframe, which is on a high side for field galaxies at comparable redshifts, but not extraordinary (Hogg et al. (1998)). The \[Ne III\] 3869 line has a flux of about 10% of the \[O II\] 3727 line. The continuum flux at $`\lambda _{obs}=5871`$Å, corresponding to $`\lambda _{rest}=2800`$Å, is $`F_\nu =0.84\mu `$Jy, uncertain by a few percent (statistical). Additional spectrophotometric flux zero-point uncertainty is estimated to be $`1020`$%.

The continuum flux from our direct photometry in the $`I`$ band, which corresponds roughly to the restframe $`B`$ band, is $`F_\nu =1.15\mu `$Jy, uncertain by $`10`$%.

## 3. Discussion

We will assume a flat cosmology with $`H_0=65`$ km s-1 Mpc-1, $`\mathrm{\Omega }_M=0.3`$, and $`\mathrm{\Lambda }_0=0.7`$. For $`z=1.0969`$, the luminosity distance is $`2.461\times 10^{28}`$ cm, and 1 arcsec corresponds to 8.8 proper kpc in projection.

The gamma-ray fluence from this burst was $`(1.71\pm 0.25)\times 10^6`$ erg cm$`2`$ (Woods et al. (1998)). The corresponding isotropic $`\gamma `$-ray energy is $`E_{\gamma ,iso}=6.2\times 10^{51}`$ erg.

From the \[O II\] 3727 line flux, we derive the line luminosity $`L_{3727}=4.0\times 10^{41}`$ erg s-1. Using the star formation rate estimator from Kennicutt (1998), we derive the SFR $`5.6M_{}`$ yr-1. From the UV continuum luminosity at $`\lambda _{rest}=2800`$Å, following Madau, Pozzetti & Dickinson (1998), we derive SFR $`3.1M_{}`$ yr-1. The difference may be due to the internal reddening within the host galaxy, but we also note that the \[O II\] line estimator is more sensitive to the current or recent massive star formation than the UV continuum estimator. This is consistent with the presence of the \[Ne III\] 3869 line, also seen in spectra of some other GRB hosts (Bloom et al. (1998); Bloom, Djorgovski & Kulkarni (2001)) which may be indicative of a recent, very massive star formation. We are completely insensitive to any fully obscured star formation component, so that these numbers represent lower limits. The only GRB host galaxy with a higher unobscured SFR measured to date is the host of GRB 980703 (Djorgovski et al. (1998)).

From the observed continuum flux in the restframe $`B`$ band, we derive $`M_B19.85`$ mag, which is about 1 mag fainter than the present-day $`L_{}`$ galaxy. Considering that an average galaxy at this redshift may be $`0.51`$ mag brighter than today due to normal evolution effects, we conclude that this galaxy (A) is moderately underluminous. The star formation rate per unit mass is thus fairly high, consistent with the large equivalent width of the \[O II\] line. We thus conclude that this galaxy is undergoing a starburst, albeit mild.

However, the most interesting feature may be the morphology of the entire system (A–E), which is highly suggestive of an interaction or early stages of a merger of at least two galaxies, some of which may be partly obscured. Redshifts of all 5 components are necessary to really test this interpretation; we note, however, that we detect a weak \[O II\] line emission coincident with the component C in our long-slit spectra, at the same redshift as the component A. Another possibility is that we are seeing a chance superposition of galaxies at very different redshifts, e.g., that the IR-bright components C and D may be unrelated to A (and perhaps also B and E), and if they are in the foreground, some gravitational lensing may be involved. However, we consider this less likely than the interaction/merger hypothesis.

The rate of galaxy interactions and mergers increases sharply with redshift (cf., Le Fèvre et al. (2000)), and finding a strongly interacting system at $`z1.1`$ is not surprising. The projected separations of components A–E (see fig. 2) are typical of intergalactic separations where dark halos are overlapping and the tidal friction inevitably leads to a merger.

The very red colors of the components C and D, $`(RK)4.4`$ and $`>5.8`$, respectively, is naturally explained due to obscuration by dust, which is consistent with the merger hypothesis. Alternatively, they could be passively evolving ellipticals which formed at a very high redshift, which requires some fine-tuning of model parameters.

It is now well established that mergers of gas-rich galaxies lead to bursts of star formation. This may be the strongest case so far for a GRB from such a system. The only other published case to date, showing morphology consistent with a mild tidal interaction is the host galaxy of GRB 990123 (e.g., Bloom et al. (1999); Fruchter et al. (1999); Holland & Hjorth (1999)). This provides a further evidence in support of the connection of GRBs with massive star formation.

We want to thank the staff of the W. M. Keck Observatory for their expert assistance, and to Dr. Harald Ebeling for obtaining the early images of the field. We thank Dana Sadava for help in $`K`$ band imaging reductions and Titus Galama for comments. JSB gratefully acknowledges the fellowship from the Fannie and John Hertz Foundation. SGD acknowledges partial funding from the Bressler Foundation. This work was supported in part by grants from the NSF and NASA to SRK.


text:
# Results from pQCD for A+A collisions at RHIC & LHC energies

## Abstract

This talk will discuss how to compute initial quantites in heavy ion collisions at RHIC (200 AGeV) and at LHC (5500 AGeV) using perturbative QCD (pQCD) by including the next-to-leading order (NLO) corrections and a dynamical determination of the dominant physical scale. The initial numbers are converted into final ones by assuming kinetic thermalization and adiabatic expansion.

The whole heavy ion physics community has entered into exciting and intense period as the first results from RHIC are expected to appear, and the era of LHC seems no more very distant. However, despite the solid status of QCD as the theory of strongly interacting matter, many uncertainties remain in the predictions for the outcome of the current experiments and various methods have been applied. Further experimental data will (hopefully) single out the best candidates for the correct approach.

As the collision energy is increased from that of the SPS, larger intrinsic scales are generated and the applicability of perturbative QCD (pQCD) becomes possible. The initial particle production is expected to be dominated by minijets, i.e. partons with $`p_T1\mathrm{}2`$GeV$`\mathrm{\Lambda }_{QCD}`$ BM87 . By assuming independent multiple semi-hard parton-parton collisions the average energy carried by minijets with $`p_Tp_0`$ at the rapidity interval $`\mathrm{\Delta }Y`$ in a central $`𝐛=0`$ $`AA`$ collision is given in leading-order (LO) by EKL89

$`\overline{E}_{AA}^T`$ $`=`$ $`T_{AA}(\mathrm{𝟎})\sigma E_T`$ (1)
$`=`$ $`T_{AA}(\mathrm{𝟎}){\displaystyle \underset{q,\overline{q},g}{}}{\displaystyle _{p_0,\mathrm{\Delta }Y}}𝑑p_t𝑑y_1𝑑y_2x_1f_{i/p}(x_1,Q^2)x_2f_{j/p}(x_2,Q^2){\displaystyle \frac{d\widehat{\sigma }^{ijjk}}{d\widehat{t}}}p_T,`$
where $`T_{AA}(𝐛)`$ is the standard nuclear overlap function and $`p_0`$ is the smallest transverse momentum scale to be considered. Collinear factorization is assumed to hold and any effects beyond it are neglected. This perturbative minijet approach suffers from two major sources of uncertainties:
(i) The next-to-leading order (NLO) corrections, which have not been known prior to LO89 ; ET00 , but rather have been simulated by ad hoc $`K`$-factors.
(ii) The determination of which value of $`p_0`$ to use at, say, RHIC or LHC energies. Should one stick to some constant universal value of $`p_0`$ or will this parameter possess some nontrivial $`\sqrt{s}`$\- and $`A`$-dependence?

In the following we will provide answers to both of these questions and combine them to obtain numerical estimates of average transverse energies and charged particle multiplicities at RHIC and LHC energies.

As is evident from formula (1), to deal with the first uncertainty, the relevant quantity to compute in NLO pQCD is the $`\sigma E_T`$, the first moment of the perturbative $`E_T`$ distribution in a pp-collision. The infrared safe NLO computation of this quantity has been presented in ET00 , where the computation was formulated via the subtraction algorithm of S. Ellis, Z. Kunszt and D. Soper KS92 .

The $`E_T`$ in central rapidity region is defined to be the total $`p_T`$ entering this region and originating from hard subprocessess in which at least an amount of $`2p_0`$ of transverse momentum is released. The numerical results are shown in fig. 1(a).

To analyze the implications of these numbers, let us define two $`K`$-factors: $`K=\text{(full NLO)}/\text{LO}`$ and $`K^{}=\text{(full NLO)}/\text{LO}^{}`$, the first of these measuring the deviation of the NLO results from the consistent LO calculation and the latter measuring rather the relative difference between two subsequent terms of the perturbation series. The magnitude of these $`K`$-factors is due to a new kinematical region which manifests itself only at NLO. If one rejects this new domain, the resulting $`K`$-factors would be a factor of two smaller than those obtained from fig. 1(a), but then a significant amount of perturbatively calculable $`E_T`$ would be neglected.

For a detailed description of the calculation and issuses such as the scale choices, see ref. ET00  and references therein. For the purposes to be considered here, the sufficient observation is that we now have control over the magnitude of NLO corrections which are stable relative to LO results even at few GeV scales, thus signalling the applicability of pQCD in this domain.

Turning to the uncertainty (ii), then, it is clear that as $`p_0`$ is decreased the cross section, as well as the uncertainty, grows. At certain value of saturation, $`p_0=p_{\mathrm{sat}}`$ the system becomes very dense and new physics enters GRL83 . For large nuclei and large collision energies this may happen already in the perturbative domain, which we concluded on the basis of the NLO analysis to include also the few GeV region. Then the corresponding values for the number of particles as well as for the amount of $`E_T`$ are easily produced via the perturbative computation at $`p_0=p_{\mathrm{sat}}`$, which effectively accounts for the contributions of all scales, since the partons with $`p_Tp_{\mathrm{sat}}`$ are rare and those with $`p_Tp_{\mathrm{sat}}`$, altough numerous, contribute negligibly to total $`E_T`$.

Various ways to determine the actual magnitude of $`p_{\mathrm{sat}}`$ can be conjectured. A simple geometric criterion has been presented in EKRT99 . This is based on the idea that if one assigns an effective area $`\pi /p_0^2`$ to each gluon produced, then at certain value of $`p_0`$ the total area of $`N_{AA}(p_0,\sqrt{s}\mathrm{\Delta }Y)`$ gluons produced will exceed the effective transverse area $`\pi R_A^2`$ of the nucleus. Therefore one can iterate the equation $`N_{AA}=p_{\mathrm{sat}}^2R_A^2`$ to determine $`p_{\mathrm{sat}}`$ for given $`A`$ and $`\sqrt{s}`$, see fig.1(b). On the basis of the NLO analysis, we take here $`K=2`$ to account for the NLO corrections and also implement nuclear shadowing via the EKS98 parametrization EKS98 .

All the initial quantities then computed can well be fitted by a scaling law of a type $`CA^b(\sqrt{s})^b`$. In particular one finds that

$`p_{\mathrm{sat}}/\mathrm{GeV}`$ $`=`$ $`0.208A^{0.128}(\sqrt{s})^{0.191}`$
$`ϵ_i/(\mathrm{GeV}/\mathrm{fm}^3)`$ $`=`$ $`0.103A^{0.504}(\sqrt{s})^{0.786}`$ (2)
$`n_i\mathrm{fm}^3`$ $`=`$ $`0.370A^{0.383}(\sqrt{s})^{0.574}`$
where the particle and energy densities are evaluated at $`p_{\mathrm{sat}}`$, and the whole production process is then considered to take place at $`\tau _0=1/p_{\mathrm{sat}}`$.

These, however, are just the initial numbers at 0.2 (0.1) fm/c at RHIC (LHC), and the major problem is how to get from these to the ones at later instants and to finally arrive at experimentally visible quantities. Let us therefore assume, not completely without reason (see EKRT99 ), that the system is initially thermalized in the sense that it possesses a correct ratio of energy per particle as far as the dominant gluonic particle content is considered, and expands conserving the total entropy $`S3.6N_i`$. As the final particles consist dominantly of pions, for which $`S4N_f`$, we find that $`N_f=0.9N_i`$. As the initial volume is $`V_i=\pi R_A^2\mathrm{\Delta }Y/p_{\mathrm{sat}}`$, formulae (2) give

$$N_f=1.245A^{0.92}(\sqrt{s})^{0.383}.$$
(3)
After the conference the very first measurements by PHOBOS collaboration have been announced PHOBOS . According to them, the charged particle multiplicity at midrapidity is $`dN/d\eta =408\pm 12\pm 30`$ at 56 AGeV and $`555\pm 12\pm 35`$ at 130 AGeV. The framework described here gives $`N_{ch}=2/3N_f=370`$ and 530 per unit $`\eta `$ respectively, when taking into account that number of participants was reported to be 330 for $`\sqrt{s}=56`$ AGeV and 343 for $`\sqrt{s}=130`$ AGeV, and that $`dN/dy=1.15dN/d\eta `$.

The final $`E_T`$ in this scenario is obtained by means of hydrodynamics EKRT99 ; KRMcLG  as

$$E_T=N_f\times [0.39+0.061\mathrm{ln}(N_f/A)].$$
(4)
Numerical values then per unit $`\eta `$ are $`E_T=260`$ GeV for $`\sqrt{s}=56`$ AGeV and $`E_T=390`$ GeV for $`\sqrt{s}=130`$ AGeV using the multiplicites computed above and the quoted participant numbers. Measurements of these final transverse energies are awaited to appear soon. These measurements will then allow us to draw conclusions on the issues such as the true degree of thermalization in the system.

Acknowledgements: I wish to thank K.J. Eskola, K. Kajantie and P.V. Ruuskanen for collaboration. The financial support from the Academy of Finland as well as from the organizers of the CIPANP2000 conference is gratefully acknowledged.


text:
# Toward the multiphoton parametric oscillators

## I The problem

Optical parametric oscillators (OPO’s) based on processes of down-conversion in a cavity have proved to be efficient sources of light with a range of unique properties. Among them it should be noted the nonclassical nature of light, including quantum correlations, sub-Poissonian statistics and a larger amount of squeezing . OPO provides a method for producing tunable radiation from a coherent fixed-frequency pump source in the intracavity down-conversion, where a pump-photon splits into a pair of subharmonic photons. So far OPO has only been realized experimentally for two-photon down-conversion in nonlinear media with second-order $`\chi ^{(2)}`$ susceptibility. Experimental and theoretical results can be found in .

In this letter we propose a new conception for the realization of multiphoton parametric oscillators using the idea of composite multiphoton interaction. Two schemes of OPO based accordingly on a three-photon and four-photon intracavity down-conversion are presented.

One of the principal motivations for the investigation of multiphoton OPO concerns to a production of nonclassical states of light with the novel features. Other motivation is their possible exploitation for quantum information technologies, because sources offering a great variety of multiparticle entangled states, first of all meaning photons, are required for the implementation of many quantum communication processes such as direct transmission of the state, quantum teleportation, quantum cryptography, and entanglement swapping .

It is known that two-photon down-conversion so far is the most standard source of pair of polarization-entangled photons which is relevant to many problems of quantum information (see, previous experiments on photon splitting in , and for polarization-entangled photon pairs ). Although two-particle entanglement has long been demonstrated experimentally the realization of entanglement between three or more particles has been considered as a difficult problem. Among the possible sources of three-photon states (so-called Greenberger-Horne-Zeilinger states ) we mention the process of degenerate three-photon down-conversion $`\omega \omega /3+\omega /3+\omega /3`$ in a $`\chi ^{(3)}`$ medium. This process has been studied for running waves . Note, that most recently the polarization entanglement for three spatially separated photons has been observed . Three-photon down-conversion in an optical cavity has also been proposed and investigated in , but too weak interaction between photons in $`\chi ^{(3)}`$ medium makes experimental realization of this scheme a difficult problem. As regards to a four-photon splitting $`\omega \omega /4+\omega /4+\omega /4+\omega /4`$ in a $`\chi ^{(4)}`$ medium, as we know, this process was not considered yet.

This letter reports the analysis of two schemes of multiphoton OPO in a double resonant-mode cavity involving two crystals of second-order susceptibilities. They are the following:

### A  Four-photon splitting in cascaded down-conversion

The scheme proposed contains two sequential nonlinear processes of degenerate two-photon splitting in $`\chi ^{(2)}`$ media taking place inside the same cavity. This scheme for the special configuration when the coupling in- and out-fields occurs at one of the ring cavity mirrors is shown in Fig.1. The pump frequency $`\omega `$ is converted to the subharmonic frequency $`\omega _2=\omega /2`$ and then into the frequency $`\omega _1=\omega /4`$ in the cascaded parametric processes $`\omega \omega _2+\omega _2`$ and $`\omega _2\omega _1+\omega _1.`$ We assume that the each of these processes takes place effectively only in the definite nonlinear crystal. The cavity is also assumed transparent at the pump frequency. So, we neglect the pump depletion effects considering the amplitude of driving field as a classical constant. We model the combination of these processes by the interaction Hamiltonian:

$$H_1=i\mathrm{}\chi _1(Ea_2^2E^{}a_2^2)+i\mathrm{}k_1(a_1^2a_2a_1^2a_2^{}),$$
(1)
where $`a_1`$ and $`a_2`$ are the operators of the modes $`\omega _1`$ and $`\omega _2,`$ $`\chi _1`$ and $`k_1`$ are the nonlinear coupling constants for the processes $`\omega \omega _2+\omega _2`$ and $`\omega _2\omega _1+\omega _1`$ respectively which are related to second-order nonlinear susceptibilities $`\chi ^{(2)}`$. $`E`$ is proportional to the coherent driving field amplitude.

### B Three-photon splitting in cascaded down-conversion

This system mainly is similar to the previous one. The only difference is that the pump frequency $`\omega `$ converts to two different frequencies$`\omega _1=\omega /3`$ and $`\omega _2=2\omega /3`$ in the parametric interaction $`\omega =\omega _1+\omega _2.`$ The subharmonic $`\omega _1`$ is also produced in the other process of down conversion $`\omega _2=\omega _1+\omega _1.`$ On the whole two subharmonic modes $`\omega _1=\omega /3`$ and $`\omega _2=2\omega /3`$ are created in the cavity due to the cascaded processes. In the undepleted pump field configuration this model is described by the following interaction Hamiltonian:

$$H_1=i\mathrm{}\chi _2(Ea_1^{}a_2^{}E^{}a_1a_2)+i\mathrm{}k_2(a_1^2a_2a_1^2a_2^{}),$$
(2)
where $`a_1`$ and $`a_2`$ are the operators of $`\omega _1`$ and $`\omega _2`$ subharmonics, and $`\chi _2`$ and $`k_2`$ are the coupling constants of the processes $`\omega \omega _1+\omega _2`$ and $`\omega _2\omega _1+\omega _1,`$ respectively.

We note, that optical schemes involving cascaded nonlinearities in a cavity are currently of considerable interest and hold promise for applications including, in particular, frequency tunable sources of light, large third-order nonlinear effects via cascaded second-order nonlinearities and squeezing phenomena. The references and experimental observation can be found in . Possible operational regimes of a combined system, as a rule, essentially differ from those in pure processes in each of nonlinear media. Some results in this direction are devoted to the cascaded frequency doubler . Concerning our schemes, it is interesting to note, that the stable stationary above-threshold regime of oscillation is realized even in undepleted pump regime. The other novelty here is comparatively low values of the threshold of generation due to the schemes proposed can be accessible for experiments. Contrary to polarization-entangled state generation in two- and three-photon devices  the basic multiphoton states of our systems are neither based on polarization nor on momentum, but on phase-space entanglement of modes.

## II Four-photon OPO

The system of interest is dissipative because the subharmonic modes suffer losses due to partially transmission of light through the mirrors of the cavity. So that the reduced density operator $`\rho `$ of both modes obeys a master equation

$$\frac{\rho }{t}=\frac{1}{i\mathrm{}}[H_1,\rho ]+\underset{i=1,2}{}\gamma _i(2a_i\rho a_i^{}a_i^{}a_i\rho \rho a_i^{}a_i),$$
(3)
where $`\gamma _1`$ and $`\gamma _2`$ are the cavity damping rates for the modes $`\omega _1`$ and $`\omega _2.`$

To proceed with further analysis we present now the semiclassical steady-state solutions and stability properties of the system proposed. In the classical limit the system may be described by the complex-field amplitudes $`\alpha _1`$ and $`\alpha _2`$ of the resonant modes $`\omega _1`$ and $`\omega _2`$ whose real and imaginary parts respectively represent the dimensionless position and momentum. They obey the equation of motion:

$`{\displaystyle \frac{\alpha _1}{t}}`$ $`=`$ $`\gamma _1\alpha _1+2k_1\alpha _2\alpha _1^{},`$ (4)
$`{\displaystyle \frac{\alpha _2}{t}}`$ $`=`$ $`\gamma _2\alpha _2+2\chi _1\alpha _2^{}Ek_1\alpha _1^2.`$ (5)
Solving these equations for the steady-state regime and carrying out the standard linearized stability analysis one can arrive at the following results.

The trivial zero-amplitude solution $`\alpha _1=\alpha _2=0`$ is stable in the region $`E<E_{th},`$ $`E_{th}=\gamma _2/2\chi _1,`$ with $`E_{th}`$ being the threshold value of the pump field amplitude, and it describes the below-threshold regime of oscillation of both subharmonics. The stable above-threshold solutions, for $`E>E_{th}`$, expressed in terms of photon numbers $`n_j`$ and phases $`\phi _j`$ of subharmonic modes ($`\alpha _j=\mathrm{exp}(i\phi _j)\sqrt{n_j},j=1,2`$) are calculated as:

$`n_1`$ $`=`$ $`\gamma _1\gamma _2(E/E_{th}1)/2k_1^2,n_2=\gamma _1^2/4k_1^2,`$ (6)
$`\phi _1`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Phi }}{4}}\pm {\displaystyle \frac{\pi }{2}}m,\phi _2={\displaystyle \frac{\mathrm{\Phi }}{2}}\pm \pi n,n,m=0,1,2\mathrm{},`$ (7)
where $`\mathrm{\Phi }`$ is the phase of the driving field.

So, this scheme of multiphoton OPO allows us to reach stable generation of subharmonics at frequencies $`\omega _1=\omega /4`$ and $`\omega _2=\omega /2`$ at $`E>E_{th}`$ in undepleted pump-field regime. It should be pointed out that the occurrence of such a regime of oscillation becomes possible just due to the combination of two kinds of down-conversion processes in a cavity. One of the peculiarities of cascading dynamics is that the threshold pump power $`P_{th}=c\mathrm{}\omega E_{th}^2/2L`$ ($`L`$ is the optical way of pump field) is only related to the coupling constant $`\chi _1`$ of $`\omega \omega /2+\omega /2`$ parametric process. This coupling constant is proportional to the second-order susceptibility that leads to the comparatively low value of the threshold. We illustrate this fact for the system depicted in Fig.1, choosing realistic losses and typical values of $`\chi _1`$ and $`k_1`$ (see detailed calculations for a similar cavity configuration in ). Assuming $`\chi _1=`$ $`k_1=3.4\times 10^4s^1`$ , $`\gamma _1=2.4\times 10^8s^1`$ , $`\gamma _2=3.6\times 10^8s^1`$, we find for the threshold pump power $`P_{th}=74mW`$. For these parameters the power of the cavity output field of $`\omega /4`$ mode increases as $`P_1^{out}=2\mathrm{}\frac{\omega }{4}\gamma _1n_1=3.4\times 10^3(E/E_{th}1)W`$, while the output power of the $`\omega /2`$ subharmonic remains constant $`P_2^{out}=2\mathrm{}\frac{\omega }{2}\gamma _2n_2=3.4\times 10^3W.`$

A very important characteristic of multiphoton OPO concerning phase information of the generated modes, is displayed in phase-space. As we see from (7), the stationary states of $`\omega /4`$ subharmonic have four-fold symmetry. There are four stable states of $`\omega _1`$ mode with equal photon numbers, but with four different phases which are: $`\frac{\mathrm{\Phi }}{4},\frac{\mathrm{\Phi }}{4}+\frac{\pi }{2},\frac{\mathrm{\Phi }}{4}+\pi ,\frac{\mathrm{\Phi }}{4}+\frac{3\pi }{2}.`$ As regards to $`\omega /2`$ subharmonic mode, it displays two-fold symmetry in phase-space, i.e. there exist two states with equal intensities but with different phases $`\frac{\mathrm{\Phi }}{2},\frac{\mathrm{\Phi }}{2}+\pi .`$

Let us now show the general character of this symmetry in the frame of a full quantum treatment of interaction. It is easy to verify that the Hamiltonian (1) as well as the density operator $`\rho `$ obeying Eq.(3) satisfy the commutation relations

$$[H_1,U]=[\rho (t),U]=0$$
(8)
with the operator

$`U=\mathrm{exp}\left(i{\displaystyle \frac{\pi }{2}}a_1^{}a_1+i\pi a_2^{}a_2\right).`$

For the reduced density operators of each of the modes $`\rho _{1(2)}=Tr_{2(1)}(\rho ),`$ (which are constructed by tracing over one of the modes), relations (8) give

$$[\rho _1(t),U_1(\frac{\pi }{2})]=[\rho _2(t),U_2(\pi )]=0,$$
(9)
where operators$`U_i(\phi )=\mathrm{exp}\left(i\phi a_i^{}a_i\right),i=1,2`$ perform rotation by the angle $`\phi `$ around the origin in phase-space of complex stochastic variables $`\alpha _1`$ and $`\alpha _2`$, respectively.

One of the most important conclusions of such symmetries is related to the Wigner functions $`W_1`$ and $`W_2`$ of the subharmonic modes $`\omega _1`$ and $`\omega _2`$ which provide a large amount of information about the states of the mode and also provide a pictorial view. To see this we address to the formula for Wigner function expressed through the density operator. Using also formulas (9) we arrive to the symmetries of Wigner functions written through the stochastic mode amplitudes $`\alpha _1`$ and $`\alpha _2`$ corresponding the operators $`a_1`$ and $`a_2.`$ In the polar coordinates $`r,\theta `$ of the complex phase-space plane $`X=\frac{\alpha +\alpha ^{}}{2}=r\mathrm{cos}\theta ,`$ $`Y=\frac{\alpha \alpha ^{}}{2i}=r\mathrm{sin}\theta `$ these symmetries take the forms:

$$W_1(r,\theta +\pi /2)=W_1(r,\theta ),W_2(r,\theta +\pi )=W_2(r,\theta ).$$
(10)
We have obtained this property of Wigner functions using the symmetry of Hamiltonian (1) and dissipation terms in master equation (3), assuming that the initial state of the modes is a vacuum state. Below, we will present the other confirmation to phase-space symmetry in quantum trajectories.

## III Three-photon OPO

The dynamics of subharmonic modes $`\omega _1=\omega /3`$ and $`\omega _2=2\omega /3`$ for this system is governed by master equation (3) involving Hamiltonian (2) instead of Hamiltonian (1). Accordingly, the classical equations of motion for the complex amplitudes are

$`{\displaystyle \frac{\alpha _1}{t}}`$ $`=`$ $`\gamma _1\alpha _1+\chi _2\alpha _2^{}E+2k_2\alpha _2\alpha _1^{},`$ (11)
$`{\displaystyle \frac{\alpha _2}{t}}`$ $`=`$ $`\gamma _2\alpha _2+\chi _2\alpha _1^{}Ek_2\alpha _1^2.`$ (12)
We list below the semiclassical results for stable values of photon numbers and phases of the modes. In the above-threshold generation regime, at $`E>E_{th}`$, where

$$E_{th}=\frac{2\sqrt{2\gamma _1\gamma _2}}{3\chi _2}$$
(13)
is the threshold value of E, the results take the form:

$`n_1`$ $`=`$ $`{\displaystyle \frac{\gamma _1\gamma _2}{18k_2^2}}\left(\epsilon +3\sqrt{\epsilon ^21}\right)^2,`$ (14)
$`n_2`$ $`=`$ $`{\displaystyle \frac{\gamma _1^2}{36k_2^2}}\left({\displaystyle \frac{\epsilon +3\sqrt{\epsilon ^21}}{\epsilon +\sqrt{\epsilon ^21}}}\right)^2,`$ (15)
where $`\epsilon =E/E_{th},`$

$$\varphi _1=\frac{\mathrm{\Phi }}{3}+\frac{2\pi }{3}n,\varphi _2=\frac{2\mathrm{\Phi }}{3}\frac{2\pi }{3}n,n=0,1,2.$$
(16)
In the region $`\epsilon <\frac{3}{2\sqrt{2}}`$ the stability condition is fulfilled only for a zero amplitude steady-state solution $`\alpha _1=\alpha _2=0,`$ i.e. the subharmonic field excitation exists at the spontaneous noise level. The photon numbers are plotted against the pump amplitude in Fig.2. They show bistable hysteresis-cycle behavior in a small domain $`1<\epsilon <\frac{3}{2\sqrt{2}}.`$ The dashed parts of the curves correspond to unstable solutions. As we see, with increasing $`E/E_{th}`$ the photon number of the $`\omega /3`$ subharmonic increases as $`n_1\frac{8\gamma _1\gamma _2}{9k_2^2}\epsilon ^2`$, while the photon number of the $`2\omega /3`$ subharmonic goes to the saturated value $`n_2=\gamma _1^2/9k_2^2.`$

The very important peculiarity of the system proposed is that the threshold value $`E_{th}`$ depends on the coupling constant $`\chi _2`$ which is related to the second-order $`\chi ^{(2)}`$ susceptibility in contrary to the scheme cited in , where the threshold of the driving field amplitude is determined by $`\chi ^{(3)}`$ susceptibility. However, the photon numbers of the generated modes depend on both nonlinear coupling constants $`\chi _2`$ and $`k_2`$. For the realistic parameters: $`k_2=2.4\times 10^3s^1,`$ $`\chi _2=2.3\times 10^3s^1,`$ $`\gamma _1=0.9\times 10^8s^1`$, $`\gamma _2=1.2\times 10^8s^1`$ parametric oscillations occur above the threshold pump power $`P_{th}=3W`$ and the threshold values of cavity-output field powers are $`P_1^{th}4mW,`$ $`P_2^{th}4mW`$. Eqs.(15),(16) imply that in the above-bistability regime there exist three states for each of the modes which have equal intensities $`n_1`$ and $`n_2`$ but different phases: $`\frac{\mathrm{\Phi }}{3},\frac{\mathrm{\Phi }}{3}+\frac{2\pi }{3},\frac{\mathrm{\Phi }}{3}+\frac{4\pi }{3}`$ and $`\frac{2\mathrm{\Phi }}{3},\frac{2\mathrm{\Phi }}{3}\frac{2\pi }{3},\frac{2\mathrm{\Phi }}{3}\frac{4\pi }{3},`$ respectively. Surprisingly, both the modes $`\omega _1=\omega /3`$ and $`\omega _2=2\omega /3`$ have the identical threefold symmetry which also follows from the specific three-photon form of Hamiltonian (2). In the full quantum treatment of cascaded three-photon down-conversion this symmetry, in particular, is displayed in the Wigner functions. Following the method presented above, we have obtained:

$$W_i(r,\theta +2\pi /3)=W_i(r,\theta ),i=1,2.$$
(17)
## IV Quantum simulation: multistabilities and threshold behavior

The second part of this paper is devoted to the study of quantum-statistical properties of cascaded multiphoton OPO in a dissipative cavity. Our aim is to analyze the phase-space properties of the subharmonic modes in a quantum regime in the presence of quantum noise on the basis of the Wigner functions.

Note that an exact analysis of quantum optical nonlinear systems including dissipation effects is a difficult problem which was solved for a few single models (see  and references therein). In this field of research, analytical expression for the Wigner function of two-photon OPO was derived in . However the quantum-statistical properties of cascaded OPO’s can not be analyzed using single analytical formulas. Our approach is based on the quantum-jump simulation method also known as the state-vector Monte-Carlo method . This method considers not the density matrix but deals with the state vector $`\mathrm{\Psi }^\alpha (t)`$ which is a member of the ensemble of state vectors and index $`\alpha `$ indicates the realizations. The procedure adopted in the quantum-jump simulation for the system Hamiltonians (1) and (2) consists in the following. The evolution of a member of the ensemble of pure states $`\mathrm{\Psi }^\alpha (t)`$ over a short time $`\delta t`$ is governed by the non-Hermitian effective Hamiltonian if there is no quantum jump

$$H_{1(2)eff}=H_{1(2)}i\mathrm{}(\gamma _1a_1^{}a_1+\gamma _2a_2^{}a_2).$$
(18)
Such evolution must be completed by the possibility of quantum jumps which change the state vectors of both modes. On the whole the state transforms in the following way

$$\mathrm{\Psi }^\alpha (t+\delta t)\{\begin{array}{cc}\left(1\frac{i}{\mathrm{}}\delta tH_{1(2)eff}\right)\mathrm{\Psi }^\alpha (t)/(1\delta p)^{1/2},\hfill & \text{with probability }1\delta p,\hfill \\ \sqrt{2\gamma _i}a_i\mathrm{\Psi }^\alpha (t)/(\delta p_i/\delta t)^{1/2},\hfill & \text{with probability }\delta p_i,i=1,2,\hfill \end{array}$$
(19)
where $`\delta p=\delta p_1+\delta p_2.`$

We generalize the method presented in  for the quantum-jump simulation of the Wigner functions. Details of analogous calculations for the second harmonic generation can be found in .The numerical simulations are performed in the truncated Fock basis of both subharmonic modes $`\omega _1`$ and $`\omega _2`$ in the regime of strong nonlinear coupling between the subharmonic modes, i.e. $`\chi _{1,2}\gamma _1,\gamma _2`$ and $`k_{1,2}\gamma _1,\gamma _2`$. We note that the OPO with such extremely large nonlinearities are not realized in practice and in this part of the Letter we do not intend to give results close to an experimental situation but discuss the fundamental problems of multiphoton entangled states.

To be short, we shall give the results for three-photon OPO, where one kind of phase-space symmetries is realized. The steady-state Wigner functions for $`\omega _1=\omega /3`$ and $`\omega _2=2\omega /3`$ modes are calculated in the basis of $`n=45`$ photon number states for each mode. The results may be written as depending on the dimensionless parameters: $`\epsilon =E/E_{th},k_2/\gamma _1,\gamma _2/\gamma _1.`$ We calculate the Wigner functions of $`\omega _1=\omega /3`$ subharmonic averaged over 1000 trajectories in the complex phase–space plane $`X,Y.`$ It is obvious that below the threshold the Wigner function is single-humped and centered at $`X=Y=0`$. This state indicates appearance of three arms as the consequence of small entanglement between the modes. Increasing the driving field $`E`$ we enter the bistability domain and observe the occurrence of three additional side-humps. They correspond to the above-threshold steady-states (12), (13) with equal photon numbers and threefold symmetric phases. We see explicitly that in this domain the behavior of our system is actually four-stable. With further increase of $`E/E_{th}`$ the central hump disappears, while the side-humps increase and we turn to a manifestly above threshold oscillation regime when we observe the phase space three-stability. This case is shown in Fig.3 for the parameters $`k_2/\gamma _1=0.2,\gamma _2/\gamma _1=0.4,\epsilon =1.59`$. Analogous behavior has the Wigner function of $`2\omega /3`$-subharmonic.. The difference between the results of two modes is displayed in the form of squeezing the humps.

## V CONCLUSION

Thus we have shown that the idea of using the cascaded photon splitting processes in a cavity opens a novel route for the realization of a multiphoton parametric oscillator. Two devices are proposed which transform the pump coherent field into two subharmonics at frequencies $`\omega /2`$ and $`\omega /4`$ (configuration A) and $`\omega /3`$ and $`2\omega /3`$ (configuration B). Furthermore, these models of OPO are to be subject to an experiment. We have been convinced that the pump threshold of each of the OPO is only expressed through the second-order susceptibility. As a result, comparatively low values of thresholds have been obtained for typical parameters. The important features of multiphoton OPO obtained are phase-space symmetries and multistabilities of subharmonics. In addition, it is shown that subharmonics ”remember” the phase of the pump field. Our numerical quantum Monte-Carlo simulations display phase-space multistabilities and threshold behavior of subharmonics in the presence of dissipation.

Acknowledgments

We acknowledge the helpful discussions with J.Bergou, H.Carmichael, S.Fauve and L.Manukyan. The work was supported in part by INTAS grant No 97–1672, and by grant No 00375 awarded by the Armenian Science Foundation.

Figure Captions

Fig.1. Principal scheme of the cascaded parametric oscillator. The cavity is resonant at two frequencies $`\omega _1=\omega /4`$ and $`\omega _2=\omega /2`$ of the subharmonics. The phase matching condition for the process $`\omega \omega /2+\omega /2`$ is satisfied in the first medium, and for the process $`\omega /2\omega /4+\omega /4`$ in the second medium.

Fig.2. The normalized mean photon numbers $`n_1/n_1^{th}`$ (curve1) and $`n_2/n_2^{th}`$ (curve 2) of $`\omega /3`$ and $`2\omega /3`$ modes versus the scaled pump intensity parameter $`\epsilon ^2=E^2/E_{th}^2`$. The threshold values of the photon numbers are $`n_1^{th}=\frac{\gamma _1\gamma _2}{18k_2^2}`$ and $`n_1^{th}=\frac{\gamma _1^2}{16k_2^2}`$. The dashed parts of the curves 1,2, and zero-amplitude (curve 3) solutions for both modes describe the unstable steady-state solutions.

Fig.3. The steady-state Wigner function of $`\omega /3`$ mode which demonstrate the phase-space multistability in the above-threshold range.