## Abstract

A novel kinetic beam scheme for the ideal quantum gas is presented for the computation of quantum gas dynamical flows. The quantum Boltzmann equation approach is adopted and the local thermodynamic equilibrium quantum distribution is assumed. Both Bose–Einstein and Fermi–Dirac gases are considered. Formulae for one spatial dimension is first derived and the resulting beam scheme is tested for shock tube flows. Implementation of high-order methods is also outlined. We only consider the system in the normal phase consisting of particles in excited states and both the classical limit and the nearly degenerate limit are computed. The flow structures can all be accurately captured by the present beam scheme. Formulations for multiple spatial dimensions are also included.

## 1. Introduction

The ideal classical gas dynamics can be described by the Maxwell–Boltzmann distribution, which corresponding to the lowest order solution of the classical Boltzmann equation when Chapman–Enskog procedure is employed (Chapman & Cowling 1970). The conservation laws based on the Maxwell–Boltzmann distribution is the well-known Euler equations of gas dynamics. In classical mechanics, a system of *N* identical particles which are distinguishable because it is possible, at least in principle, to label them according to their trajectories in phase space and they satisfy the Maxwell–Boltzmann statistics. In quantum mechanics, identical particles are absolutely indistinguishable from one another and an *N*-particle system can be described by a wave function with permutation symmetry. In nature, it is found that particles with antisymmetric wave functions are called fermions, which obey Fermi–Dirac statistics and particles with symmetric wave functions are called bosons, which obey Bose–Einstein statistics. The statistical properties of fermion and boson systems are profoundly different at low temperature. However, in the classical limit, both quantum distributions reduce to the Maxwell–Boltzmann distribution. Similar to the classical Boltzmann equation, a quantum Boltzmann equation for transport phenomenon can be developed for fermions and bosons, see Kadanoff & Baym (1962). This work considers numerical methods for solving the equilibrium limit solution of the quantum Boltzmann equation.

An interesting explicit scheme, which is called the beam scheme has been presented by Sanders & Prendergast (1974) for solving the equilibrium limit of the classical Boltzmann equation. The derivation was based on the local thermodynamic equilibrium Maxwell–Boltzmann distribution and resulted in a novel method for solving the transport processes governed by the Euler equations of Newtonian gas dynamics. Later, the beam scheme was successfully extended to relativistic Boltzmann transport equation based on the Jüttner distribution by Yang *et al*. (1997). In the beam scheme, a presumed local thermodynamic equilibrium distribution function is approximated by a combination of several distinct Dirac delta functions or discrete beams of particles in each cell. These beams are permitted to move over a time-step transporting mass, momentum and energy into adjacent cells. The motion of each beam is followed to first-order accuracy. The transport is taken into account to determine the new mass, momentum and energy in each cell; and these macroscopic moments are used to describe the new local equilibrium distribution for each cell. The entire process is then repeated and advanced to the next time-step. The choice of the size of the time-step is dictated by the Courant–Friedrichs–Lewy stability condition that physically no beam of gas particles travels farther than one cell spacing in one time-step. The beam scheme, although it is a particle scheme, has all the desirable features of modern characteristics-based wave propagating numerical methods for hyperbolic conservation laws of gas dynamics. For further details and discussion of the beam scheme, see Sanders & Prendergast (1974) and Yang *et al*. (1997). Both the beam scheme and the lattice Boltzmann methods (Benzi *et al*. 1992; Succi 2001) are gas kinetics based schemes for solving hydrodynamics. A study on the connection between the beam scheme and the lattice Boltzmann method has been given by Xu & Luo (1998), which analysed and compared both methods in detail. The distinct difference between the two methods lies in their equilibrium distribution function. The lattice Boltzmann equation expands the equilibrium at zero mean velocity and uses a polynomial (of mean velocity) to approximate the Maxwellian, while the beam scheme obtains particle beams around the average velocity of the Maxwellian distribution thus avoiding the low Mach number expansion. The basic beam scheme is first-order accurate and thus numerically diffusive like most of the exact or approximate Riemann solvers for the classical gas dynamics (Hirsch 1988; Toro 1999). However, this numerical diffusion of first-order upwind methods can be significantly reduced by the implementation of higher-order methods. A class of high-order methods using essentially non-oscillatory methods were also presented for the beam scheme (Yang *et al*. 1997).

In this work, we shall adopt the concept of beam scheme to devise a numerical method for the computation of ideal quantum gas dynamics. We consider the equilibrium limit of the quantum Boltzmann transport equation, in which the particles obey the quantum statistics and the particles of the system are in the excited states. It is well known that the boson or fermion particles in the excited states are in the gas phase and in the low-density (or high-temperature) limit (*βμ*≈−∞), they behave like classical ideal gas governed by the Maxwell–Boltzmann distribution. Here, *μ* is the chemical potential and *β*=1/*k*_{B}*T*, with *k*_{B} the Boltzmann constant and *T* the temperature. At low temperatures, in the neighbourhood of degeneracy temperature, *T*_{0}, corresponding to , (here, *n* is the number density and *λ*_{B} is the thermal de Broglie wavelength), the quantum effects become important (onset of quantum effects). The present method considers both the Bose–Einstein and Fermi–Dirac gas behaviours in the regime between the nearly degenerate (onset of quantum effects) limit and the classical limit. In Bose–Einstein statistics, the condensed phase of boson gas is not considered. Thus, similar to the ideal classical gas dynamics based on Maxwell–Boltzmann distribution, we can claim that the present method is for ideal quantum gas dynamics based on the Bose–Einstein and Fermi–Dirac distributions. Considering the existing extensive works on experimental and computational classical gas dynamics, the present quantum beam scheme can be potentially useful for revealing various dynamical aspects of ideal quantum gas through mathematical and physical analogies. We are also alert to the limitation of the present beam scheme, which considers only ideal non-interacting quantum particles. Interacting quantum particle systems which include Bose–Einstein condensates and ultracold fermions are beyond the scope of the present study. An excellent review on the development of numerical techniques for solving the mean-field problems and for studying the strongly correlated regimes of atomic quantum gases has been given by Minguzzi *et al*. (2004).

The paper is organized as following. We first briefly describe the elements of quantum Boltzmann transport equation in §2. In §3, the basic kinetic beam scheme in one spatial dimension based on the local thermodynamic equilibrium Bose–Einstein distribution and Fermi–Dirac distribution are derived. Implementations of high-order methods using weighted essentially non-oscillatory (WENO) schemes are outlined in §4. In §5, numerical experiments with quantum shock tube flows in both classical and nearly degenerate limits are given to illustrate the method. Lastly, discussion and concluding remarks are given in §6. Formulations for multiple spatial dimensions are given in appendix A.

## 2. Elements of quantum Boltzmann equation

In this section, we briefly describe the elements of quantum Boltzmann transport equation appropriate for the development of the present work. Following Kadanoff & Baym (1962), we consider the Boltzmann equation,(2.1)where *m* is the particle mass, *U* is the externally applied field and *f*(, , *t*) is the distribution function which represents the average density of particles with momentum at the space-time point , *t*. The (*δf*/*δt*)_{collision} denotes the collision term. A formal solution procedure which generalizes the Chapman–Enskog method to solve equation (2.1) was given by Uehling & Uhlenbeck (1933), where the first and second approximations of the distribution function and expressions for the viscosity and heat conductivity coefficients were given. In recent years, the development of numerical methods for solving the quantum Boltzmann equation has become an active research subject (Markowich & Pareschi 2002; Garcia & Wagner 2003). A recent review of the numerical methods for the dynamics and transport of atomic quantum gases has been given (see Minguzzi *et al*. 2004 and references therein). Here, we consider only the lowest order (first approximation) of solution of the above Boltzmann equation, requiring that(2.2)The solution to equation (2.1) is given by:(2.3)where *θ*=+1 denotes the Fermi–Dirac statistics and *θ*=−1 the Bose–Einstein statistics. To complete the equilibrium solution we have to determine the five unknown functions *T*(, *t*), *μ*(, *t*) and (, *t*), which appear in equation (2.3). These five flow parameters can be determined by making use of the conservation laws for number of particles, momentum and energy. These conservation laws are obtained by multiplying equation (2.1) by 1, or ^{2}/2*m*, and then integrating the resulting equations over all . The integrals of the collision terms in all three cases vanish automatically and we have the differential conservation laws for the conserved macroscopic quantities, i.e. the particle number density *n*(, *t*), the momentum density, =*m*, and the energy density, as follows:(2.4)(2.5)(2.6)where(2.7)(, *t*) is the number density flux,(2.8)and(2.9)Other higher-order moments can also be defined such as stress tensor and the heat flux vector. For the local equilibrium solution, one can obtain these macroscopic quantities in closed form in terms of the Bose or Fermi function (Huang 1987; Pathria 1996). In this work, to derive the basic quantum beam scheme and without loss of generality, we shall first neglect the effect of the externally applied field *U*(, *t*). To illustrate the method, we first consider the following local equilibrium Bose–Einstein distribution in one spatial dimension:(2.10)where is the fugacity and *u*_{x}(*x*, *t*) is the mean velocity. Then the number density *n*(*x*, *t*) is given by(2.11)the momentum *j*(*x*, *t*),(2.12)and the energy density *ϵ*(*x*, *t*),(2.13)where is the thermal wavelength and *β*=1/*k*_{B}*T*(*x*, *t*).

In the above, *g*_{ν}(*z*) denotes the Bose function of order *ν* which is defined by(2.14)We note that Bose functions with *ν*≤1 diverge as *z*→1 and *g*_{3/2}(*z*) remains finite at *z*=1 but with infinite slope while *g*_{ν}(*z*) with *ν*>3/2 have finite values and slopes at *z*=1.

In the case of Fermi gas, we start with the Fermi–Dirac distribution. The definitions for the number density, momentum and energy densities for the Fermi–Dirac gas are completely identical to the above procedures except that we only need to replace the Bose function with Fermi function which is defined by(2.15)We note that fermi functions *f*_{ν}(*z*)→*z* for small *z*. Since the chemical potential for a dilute system is large and negative, the classical limit corresponds to small *z*. Conversely, at low temperature the chemical potential approaches the Fermi energy so that the fugacity for a Fermi system is large. An asymptotic expansion of the Fermi function for large fugacity *z* is due to Sommerfeld.

Since we will derive beam schemes from one to three spatial dimensions and they involve the Bose (or Fermi) functions *g*_{d/2}(*z*), *g*_{d/2+1}(*z*), and *g*_{d/2+2}(*z*) (or *f*_{d/2}(*z*), *f*_{d/2+1}(*z*), and *f*_{d/2+2}(*z*)) for the beam scheme in *d* spatial dimensions (*d*=1, 2, 3). It is also noted that thermodynamic functions and the Bose–Einstein condensation of an ideal (non-interacting) gas of *N* bosons are peculiarly related to the space dimensionality (de Groot *et al*. 1950; Aguilera-Navarro *et al*. 1999). Correspondingly, it is also noted that thermodynamic functions and the Fermi–Dirac degeneracy of an ideal gas of *N* fermions are peculiarly related to the space dimensionality (Grether *et al*. 2003). Although these results look quite similar, the two functions behave differently in the nearly degenerate limit. For example, in boson gas, the fugacity *z* cannot be larger than 1 because of non-negative density while in Fermi gas such a restriction does not exist.

## 3. The kinetic beam scheme

In this section, the basic first-order beam scheme is derived. We consider the problems in one spatial dimension in which the macroscopic quantities are assumed to be functions of *x* and *t*. Divide the space into a number of cells of size Δ*x*_{i}=*x*_{i+1}−*x*_{i}. Without loss of generality, we assume that the cells are of equal size Δ*x*. The local state of gas in each cell *i* at any particular time *t* is specified by the three conserved macroscopic properties , which are the mass density, momentum density and energy density in cell *i*, respectively. The fundamental approximation of the beam scheme of Sanders & Prendergast (1974) for classical gas dynamics is to approximate the Maxwell–Boltzmann distribution by a combination of a finite number of Dirac delta functions of discrete velocities (or discrete beams of particles) which will reproduce the appropriate moments of the original continuous distribution function. Here in the ideal quantum gas dynamics, we first approximate the distribution function *f*_{1}(*p*_{x}, *x*, *t*) in each cell *i* at a given time *t* by(3.1)The unknown parameters in each cell *i* in equation (3.1) are *a*_{i}, *b*_{i}, *p*_{0} and Δ*p*, and they are determined in such a way that the following macroscopic moments are preserved:(3.2)(3.3)(3.4)We have four unknowns, *a*_{i}, *b*_{i}, Δ*p* and *p*_{0}, but with only three equations, an auxiliary condition is needed to close the equations. Here, we employ the fourth-order moment,(3.5)With the four equations we can solve the coefficients *a*, *b*, *p*_{0} and Δ*p* as follows:(3.6)(3.7)(3.8)(3.9)After obtaining the parameters *a*, *b*, Δ*p* and *p*_{0}, we can determine the conservative quantities carried by each beam in cell *i* as *Q*_{σ,i}=(*R*_{σ,i}, *M*_{σ,i}, *E*_{σ,i}) with(3.10)(3.11)(3.12)where, denote *p*_{0,i}, *p*_{0,i}+Δ*p* and *p*_{0,i}−Δ*p*, respectively, and *c*_{σ,i}=*a*_{i} for *σ*=1 and *c*_{σ,i}=*b*_{i}, if *σ*=2, 3. We can write the macroscopic quantity of each beam as follows:It is noted that .

Denote Δ*u*=Δ*p*/*m* and define the velocity of each beam, *V*_{σ,i}, (*σ*=1, 2, 3) in cell *i* as(3.13)During an interval of time Δ*t*, these discrete beams will move and transfer mass, momentum and energy into adjacent cells. The time-step Δ*t* is subjected to the condition that no single beam can move farther than a cell width Δ*x* during Δ*t*, i.e.(3.14)which is the well-known Courant–Friedrichs–Lewy (CFL) stability condition.

During time-step Δ*t*, particle beam *σ* in cell *i* moves a distance either into the left cell *i*−1 or into the right cell *i*+1or remains in cell *i* depending on the sign of the beam velocity. Define the fraction that transfers from cell *i* into cells *k*=*i*±1 as(3.15)and the fraction that remains in cell *i* is(3.16)After the time-step Δ*t*, the new values of density, momentum and energy in each cell, , taking into account transfer of these quantities from adjacent cells are given by(3.17)Now the three conservation laws in one space dimension become(3.18)where(3.19)Here, *Q* is the conserved state vector, *F*(*Q*) is the flux vector, *ρ*(*x*, *t*)=*mn*(*x*, *t*) is the mass density and *p* is the gas pressure. Denoting *U* as the internal energy density, we have(3.20)At this stage, we have completed the description of the beam scheme for the Bose–Einstein gas. In the computation, we employ dimensionless quantities and they are defined by the following:where , and .

We get new macroscopic properties (*n*, *u*_{x}, *ϵ*) after one-step beam transfer. We also need new values of fugacity *z* and temperature (or *β*) to determine the beam parameters (involving Bose functions *g*_{ν}(*z*)) in each cell as required by equations (3.6)–(3.8). Here, we can solve the equation for fugacity *z* through the combination of the three macroscopic equations (in dimensionless form) which define , i.e. equations (2.11)–(2.13), in the following form:(3.21)The above equation can be solved by the numerical method. After obtaining fugacity *z*, we can get temperature *T* easily and we can advance the solution for the next time-step.

The derivations of the kinetic beam scheme for the Fermi–Dirac gas are completely identical to the above procedures except that we only need to replace the Bose functions by Fermi functions correspondingly.

We note that the external potential terms on the right-hand side of equations (2.5) and (2.6) can be treated as source terms and directly added without much modification to the basic beam scheme.

## 4. Implementation of WENO schemes

For the purpose of implementing high-order methods we can further express the basic quantum beam scheme defined by equation (3.17) in terms of numerical flux as follows:(4.1)where the numerical flux is given by(4.2)Here, the split fluxes are defined by(4.3)and(4.4)It is noted that the basic beam scheme turns out to be in the form of a flux split method.

The above scheme is of first-order accuracy and thus is very diffusive. In practical applications we need high-order methods. In this section, we adopt the WENO interpolation method (Liu *et al*. 1994; Jiang & Shu 1996) to the basic first-order quantum beam scheme to result in a class of high-resolution methods for the computation of quantum ideal gas dynamical flows. The numerical flux defined in equation (4.2) is further expressed as(4.5)Here, we consider a fifth-order accurate (*r*=3) WENO scheme for the spatial difference of numerical fluxes. The WENO scheme for *r*=3, denoted as WENO3, can be expressed as(4.6)(4.7)(4.8)(4.9)The numerical flux for the negative part is given by(4.10)(4.11)(4.12)(4.13)Similar expressions for the case of third-order accurate WENO scheme (*r*=2), denoted as WENO2, can be given.

## 5. Numerical experiments

In this section, we report some numerical experiments to illustrate the present quantum beam scheme. We consider ideal quantum gas flows in a shock tube. Both Bose–Einstein and Fermi–Dirac gases are considered. In this problem a diaphragm, which is located at *x*=0.5, separating two regions, each remains in a constant equilibrium state at initial time *t*=0. The macroscopic properties on both sides of the diaphragm are different. Then the diaphragm is removed and flow structures consisting of moving shock wave, contact line and expansion fan, are generated which are similar to that in the classical gas dynamics. We consider two cases for each quantum gas. The first two cases are set up for boson gas in such a way that one corresponds to the classical limit and the other is set up with conditions that *z* is very close to unity which corresponds to the nearly degenerate limit. The third and fourth cases are set up for fermion gas. One corresponds to the classical limit and the other is set up for large *z*, close to the nearly degenerate limit.

## Example 1 Boson gas in near classical limit

The initial condition is specified as (*m*=1 assumed) =(2.33756, 0, 0.62985) for 0<*x*≤0.5 and =(1.24938, 0, 0.248107) for 0.5<*x*<1. This corresponds to *z*_{L}=0.8 and *z*_{R}=0.7 and we assume that *μ*_{L}=*μ*_{R}. The solution behaviours are similar to that in the classical limit. The calculations were done using several different grid cells with cell size Δ*x*=0.005, 0.0025 and 0.00125, respectively. The results are output at time *t*=0.3 (dimensionless time) and the CFL number used was 0.7. The computed number density, velocity, energy, temperature, fugacity and chemical potential profiles are depicted in figure 1 for several grid systems. Open square denotes the solution using 200 cells, open triangle that using 400 cells, and the solid line the solution using the finest 800 cells. The flow structures of moving shock, contact line and expansion fan can be identified and captured by the present beam scheme. It can be seen from the figures that when the grid is getting finer, the profiles becoming crisper indicates the grid convergence of the solution. A plot of the Bose–Einstein distribution functions *f*_{1}(*p*_{x}, *x*, *t*) at several *x* stations at time *t*=0.3 is also shown in figure 2. Here, 20 lines are plotted for 0≤*x*≤1 with equal spacing Δ*x*=0.05. The profiles at most stations are close to Gaussian (Maxwell–Boltzmann) distribution in this near classical limit condition.

## Example 2 Boson gas in nearly degenerate limit

The initial condition in this case is set as =(16.2218, 0, 1.13583) for 0<*x*≤0.5, and =(7.76706, 0, 0.374899) for 0.5<*x*<1. We assume *μ*_{L}=*μ*_{R} and set the values of fugacity as *z*_{L}=0.99 and *z*_{R}=0.98 which are quite close to unity. Although in one spatial dimension Bose–Einstein condensation does not occur because the excited states can accommodate essentially infinite particle number due to that boson functions may diverge at *z*=1. This example was selected to simulate the situation close to degenerate limit. Again, the calculations were done using several different grid cells with cell size Δ*x*=0.005, 0.0025 and 0.00125, respectively. The results are output at time *t*=0.3 and the CFL number used was 0.7. The same macroscopic properties as in figure 1 are depicted in figure 3 for three grid systems (200, 400 and 800 cells). The finest grid solution is represented by solid lines and used for comparison. It can be seen from the figures that when the grid is getting finer, the profiles are becoming sharper which indicates the grid convergence of the solution. The general flow structures of shock wave, contact line and expansion wave can be captured quite well. Finally, a plot of the Bose–Einstein distributions of this nearly degenerate limit case at several stations at time *t*=0.3 is depicted in figure 4. It can be seen that the distributions steepen markedly with high peaks and narrower spreads which reflect the fact that the fugacity value *z* is close to unity in this case.

## Example 3 Fermion gas in near classical limit

The initial condition is specified as =(0.521880, 0, 0.319701) for 0<*x*≤0.5 and for 0.5<*x*<1. This corresponds to *z*_{L}=0.8 and *z*_{R}=0.7 and we assume that *μ*_{L}=*μ*_{R} and *m*=1. The solution behaviours are similar to that in the classical limit. The calculations were done using two grid cells with uniform cell size Δ*x*=0.005 and 0.0025, respectively. The results are output at time *t*=0.06 (dimensionless time) and the CFL number used was 0.7. The computed number density, energy density, temperature and chemical potential profiles are depicted in figure 5 for these two grid systems. Open square denotes the solution using 200 cells and solid line using 400 cells. The flow structures of moving shock, contact line and expansion fan can be identified and captured by the present beam scheme. It can be seen from the figures that the grid convergence of the solution is ensured. Computations using WENO2 and WENO3 schemes with 200 cells for example 3 are also given in figure 6. The improvement of accuracy and resolution of the flow structures as compared to the first-order method can be easily seen in every aspect. A plot of the Fermi–Dirac distribution functions *f*_{1}(*p*_{x}, *x*, *t*) at several *x* stations at time *t*=0.06 is also shown in figure 6*d*. Here, 20 lines are plotted for 0≤*x*≤1 with equal spacing Δ*x*=0.05 and the range of *p*_{x}, −5≤*p*_{x}≤+5. The profiles at most stations are close to Gaussian (Maxwell–Boltzmann) distribution in this near classical limit condition.

## Example 4 Fermion gas in nearly degenerate limit

The initial condition in this case is set as =(3.53544, 0, 23.7423) for 0<*x*≤0.5, and =(3.52148, 0, 23.9943) for 0.5<*x*<1. We assume *μ*_{L}=*μ*_{R} and set the values of fugacity as *z*_{L}=20 000 and *z*_{R}=1500, using Sommerfeld expansion to get Fermi function.

Again, the calculations were done using two grid systems with cell size Δ*x*=0.05 and Δ*x*=0.0025, respectively. The results are output at time *t*=0.05 (dimensionless time) and the CFL number used was 0.5. The same macroscopic properties as in figure 5 are depicted in figure 7 for two grid systems (200 and 400 cells). The fine grid solution is represented by solid lines. Again, the grid convergence of the solutions can be observed. The general flow structures of shock wave, contact line and expansion wave can be captured quite well. Computations using WENO2 and WENO3 schemes with 400 cells for example 4 are also given in figure 8. The improvement of accuracy as compared to the basic first-order method can be clearly observed. A plot of the Fermi–Dirac distributions of this nearly degenerate limit case at several stations at time *t*=0.05 is also depicted in figure 8*d*. Twenty lines with Δ*x*=0.05 are plotted and the range of *p*_{x} is −8≤*p*_{x}≤+8. The approach to the typical double step function shape distribution is quite obvious in this large fugacity limit.

Although only one-dimensional method and results are presented here, however, the extension of the present numerical scheme to multiple dimensions can be done directly using Strang-type dimensional splitting (see Strang 1968).

## 6. Concluding remarks

In this work, the basic first-order beam schemes of Sanders & Prendergast for Newtonian gas dynamics and that of Yang *et al*. for relativistic gas dynamics have been successfully generalized to the ideal quantum gas dynamics based on Bose–Einstein statistics and Fermi–Dirac statistics. The beam scheme is first derived for one spatial dimension and then tested and illustrated by computing shock tube flows in ideal quantum gases. The regimes in the classical limit and nearly degenerate limit are tested. Numerical experiments indicate that the flow structures such as shock wave, contact line and expansion fan can be properly captured using the present quantum beam scheme. The derivation of the quantum beam schemes in multiple dimensions is also given in appendix A. Higher-order methods using WENO scheme are implemented to yield a class of efficient and accurate quantum Euler solvers. The present method applies for Bose–Einstein gases and Fermi–Dirac gases at temperatures above the critical temperature, provided that the hydrodynamic limit is applicable. Our purpose here is simply to give all the elements of the underlying method of beam scheme, which will allow us to build more complex and sophisticated (say parallel) codes to treat more general problems later. For example, together with general coordinates system, one will be able to tackle practical flows in ideal quantum gas dynamics just as one does for the Euler equations in classical gas dynamics. We also note that the discrete beam transfer provides a natural way to implement parallel computing similar to that of the lattice Boltzmann method. Although we did not consider the external potential due to magnetic or optical fields that are used to cool and hold the atomic gas in experiments, we note that it can be directly added into the present beam scheme without much modification. This will be the subject of a future study.

## Acknowledgments

The authors thank Dr J. C. Huang for many fruitful discussions on high-order beam schemes. This work is done under the auspices of National Science Council, TAIWAN through grants NSC 92-2210-E002-039 and 93-2210-E002-040.

## Footnotes

- Received August 18, 2005.
- Accepted November 17, 2005.

- © 2006 The Royal Society