text:
# On anharmonicities of giant dipole excitations

## Abstract

The role of anharmonic effects on the excitation of the double giant dipole resonance is investigated in a simple macroscopic model. Perturbation theory is used to find energies and wavefunctions of the anharmonic oscillator. The cross sections for the electromagnetic excitation of the one- and two-phonon giant dipole resonances in energetic heavy ion collisions are then evaluated through a semiclassical coupled-channel calculation. It is argued that the variations of the strength of the anharmonic potential should be combined with appropriate changes in the oscillator frequency, in order to keep the giant dipole resonance energy consistent with the experimental value. When this is taken into account, the effects of anharmonicities on the double giant dipole resonance excitation probabilities are small and cannot account for the well known discrepancy between theory and experiment.

The double giant dipole resonance (DGDR) has attracted considerable interest in the last decade. Several experiments to measure the DGDR cross section using relativistic heavy ion beams have been performed. Comparison with the predictions of the harmonic oscillator model has clearly demonstrated a systematic discrepancy. The experimental values for the DGDR cross sections exceed the theoretical predictions by a considerable amount. One of the attempts to explain these differences was made by Bortignon and Dasso, using a macroscopic anharmonic oscillator model. These authors found that with a small anharmonic perturbation of the $`r^4`$-type one can reproduce both the experimentally observed DGDR excitation energy (which only marginally differs from that obtained in the harmonic approximation) and the DGDR cross section for the $`{}_{}{}^{208}Pb+^{208}Pb`$ collision at $`640`$A$``$MeV. They reached a similar conclusion for the $`{}_{}{}^{136}Xe+^{208}Pb`$ collision at $`700`$A$``$MeV, where a much greater discrepancy from the harmonic model appears. The purpose of this paper is to point out that this model does not lead to the enhancement found in Ref. , if proper renormalization of the oscillator frequency is performed in order to guarantee that the theoretical giant dipole resonance (GDR) excitation energy is kept at the experimental value.

The model of Refs.  is based on the following Hamiltonian

$$H=H_0+F(x,y,z;t),$$
(1)
where $`H_0`$ is the anharmonic oscillator describing the intrinsic motion of the projectile,

$$H_0=\frac{1}{2D}(p_x^2+p_y^2+p_z^2)+\frac{C}{2}(x^2+y^2+z^2)+\frac{B}{4}(x^2+y^2+z^2)^2,$$
(2)
where $`D`$ is the mass parameter, $`C`$ is the oscillator strength and $`B`$ is the strength of the anharmonicity. Here, we take the mass parameter to be the reduced mass for the motion of the protons against the neutrons,

$`D={\displaystyle \frac{NZ}{A}}m_0,`$

where $`m_0`$ is the average nucleon mass. The beam is assumed to be parallel to the $`x`$axis and the coupling interaction $`F`$ is derived from the Lienard-Wiechert potential  in the projectile frame

$$\varphi (x,y,z,t)=\frac{Z_Te\gamma }{[\gamma ^2(xvt)^2+(yb)^2+z^2]^{1/2}},$$
(3)
were $`Z_Te`$ is the charge of the target, $`b`$ is the impact parameter, and $`\gamma `$ is the Lorentz factor, $`\gamma =1/\sqrt{1(v/c)^2}`$.

To be specific, we study the $`{}_{}{}^{208}Pb+^{208}Pb`$ collision at $`640`$A$``$MeV. We first solve the Schrödinger equation for the intrinsic motion, described by $`H_0`$. For this purpose it is convenient to recast the intrinsic Hamiltonian into the following equivalent form

$$H_0=\mathrm{}\omega \left[\frac{1}{2}\left(\pi ^2+\rho ^2\right)+\beta \rho ^4\right].$$
(4)
In the above, the commonly used variable transformations

$$\rho _i=\sqrt{\frac{D\omega }{\mathrm{}}}r_i;\pi _i=\frac{p_i}{\sqrt{D\mathrm{}\omega }}$$
(5)
have been made, where $`r_i`$ and $`p_i`$ stand for the components of the position and momentum operators respectively. The oscillator frequency is given by

$$\mathrm{}\omega =\mathrm{}\sqrt{\frac{C}{D}},$$
(6)
and the dimensionless strength $`\beta `$ is related to $`B`$ as

$$B=\left[\frac{4\left(\mathrm{}\omega \right)^3D^2}{\mathrm{}^4}\right]\beta .$$
(7)
In Fig. 1, we show the ratios $`E_{DGDR}^{l=0}/(2E_{GDR})`$ and $`E_{DGDR}^{l=2}/(2E_{GDR})`$ as a function of $`B,`$ in the same range as chosen in Ref. . In this range, the anharmonicity can be treated using first order perturbation theory to great accuracy ($``$ 2%). The GDR and DGDR energies, to first order in $`\beta `$, are given by

$`E_{GDR}(\beta )`$ $`=`$ $`\mathrm{}\omega \left(1+5\beta \right),`$ (8)
$`E_{DGDR}^{l=0}(\beta )`$ $`=`$ $`2\mathrm{}\omega \left(1+7.5\beta \right),`$ (9)
$`E_{DGDR}^{l=2}(\beta )`$ $`=`$ $`2\mathrm{}\omega \left(1+6\beta \right).`$ (10)
Fig. 1 is equivalent to that shown in Ref.  and our results are essentially identical to theirs.

The reduced transition matrix elements can also be easily calculated to first order in the parameter $`\beta `$. We find

$`GDRE1GS`$ $`=`$ $`e\left({\displaystyle \frac{S_1}{\mathrm{}\omega }}\right)^{1/2}\left(12.5\beta \right),`$
$`DGDR,l=0E1GDR`$ $`=`$ $`e\left({\displaystyle \frac{S_1}{\mathrm{}\omega }}\right)^{1/2}\sqrt{{\displaystyle \frac{2}{3}}}(15\beta ),`$
$`DGDR,l=2E1GDR`$ $`=`$ $`e\left({\displaystyle \frac{S_1}{\mathrm{}\omega }}\right)^{1/2}\sqrt{{\displaystyle \frac{10}{3}}}\left(13.5\beta \right),`$
where $`e`$ is the absolute value of the electron charge and $`S_1`$ is given by the energy-weighted sum rule,

$`S_1={\displaystyle \frac{9}{4\pi }}{\displaystyle \frac{\mathrm{}^2}{2m_0}}{\displaystyle \frac{NZ}{A}}.`$

The energy-weighted sum rule for transitions from the ground state and from the GDR are satisfied to first order in the parameter $`\beta `$, using the above energies and reduced matrix elements.

In order to maintain $`E_{GDR}(\beta )`$ at the experimental value, namely $`E_{GDR}(\beta )=E_{GDR}^{\mathrm{exp}}`$ (13.4 MeV, in the present case), the oscillator frequency must be renormalized as $`\beta `$ is changed. The resulting renormalized frequency, from Eq.(8), is

$$\mathrm{}\omega (\beta )=\frac{E_{GDR}^{\mathrm{exp}}}{\left(1+5\beta \right)}.$$
(11)
Note that in the $`B`$-range of Fig. 1, the dimensionless parameter varies in the range $`0.014<\beta <0.014`$ which yields $`1.08E_{GDR}^{\mathrm{exp}}>\mathrm{}\omega (\beta )>0.93`$ $`E_{GDR}^{\mathrm{exp}}`$. Whereas our oscillator frequency is a function of the anharmonicity parameter, in Ref.  it is kept constant at the harmonic value, $`\mathrm{}\omega (\beta =0)=E_{GDR}^{\mathrm{exp}}`$. This difference does not affect the ratio $`E_{DGDR}^l/(2E_{GDR})`$ shown in Fig. 1, since the oscillator frequency cancels out in this case (see eqs. (8) to (10)). When the renormalized frequency is used in both the GDR and DGDR energies and matrix elements, the sum rules for transitions from the ground state and from the GDR are still satisfied. However, use of the renormalized frequency substantially changes the excitation probability of the DGDR, as will be shown below.

The calculation of electromagnetic excitation probabilities and cross sections is performed with the code RELEX, based on the Winther and Alder theory. With this code, we perform a full coupled-channels calculation of the electromagnetic excitation of the GDR and DGDR. Similar results (about 10% larger) would be obtained when perturbation theory is used for the collision dynamics. In Fig. 2, we show the enhancement of the DGDR excitation probability relative to its harmonic value as a function of $`B`$ for the impact parameter $`b=30`$ fm. We find that for $`B100`$ MeV/fm4 (which in this case corresponds to $`\beta 0.7`$ $`\times 10^2`$) the overall enhancement is 6%. For purposes of comparison, we have also performed calculations using a constant frequency ($`\mathrm{}\omega =13.4`$ MeV in this case). We then obtain an enhancement of 35%, as shown by the dashed in line in Fig. 2, in agreement with Ref.  (see their Fig. 1).

In Fig. 3a, we show the enhancement in the impact-parameter integrated DGDR cross section (solid line) vs. $`B`$, for the same system. In the cross section calculations, impact parameters up to 200 fm are taken into account and a lower cut-off at 15 fm is used to eliminate nuclear effects. The full line in Fig. 3a represents the result of the present work, in which an enhancement of only $`4\%`$ is obtained for $`B=100`$ MeV/fm4. The dashed line, obtained using a fixed value of the oscillator frequency, yields an enhancement of the DGDR cross section of 22% for the same value of $`B`$. The $`GDR`$ cross section ratio $`\sigma _{GDR}(B)/\sigma _{GDR}(B=0)`$ obtained with fixed GDR energy, shown as a solid line in Fig. 3b, is close to one over the entire range of $`B`$ values but is slightly less than one for large, negative anharmonicities, (about $`0.5\%`$ at $`B=100`$ MeV/fm4). This small deviation is due to the increase in the population of the DGDR at these values of $`B`$ and the correponding depopulation of the GDR. The GDR cross section ratio obtained with fixed oscillator frequency is shown as a dashed line in Fig. 3b. In this case, we find the GDR cross section to be enhanced by about 10% at $`B=100`$ MeV/fm4. The enhancement of 10% in the GDR cross section of Fig. 3b is clearly responsible for the large enhancement of $`22\%`$ in the DGDR cross section of Fig. 3a at $`B=100`$ MeV/fm4 .

The above conclusions do not change noticeably when the calculations are extended to other systems, such as 136Xe + 208Pb at $`700`$A$``$MeV. The microscopic study of Ref.  established that the anharmonicity parameter scales as $`A^1`$ with the mass number. Thus, if $`B=100`$ MeV/fm4 represents a reasonable value for 208Pb, then for 136Xe a corresponding value would be $`B=150`$ MeV/fm$`^4.`$ In Fig. 4, we display the results of calculations for this system as a function of the anharmonicity parameter $`B`$ in Fig. 4. The solid line in the figure again shows the results of calculations in which the oscillator frequency is varied to maintain the GDR energy constant, while the dashed line represents the results of calculations in which the oscillator frequency is maintained fixed. Similar to the previous case, we find the enhancement of the DGDR cross section to be greatly reduced when the GDR resonance energy is maintained at a fixed value. As can be seen in Fig. 4, at $`B=150`$ MeV/fm4, the DGDR cross section is enhanced by 62% when the oscillator frequency is maintained constant, but is enhanced by less than 10% when the GDR energy is maintained at its physical value.

Before ending we comment briefly on the connection between the Bortignon-Dasso model used in this paper and microscopic models  that aim to assess the importance of the anharmonic effects both on the spectrum and on the transition operator. Ref.  finds, within the Lipkin model, small effects on the spectrum (which scale roughly as 1/A). Hamamoto finds, within nuclear field theory, that the nonlinear effects in the 1-phonon to 2-phonon transition operator are also quite small and scale as 1/A. As mentioned above, Ref. , through detailed microscopic calculations, finds that the anharmonic effects are indeed small and scale as 1/A. The values of the parameter $`B`$ in both the Bortignon-Dasso and present calculations are taken to be small enough to be in line with the microscopic findings but also with the experimentally observed DGDR excitation energies (although the enhancement of the DGDR cross section could be increased thorough an artificially large $`B`$, there is no choice for this parameter that whould simultaneously explain the observed cross section enhancement and the only very small deviations of the DGDR excitation energy from the harmonic limit).

Another interesting point to mention is that the GDR has a width, which is considered neither by Bortignon and Dasso nor in the present calculation. The effect of the width of the GDR on the excitation of the DGDR has been recently studied within a harmonic picture. The overall effect of the width, at the energies considered here is, to produce a slight increase in the DGDR cross section, although not enough to explain all the available data. It would certainly be of interest to extend the present calculation within the anharmonic model by coupling the oscillator to other degrees of freedom (which would generate the damping width).

In conclusion, we have investigated the effect of anharmonicities in the excitation of the DGDR in relativistic heavy ion collisions, with the same macroscopic model used by Bortignon and Dasso. We point out that variations of the anharmonicity strength must be accompanied by a renormalization of the oscillator frequency, in order to maintain the GDR energy at a value consistent with the experimental one. We have found that this condition strongly reduces the enhancement in the DGDR excitation probabilities and corresponding cross sections, so that they remain much below the experimental results.

This work was supported in part by DAAD/CAPES cooperative agreement no. 415-bra-probral/bu, CNPq and the MCT/FINEP/CNPq(PRONEX) under contract no. 41.96.0886.00. D.T.P and L.F.C. acknowledge partial support from the Fundação Universitária José Bonifácio, and M.S.H. and B.V.C. acknowledge support from the FAPESP.

Figure Captions

* Figure 1: The ratio $`E_{DGDR}^l/(2E_{GDR})`$ vs the anharmonicity parameter $`B,`$ for $`{}_{}{}^{208}Pb`$. The solid line is for $`l=2`$ and the dashed line for $`l=0.`$ The reduced mass for the oscillation of protons against neutrons is used for the mass parameter $`D`$.
* Figure 2: The enhancement in the excitation of the DGDR in the collision of $`{}_{}{}^{208}Pb+^{208}Pb`$ at $`640`$A$``$MeV for the impact parameter $`b=30`$ fm. The solid line represents the results of the present calculation while the dashed line corresponds to a constant oscillator frequency.
* Figure 3: Enhancement factor of the (a) DGDR and (b) GDR cross sections in the collision of $`{}_{}{}^{208}Pb+^{208}Pb`$ at $`640`$A$``$MeV. The dashed lines correspond to the results obtained with fixed oscillator frequency, while the full lines correspond to a fixed $`E_{GDR}`$.
* Figure 4: Enhancement factor of the DGDR cross section in the collision of $`{}_{}{}^{136}Xe+^{208}Pb`$ at $`700`$A$``$MeV. The dashed line corresponds to the results obtained with fixed oscillator frequency, while the full line corresponds to a fixed $`E_{GDR}`$.


format: text

text:
# Abstract

## Abstract

The NASA program in X-ray astronomy has two long term goals: 1) to achieve sufficient angular resolution to image the event horizon of a black hole (0.1 micro arc sec) and 2) to achieve sufficient collecting area (50-150 sq m) and angular resolution (0.1–1.0 arc sec) to observe in detail the first black holes and galaxies at high redshift. These ambitous goals can be used to map out a series of missions and a technology program. The next major mission will be Constellation-X which will be dedicated to high resolution X-ray spectroscopy for launch in $``$ 2010. This mission is a critical step in the roadmap to achieve these goals. Following Constellation-X NASA is considering two very ambitious vision missions: MAXIM and Generation-X that will achieve the ultimate capabilities. The modest missions Astro-E2 and Swift address more focussed science goals on a rapid development cycle and provide important pathfinders to the larger missions.

Key-words X-ray Astronomy, Missions

## 1 Introduction

The current generation of large X-ray observatories: the Chandra X-ray Observatory from NASA and the XMM-Newton Observatory from ESA are delivering an impressive array of results. As we look to the future the requirements for future X-ray observatories are driven by some basic questions: What is the nature of space and time? When did the first black holes appear? What is the Universe made of? X-ray observations are central to addressing these questions. Using these questions we can plan a roadmap for the future of X-ray astronomy.

### 1.1 What is the nature of Space and Time?

Black holes represent the most extreme environments known. They provide the only known regime where the effects of General Relativity (GR) can be tested in the strong gravity regime. Material falling into the black hole is heated to high temperatures and the bulk of the energy is emitted in X-rays. Results from the Japanese-US ASCA mission have demonstrated that X-ray spectral features are coming from close to the event horizon of the black hole. These results suggest that X-ray observations hold the promise of using black holes to test GR. The most direct way to do this is to resolve the event horizon of a black hole and directly observe the predicted distortions of space time. The best candidates for this are the super massive black holes in nearby active galactic nuclei, such as M87. To resolve the event horizon will require a resolution of 0.1 to 1.0 micro arc sec. This can only be achieved with X-ray interferometry and is a similar challenge in scope to the plans for the Terrestrial Planet Finder (TPF) and Darwin missions. The MAXIM mission will address this science goal. In the nearer term, the study of spectroscopic line signatures from close to the event horizon requires an 100 times increased sensitivity compared to Chandra and XMM-Newton for high resolution spectroscopy with a broad 0.25 to 40 keV energy range. Such a capability is more accessible with technology that is an extension of the current generation of missions and will be addressed by the Constellation-X mission.

### 1.2 When did the first black holes appear?

Current theories predict that the first black holes formed at redshifts of 10-20 with masses of 103 to 106 solar masses and that they play a crucial role in galaxy formation. X-ray observations are essential to observe these black holes turn on as proto-quasars, and evolve with cosmic time. The Chandra and XMM-Newton deep surveys are revealing that active super-massive ($`>`$ 1 million solar mass) black holes are present in up to 10% of galaxies at redshifts of 1-3 (compared to 1% in the local universe). These results are demonstrating the power of X-ray surveys to find accreting black holes in otherwise normal galaxies. X-ray spectroscopic observations can be used to determine the evolution of the black hole spin with cosmic time, and infer black hole mass from reverberation analysis and X-ray luminosity scaling arguments. To observe the first black holes requires a collecting area 100 to 1,000 times Chandra and XMM-Newton, with $``$ 0.1–1 arc sec angular resolution. This capability would also allow the X-ray detection of normal galaxies essentially throughout the Universe, which will be important to determine whether black holes formed, before, during, or after, the era of galaxy formation. The Generation-X mission is tuned to directly achieve this goal. The nearer term Constellation-X mission will provide critical measurements of the evolution of active galactic nuclei with redshift.

### 1.3 What is the Universe made of?

A pressing challenge in astronomy today is determining the composition of the Universe. Using clusters of galaxies as cosmological probes, X-ray observations have constrained the matter content (dark and baryonic) to be about 35% of the closure density. This result was a precursor to the recent revelation from observations of the cosmic microwave background and type 1a supernovae that the remaining 65% of the universe is made up of a mysterious dark energy. X-ray observations are central to this problem because the soft and medium X-ray band contains spectral diagnostics from all the abundant elements formed after the big bang. X-ray observations of supernovae, the hot interstellar and intergalactic medium, galaxies and clusters of galaxies provide a direct probe of when the elements from which we are made were created and the formation of large scale structure. Current cosmological models predict that most of the baryons in the Universe are in a hot intergalactic medium, and that sufficiently sensitive X-ray observations will detect this gas in the absorption spectra of high redshift quasars, or other bright distant objects. This requires an increase in sensitivity by a factor of 100 relative to Chandra and XMM-Newton for high resolution (R$`>`$300) spectroscopy. The same increase in capability will also allow precise measurements of the first clusters of galaxies at high redshift which, in concert with future observations of the microwave background and of type Ia supernovae, will tightly constrain the overall energy and mass content of the universe. In addition, an extension of imaging observations to the hard X-ray band will be very important to investigate the contribution of non-thermal processes. Constellation-X will directly address this question via high resolution spectroscopy.

## 2 Constellation-X

The Constellation-X mission is a large collecting area X-ray observatory emphasizing high spectral resolution (R $``$ 300 to 3000) and a broad energy band (0.25 to 60 keV). It is the next major mission being planned by NASA for X-ray astronomy. Constellation-X is a high priority in the US National Academy of Sciences McKee-Taylor 2000-2010 Survey for a new start and flight in the first decade of the 21st century.

The 0.25–10 keV X-ray band contains the K-shell lines for all of the abundant metals (carbon through zinc), as well as many of the L-shell lines. The detailed X-ray line spectra are rich in plasma diagnostics which also provide unambiguous constraints on physical conditions in the sources. A spectral resolving power of at least 300 is required to separate the He-like, density sensitive, triplet from ions such as O, Si, and S. In the region near the iron K complex a resolving power exceeding 2000 is necessary to distinguish the lithium-like satellite lines from the overlapping helium-like transitions. A resolution of 300–3000 provides absolute velocity measurements of 10–100 km/s, commonly found in many astrophysical settings.

By increasing the telescope aperture to 3 m2 with 5–15 arc sec resolution and utilizing efficient spectrometers the mission will achieve the required 25 to 100 increased sensitivity over current high resolution X-ray spectroscopy missions. The use of focusing optics across the 10 to 60 keV band will provide a similar factor of 100 increased sensitivity in this band, essential for constraining the nature of the underlying continuum. The effective area curves of Constellation-X are shown in Figure 1, and compared with the equivalent curves from the high resolution spectrometers on Chandra and XMM-Newton. The peak telescope area of 30,000 cm2 when combined with the efficiency of the spectrometers and detectors gives the required 15,000 cm2 at 1 keV and 6,000 cm2 at 6 keV. The effective area of 1,500 cm2 at 40 keV is set by using a typical AGN spectrum to match the sensitivity in the 0.25–10 keV band. A $`10^5`$ s exposure of a typical AGN with the Constellation-X will obtain $``$ 1000 counts in the calorimeter at a flux of $``$ $`2\times 10^{15}`$ erg cm-2 s-1 (0.2–2 keV band).

A spectroscopy X-ray telescope (SXT) covers the 0.25–10 keV band. X-ray spectroscopy requires that a maximum number of photons be collected. To maximize the collecting area per unit mass an angular resolution requirement of $`15^{^{\prime \prime }}`$ half power diameter (HPD) is set for the SXT as the baseline requirement, with a goal of $`5^{^{\prime \prime }}`$. The $`15^{^{\prime \prime }}`$ HPD requirement is consistent with the confusion limit for the limiting flux of the sources to be studied. The SXT uses two complementary spectrometer systems to achieve the desired energy resolution of 300–3000: an array of high efficiency quantum micro-calorimeters with energy resolution of 2eV, and a set of reflection gratings with a resolution of 0.05$`\AA `$ in the first order. The gratings deflect part of the telescope beam away from the calorimeter array to a CCD array in a design similar to XMM-Newton, except that the Constellation-X direct beam falls on a high spectral resolution quantum calorimeter instead of on a CCD. The baseline field of view of the calorimeter is $`2.5^{^{}}`$ square, with a 30 x 30 array. This gives $`5^{^{\prime \prime }}`$ pixels, sized to adequately sample the telescope point spread function.

The Constellation-X hard X-ray telescope (HXT) baseline design uses multi-layers to provide a focusing optics system that operates in the band above 10 keV band. The baseline upper energy is 40 keV, with a goal of reaching 60 keV. The improvement in the signal to noise results in a factor of 100 or more increased sensitivity over non-focussing methods used in this band. The requirement for the HXT angular resolution is $`1^{^{}}`$ HPD (with a goal of $`30^{^{\prime \prime }}`$) which is quite sufficient to resolve sources in this energy band. The FOV is $`8^{^{}}`$, or larger. There are no strong atomic lines expected above 10 keV, so there is a relatively modest HXT spectral resolution requirement of R $`>`$10 across the band.

The required large collecting area is achieved utilizing several mirror modules, each with its own spectrometer/detector system. The Constellation-X design takes advantage of the fact that several smaller spacecraft and more modest launch vehicles each carrying one “science unit” can cost less than one very large spacecraft and launcher. The program is then very robust in that risks are distributed over several launches and spacecraft with no single failure leading to loss of mission. The current baseline mission is four satellites that are carried in pairs on either two Atlas V or Delta IV launchers. The interval between the two launches will be $``$ 1 year, with the first launch in the 2010 timeframe. A L2 orbit will facilitate high observing efficiency, provide an environment optimal for cyrogenic cooling, and simplify the spacecraft design. The telescope collecting area is sized assuming $``$ 95% observing efficiency, which maximizes the observing program for a minimum 5 year mission lifetime (with all satellites operating).

## 3 Vision Missions

To set the long term goals, NASA has put on the road map two ”vision missions” for the $`>`$2015 timeframe

### 3.1 Generation-X

Generation-X is designed to observe the high redshift (5-20) Universe to observe X-ray emission from the first black holes and starburst galaxies at redshift 5-20, and to observe the X-ray evolution of galaxies over all redshifts. These science objectives require that the mission has 50-150 m2 collecting area at 1 keV combined with Chandra-like angular resolution of 0.1 to 1.0 arc sec resolution. This resolution is necessary to provide the necessary detection sensitivity and avoid source confusion. The large collecting area would be achieved by developing precision light weight X-ray optics, with a mass per unit collecting area similar to the Astro-E optics, but with a factor of 100 improved angular resolution. As currently envisioned the mission would consist of 3-6 identical satellites, each carrying a fraction of the collecting area. Deployable optics would be used to maximize the collecting area per launch vehicle. The focal length of 25 to 50 m would be achieved by using an extendible mast. The mission would be placed at L2 to provide a stable thermal environment and to maximize the viewing efficiency. The enabling technology for this mission is lightweight, precision, grazing incidence optics and the mission requirements feed into NASA’s long term technology program.

### 3.2 MAXIM

The Micro-Arcsecond X-ray Imaging Mission (MAXIM) will have the angular resolution of 0.1 to 1 micro arc second required to resolve the event horizon of accreting black holes at the center of nearby galaxies (e.g. M87). MAXIM will employ for the first time X-ray interferometry and in doing so will achieve a spectacular 10 million times increase in angular resolution compared to Chandra. The basic technique has recently been demonstrated in the laboratory. To obtain 0.1 micro arcsec requires a 100 m to 1000 m baseline. The mission is being studied by NASA’s Institute for Advanced Concepts (NIAC). It involves up to 34 spacecraft flying in formation, with precision metrology and tolerances at the nano-meter level. On orbit metrology and stability requirements are similar to TPF/Darwin. A two spacecraft Pathfinder mission is being studied that will have a 1 m baseline giving 0.1 milli-arc second imaging. This will provide the required demonstration of X-ray interferometer in space and at the same time give a 1,000 times improvement compared to Chandra, which in itself will allow exciting new science e.g. imaging coronal structures on nearby stars and the base of jets in active galactic nuclei. The full MAXIM is still very early in a concept study phase and the technology requirements are very challenging.

## 4 Modest Missions

It is also important that there is an adequate number of more modest missions that focus on more specialized scientific goals. NASA pursues these small to medium sized opportunities via its Explorer program, which allows rapid response to topical science challenges.

### 4.1 Swift

Swift is a multiwavelength observatory for gamma-ray burst (GRB) astronomy, that is in development for launch in 2003. It is the first-of-its-kind multi-wavelength, autonomous rapid slewing satellite for studying explosive astronomical events. It will be far more powerful than any previous GRB mission, observing 150-300 bursts per year and performing detailed X-ray and UV/optical afterglow observations spanning timescales from 1 minute to several days after the burst. The mission is being developed by an international collaboration led by NASA, with Italy and the UK providing significant parts of the science payload and operations. The mission will carry three instruments: a new-generation wide-field gamma-ray detector and narrow-field X-ray and optical telescopes. Innovations from the Swift program include: 1) a large-area gamma-ray detector using CdZnTe detectors; 2) an autonomous rapid slewing spacecraft; 3) a multi-wavelength payload combining optical, X-ray, and gamma-ray instruments; 4) immediate data flow to the community.

### 4.2 Astro-E2

The Astro-E2 mission is a Japanese mission, with substantial US participation. It is a re-flight of Astro-E, which suffered a launcher failure on February 10, 2000. Astro-E2 is approved to be launched with an ISAS M-V rocket in early 2005. Astro-E2 will basically be identical to Astro-E. The spacecraft will carry five X-Ray Telescopes (XRTs). One telescope module will feed the X-ray Spectrometer (XRS) and the other four telescopes will have X-ray Imaging Spectrometer (XIS) units at their focal planes. At the same time, the Hard X-ray Detector (HXD) will be able to measure high-energy X-rays. The development of the XRT, XRS and XIS are being performed in collaboration with NASA. Even though Astro-E2 will be realized about 5 years after Astro-E was expected these capabilities remain unsurpassed. The lifetime of the XRS will be limited to 2–3 yr by expendable cryogenics. The broad band-pass and spectroscopic response will provide a unique capability for addressing a wide variety of fundamental problems in astrophysics including the origin of the elements and structures in the Universe as well as the evolution and dynamics of these structures.


text:
# Stochastic volatility and leverage effect

## Abstract

We prove that a wide class of correlated stochastic volatility models exactly measure an empirical fact in which past returns are anticorrelated with future volatilities: the so-called “leverage effect”. This quantitative measure allows us to fully estimate all parameters involved and it will entail a deeper study on correlated stochastic volatility models with practical applications on option pricing and risk management.

The multiplicative diffusion process known as the geometric Brownian motion (GBM) has been widely accepted as one of the most universal models for speculative markets. The model, started out by Bachelier in 1900 as an ordinary random walk and redefined in its final version by Osborne in 1959 cootner , presupposes a constant “volatility” $`\sigma `$ which is equivalent to a constant diffusion coefficient $`D=\sigma ^2`$. However, and especially after the 1987 crash, there seems to be ample empirical evidence, given by the so-called “stylized facts”, that the assumption of constant volatility does not properly account for important features of markets. Some relevant examples of stylized facts are the “leverage” and “smile” effects and the skewness and “fat tails” in probability distributions cont .

One of the main ideas that has come out to explain those features is that the market dynamics has intrinsically changed in the sense that volatility is no longer constant. It is not a function of time either (as might be inferred by the evidence of nonstationarity in financial time series lo ) but a random variable. In its most general form one therefore assumes that the volatility $`\sigma `$ is a function of a random process $`Y(t)`$, i.e., $`\sigma (t)=\varphi (Y(t))`$. We could make an analogy from physics saying that speculative prices $`S(t)`$ evolve in a “random medium” determined by a random diffusion coefficient. Most of the stochastic volatility (SV) models presented up to date suppose that $`Y(t)`$ is itself a diffusion process that may or may not be correlated with price and different models differ from each other basically in the form of the function $`\varphi `$ ghysels ; fouquebook .

As mentioned above, the hypothesis of stochastic volatility (SV) has been suggested to explain significant peculiarities observed in real markets. Thus, the smile effect, related to the implicit volatility in option prices, has been thoroughly studied both qualitatively and quantitatively smile ; heston . Other features such as leverage and long tails are less studied and some works on SV only address to them from a qualitative point of view fouque  while others measure them by giving numerical coefficients, based on ARCH-GARCH models, for kurtosis and skewness englepat . However, the time evolution and structure of the leverage and tails have never been investigated. Our main objective here is to prove that the leverage effect can be quantitatively explained by a wide class of correlated SV models. This will allow us to overturn the main objection against SV models: the impossibility of fitting all parameters appearing in these models fouquebook ; fouque  which, in turn, opens the door to simulations of real markets with far reaching practical consequences on option pricing and risk management.

We recall that the leverage effect has its origin in the observation that volatility seems to be negatively correlated with stock returns which, in continuous time finance and in terms of speculative prices $`S(t)`$, are defined by $`R(t)=\mathrm{ln}[S(t)/S(0)]`$. The first explanation to this empirical fact was given by Black black  and Christie christie  in the sense that negative returns increase financial leverage which extend the risk of the company and therefore its volatility. Hence the name of “leverage effect”. Nevertheless, the cause of this effect is still unclear since another explanation is just the contrary, that is, an increase of volatility makes the stock riskier which produces a fall of demand and the price drops ghysels .

In a very recent paper, Bouchaud et al. bouchaud  have performed a complete empirical study of the leverage effect for both individual stocks and indices using daily data. The volatility-return (negative) correlation is clearly shown to have a definite direction in time – a very confusing fact in the literature – since correlations are shown to be between future volatilities and past returns. Bouchaud et al. conclusively prove from data that the negative correlation decays exponentially in time, faster for indices than for individual stocks bouchaud .

In this letter, we present a theoretical study on these correlations and show that a wide class of SV models completely explain the leverage effect, both qualitatively and quantitatively, in complete agreement with experimental observations. The starting point is the GBM model:

$$dR=\mu dt+\sigma (t)dW_1(t),$$
(1)
where $`\mu `$ is the drift and $`\sigma (t)=\sigma (Y(t))`$ is a random volatility and $`Y(t)`$ is a diffusion process:

$$dY=f(Y)dt+g(Y)dW_2(t).$$
(2)
In these equations $`W_i(t)`$ $`(i=1,2)`$ are Wiener processes, i.e., $`dW_i(t)=\xi _i(t)dt`$ where $`\xi _i(t)`$ are zero-mean Gaussian white noises with cross-correlation given by

$$\xi _1(t)\xi _2(t^{})=\rho \delta (tt^{}),$$
(3)
$`(1\rho 1)`$. As is common in finance, Eqs. (1)-(2) are interpreted in the sense of Ito and for the rest of the paper we will follow Ito convention perello1 .

Bouchaud et al. bouchaud , quantify the leverage effect by means of a leverage correlation function defined by

$$(\tau )\frac{1}{Z}[dX(t+\tau )]^2dX(t)$$
(4)
where

$$X(t)R(t)\mu t$$
(5)
is the zero-mean return and $`Z=[dX(t)^2]^2`$ is a convenient normalization coefficient. Bouchaud et al. have analyzed a large amount of daily relative changes for either market indices and stock share prices and find that bouchaud

$$(\tau )=\{\begin{array}{cc}Ae^{b\tau },\hfill & \text{if }\tau >0\text{;}\hfill \\ 0,\hfill & \text{if }\tau <0\text{;}\hfill \end{array}$$
(6)
$`(A,b>0)`$. Hence, there is a negative correlation with an exponential time decay between future volatility and past returns changes but no correlation is found between past volatility and future price changes. In this way, they provide a sort of causality to the leverage effect which, to our knowledge, has never been previously mentioned in the literature ghysels ; fouquebook .

Let us sketch how correlated SV models are able to exactly reproduce this result. We first combine Eqs. (1) and (5) to produce $`(\tau )=\sigma (t)dW_1(t)\sigma (t+\tau )^2dW_1(t+\tau )^2/Z`$. Note that when $`\tau <0`$ Ito’s rules tell us that $`dW_1(t)`$ is uncorrelated with the rest of terms then, recalling that $`dW_1(t)=0`$, we have $`(\tau )=0`$ if $`\tau <0`$. On the other hand, when $`\tau >0`$, $`dW_1(t+\tau )`$ is uncorrelated with the rest and, since $`dW_1(t+\tau )^2=dt`$, we conclude that

$$(\tau )=\theta (\tau )\sigma (t)dW_1(t)\sigma (t+\tau )^2dt/Z,$$
(7)
where $`\theta (\tau )`$ is the Heaviside step function and $`Z=[\sigma ^2(t)]^2dt^2`$. Note that we have proved the existence of correlations between future volatilities and past returns but not vice-versa. Note also that this have been proved independently of the underlying volatility process $`Y(t)`$ (which needs not to be a diffusion process) and of the specific form of $`\sigma `$ in terms of $`Y`$.

Suppose now that $`Y(t)`$ is a diffusion process given by Eq. (2). As is well known any pair of correlated Wiener process, such as $`W_1(t)`$ and $`W_2(t)`$, satisfy the identity $`dW_1(t)=\rho dW_2(t)+\sqrt{1\rho ^2}dW(t)`$, where $`W(t)`$ is a Wiener process independent of $`W_2(t)`$ (therefore, $`W(t)`$ is independent of $`\sigma `$). Substituting this identity into Eq. (7) we get $`(\tau )`$ in terms of the average $`\sigma (t)\sigma (t+\tau )^2\xi _2(t)`$. This average can be calculated by means of Novikov’s theorem novikov  with the result pm

$$(\tau )=\frac{2\rho \theta (\tau )}{[\sigma ^2(t)]^2}\sigma (t)\sigma (t+\tau )\sigma ^{}(t+\tau )\frac{\delta Y(t+\tau )}{\delta \xi _2(t)},$$
(8)
where $`\sigma ^{}=\sigma (Y)/Y`$ and $`\delta Y(t+\tau )/\delta \xi _2(t)`$ is the functional derivative of $`Y(t+\tau ,[\xi _2])`$ with respect to $`\xi _2(t)`$ novikov .

There is a wide consensus on volatility being “mean reverting”. This means that there exists a normal level of volatility footnote  to which volatility will eventually return englepat . For a general SV model such as (1)-(2), the presence of mean-reversion implies restrictions on the form of the drift coefficient $`f(Y)`$. In order to include this experimental fact in the model, the simplest choice is to assume that $`f(Y)`$ is a linear function. That is,

$$\dot{Y}=\alpha (Ym)+g(Y)\xi _2(t),$$
(9)
where $`\alpha <0`$. The formal solution to this equation in the stationary state reads

$$Y(t)=m+_{\mathrm{}}^te^{\alpha (tt^{})}g(Y(t^{}))\xi _2(t^{})𝑑t^{},$$

from which we get pm

$$\frac{\delta Y(t+\tau )}{\delta \xi _2(t)}=\theta (\tau )e^{\alpha \tau }g(Y(t))\mathrm{exp}\left[_t^{t+\tau }g^{}(Y(s))\xi _2(s)𝑑s\right].$$

Substituting this expression into Eq. (8) yields

$$(\tau )=\rho \theta (\tau )B(\tau )e^{\alpha \tau },$$
(10)
where

$$B(\tau )=\frac{2\sigma (t)\sigma (t+\tau )\sigma ^{}(t+\tau )G(t,t+\tau )}{[\sigma ^2(t)]^2},$$
(11)
and

$$G(t,t+\tau )=g(Y(t))\mathrm{exp}\left[_t^{t+\tau }g^{}(Y(s))\xi _2(s)𝑑s\right].$$
(12)
In consequence, any SV model of the form given by Eqs. (1) and (9) and whose function $`B(\tau )`$ does not increase faster than $`e^{\alpha \tau }`$ as $`\tau \mathrm{}`$, satisfies an exponentially decaying leverage as expressed by Eq. (6). Moreover, if $`\sigma (Y)`$ is an increasing function of $`Y`$ with definite sign and $`g(Y)`$ is positive definite (or $`\sigma `$ is decreasing and $`g`$ is negative) we see from Eq. (10) that the correlation coefficient $`\rho `$ must be negative and driving noises $`W_1(t)`$ and $`W_2(t)`$ are anticorrelated. Eqs. (10)-(12) constitute the main result of the paper.

The exact form of $`(\tau )`$ will depend on the expression of $`B(\tau )`$ which in turn will depend on the SV model chosen. Within diffusion theory, as is the case of Eq. (9), there are basically three different SV models fouquebook : 1) The Ornstein-Uhlenbeck (OU) model where $`\sigma =Y`$ and $`g(Y)=k`$ (a positive constant) stein , 2) the exponential Ornstein-Uhlenbeck (expOU) model where $`\sigma =e^Y`$ and $`g(Y)=k`$ fouque  and 3) the Cox-Ingersoll-Ross (CIR) model where $`\sigma =\sqrt{Y}`$ and $`g(Y)=k\sqrt{Y}`$ heston . For all these models the leverage function has the form given by Eq. (10). In the OU model and in the CIR model (the latter with zero mean-reversion, i.e., $`m=k^2/4\alpha `$) the leverage function $`(\tau )`$ is respectively given by pm

$$_1(\tau )=2k\rho \theta (\tau )\left[\frac{m^2+(k^2/2\alpha )e^{\alpha \tau }}{(m^2+k^2/2\alpha )^2}\right]e^{\alpha \tau },$$
(13)
$$_3(\tau )=4\frac{\rho \alpha }{k}\theta (\tau )e^{\alpha \tau },$$
(14)
while for the expOU model we have

$$_2(\tau )=2\rho k\theta (\tau )\mathrm{exp}\left[m+k^2(e^{\alpha \tau }3/4)/\alpha \right]e^{\alpha \tau }.$$
(15)
In Fig. 1 we show the leverage effect for the Dow-Jones daily index (1900-2000) and plot the leverage function for the OU model.

Any market model, besides being able to reproduce the market dynamics, must provide a systematic way of evaluating its parameters. Almost all current SV models have four parameters to estimate: $`\rho `$, $`k`$, $`m`$, and $`\alpha `$. To our knowledge, all works on SV models presented up till now are able to evaluate only two of them. Thus, for instance, Fouque et al. fouque  estimate $`k`$ and $`m`$ from the empirical second and fourth moment of daily data but they can’t give a clear estimation of $`\alpha `$ and $`\rho `$. This constitutes the main criticism to SV models, say, their inability to estimate all parameters involved. This situation changes completely when one measures leverage. Indeed, $`k`$ and $`m`$ are obtained, as usual, from the empirical second and fourth moment. Next, by adjusting $`e^{\alpha \tau }`$ to leverage empirical data we estimate $`\alpha `$. Moreover, comparing the theoretical and empirical leverage at $`\tau =0^+`$, $`(0^+)`$, we finally obtain $`\rho `$. Following this procedure, for the Dow-Jones daily index and using the OU model, we have estimated $`k=1.4\times 10^3\text{ days}^1`$, $`m=18.9\%\text{ year}^{1/2}`$, $`\alpha =0.05\text{ day}^1`$, and $`\rho =0.58`$.

As an illustration we have simulated, using Eqs. (1) and (9), the OU resulting process with the parameters estimated above. We follow the random dynamics of the daily changes of the zero-mean return $`X(t)`$ and compare it with the empirical Dow-Jones time series during one trading year. We have also simulated the geometric Brownian motion assuming a constant volatility $`\sigma `$ whose value is directly estimated from Dow-Jones one-century data. We present these results in Fig. 2 and observe that GBM cannot describe either the largest or the smallest fluctuations of daily returns. We nonetheless see in the figure that the SV model chosen describes periods of high volatility together with periods of very low volatility, resulting in a more similar trajectory to the Dow-Jones index than that of the GBM. This is quite remarkable, because we have simulated last year trajectory using all past data (one century) of the Dow-Jones index thus showing the stability of parameters.

These results may encourage to deepen in the study of statistical properties of SV models showing leverage. Several models have been presented in the literature without being able to discern which is the more realistic one. Now, thanks to leverage correlation, it is possible to estimate all parameters involved in SV models. This will allow us to confront in detail different models with the empirical statistical properties of markets. Finally, a better knowledge of SV models has non trivial consequences on option pricing (since classical Black-Scholes method is still suitable in SV models smile ; heston ) and, more generally, on risk management.

###### Acknowledgements.

This work has been supported in part by Dirección General de Investigación under contract No. BFM2000-0795, by Generalitat de Catalunya under contract No. 2000 SGR-00023. We thank Marian Bogunya for useful suggestions and careful reading of the manuscript.


text:
# 1 Introduction

## 1 Introduction

The QCD vacuum is characterized by two striking phenomena, the breaking of chiral symmetry and the confinement of color charge. Chiral symmetry breaking may be understood in terms of localized topological excitations of the gluon field and their associated quark zero-modes that produce a non-vanishing value of the chiral condensate. Classical instanton  solutions of the Yang Mills equations with topological charge $`Q=1`$ and their quantum fluctuations provide a physical foundation for these topological excitations and thus a natural understanding of chiral symmetry breaking.

In contrast, the mechanism for confinement is not presently well understood, and various pictures have been investigated to try to explain it in terms of relevant structures in the QCD vacuum. Various point-like solutions to the Yang Mills equations, which fall off at large distances in all space-time dimensions, have been considered. Although Q=1 instanton solutions provide an understanding of chiral symmetry breaking, in the dilute gas and instanton liquid approximations they do not lead to confinement. Merons, topological charge $`\frac{1}{2}`$ solutions found by De Alfaro, Fubini and Furlan , are more strongly disordering objects than instantons and were proposed as a mechanism for confinement by Callan, Dashen and Gross . Fractons, also solutions of the Yang Mills equations of motion with fractional topological charge, appear on the four-dimensional torus, $`T^4`$, when twisted boundary conditions are imposed . The possible relevance of these objects to confinement was pointed out in , and a scenario for confinement based on the fractional charge solution found in reference , was proposed by González-Arroyo and Martínez .

One and two-dimensional structures in the QCD vacuum have also been considered as mechanisms for confinement. In the dual superconductor picture , the condensation of monopoles in the QCD vacuum leads to confinement. Monopoles are one-dimensional curves in space-time that appear in QCD as defects in the abelian gauges proposed by ’t Hooft . The gauge is fixed up to the Cartan subgroup of the gauge group and monopoles appear at points in space where this gauge fixing is ill-defined, leaving a gauge freedom larger than the abelian subgroup. In the vortex theory, confinement is due to the condensation of vortices. Vortices are two-dimensional surfaces carrying flux in the center of the SU(N) group, which means that a Wilson loop intersecting the surface of the vortex takes the value of one of the elements of the center of the group. Classical vortex solutions to the SU(N) Yang Mills equations have been found numerically .

The mechanism for chiral symmetry breaking and the alternative descriptions of confinement are not mutually exclusive - rather they are highly interrelated. The fact that the intersection of two vortices has topological charge $`\frac{1}{2}`$  provides a provocative connection between chiral symmetry breaking and confinement and suggests that the confinement properties of charge $`\frac{1}{2}`$ merons may also be understood in terms of the intersections of vortices. In addition, as elaborated below, monopole lines lie on vortex surfaces, so that both structures coexist and may be studied simultaneously. In this picture, a meron pair corresponds to the intersection of two closed vortex sheets containing closed monopole loops and provides the simplest system in which one could explore this structure quantitatively. As the separation between the merons decreases to zero and they merge into an instanton, one would expect a vortex sheet and a monopole loop on it to shrink to a point at the center of the instanton . A similar picture of the separation of an instanton into two fractionally charged objects connected by hedgehog world lines is given in reference .

In this article we investigate numerically the monopole and vortex content of a meron pair in SU(2) Yang Mills theory by calculating the points at which Laplacian Center Gauge fixing is ill-defined and by calculating the behavior of Wilson loops. The monopole and vortex content of an isolated meron has already been studied analytically by Reinhardt and Tok  using Laplacian Center Gauge fixing and Wilson loops, and provides an essential foundation for the present work. Since their work, as well as that of others, has shown Laplacian Center Gauge fixing to be an imperfect tool, in this study we also explore the limitations of this tool as well as the physics of the QCD vacuum.

The outline of this letter is the following. In section 2 we describe the meron pairs that we study and in section 3 we use Wilson loops to explore their vortex content. Section 4 presents the monopole and vortex content of these configurations determined from Laplacian Center Gauge defects and section 5 summarizes our conclusions.

## 2 The meron pair

Merons  are solutions to the classical Yang-Mills equations of motion in four Euclidean dimensions, which can be written as

$$A_\mu ^a(x)=\eta _{a\mu \nu }\frac{x^\nu }{x^2},$$
(1)
where $`\eta _{a\mu \nu }`$ is the ’t Hooft symbol. Using the conformal symmetry of the classical Yang-Mills action, it can be shown that in addition to a meron at the origin, there is a second meron at infinity, and these two merons may be mapped to arbitrary positions. The gauge field for the two merons  is

$$A_\mu ^a(x)=\eta _{a\mu \nu }\left[\frac{x^\nu }{x^2}+\frac{(xd)^\nu }{(xd)^2}\right].$$
(2)
This gauge field for the meron pair has infinite action density at points $`x_\mu =\{0,d_\mu \}`$. To avoid the problem of these singularities, we use the following expression

$$A_\mu ^a(x)=\eta _{a\mu \nu }x^\nu \{\begin{array}{ccc}\frac{2}{x^2+r^2},& \sqrt{x^2}<r,\hfill & \\ \frac{1}{x^2},& r<\sqrt{x^2}<R,\hfill & \\ \frac{2}{x^2+R^2},& R<\sqrt{x^2}.\hfill & \end{array}$$
(3)
Here, the singular meron fields for $`\sqrt{x^2}<r`$ and $`\sqrt{x^2}>R`$ are replaced by instanton caps, each containing topological charge $`\frac{1}{2}`$ to agree with the topological charge carried by each meron. We study the monopole and vortex content of this configuration by putting the gauge field on a lattice of size $`N_t\times N_s^3`$. For details of the procedure for putting the meron pair on the lattice and relaxing it to a solution of the field equations, see reference .

In this article, we analyze four meron pair configurations obtained on $`N_t\times N_s^3`$ lattices with $`N_s=16,24`$ and $`N_t=2N_s`$. We study configurations with different cap sizes, $`c`$, distances between merons, $`d`$, and sizes of the lattice, $`N_s`$. We used a configuration with $`N_s=16`$, $`c=4`$ and $`d=10`$ (configuration I), and three configurations with $`N_s=24`$: one with $`c=1`$ and $`d=12`$ (configuration II), one with $`c=5`$ and $`d=12`$ (configuration III), and one with $`c=1`$ and $`d=16`$ (configuration IV). We have checked that the field strength from each of the lattice configurations has the essential properties described in reference  for the continuum field strength. We have also applied up to five cooling sweeps to the meron pair configurations in order to relax them close to lattice solutions, and checked that the monopole and vortex content for these meron pair configurations are independent of this cooling. Although we do not explicitly address Dirac zero modes in this work, note that the zero mode for a meron pair configuration has been calculated for a range of separations in reference  and displays two peaks at the positions of the merons.

## 3 Vortex content from Wilson loops

Before considering Laplacian Center Gauge fixing, it is useful to describe the vortex content obtained from calculating Wilson loops. For a single, singular meron at the origin, it has been shown that a circular Wilson loop around the origin in any of the six planes defined by a pair of coordinate axes ($`x,y,z,t`$) has the value $`1`$ for any size of the circle . Hence, Wilson loops indicate the presence of a vortex surface on all the planes defined by pairs of coordinate axes. For our configurations, regularized meron pairs, we studied two sets of Wilson loops. For the $`xy`$, $`xz`$ or $`yz`$ plane, planes orthogonal to the line joining both merons, we calculated a square Wilson loop of size $`r\times r`$ with one of the merons in the center of the loop. The results for configurations II and III (distance between merons $`d=12`$ and cap sizes $`c=1,5`$), and for the $`xy`$ plane, are shown in figure 1A. We see that at short distance, the value of the Wilson loop goes from $`+1`$ towards the value $`1`$, as for a single meron, and only changes this behavior at large distance where the contribution of the second meron starts to be significant, approaching the value $`+1`$ when the loop is bigger than the distance between the merons. We also see in figure 1A the effect of the cap size. The cap gives a characteristic size $`c`$ to the meron, which is reflected in the distance one must go for the value of the Wilson loop to start to be approach $`1`$ and thus enclose the vortex flux. Note that for the original singular meron, this size would be zero. The results obtained for the other two planes, $`xz`$ and $`yz`$, are the same, showing the underlying spherical symmetry in the spatial directions of the meron pair. Results for the other two configurations were completely analogous, with the curves simply reflecting the corresponding cap sizes and separation between merons.

For the $`xt`$, $`yt`$ or $`zt`$ planes, which include the line joining both merons, we calculated a rectangular Wilson loop of size $`r\times 2r`$ with one of the merons in the center of the loop. The results in the $`xt`$ plane for the same configurations as in figure 1A, are shown in figure 1B. We see that again at short distances, the Wilson loop goes from $`+1`$ to $`1`$ as the size of the loop increases, and as $`r`$ exceeds half the separation between merons, the Wilson loop begins to approach +1, which it will reach when both merons are included. Again, the loop must be larger than the cap size, $`c`$, to enclose all the vortex flux. As before, the results for the other two planes, $`yt`$ and $`zt`$, are the same, and the other configurations show analogous behavior reflecting the other cap sizes and separations.

The conclusion from this study of Wilson loops in a meron pair is that, like an isolated meron, a meron in a pair behaves like a source or sink for flux in non-trivial elements of the center of the group for all six planes defined by the Cartesian axes, and the size of the source or sink is of order of the cap size, $`c`$. Thus, each meron corresponds to the intersection of orthogonal pairs of vortices.

## 4 Monopole and vortex content from <br>Laplacian Center Gauge defects

In this section, we present the monopole and vortex content of the meron pair configurations described in the previous section, as inferred from the points at which gauge fixing to the Laplacian Center Gauge is ill-defined.

Fixing the gauge to Laplacian Center Gauge  involves the use of the two eigenvectors with lowest eigenvalues, $`\psi _1^a(n)`$ and $`\psi _2^a(n)`$, of the Laplacian operator,

$$_{nm}^{ab}(R)=\underset{\mu }{}\left(2\delta _{nm}\delta ^{ab}R^{ab}(n,\mu )\delta _{m,n+\widehat{\mu }}R^{ba}(m,\mu )\delta _{n,m+\widehat{\mu }}\right)$$
(4)
in the presence of a gauge field $`R^{ab}(n,\mu )`$ in the adjoint representation of the gauge group. The lowest eigenvector, $`\psi _1^a(n)`$, is rotated to the ($`\sigma _3`$) direction in color space. This step fixes the gauge up to the abelian subgroup of the SU(2) group. The U(1) abelian freedom is fixed by imposing the additional condition that the $`\psi _2^a(n)`$ eigenvector is rotated to lie in the positive ($`\sigma _1,\sigma _3`$) half-plane. After these two steps, the gauge is completely fixed up to the center degrees of freedom.

Monopoles and vortices are found in Laplacian Center Gauge as defects of the gauge fixing procedure, which means we have to look at the points at which the gauge fixing prescription is ill-defined. The first step, rotation of the first eigenvector to the third direction in color space, is ill-defined if $`\psi _1^a(t,x,y,z)=0`$. This defines lines in four-dimensional space and these lines are identified as monopole lines. The second step, rotation of the second eigenvector to the positive $`(\sigma _1,\sigma _3)`$ half-plane, is ill-defined at points at which the first and second eigenvectors are parallel. This condition defines surfaces in four dimensional space and these surfaces are identified as vortex sheets.

To fix to the Laplacian Center Gauge we use the algorithm presented in  to calculate the lowest eigenvectors of the Laplacian operator. We calculate the four eigenvectors with lowest eigenvalues, and find that the three lowest eigenvalues are degenerate. With two vectors chosen from these three, or from linear combinations of these three, we can fix the gauge to Laplacian Center Gauge. Note that because of the degeneracy in the lowest eigenvalues, the monopole and vortex content is ambiguously defined, and in this work we will consider all the different monopole and vortex patterns that may be obtained from the lowest eigenvectors.

Before considering the monopole and vortex content of our meron pair configurations, it is useful to review the monopole and vortex content of two limiting cases, an instanton and a single meron. The eigenfunctions of the lowest state of the Laplacian for these two cases are known analytically. For an isolated instanton there are three degenerate eigenfunctions, a monopole is only be located at the origin since it is the only point at which the eigenfunctions vanish, and there are no vortices because the three eigenfunctions are always mutually orthogonal. For the single meron, there are four degenerate eigenfunctions and the monopole content depends on the choice of lowest eigenvector. For the functions given in  it is easy to see that the $`i^{th}`$ eigenvector has a monopole line along the $`i^{th}`$ coordinate axis. We may think of this monopole line as joining the meron at the origin with the second meron at infinity, and the different lines arising from the different eigenvectors or combinations of them simply reflects the fact that the second meron may be reached at any position on the sphere at infinity. The vortex content also depends on the choice of the two lowest eigenvectors. It is easy to see that if one takes the $`i^{th}`$ and $`j^{th}`$ eigenvectors, the vortex content is given by the plane generated by the $`i^{th}`$ and $`j^{th}`$ coordinate axis. Taking other combinations of these four eigenvectors produces more complicated results, like one of the examples presented in , in which one obtains three vortex sheets given by three planes intersecting at the origin. It is noteworthy that this construction never generates the expected geometry of two intersecting vortices which in turn contain monopole lines, revealing that the Laplacian Center Gauge defects do not provide a completely satisfactory picture even in this analytically solvable case.

We now consider the monopole and vortex content of our meron pair configurations extracted from the three degenerate eigenvectors of the Laplacian. We will show that there is a vortex surface that looks like an ellipsoid of revolution touching the center of each meron at each tip, and that the monopole loops lie on this surface. The zeros in figures 2B and 2C correspond to longitudinal and axial sections of this surface respectively. The first quantity we examine is the modulus $`\mathrm{\Psi }(n)=\sqrt{_{a=1}^3(\psi ^a(n))^2}`$ and we look for points at which $`\mathrm{\Psi }(n)=0`$. One of the problems we have to face is the different monopole pictures obtained by choosing different eigenvectors or a linear combinations thereof. We have looked for a combination in which the monopole content is particularly clear and found that it was useful to minimize the sum of the modulus in a specified plane by taking linear combinations of the three eigenvectors. Doing this, we have found a combination in which the first eigenvector, I, has a monopole loop in the $`xt`$ plane joining the two merons, the second eigenvector, II, has the same monopole loop but in the $`yt`$ plane and the third eigenvector, III, has the same monopole loop in the $`zt`$ plane. A picture of the loop coming from eigenvector I, the loop in the $`xt`$ plane, can be seen in figure 2B, which is explained in more detail below. The vortex content obtained from eigenvectors I, II and III is also quite clear. If we choose eigenvectors I and II as the two lowest eigenvectors, the vortex surface looks like an ellipsoid of revolution around the $`t`$ axis in the $`xyt`$ coordinates, with this surface including the monopole loop in the $`xt`$ plane coming from eigenvector I and the monopole loop in the $`yt`$ plane coming from eigenvector II. If we choose eigenvectors I and III as the two lowest eigenvectors, the vortex surface is the same but in this case in the $`xzt`$ coordinates and if we choose eigenvectors II and III we see the same surface in the $`yzt`$ coordinates.

To obtain the locus of all the points which can be monopoles or vortices for our meron configurations, we calculate the determinant of these three vectors, I, II and III, at each lattice point. First, note that if any of the vectors is zero, the condition to find monopoles, the determinant is zero. Second, note that if there is a linear combination between them giving a zero vector, the condition to find vortices, the determinant is also zero. Finally, note that the determinant is independent under linear combinations of the three vectors. Hence, all the points for our meron configurations that can be monopoles or vortices are determined by the condition that the determinant vanishes.

The result we obtain is the following. We find a region on the lattice in which the determinant is always positive and another one in which it is always negative, both regions separated by a three-dimensional volume in which the determinant vanishes and defines all positions that can be monopoles or vortices. We describe the shape of this volume by showing some of its two-dimensional sections. First, we show its temporal dependence. Consider the determinant as a function of $`x`$ and $`t`$, for values of $`y`$ and $`z`$ fixed to the values that maximize the action density for these coordinates. Figure 2B shows the absolute value of the determinant (raised to the $`1/4`$ power to see the curve more clearly). The curve on which the determinant vanishes is similar to an ellipse, and joins the two maxima of the action density in the $`x`$, $`t`$ coordinates. To see that this loop joins the maxima in the action density to within a fraction of a lattice spacing, we show in figure 2A the action density as a function of $`x`$ and $`t`$, and for the same fixed values of $`y`$ and $`z`$. This loop is the same as the monopole loop described above defined by the points at which the modulus of eigenvector I vanishes. If we look at the determinant as a function of $`y`$ and $`t`$ (or $`z`$ and $`t`$), and for values of $`x`$ and $`z`$ (or $`x`$ and $`y`$) fixed to the maxima in the action density for these coordinates, we again obtain the same curves shown in figures 2A and 2B with $`x`$ interchanged with $`y`$ (or $`z`$). This curve in the $`y`$ and $`t`$ coordinates is also the monopole loop defined by the points at which the modulus of eigenvector II vanishes (and the curve in the $`z`$ and $`t`$ coordinates is the monopole loop defined by eigenvector III).

Second, we show the spatial dependence for fixed time positions. If we label the temporal positions for the maxima in the action density for each meron as $`t_1`$ and $`t_2`$, we find that for values of the lattice position $`t`$ belonging to $`[0,t_1]`$ and $`[t_2,N_t]`$, the determinant is always positive for all points in the three-dimensional lattice defined at each temporal point. For values of $`t`$ between the two merons, $`t_1<t<t_2`$, we find a spherical surface in which the determinant changes sign, and inside the sphere the determinant is negative. The radius of this sphere at each temporal point may be seen in figure 2B, where it corresponds to half the width of the monopole loop at each temporal point. A section of this sphere is plotted in figure 2C, which shows the absolute value of the determinant $`D(t,x,y,z)`$ to the $`1/4`$ power as a function of $`x`$ and y, with $`z`$ fixed to the value that maximizes the action density and $`t`$ fixed midway between the two merons. As claimed, this section is clearly observed to be circular.

We have obtained analogous results for all four configurations we have studied. The locus of all the points that can be monopoles or vortices is a three-dimensional volume as described above, joining the two meron components. The width of this volume joining the merons increases with increasing separation between the merons, and the maximum width at the midpoint is approximately 4, 4.5, and 6 for $`d`$ = 10, 12, and 16 respectively. The only effect of a few cooling sweeps applied to these configurations is a small change in the positions and widths of the merons, and the resulting monopole and vortex content is the same as described above relative to the new positions of the merons.

Finally, it is interesting to consider how the monopole and vortex content we found for the meron pair connects with the two limiting cases discussed before, an isolated instanton and a single meron. If we think of an instanton as a meron pair and dissociate it into two separated merons (keeping one of them fixed at the origin), we see that the monopole content of the instanton, a point at the maximum, becomes a loop of the form shown in figure 2B joining both merons. As we continue to separate both components, this loop becomes larger, and when one of the components approaches infinity, this loop joins the meron at the origin with the meron at infinity, the result we have already discussed for the single meron solution. The same argument is valid for the vortex content. As we separate the merons, a vortex surface joins both components, and as the separation approaches infinity, in the vicinity of one meron the surface locally looks like the planar vortex surface of a single meron.

## 5 Conclusions

We have investigated the monopole and vortex content of a meron pair by calculating the points at which the gauge transformation fixing the gauge to the Laplacian Center Gauge is ill-defined. Threefold degeneracy of the lowest eigenvalues of the Laplacian allows the choice of different pairs of vectors to define the Laplacian Center Gauge, giving rise to different pictures for the monopole and vortex content of the configuration. The determinant of these three eigenfunctions at each lattice point defines the locus of all the points that can be monopoles or vortices. This locus is a three-dimensional volume joining both merons, and at each time plane between the merons, the locus is the surface of a sphere with its center on the line connecting the merons, and with a diameter given by the width of the curve shown in figure 2B. One particular choice of degenerate eigenvectors has a monopole line joining both merons in the $`xt`$ plane for the first vector and in the $`yt`$ plane for the second vector. The corresponding vortex surface looks like an ellipsoid of revolution around the $`t`$ axis in the $`x,y,t`$ coordinates. Many other choices of two combinations of the three degenerate eigenvectors are possible, but all monopole lines and vortex surfaces must lie in the volume where the determinant vanishes. In particular, this implies that at the position of the meron, the vortex must always be in a purely spatial plane and can never be in a space-time plane.

We have also investigated the vortex content of the meron pair by calculating Wilson loops in all Cartesian planes containing the merons. This calculation showed that as in the case of an isolated meron, each meron in a pair behaves like a source or sink for flux in non-trivial elements of the center of the group for all six planes defined by the Cartesian axes and thus corresponds to the intersection of orthogonal pairs of vortices. Thus the Wilson loops imply that in addition to a vortex at the position of a meron in a purely spatial plane, there must also be a second vortex in a space-time plane.

Although these two complementary investigations have provided interesting insight into the vortex structure of a meron pair, it is clear from comparing the results that Laplacian Center Gauge fixing is not a sufficiently powerful tool to reveal the full structure of intersecting vortices. Whereas Wilson loops clearly imply the intersection of both spatial and space-time vortices at the merons, Laplacian Center Gauge fixing only finds vortex surfaces joining the two merons that intersect the merons in spatial planes. We note that the high symmetry of the background field produces a highly atypical situation including, for example, the intersection of monopole loops and a high degeneracy of equivalent solutions. It is possible that the introduction of a small perturbation would not only remove the intersections and degeneracy, but also produce a more generic situation of intersecting vortices. If this is not the case, more powerful techniques will be required to fully analyze the vortex structure.

Finally, looking at the combination of the results we obtain for the meron pair from Wilson loops and Laplacian Center Gauge fixing, it is reasonable to conclude that in a pair as well as in isolation, a meron is a localized source of monopole trajectories and a localized object with topological charge $`\frac{1}{2}`$ carrying center flux in six orthogonal space-space and space-time planes.


text:
# (Quantum) Chaos in BECs

The achievement of Bose-Einstein condensation of a dilute vapour of Rubidium atoms in 1995  heralded the beginning of a rapidly growing field of experimental and theoretical endeavour, each of which has further stimulated the other. In the intervening years, the experimental production of Bose-Einstein condensates (BECs), particularly of Rubidium and Sodium vapours, has become almost routine, and attention has shifted to either the condensation of more exotic species (of particluar note is the recent successful condensation of metastable helium ); or to deeper experimental investigations of the properties of an atomic BEC, up to the point of viewing a BEC more as a tool rather than an end in itself. The purpose here is to show links between chaotic dynamics, both in a classical and a quantum mechanical sense, and the possibly “chaotic” dynamical and stability properties of BECs.

To begin: considering only the dynamics of a single classical point particle of mass $`m`$, one can reasonably speak of a class of “chaos-inducing” potentials $`V_c(𝐱)`$, such that a particle can exhibit chaotic dynamics under the influence of such a potential. One can qualitatively visualize chaotic dynamics as being highly irregular motion in phase space. For a particle moving in one spatial dimension this is a two-dimensional space, where (subject to canonical changes of variables) the axes correspond to the particle’s position $`x`$ and momentum $`p`$. Such highly irregular motion is generally associated with exponential sensitivity to initial conditions, i.e. small differences in a point particle’s initial position and momentum can make a large difference in its long-time behaviour. The appropriate classical (Hamilton’s) equations of motion can then be derived from a Hamiltonian function of the general form:

$$H=\frac{𝐩^2}{2m}+V_c(𝐱).$$
(1)
One can readily make the jump to quantum mechanics by replacing the number-valued position and momentum variables with operator quantities. For a single particle, the dynamics are then generally most conveniently described by the following Schrödinger equation:

$$i\mathrm{}\frac{\psi (𝐱)}{t}=\left[\frac{\mathrm{}^2^2}{2m}+V_c(𝐱)\right]\psi (𝐱).$$
(2)
This is clearly a linear partial differential equation, where the time-derivative of the wavefunction $`\psi (𝐱)`$ is proportional to the Hamiltonian operator \[the momenta $`𝐩=(p_x,p_y,p_z)`$ have been replaced by differential operators\] acting on $`\psi (𝐱)`$. As chaos is generally associated with nonlinear differential equations, this can cause some to question how one can speak of quantum chaos at all. In classical mechanics we are interested in the motion of the particle in phase space, and thus in the changes of its position and momentum with time. The classical dynamics of a point particle’s position and momentum, and the dynamics of the quantum wavefunction $`\psi (𝐱)`$ are thus not really directly comparable: better analogies exist between Heisenberg’s equations and Hamilton’s equations of motion, or alternatively between Schrödinger’s equation and Liouville’s equation. Here quantum chaos is considered to be the study of quantum dynamical systems, the classical limits of which are capable of exhibiting chaotic dynamics. That Schrödinger’s equation is a linear partial differential equation is in this context an irrelevancy; nevertheless it is interesting to consider the effect of introducing an explicit nonlinearity into a quantum-chaotic Schrödinger equation:

$$i\mathrm{}\frac{\psi (𝐱)}{t}=\left[\frac{\mathrm{}^2^2}{2m}+V_c(𝐱)+u|\psi (𝐱)|^2\right]\psi (𝐱).$$
(3)
One could regard such an equation as describing “nonlinear quantum chaos” or “nonlinear wave chaos.” Such nonlinear Schrödinger equations are relevant for nonlinear optics, but more pertinently are extensively used to describe BECs, when they are called Gross-Pitaevskii equations . A Gross-Pitaevskii equation is used as an approximate description of the dynamics of a large number $`N`$ of (bosonic) interacting quantum mechanical particles, which is more fully described by a quantum field; if nearly all of these particles can be considered to be in the same motional state $`\psi (𝐱)`$ (i.e. if we have a BEC), then it is possible to describe most of the underlying quantum field by a classical field, and to regard what is left over as being “small.” The Gross-Pitaevskii equation is then the equation of motion of this classical field. It is nevertheless sometimes necessary to consider the dynamics of the full quantum field; in second-quantized form, the Hamiltonian operator leading to such a Gross-Pitaevskii equation (with $`u=gN`$ and $`g`$ being a parameter essentially describing how strongly individual particles interact) is given by

$$\widehat{H}=d^3𝐱\widehat{\psi }^{}(𝐱)\left[\frac{\mathrm{}^2^2}{2m}+V_c(𝐱)+\frac{g}{2}\widehat{\psi }^{}(𝐱)\widehat{\psi }(𝐱)\right]\widehat{\psi }(𝐱),$$
(4)
where $`\widehat{\psi }^{}(𝐱)`$ and $`\widehat{\psi }(𝐱)`$ are bosonic field operators, which create and annihilate, respectively, a particle at position $`𝐱`$. In this context, one could thus possibly speak of “quantum field chaos” when considering the dynamics induced by such a Hamiltonian.

A brief discussion of integrability, a concept closely tied to the (non) chaoticity of dynamical systems, is now necessary . Recalling the classical dynamics of a single particle, such a system is considered integrable if there exist as many independent conserved quantities as there are motional degrees of freedom (i.e. the number of relevant spatial dimensions). If the system is integrable, then the phase-space dynamics of the particle are so restricted by the necessity of observing these conservation laws that chaos is impossible. To take the rather trivial example of the one-dimensional harmonic oscillator \[potential $`V(x)=m\omega ^2x^2/2`$\], the solution of Hamilton’s equations of motion (in slightly unconventional form) is given by:

$$\alpha (t)=e^{i\omega t}\alpha (0),$$
(5)
where $`\alpha `$ is a complex variable, defined as

$$\alpha =\frac{m\omega x+ip}{\sqrt{2m}}.$$
(6)
It is clear that the particle motion is restricted to closed circles in the two-dimensional phase space, making both “irregular” motion in phase space, and sensitivity to initial conditions completely impossible. In the language of Hamilton-Jacobi theory, $`|\alpha (t)|`$ is an action variable, which is conserved (and equal to the square root of the energy), and the dynamics are fully described by the time-evolution of the phase angle of $`\alpha (t)`$, which is then the canonically conjugate angle variable $`\theta (t)`$ to $`|\alpha (t)|`$, defined by

$$\theta =\mathrm{arctan}\left(\frac{p}{m\omega x}\right).$$
(7)
Similar analyses can be made for any one-dimensional, time independent system; due to the time-independence of the Hamiltonian, the energy is always a conserved quantity, which is enough to make a one-dimensional system integrable. Two- and three-dimensional systems will require additional independent conserved quantities to be integrable, and time dependence in a one-dimensional system generally lifts the property of integrability. Concepts of quantum integrability generally, by analogy, consider conserved observable quantities, which in quantum mechanics are described by operators ; it is possible to consider something different however.

A wavefunction $`\psi (𝐱)`$ can generally be decomposed in terms of some complete orthonormal basis :

$$\psi (𝐱)=\underset{k=0}{\overset{\mathrm{}}{}}c_k\phi _k(𝐱).$$
(8)
If additionally the $`\phi _k(x)`$ are eigenfunctions of the (assumed time-independent) Hamiltonian, then any time dependence of $`\psi (x)`$ is contained completely within the coefficients $`c_k`$: thus

$$c_k(t)=e^{i\omega _kt}c_k(0).$$
(9)
For example, in the case of a one-dimensional quantum harmonic oscillator, the $`\phi _k(x)`$ are Gauss-Hermite polynomials, and $`\omega _k=\omega (k+1/2)`$. There is an now an obvious similarity between equations (5) and (9). Extending the analogy, we can consider motion in an infinite-dimensional pseudo phase space of a “point particle” having “position” values

$$\chi _k=c_k+c_k^{},$$
(10)
and “momentum” values of

$$\rho _k=i(c_kc_k^{}).$$
(11)
Within this pseudo phase space however, as in the case of the classical harmonic oscillator, motion is severely restricted by the fact that all the $`|c_k|`$ are conserved . The $`|c_k|`$ (effectively the “action variables” in the pseudo phase space) in fact form a complete set of conserved quantities (i.e. there are as many conserved quantities as degrees of freedom), meaning that, in this sense, the Schrödinger equation is always integrable and thus never chaotic, a fact which is indeed a direct consequence of the Schrödinger equation’s linearity.

For a general Gross-Pitaevskii equation, the additional nonlinear term means that integrability clearly cannot be taken for granted. It is, for example, known that the one dimensional, homogeneous case

$$i\mathrm{}\frac{\psi (x)}{t}=\left[\frac{\mathrm{}^2}{2m}\frac{^2}{x^2}+u|\psi (x)|^2\right]\psi (x)$$
(12)
is integrable , and that the conserved quantities can be determined using inverse scattering techniques. However, even for the seemingly trivial extension of simply adding a harmonic potential (generally more relevant to current experiments), it is unknown whether the resulting Gross-Pitaevskii equation is integrable or not; the application of a $`\delta `$-kicking potential (as is widely considered in the study of classical and quantum chaos)

$$V_c(x)=K\mathrm{cos}(kx)\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\delta (tnT),$$
(13)
will produce a Gross-Pitaevskii equation which is unlikely to be integrable . This general lack of integrability therefore raises the possibility of chaotic dynamics in the pseudo phase space, including extreme sensitivity to initial conditions.

In classical machanics, sensitivity to initial conditions is generally gauged by the Lyapunov exponent, so that a positive Lyapunov is considered a strong indicator of chaotic dynamics. The Lyapunov exponent is concerned with the growth with time of the euclidean “distance” $`d`$ in phase space between two initially close point particles. In one dimension, $`d`$ is thus given by

$$d(t)=\sqrt{[x_1(t)x_2(t)]^2+[p_1(t)p_2(t)]^2}.$$
(14)
In the pseudo phase space produced by some orthonormal representation of two initial wavefunctions $`\psi _1(x)`$ and $`\psi _2(x)`$, the equivalent $`d`$ is given by

$`d(t)`$ $`=`$ $`\sqrt{{\displaystyle \underset{k}{}}[\chi _{k1}(t)\chi _{k2}(t)]^2+[\rho _{k1}(t)\rho _{k2}(t)]^2}`$ (15)
$`=`$ $`\sqrt{2{\displaystyle d^3𝐱\left[\psi _1^{}(𝐱,t)\psi _2(𝐱,t)+\psi _2^{}(𝐱,t)\psi _1(𝐱,t)\right]}},`$
where in the case of the linear Schrödinger equation, by unitarity $`d^3𝐱\psi _1^{}(𝐱)\psi _2(𝐱)`$ cannot change, and $`d`$ will therefore be constant . Upon the addition of a nonlinearity this is of course no longer guaranteed, which in the case of the Gross-Pitaevskii equation can have important consequences for the stability of the condensate, due to a correspondance between the propagation of linearized perturbations around the condensate wavefuntion, and propagation (and growth) of the non-condensate fraction.

It is now necessary to return briefly to a fuller quantum mechanical description of the BEC. If the number of particles in the condensate mode is nearly $`N`$, the field operator $`\widehat{\psi }(𝐱)`$ can be written as

$$\widehat{\psi }(𝐱)\widehat{a}\left[\psi (𝐱)+\frac{1}{\sqrt{N}}\underset{k}{}\widehat{b}_ku_k(𝐱)+\widehat{b}_k^{}v_k^{}(𝐱)\right].$$
(16)
The $`\widehat{a}`$ operator annihilates a particle in the condensate mode $`\psi (𝐱)`$ (it is approximately this mode which is propagated by the Gross-Pitaevskii equation if the number of particles in the condensate mode $`\widehat{a}^{}\widehat{a}N`$), and the $`\widehat{b}_k^{}`$ and $`\widehat{b}_k`$ are creation and annihilation operators of excitations above the condensate, the $`u_k(𝐱)`$ and $`v_k^{}(𝐱)`$ describing the spatial dependence of these excitations. These excitation modes clearly describe non-condensate particles, and for the idealized case of zero temperature, the number of non-condensate particles is then given by

$$\underset{k}{}d^3𝐱|v_k(𝐱,t)|^2.$$
(17)
The excitation modes are propagated by slightly modified Bogoliubov equations, which are identical to the equations derived by considering perturbations around and orthogonal to $`\psi (𝐱)`$ propagated by the Gross-Pitaevskii equation up to linear order in the perturbation only. If one then considers

$$\psi _1(𝐱,0)=\psi (𝐱,0),$$
(18)
and

$$\psi _2(𝐱,0)=\psi (𝐱,0)+\frac{1}{\sqrt{N}}[u_k(𝐱,0)+v_k^{}(𝐱,0)],$$
(19)
it is clear that a reduction with time of $`|d^3𝐱\psi _1^{}(𝐱,t)\psi _2(𝐱,t)|`$ directly implies growth in the number of non-condensate particles, and, more generally, growth in the number of non-condensate particles is implied by any kind of linear instability in a solution of the Gross-Pitaevskii equation. Chaotic dynamics in the pseudo phase space thus imply the rapid increase in number of non-condensate particles, which in turn implies rapid depletion of the condensate. There is a caveat to this: the derivation of the equations used in the analysis formally take the limit $`N\mathrm{}`$, $`g0`$, while keeping $`gN`$ constant. There is therefore formally an infinite number of condensate particles, which stays infinite as the number of non-condensate particles increases, and the Gross-Pitaevskii equation is always valid.

In any real situation the total number of particles must of course be finite. Recalling that it is assumed that nearly all the particles are in the condensate mode $`\psi (𝐱)`$, it is clear that once depletion starts to become significant, this approach can be problematical. In the systematic (and equivalent) expansions of C.W. Gardiner , and of Castin and Dum , there are in principle further terms which can be incorporated to account for higher order effects. Morgan  has taken a slightly different approach by only requiring that the occupation of the condensate mode be large, i.e. it may be significantly different from the total particle number, and adding in perturbative corrections to the Gross-Pitaevskii and Bogoliubov equations. This treatment, however, is designed for studying essentially steady-state BECs at finite temperature, rather than the explicitly dynamical sorts of situations described here.

The theoretical treatment of such situations is clearly a difficult problem, it is however likely relevant to coming generations of BEC experiments. If one actually wants to do something with a BEC, this will involve dynamics. As generic Gross-Pitaevskii equations are presumably non-integrable, there exists the definite possibility of chaos and instability in some circumstances; this in turn implies the possibility of rapid depletion of the condensate mode, and concomitant loss of the desirable property of coherence of the macroscopic matter wave. Given this, it is equally desirable that such processes be better understood; this thus seems a suitable “quantum challenge” for the 21st century!