text:
# Self-consistent Coulomb effects and charge distribution of quantum dot arrays

## I Introduction

The strong localization of few electrons in three dimensions is one of the key features of quantum dots bim99  (QDs). This results in a variety of properties making them potential candidates for new semiconductor applications such as laser and memory devices gru00 ; fin98 ; yus97 . One of those properties is the singular density of states, which relies on the alignment of the energy levels in the whole QD ensemble. On the one hand this requires a fairly homogeneous QD size distribution, on the other hand the statistical occupation of electronic QD levels gru97  changes the energy of levels of the same dot and all neighboring QDs due to Coulomb interaction.

The amount of Coulomb repulsion between QD electrons, especially at different QDs, is expected to depend on the details of the complicated bulk charge distributions of the QD device. Additionally, the microscopic charging and discharging processes of the QDs depend on the energy differences between the QD levels and the non-equilibrium quasi-Fermi level, resulting in a very complicated self-consistent problem. The microscopic charge distribution in the QD array is never stationary due to the statistical character of microscopic processes, but for large ensembles of QDs or in a long time average a stationary mean QD sheet charge density is attained. Standard approaches use quasi one-dimensional models where the charge in the QDs is approximated by a homogeneous sheet charge bro98b ; bru99 . This is suitable to obtain the QD energy levels by a fit to experimental capacitance-voltage (CV) data wet00 , where a Gaussian distribution of QD energy levels was assumed. Nevertheless the interpretation of the fitted broadenings remains controversial, since they may result both from size fluctuations and inhomogeneous charging effects of the single QDs.

The model used in this work calculates self-consistently the charge distribution in an array of QDs, with respect to the three-dimensional potential distribution under bias conditions. From an ensemble average of many QD charge configurations, the mean stationary QD sheet charge is calculated, to obtain macroscopic quantities like the depletion capacitance and the CV characteristics. The use of the electrostatic Green’s function method makes this expensive self-consistent problem solvable and charging effects in the QD array can be investigated separately from structural fluctuations.

## II Theory

Since only a small but essential part of the structure is not translationally invariant in the lateral ($`x`$,$`y`$) directions, but most of the device is inhomogeneous only in the growth direction ($`z`$) the problem of solving the three-dimensional non-linear Poisson’s equation is divided up into calculating two one-dimensional potentials $`\mathrm{\Phi }`$ and $`\mathrm{\Phi }^{\text{BQ}}`$ and a suitable 3D Green’s function. The potential $`\mathrm{\Phi }`$ is introduced as the self-consistent solution of the non-linear one-dimensional Poisson’s equation

$`ϵ_0_z[ϵ(z)_z\mathrm{\Phi }(z)]=\rho (\mathrm{\Phi },E_{\text{Fn}},E_{\text{Fp}})\sigma ^{\text{QD}}\chi _{\text{QD}}(z)`$ (1)
and the current equations neglecting generation and recombination in the active region

$$\begin{array}{cc}& _zj_n[\mathrm{\Phi },E_{\text{Fn}}]=0\hfill \\ & _zj_p[\mathrm{\Phi },E_{\text{Fp}}]=0\text{,}\hfill \end{array}$$
(2)
where the QD charge is treated as a sheet charge $`\sigma ^{\text{QD}}`$ (per unit area) averaged over the lateral directions. The bulk charge density $`\rho `$ is a function of the local potential $`\mathrm{\Phi }`$, and the electron and hole quasi-Fermi levels $`E_{\text{Fn}}`$ and $`E_{\text{Fp}}`$. The electron current $`j_n`$ and the hole current $`j_p`$ are calculated within a drift-diffusion model. Ohmic contacts are used to obtain boundary conditions for $`\mathrm{\Phi }`$, $`E_{\text{Fn}}`$ and $`E_{\text{Fp}}`$. For further details see Ref. wet00 . Here, $`ϵ_0`$ and $`ϵ`$ denote the absolute and relative permittivities of the semiconductor material, and $`\chi _{\text{QD}}`$ is the characteristic function

$`\chi _{\text{QD}}(z)`$ $`=`$ $`1/h\text{for}z\text{QD layer}`$ (3)
$`\chi _{\text{QD}}(z)`$ $`=`$ $`0\text{else,}`$
where $`h`$ is the height of the QDs, accounting for the localization in $`z`$ direction.

$`\mathrm{\Phi }^{\text{BQ}}`$ is the potential which is felt by a single QD charge, originating from the bulk charge if the QD charges would not exist. $`\mathrm{\Phi }^{\text{BQ}}(z_{\text{QD}})`$, where $`z_{\text{QD}}`$ is the position of the QDs, can be conceived as a background potential for the QDs. $`\mathrm{\Phi }^{\text{BQ}}`$ is calculated as $`\mathrm{\Phi }`$ but without $`\sigma ^{\text{QD}}`$ from eqn. (1).

In this paper the intrinsic QD level energy $`E_{\alpha \mu }^{\text{intr}}`$ of a state $`\mu =(n_x,n_y,n_z)`$ in QD $`\alpha `$ (where $`n_x`$, $`n_y`$, $`n_z`$ are the respective quantum numbers) is measured from the conduction band edge. In order to specifically account for the size dependence of the level energies in InAs QDs we use a fit to the results of $`𝒌\mathbf{}𝒑`$ calculations sti99  as shown in Fig. 1.

The charging energy $`E_{\alpha \mu }^{\text{char}}`$ originating from the Coulomb interaction of the QD electrons is treated as a first order perturbation. Thus the QD electron energy can be expressed as

$$E_{\alpha \mu }^{\text{QD}}=E_{\text{C0}}(z_{\text{QD}})e\mathrm{\Phi }^{\text{BQ}}(z_{\text{QD}})+E_{\alpha \mu }^{\text{intr}}+E_{\alpha \mu }^{\text{char}}\text{.}$$
(4)
Here $`E_{\text{C0}}`$ denotes the intrinsic conduction band edge, and $`e>0`$ is the elementary charge. The Coulomb charging energy $`E_{\alpha \mu }^{\text{char}}`$ can be expressed by the use of the capacitance matrix $`C`$ wha96  as

$$E_{\alpha \mu }^{\text{char}}=e^2\underset{\beta \nu }{}p_{\beta \nu }\left(C^1\right)_{\alpha \beta }\left(1\delta _{\alpha \mu ,\beta \nu }\right)\text{.}$$
(5)
Here $`p_{\beta \nu }`$ denotes the occupation numbers of the QD state $`\nu `$ in QD $`\beta `$ ($`p_{\beta \nu }=1`$ if occupied, $`p_{\beta \nu }=0`$ if not occupied). $`\delta _{\alpha \mu ,\beta \nu }`$ ensures that an electron does not interact with itself. The charging energies obtained by summing over $`\beta =\alpha `$ will be referred to as intradot charging energy, where the sum over $`\beta \alpha `$ is referred to as interdot charging energy in this paper. The elements of the inverse of the capacitance matrix $`\left(C^1\right)_{\alpha \beta }^{\mu \nu }`$ can be expressed by the use of the Green’s function $`G`$ war98 :

$$\left(C^1\right)_{\alpha \beta }^{\mu \nu }\underset{V}{}\underset{V}{}\left|\psi _{\alpha \mu }(𝒓)\right|^2\left|\psi _{\beta \nu }(𝒓^{})\right|^2G(𝒓,𝒓^{})d^3𝒓^{}d^3𝒓$$
(6)
Here $`\psi _{\alpha \mu }`$ and $`\psi _{\beta \nu }`$ symbolize the QD electron wavefunctions. Assuming that the QD electrons are distributed homogeneously within cylinder shaped QDs with height $`h`$ and a diameter $`D_\alpha `$, neglecting the details of the QD wave functions leads to

$$\left(C^1\right)_{\alpha \beta }=\frac{1}{V_\alpha V_\beta }\underset{V_\alpha }{}\underset{V_\beta }{}G(𝒓,𝒓^{})d^3𝒓^{}d^3𝒓\text{,}$$
(7)
where $`V_\alpha `$ and $`V_\beta `$ denote the volumes of two QDs $`\alpha `$ and $`\beta `$. In our model the charging energies are the same for all states $`\mu ,\nu `$ within a particular QD. A Green’s function $`G`$ including the static screening effects by free bulk charges can be found by defining a screening parameter $`f_{\text{scr}}`$ as

$$f_{\text{scr}}(z)\frac{\rho }{\mathrm{\Phi }}(\mathrm{\Phi }^{\text{BQ}}(z))\text{.}$$
(8)
Then $`G`$ is calculated as lan91

$`\left[\mathrm{\Delta }_𝒓+_z\left(\mathrm{ln}ϵ(z)_z\right)+{\displaystyle \frac{f_{\text{scr}}(z)}{ϵ(z)ϵ_0}}\right]G(z,r,z^{})`$ (9)
$`={\displaystyle \frac{1}{ϵ(z)ϵ_0}}\delta ((0,z^{})(r,z))\text{,}`$
where $`\mathrm{\Delta }_𝒓`$ is the Laplace operator and $`r`$ is the radial distance in the $`xy`$ plane. The solution of eqn.(9) has a cylindrical symmetry, since $`f_{\text{scr}}`$ and $`ϵ`$ depend only on $`z`$. Note that $`G`$ is the potential of a unit point-like charge positioned at $`z^{}`$, $`r^{}=0`$.

The QDs are occupied statistically according to the level energies $`E_{\alpha \mu }^{\text{QD}}`$ and the contact quasi-Fermi level $`E_{\text{Fn}}^c`$ with the occupation probability

$$f=\frac{1}{1+\mathrm{exp}\left(\frac{E_{\alpha \mu }^{\text{QD}}E_{\text{Fn}}^c}{k_BT}\right)}\text{,}$$
(10)
where $`T`$ is the temperature and $`k_B`$ denotes Boltzmann’s constant. Note that every change in the occupation numbers $`p_{\alpha \mu }`$ changes all QD level energies. The stationary, mean QD sheet charge density $`\sigma ^{\text{QD}}`$ used in eqn. (1) is found as a charge configuration average

$$\sigma ^{\text{QD}}=\frac{e}{A}\frac{1}{N_i}\underset{i=1}{\overset{N_i}{}}\underset{\alpha ,\mu }{}p_{\alpha \mu }^i\text{,}$$
(11)
where $`i`$ denotes $`N_i`$ simulated realizations of configurations of occupied QD levels, and $`A`$ is the quadratic contact cross section simulated. Periodic boundary conditions are used in the four parallel and four diagonal directions in the QD plane. Since $`\sqrt{A}`$, which takes values from 1 $`\mu `$m to 3 $`\mu `$m in our simulations, is always larger than the maximum expected screening length for this structure, the influence of boundaries is negligible. In other words, a QD charge does not “see” itself in the periodically continued plane. The stationary depletion layer capacitance is determined as in Ref. wet00  by the electric field in the depletion layer given by $`\mathrm{\Phi }`$.

## III Results

The pn structure considered here consists of highly doped $`n^+`$ and $`p^+`$ GaAs contact regions, as well as a $`3\times 10^{16}cm^3`$ $`n^{}`$-doped GaAs region. A layer of self-organized InAs QDs is positioned 400 nm from the pn interface. The structure is displayed in the inset of Fig. 1. Similar pn structures have been investigated experimentally by CV spectroscopy and DLTS measurements kap99 ; kap00a . We simulate arrays of up to $`30`$ by $`30`$ QDs with a QD sheet density of $`N_{QD}=10^{10}cm^2`$ for various temperatures and ensembles of QDs. The average QD diameter is 12 nm. A QD height of 3 nm is used.

### III.1 Ideal dot distribution

When a reverse bias $`U`$ is applied to the diode, bulk carriers as well as QD states will be successively discharged wet00 . First we simulate the structure at 4 K, 35 K, 70 K, 150 K and 300 K with periodically ordered QDs and with all QDs having the same size. The resulting CV characteristics are shown in Fig. 2. The main drop of the capacitance with increasing reverse bias is due to the increasing bulk depletion region sze81 . On top of this background three peaks can clearly be observed for temperatures less than 150 K. The peak at $`U`$-2.2 V originates from the discharging of the excited states, which is split since piezoelectric fields break the degeneracy of the states (100) and (010) (see Fig. 1). At $`U`$-3.6 V the two ground state electrons are discharged successively leading to a double peak. Here, the energy degeneracy is broken by the Coulomb charging effects. The Coulomb splitting has also been observed experimentally for different structures med97 ; med97a ; luy99 . The peak structure is broadened with rising temperature, resulting in a plateau in the CV characteristics at 300 K. Our calculations use a charging energy of about 20 meV. From the ground state peak splitting in the 4 K curve of about 100 mV we can estimate a lever-arm factor of 200 meV/V which, when applied to the width of the ground state peaks, results in an energy broadening of 2.5 meV. This value coincides with the unscreened Coulomb potential (in units of energy) of two electrons at a distance of 100 nm from a QD (2 $`\times `$ 1.25 meV). This means that the broadening at 4 K originates from the charge fluctuations within the four nearest neighbor QDs, in particular if every second one is charged. At 300 K, where occupation is rather insensitive to the QD level spacing, the full plateau width of about 500 mV corresponds to an energy broadening of 100 meV. This value can be attributed to the broadening in InAs QDs of the Fermi distribution function of about $`2`$ $`k_BT`$ at 300 K of the two ground state levels. Consequently, at large temperatures, the charging energy is given by the averaged, statistically distributed charge stored in the QD array. This is the reason why the CV characteristics shows a flat plateau at 300 K and is hardly distinguishable from that of a quantum well wan96 .

Additionally, the ground state peaks show a shift to lower reverse bias for rising temperatures, where the excited states peaks do not show this feature. This shift can be explained as follows: When the ground states electrons are discharged, interdot charging energies become more important than for the excited states (explained later in detail). The interdot interaction depends on how the free electrons can screen the interdot Coulomb potential. This ability rises with lower temperatures leading to shorter screening lengths. Therefore the charging energies are lower at low temperatures, and a higher reverse bias has to be applied to discharge the ground state electrons.

### III.2 Impact of structural fluctuations

Now we investigate the impact of structural fluctuations of the QD ensemble such as position disorder and size fluctuations. At first we distribute the equally sized QDs randomly instead of periodically in the lateral plane leaving the QD sheet density constant. Compared to the characteristics without fluctuations (solid lines in Fig. 3) at 35 K and 70 K the ground state peaks for randomly distributed QDs are only slightly smeared out (dashed lines in Fig. 3). Since the electrons in the excited states can be thermally activated more easily due to fluctuations in the interdot Coulomb interaction, resulting in an inhomogeneous charging of these levels, the respective peaks are more strongly broadened.

Next we simulate an ensemble of periodically ordered QDs with random QD diameters between 10 and 14 nm at 70 K, corresponding to a fluctuation of $`\pm `$ 17 %. The resulting CV characteristic (see dash-dotted line in Fig. 3) shows a very broad ground state peak. Again, the excited states are more sensitive to the fluctuations, i.e., the respective peak disappears almost completely. The increase of the fluctuations of the diameter to $`\pm `$33 % (8 - 16 nm) leads to a broad ground state plateau in the characteristics (dotted lines in Fig. 3). Such a plateau has been observed experimentally for similar pn structures kap99 ; kap00a , and a comparison with our results indicates that the size fluctuations in self-organized QDs is of this order of magnitude. This size fluctuation corresponds to an energy fluctuation of the ground state of about 100 meV (Fig. 1), which coincides with the typical value for the broadening of the ground state in the photoluminescence signal of InAs QDs kap99 . For room temperature the CV characteristics with $`\pm `$33 % size fluctuation shows almost no deviation from the characteristics without fluctuations. Here, thermal charge fluctuations dominate the behavior.

As already mentioned, the QDs are discharged by applying a reverse bias. In Fig. 4 the mean number of electrons per QD is displayed as a function of bias for 35 K, 70 K and 300 K. At 35 K and 70 K the QDs are charged with three electrons at zero bias and are discharged monotonically in a step-like manner at biases where peaks are observed in the CV characteristics (solid lines). Structural fluctuations lead to a broadening of these curves (dashed, dash-dotted and dotted lines). Additionally, with size fluctuations of $`\pm `$33 % some QDs can store up to four electrons, leading to a higher average QD occupation. For 300 K the QDs are charged with less electrons, because of the lower degree of degeneracy of bulk electrons in the contact region, and consequently the lower lying quasi-Fermi level at higher temperatures.

### III.3 Variation of the charging energy

The mean charging energy per QD electron as shown in Fig. 5 shows an interesting, non-monotonic behavior. At 70 K without fluctuations (solid lines) the charging energy is at first constant for $`U<`$-0.5 V and then rises until the excited state electrons are rapidly discharged. Between $`U`$=-2.2 V and $`U`$=-3.6 V the charging energy again rises, and then drops to zero when the residual ground state electrons are discharged. The continuous increase of the charging energies can be explained by considering the bulk electrons in the vicinity of the QDs being depleted with reverse bias. This depletion increases the interdot charging energy, and therefore the total QD level energies. Fig. 5 shows that these negative and positive differential charging energies persist even for $`\pm `$33 % size fluctuations at 70 K. Also for 300 K this behavior can be observed, but it vanishes when taking into account additionally $`\pm `$33 % size fluctuations. The mean height of the potential barrier $`e[\mathrm{\Phi }(z_{\text{QD}})\mathrm{\Phi }^{\text{BQ}}(z_{\text{QD}})]`$ at the QD sheet as obtained from the one-dimensional Poisson’s equation is also displayed in Fig. 5. It shows the big difference of the charging energies between one-dimensional and three-dimensional calculations, especially if the QDs are charged with more than two electrons. But nevertheless, it shows the same differential behavior, which has an inhibitory effect on leakage and recombinations currents in QD devices, and therefore may provide an explanation for the observed bistabilities in similar pn structures kie03 .

To explain the origin of the negative differential charging energies in more detail, the charging energies are plotted in Fig. 6 as a function of bias, assuming that all QDs are charged with one electron at any bias and adding one extra electron into a single QD. We also differentiate between the intradot charging energy (solid line) originating from the second electron in a QD, and the interdot charging energy (dashed line) from the Coulomb repulsion with all neighboring QDs. The ratio of the Green’s function and the unscreened Green’s function ($`f_{\text{scr}}=0`$) is also displayed in the two insets of Fig. 6 for $`U`$=-1.5 V and $`U`$=-4.5 V as a function of the lateral position $`r`$ and the vertical position $`z`$ as contour plots. As shown in Section II, the Green’s function determines the capacitance matrix and therefore the charging energies. As the Green’s function spreads out more and more, especially in the lateral directions, the charging energies due to interdot charging effects become larger. The effect on the intradot charging energy is rather small, since the size of the QDs is small compared to the distance of the free electron gas and the respective screening length. But by depletion of the free electrons with the applied reverse bias, the interdot charging rises. The slope of the interdot charging energy can clearly be identified with the slope of the mean charging energy of Fig. 5 for biases where the QD charge is constant but non-zero.

Fig. 7 shows the QD charge distribution within the QD layer for $`U`$=-2.2 V, $`U`$=-2.8 V and $`U`$=-3.6 V and the energy distribution of the ground state (000) at $`U`$=0.0 V, $`U`$=-2.2 V, $`U`$=-2.8 V and $`U`$=-3.6 V for 70 K and without structural fluctuations. The normalizations of the distributions have been chosen differently for clarity. The QD charge distribution is inhomogeneous when the QDs are discharged at $`U`$=-2.2 V and $`U`$=-3.6 V, but it is homogeneous at $`U`$=-2.8 V. Fig. 7 shows that the initial ground state energy at $`U`$=0.0 V is lowered due to discharging one electron at $`U`$-2.2 V. Then the energy rises (up to $`U`$-2.8 V), due to the rising interdot charging energies as argued above. The discharging of the ground state electron then leads to a lowering of the QD electron energies again. Overall, the QD energy shows a zig-zag behavior when sweeping the bias.

The inhomogeneity of the QD charge leads also to a line broadening of the electron states. For $`U`$=-2.8 V and $`U`$=0.0 V the charge distribution is homogeneous, and the broadening vanishes. As shown above the magnitude of the Coulomb interaction between the QDs depends on the bias (Fig. 6). Therefore the line broadening at $`U`$=-2.2 V is still small even for an inhomogeneous QD charge, since the interdot interaction is close to zero. The broadening becomes larger (about 2 meV) at $`U`$=-3.6 V, since at this bias the charge distribution again is inhomogeneous, and the interdot interaction is large (see Fig. 6).

## IV Conclusion

In conclusion, we have proposed a model to calculate the self-consistent three-dimensional charge distributions and charging energies within a QD device. We have shown that the charging energies depend on the electronic structure, which is determined by doping concentration, bias and temperature. The effects of size fluctuations of the dots and the role of a spatially periodic or random distribution of the dots on the capacitance-voltage characteristics have been investigated. We have found that size fluctuations have a much stronger influence on the broadening of the peak structure of the CV characteristics than the positional disorder of the QDs. Our detailed simulations have shown that the QD energy can be a strongly non-monotonic function of the applied bias.

Acknowledgement

This work was supported by DFG in the framework of Sfb 296. Helpful discussions with C. Kapteyn and R. Heitz are acknowledged.


format: text

text:
# Splitting of a doubly quantized vortex through intertwining in Bose-Einstein condensates

## Abstract

The stability of doubly quantized vortices in dilute Bose-Einstein condensates of 23Na is examined at zero temperature. The eigenmode spectrum of the Bogoliubov equations for a harmonically trapped cigar-shaped condensate is computed and it is found that the doubly quantized vortex is spectrally unstable towards dissection into two singly quantized vortices. By numerically solving the full three-dimensional time-dependent Gross-Pitaevskii equation, it is found that the two singly quantized vortices intertwine before decaying. This work provides an interpretation of recent experiments \[A. E. Leanhardt *et al.* Phys. Rev. Lett. 89, 190403 (2002)\].

The stability of single-quantum vortices has been studied extensively after such vortices were first observed in dilute alkali atom Bose-Einstein condensates (BECs) vorform . Since a single-valued complex order parameter describes the state of the condensate, its phase must undergo a $`2\pi n`$ change along a loop encircling a vortex, where $`n`$ is the quantum number of the vortex. However, the creation of multiquantum vortices is impossible just by rotating the harmonic trapping potential at a high frequency, since the existence of many singly quantized vortices is energetically more favorable than a single multiquantum vortex energetics . Indeed, vortex lattices composed of single-quantum vortices were observed in such experiments vortexlattice .

Verifying the proposal of topological phase engineering by Nakahara et al. nakahara , Leanhardt et al. topological  have recently succeeded in creating vortices simply by reversing the bias magnetic field used to trap the condensate. The vortices created display winding numbers with $`n=2`$ ($`4\pi `$ phase winding) or with $`n=4`$ ($`8\pi `$ phase winding), depending on the hyperfine spin states used for 23Na condensates. It was confirmed that the axial angular momentum per particle is $`2\mathrm{}`$ ($`4\mathrm{}`$) for the doubly (quartically) quantized vortex.

In this experiment, after creating a vortex it is held for $`20\text{ ms}`$, watching the vortex core to split into single-quantum vortices. However, no splitting was observed in this time span. Since a doubly quantized vortex is expected to decay spontaneously into two singly quantized vortices owing to energetics, this observation seems to be puzzling. This motivates our investigation of the detailed dynamics of the decay process of multiply quantized vortices in view of the present experimental situation. In dilute BECs, these were the first experimentally realized double-quantum vortices which are known to exist in other superfluids, such as superfluid 3He-A nature . Therefore, we have a unique opportunity to examine the physics of multiply quantized vortices.

There are several theoretical investigations on the instability of the doubly quantized vortex: The instability due to the bound state in the vortex core was pointed out by Rokhsar rokhsar , who claimed that the decay of vortex states requires the presence of thermal atoms. On the other hand, Pu et al. pu  found that even in the absence of thermal atoms, *i.e.* in the non-dissipative system, appearance of the modes with complex eigenvalues leads to the spontaneous decay of multiply quantized vortices.

Here we first study the Bogoliubov eigenvalue problem for a doubly quantized vortex in a cylindrically symmetric system in which the existence of modes with complex eigenfrequency, *i.e.*, spectral instability, depends on the interaction strength. Then we investigate full three-dimensional cigar-shaped systems similar to the MIT experiments topological . It is found that the splitting of the vortex core is initiated in different positions along the vortex axis, resulting in entangling vortices.

Let us consider the spectral stability problem of the doubly quantized vortex by finding the spectrum of the quasiparticle excitations of the vortex state described by the order parameter $`\mathrm{\Phi }(𝐫)`$. We apply the Bogoliubov eigenvalue equations bogo

$`(\mathrm{}\omega _\alpha \mu )u_\alpha (𝐫)+g[\mathrm{\Phi }(𝐫)]^2v_\alpha (𝐫)`$ $`=`$ $`0,`$
$`g[\mathrm{\Phi }^{}(𝐫)]^2u_\alpha (𝐫)+(+\mathrm{}\omega _\alpha \mu )v_\alpha (𝐫)`$ $`=`$ $`0,`$ (1)
where $`\omega _\alpha `$ are the eigenfrequencies related to the normal-mode functions $`u_\alpha (𝐫)`$ and $`v_\alpha (𝐫)`$. The chemical potential for the system of $`N`$ bosons each having mass $`m`$ is denoted as $`\mu `$ and the interaction strength is described by $`g=4\pi \mathrm{}^2a/m`$, where the $`s`$-wave scattering length is denoted as $`a`$. The operator $`=\mathrm{}^2^2/2m+V_{\text{tr}}(𝐫)+2g|\mathrm{\Phi }(𝐫)|^2`$ consists of the kinetic energy, the trapping potential and the interaction term. In this paper, the trapping potential is of the form $`V_{\text{tr}}(𝐫)=m(\omega _x^2x^2+\omega _y^2y^2+\omega _z^2z^2)/2`$, where $`\omega _i`$ is the frequency of the trap in direction $`i`$. The doubly quantized vortex state created in the experiments topological  is obtained by finding the minimum-energy solution of the stationary version of the time-dependent Gross-Pitaevskii (GP) equation

$$i\mathrm{}\frac{\mathrm{\Phi }}{t}=\left[\frac{\mathrm{}^2}{2m}^2+V_{\text{tr}}+g|\mathrm{\Phi }|^2\mu \right]\mathrm{\Phi },$$
(2)
with the constraint $`|\mathrm{\Phi }(𝐫,t)|^2𝑑𝐫=N`$ and ansatz $`\mathrm{\Phi }(𝐫,t)=\varphi (r,z)e^{i2\theta }`$ footnote , where $`(r,\theta ,z)`$ denote cylindrical coordinates. The ansatz for the order parameter leads to constraints for the normal-mode functions $`u_\alpha (𝐫)=\stackrel{~}{u}_\alpha (r,z)e^{i(q+2)\theta }`$, $`v_\alpha (𝐫)=\stackrel{~}{v}_\alpha (r,z)e^{i(q2)\theta }`$. Since the Bogoliubov eigenvalue equation is non-Hermitian, the existence of complex eigenvalues is not ruled out. If such an eigenvalue exists, the state is said to be spectrally unstable and small perturbations of the order parameter may grow exponentially in time.

The Bogoliubov eigenvalue spectrum for double-quantum vortex states is computed for cylindrically symmetric condensates with no phase modulation along the $`z`$-direction using a harmonic trapping potential $`\omega _x=\omega _y=\omega _{}`$ and $`\omega _z=0`$. In this case, there exist complex-frequency eigenmodes with $`q=\pm 2`$ for certain values of the interaction strength. The imaginary part of the eigenfrequency is plotted in Fig. 1(a) as a function of the interaction strength $`an_z=a|\varphi (r)|^2𝑑x𝑑y`$.

The emerging complex eigenvalues imply that by adding an infinitesimal perturbation, a pair of modes with complex conjugate eigenvalues is created spontaneously and grow exponentially in time. In numerical calculations trap anisotropy $`\lambda `$ is used as the perturbation. One of these modes has $`q=2`$ and is confined into the vortex core. This mode describes the instability of the double-quantum vortex and leads to its splitting into a pair of single-quantum vortices rokhsar . The remaining mode with $`q=+2`$ is created in order to ensure energy conservation. When the strength of interactions is weak ($`an_z<3`$), this mode corresponds to the co-rotating quadrupole mode. As the interaction strength $`an_z`$ increases, however, this mode changes from a nodeless mode to the nodal mode along the $`r`$-axis and the real part of its frequency grows.

To verify the role of the complex modes, using a slightly anisotropic trap with $`\omega _y=\omega _x/\lambda =2\pi \times 250`$ Hz ($`\lambda =1.01`$), we first solved the time-dependent GP equation (2), which contains no term which can describe dissipative effects. The time evolution of the density profiles with $`an_z=1.5`$ is shown in the upper panels of Fig. 1(b). At $`t20`$ ms, the vortex splits into a pair of single-quantum vortices and the condensate simultaneously changes into quadrupolar shape. Due to the co-rotating quadrupole excitation with $`q=+2`$, the condensate rotates with a constant frequency and “ghost” vortices tsubota  nucleate near the condensate surface, which is clearly seen through the density modulation around the surface at $`t=36`$ ms. Owing to the absence of dissipation, the pair of single-quantum vortices precesses around the trap center and cannot escape from the trap.

For the interaction strength corresponding to the experiments ($`an_z=13.75`$), however, it is seen that within the time interval where the condensates are held experimentally the splitting does not occur clearly in this cylindrically symmetric case shown in the lower panels of Fig. 1(b). It should be noted that the characteristic time for the splitting is proportional to the inverse of the imaginary part of the eigenfrequency, but the mode excited with $`q=2`$ has not yet split the vortex in Fig. 1(b). The above calculations have also been carried out for the anisotropies $`\lambda =1.02\text{ and }1.03`$, but no difference was observed in the splitting of the vortex. The detailed dependence of the splitting on trap anisotropy is left for future research.

In order to investigate the dissipative dynamics, we introduce a phenomenological damping term choi ; tsubota  to the GP equation through the substitution $`t(1i\gamma )t`$, $`\gamma =0.005`$. The temporal development of the angular momentum is plotted in Fig. 1(c) for several values of the interaction strength and dynamics of the double-quantum vortex with $`an_z=13.75`$ is shown in the insets. It is found that when modes with a complex eigenvalue exist, the angular momentum decreases with time, which is related to the decay into a pair of single-quantum vortices which spiral out to the condensate surface. In contrast, when there occur no complex eigenvalues, the double-quantum vortex becomes stable for a longer time.

The experimentally realized multiquantum vortices were created in a cigar-shaped condensate confined before the vortex-generation process in a harmonic trap with the frequencies $`\omega _{}=2\pi \times 210`$ Hz and $`\omega _z=2\pi \times 6`$ Hz topological . After the creation, the radial trapping frequency was $`\omega _{}=2\pi \times 250`$ Hz and the condensate with $`N=10^6`$ atoms was unconfined along the $`z`$-direction. The unconfinement is a technical difficulty which can be removed by reversing the magnetic field trapping the atoms along the $`z`$-direction when the bias field vanishes or by adding an optical potential along the $`z`$-direction so that the trapping frequency stays constant. Since it is possible to realize non-expanding condensate with a double-quantum vortex and the Bogoliubov equations can not be solved for an unconfined case, we restrict our computation to the confined case with $`\omega _{}=2\pi \times 250`$ Hz and $`\omega _z=2\pi \times 6`$ Hz. Even though the condensate expands in the experiments topological , the qualitative behavior of the splitting of the double-quantum vortex is revealed by our results.

After finding the minimum-energy solution for the condensate with a double-quantum vortex along the $`z`$-axis footnote , we plotted the strength of the interactions $`an_z(z)=a|\varphi (z,x,y)|^2𝑑x𝑑y`$ as a function of the $`z`$-coordinate in Fig. 2(a), from which it is seen that the linear densities $`n_z`$ in the regions $`|z|<51`$ $`\mu \mathrm{m}`$ and $`|z|>142`$ $`\mu \text{m}`$ correspond to the two leftmost complex eigenfrequency regions in Fig. 1(a). Therefore, it is expected that the splitting of the vortex begins in those regions and that the Bogoliubov spectrum features many complex-frequency modes, since the cigar-shaped condensate may be visualized to consist of several non-interacting condensates in $`xy`$-planes if one neglects the kinetic energy along the $`z`$-direction.

The above argument is verified by first solving the Bogoliubov equations for the cigar-shaped rotational BEC, which show that the vortex state features many complex- frequency eigenmodes with $`q=2`$, plotted in Fig. 3(a). The frequencies possess noticeable finite imaginary parts when their real part is in either one of two distinct domains which are related to the two regions in the $`z`$-coordinate where the splitting of the vortex is expected to begin. The mode marked with a rectangle in Fig. 3(a), is localized in the region $`|z|<60\text{ }\mu \text{m}`$ (see Fig. 3(b)), while the other complex modes with Im$`(\omega _\alpha )>0.02`$ occupy the region $`|z|>100\text{ }\mu \text{m}`$ (see Fig. 3(c)).

We solve the non-dissipative temporal evolution of the doubly quantized vortex state using the GP equation with the trap anisotropy $`\lambda =1.01`$. Isosurface plots of the particle density after $`18.1`$ and $`41.4\text{ ms}`$ of time propagation are shown in Fig. 2(b), from which it is observed that the doubly quantized vortex had split into two single-quantum vortices first at the ends and then in the middle of the condensate, consistent with the above discussion. Also the particle density integrated in the $`z`$-direction $`n_{}(x,y)=|\varphi (x,y,z)|^2𝑑z`$ is displayed to illustrate the density distribution observed in the experiments topological . If a dissipation term $`\gamma =0.005`$ is used, the splitting is slightly faster, but the shape of the condensate is not deformed in time as much as in the non-dissipative system, since the surface waves seen at the ends of the condensate are damped.

Simulations performed for the cylindrically symmetric system show that an off-axis vortex tends to precess in the positive direction and the frequency of precession is an increasing function of the precession radius jukka . However, numerical simulations of the GP equation show that the repulsion between the two vortices in a vortex pair renders the frequency of precession to be a decreasing function of the precession radius. Analogously to the cylindrically symmetric system, two vortices tend to twist around each other in a cigar-shaped BEC such that the parts where the distance between the vortices is large lag in precession compared with the parts where the distance is small, resulting in an entangled vortex pair 111The coupling of the Kelvin modes with dynamical instabilities may also play a role in the entanglement of vortices. The consequences of this phenomenon is left for future research..

Following the above argument, the twisting of the vortices at $`t=41.4\text{ ms}`$ is negative in the region $`75\text{ }\mu \text{m}<z<150\text{ }\mu \text{m}`$, *i.e.*, the azimuthal angle of the position of the vortices is a decreasing function of the $`z`$-coordinate. The azimuthal angle decreases approximately by $`10\pi `$ as the $`z`$-coordinate changes from $`z=78\text{ }\mu \text{m}`$ to $`130\text{ }\mu \text{m}`$.

Even though the vortex had clearly split in the middle and at the ends of the condensate, only a single vortex core is seen in the $`z`$-integrated particle density. The reasons are that the twisting of the two singly quantized vortices causes the splitting to average out in the $`z`$-integration and that the vortex remains doubly quantized in the range $`60\text{ }\mu \text{m}<z<120\text{ }\mu \text{m}`$ which covers a significant fraction of the condensed particles. This suggests an explanation of why no splitting was observed in the experiments topological .

In conclusion, the stability and the splitting of the double-quantum vortex created in Ref. topological  has been analyzed via the Bogoliubov excitation spectrum and the time-dependent GP equation. The vortex is found to be spectrally unstable and the splitting of the vortex is in close analogy with the calculations for cylindrically symmetric BECs. The emerging two single-quantum vortex lines tend to twist around each other, which makes it difficult to study the splitting of the doubly quantized vortex via particle density graphs integrated in the direction of the vortex line. Therefore, we suggest particle density profiles integrated in a direction perpendicular to the vortex line to be extracted xdens  in future experiments. We also encourage an experiment with different total particle numbers to investigate the various combinations of instability regions shown in Fig. 1(a). Our simulations show that the splitting of the doubly quantized vortex begins only from the ends of the condensate for peak interaction strength $`an_z<3`$. A confining potential should also be used as suggested to increase the lifetime of the condensate.

The authors thank CSC-Scientific Computing Ltd (Espoo, Finland) for resources and M. Nakahara, T. Ohmi, T. P. Simula, and S. M. M. Virtanen for stimulating discussions. MM acknowledges Foundation of Technology (Helsinki, Finland) for financial support.


text:
# Spherical Relativistic Hartree theory in a Woods-Saxon basis

## I Introduction

The existence of an average field in atomic nuclei revealed by the exceptional role of the nuclear magic numbers provides the foundation of the nuclear shell model and various mean field approaches Mayer55 ; Bohr69 ; Ring80 . This average field is believed to be approximated most closely by a Woods-Saxon (WS) potential Woods54  either from analyzing the radial dependence of the nuclear force or by deriving it from a microscopic two-body force.

Since the eigenfunctions for the WS potential can not be given analytically, as good approximations for stable nuclei, one often adopts the harmonic oscillator (HO) potential or the square well, in particular the former, in shell model calculations for both spherical Mayer55  and deformed nuclei Nilsson55 . The HO eigenfunctions also often serve as a complete basis in solving equations in both non-relativistic and relativistic mean field approximations, such as the Skyrme Hartree-Fock (SHF), Hartree-Fock-Bogoliubov (HFB), relativistic Hartree (RH) and relativistic Hartree-Bogoliubov (RHB) theories. In these approaches, solution of corresponding equations is transformed to a matrix diagnolization problem which can be easily dealt with.

However, due to the incorrect asymptotic property of the HO wave functions, the expansion in the localized HO basis are not appropriate for the description of drip line nuclei Dobaczewski96 ; Meng96 ; Zhou00  which display many interesting features because of the extremely weakly bound property, e.g., the coupling between bound states and the continuum due to the pairing correlation, large spacial density distributions, possible modifications of shell structure, et al. One must improve the asymptotic behavior of HO wave functions, e.g., by performing a local scaling transformation Stoitsov98 . However, one does not know the scaling parameter beforehand thus the predictive power in this method is lost.

A proper representation to solve the HFB or RHB equations for drip line nuclei is the coordinate space Dobaczewski84 ; Meng96 ; Terasaki96 ; Poeschl97  where wave functions are approximated on a spatial lattice and the continuum is discretized by suitably large box boundary conditions. The HFB method solved in $`r`$ space can take fully into account all the mean-field effects of the coupling to the continuum Dobaczewski96 ; Meng96 ; Dobaczewski84 ; Belyaev87 . Nevertheless for deformed nuclei, working in $`r`$ space becomes much more difficult and numerically very sophisticated Zhou00 . Particularly, it become very time consuming when the pairing correlation is included. Therefore much effort is made towards a more efficient solution of HFB or RHB equations, e.g., using natural orbitals Reinhard97  or working on basis-spline Galerkin lattices Oberacker99 , etc..

A reconciler between the HO basis and the $`r`$ space may be the WS basis because (i) the WS potential represents the nuclear average field more suitably than the HO potential and (ii) in principal there is no localization restrictions in the WS potential. Although analytical wave functions can not be given for the WS potential, one may easily find numerical solutions for a spherical WS potential in the $`r`$ space by virtue of various effective methods of solving ordinary differential equations NumericalRecipe . One can still use a large box boundary condition to discretize the continuum. These WS wave functions can thus be used as a complete basis for spherical or deformed systems and one finally comes back to the familiar matrix diagnolization problem.

In the present work we restrict the application of this method to nuclei with spherical symmetry which largely facilitates the discussion of basic principles and allows presenting illustrations for the radial dependence of all relevant physical quantities like density distributions. We combine this approach with the relativistic Hartree theory Serot86  which provides a framework for describing the nuclear many body problem as a relativistic system of baryons and mesons and, together with its extensions with deformation and/or pairing included, have been successfully applied in calculations of nuclear matter and properties of finite nuclei throughout the periodic table Reinhard89 ; Ring96 .

The paper is organized as follows. In Sec. II, we give a brief reminder of the formalism of relativistic Hartree theory. The numerical details of solving it in the WS basis are given in Sec. III. In Sec. IV we present our results and compare them with those obtained in the HO basis and in the $`r`$ space. We also discuss the contribution from negative levels in the Dirac sea in the same section. Finally, the work is summarized in Sec. V.

Throughout the paper, the relativistic Hartree theories solved in the $`r`$ space, in the HO basis and in the WS basis are abbreviated as “SRHR”, “SRHHO” and “SRHWS” where the first “S” represents “spherical”. We use “SWS” and “DWS” to distinguish the WS basis which is obtained from solving the Schrödinger equation or the Dirac equation with initial WS potentials, respectively. Thus we have “SRHSWS” and “SRHDWS” theories.

## II Basic formalism of relativistic Hartree theory

The starting point of the relativistic Hartree theory is a Lagrangian density where nucleons are described as Dirac spinors which interact via the exchange of several mesons ($`\sigma `$, $`\omega `$, and $`\rho `$) and the photon Serot86 ; Reinhard89 ; Ring96 ,

$``$ $`=`$ $`\overline{\psi }_i\left(i\text{/}M\right)\psi _i+{\displaystyle \frac{1}{2}}_\mu \sigma ^\mu \sigma U(\sigma )g_\sigma \overline{\psi }_i\sigma \psi _i`$ (1)
$`{\displaystyle \frac{1}{4}}\mathrm{\Omega }_{\mu \nu }\mathrm{\Omega }^{\mu \nu }+{\displaystyle \frac{1}{2}}m_\omega ^2\omega _\mu \omega ^\mu g_\omega \overline{\psi }_i\text{/}\omega \psi _i`$
$`{\displaystyle \frac{1}{4}}\stackrel{}{R}_{\mu \nu }\stackrel{}{R}^{\mu \nu }+{\displaystyle \frac{1}{2}}m_\rho ^2\stackrel{}{\rho }_\mu \stackrel{}{\rho }^\mu g_\rho \overline{\psi }_i\text{/}\stackrel{}{\rho }\stackrel{}{\tau }\psi _i`$
$`{\displaystyle \frac{1}{4}}F_{\mu \nu }F^{\mu \nu }e\overline{\psi }_i{\displaystyle \frac{1\tau _3}{2}}\text{/}A\psi _i,`$
with the summation convention used and the summation over $`i`$ runs over all nucleons, $`\text{/}x\gamma ^\mu x_\mu =\gamma _\mu x^\mu `$, $`M`$ the nucleon mass, and $`m_\sigma `$, $`g_\sigma `$, $`m_\omega `$, $`g_\omega `$, $`m_\rho `$, $`g_\rho `$ masses and coupling constants of the respective mesons. The nonlinear self-coupling for the scalar mesons is given by Boguta77

$$U(\sigma )=\frac{1}{2}m_\sigma ^2\sigma ^2+\frac{g_2}{3}\sigma ^3+\frac{g_3}{4}\sigma ^4,$$
(2)
and field tensors for the vector mesons and the photon fields are defined as

$`\{\begin{array}{ccc}\hfill \mathrm{\Omega }_{\mu \nu }& =& _\mu \omega _\nu _\nu \omega _\mu ,\hfill \\ \hfill \stackrel{}{R}_{\mu \nu }& =& _\mu \stackrel{}{\rho }_\nu _\nu \stackrel{}{\rho }_\mu g_\rho (\stackrel{}{\rho }_\mu \times \stackrel{}{\rho }_\nu ),\hfill \\ \hfill F_{\mu \nu }& =& _\mu A_\nu _\nu A_\mu .\hfill \end{array}`$ (6)
The classical variation principle gives the equations of motion for the nucleons, mesons and the photon. As in many applications, we study the ground state properties of nuclei with time reversal symmetry, thus the nucleon spinors are the eigenvectors of the stationary Dirac equation

$$\left[𝜶𝒑+V(𝒓)+\beta (M+S(𝒓))\right]\psi _i(𝒓)=ϵ_i\psi _i(𝒓),$$
(7)
and equations of motion for the mesons and the photon are

$`\{\begin{array}{ccc}\hfill \left(\mathrm{\Delta }+_\sigma U(\sigma )\right)\sigma (𝒓)& =& g_\sigma \rho _s(𝒓),\hfill \\ \hfill \left(\mathrm{\Delta }+m_\omega ^2\right)\omega ^0(𝒓)& =& g_\omega \rho _v(𝒓),\hfill \\ \hfill \left(\mathrm{\Delta }+m_\rho ^2\right)\rho ^0(𝒓)& =& g_\rho \rho _3(𝒓),\hfill \\ \hfill \mathrm{\Delta }A^0(𝒓)& =& e\rho _p(𝒓),\hfill \end{array}`$ (12)
where $`\omega ^0`$ and $`A^0`$ are time-like components of the vector $`\omega `$ and the photon fields and $`\rho ^0`$ the 3-component of the time-like component of the iso-vector vector $`\rho `$ meson. Equations (7) and (12) are coupled by the vector and scalar potentials

$`\{\begin{array}{ccc}V(𝒓)\hfill & =\hfill & g_\omega \omega ^0(𝒓)+g_\rho \tau _3\rho ^0(𝒓)+e{\displaystyle \frac{1\tau _3}{2}}A^0(𝒓),\hfill \\ S(𝒓)\hfill & =\hfill & g_\sigma \sigma (𝒓),\hfill \end{array}`$ (15)
and various densities

$`\{\begin{array}{ccc}\hfill \rho _s(𝒓)& =& _{i=1}^A\overline{\psi }_i(𝒓)\psi _i(𝒓),\hfill \\ \hfill \rho _v(𝒓)& =& _{i=1}^A\psi _i^{}(𝒓)\psi _i(𝒓),\hfill \\ \hfill \rho _3(𝒓)& =& _{i=1}^A\psi _i^{}(𝒓)\tau _3\psi _i(𝒓),\hfill \\ \hfill \rho _c(𝒓)& =& _{i=1}^A\psi _i^{}(𝒓){\displaystyle \frac{1\tau _3}{2}}\psi _i(𝒓).\hfill \end{array}`$ (20)
For spherical nuclei, meson fields and densities depend only on the radial coordinate $`r`$. The spinor is characterized by the angular momentum quantum numbers ($`l`$,$`j`$), $`m`$, the parity, the isospin $`t=\pm 1/2`$ (“+” for neutrons and “$``$” for protons), and the radial quantum number $`\alpha `$. The Dirac spinor has the form

$$\psi _{\alpha \kappa m}(𝒓,s,t)=\left(\begin{array}{c}i\frac{G_\alpha ^\kappa (r)}{r}Y_{jm}^l(\theta ,\varphi )\\ \frac{F_\alpha ^\kappa (r)}{r}Y_{jm}^{\stackrel{~}{l}}(\theta ,\varphi )\end{array}\right)\chi _{t_\alpha }(t),j=l\pm \frac{1}{2},$$
(21)
with $`G_\alpha ^\kappa (r)/r`$ and $`F_\alpha ^\kappa (r)/r`$ the radial wave functions for the upper and lower components and $`Y_{jm}^l(\theta ,\varphi )`$ the spin spherical harmonics where $`\kappa =(1)^{j+l+1/2}(j+1/2)`$ and $`\stackrel{~}{l}=l+(1)^{j+l1/2}`$. The value of $`\kappa `$ of the upper component is used to label a state both for normal levels in the Fermi sea and for negative ones in the Dirac sea. States with the same $`\kappa `$ form a “block”. The radial equation of the Dirac spinor, Eq. (7), is reduced as

$$\{\begin{array}{ccc}ϵ_\alpha G_\alpha ^\kappa \hfill & =\hfill & \left(\frac{}{r}+\frac{\kappa }{r}\right)F_\alpha ^\kappa +\left(M+S(r)+V(r)\right)G_\alpha ^\kappa ,\hfill \\ ϵ_\alpha F_\alpha ^\kappa \hfill & =\hfill & \left(+\frac{}{r}+\frac{\kappa }{r}\right)G_\alpha ^\kappa \left(M+S(r)V(r)\right)F_\alpha ^\kappa .\hfill \end{array}$$
(22)
The meson field equations become simply radial Laplace equations of the form

$$\left(\frac{^2}{r^2}\frac{2}{r}\frac{}{r}+m_\varphi ^2\right)\varphi (r)=s_\varphi (r).$$
(23)
$`m_\varphi `$ are the meson masses for $`\varphi =\sigma ,\omega ,\rho `$ and zero for the photon. The source terms are

$`s_\varphi (r)=\{\begin{array}{cc}g_\sigma \rho _s(r)g_2\sigma ^2(r)g_3\sigma ^3(r),\hfill & \text{for}\sigma ,\hfill \\ g_\omega \rho _v(r),\hfill & \text{for}\omega ,\hfill \\ g_\rho \rho _3(r),\hfill & \text{for}\rho ,\hfill \\ e\rho _c(r),\hfill & \text{for}A,\hfill \end{array}`$ (28)
with

$`\{\begin{array}{ccc}4\pi r^2\rho _s(r)\hfill & =\hfill & _{i=1}^A(|G_i(r)|^2|F_i(r)|^2),\hfill \\ 4\pi r^2\rho _v(r)\hfill & =\hfill & _{i=1}^A(|G_i(r)|^2+|F_i(r)|^2),\hfill \\ 4\pi r^2\rho _3(r)\hfill & =\hfill & _{i=1}^A2t_i(|G_i(r)|^2+|F_i(r)|^2),\hfill \\ 4\pi r^2\rho _c(r)\hfill & =\hfill & _{i=1}^A\left(\frac{1}{2}t_i\right)(|G_i(r)|^2+|F_i(r)|^2).\hfill \end{array}`$ (33)
The above coupled equations have been solved in the $`r`$ space Horowitz81  and in the HO basis Gambhir90  with the no sea and the mean field approximation. Here we depict briefly the procedure of solving these coupled equations. With a set of estimated meson and photon fields, the scalar and vector potentials are calculated and the radial Dirac equation is solved. Thus obtained nucleon wave functions are used to calculate the source term of each radial Laplace equation for mesons and the photon. New meson and photon fields are calculated from solving these Laplace equations. This procedure is iterated until a demanded accuracy is achieved. Laplace equations are usually solved using the Green’s function method Horowitz81 ; Gambhir90  though in Ref. Gambhir90  Laplace equations for mesons are solved in the HO basis. SRHR, SRHHO and SRHWS differ from each other mainly in how to solve the Dirac equation. In the following the numerical solution of the Dirac equation in the WS basis will be presented.

## III  Solving the Dirac equation in a Woods-Saxon basis and numerical details

### III.1 Woods-Saxon basis from solving a Schrödinger equation (the SWS basis)

For the Schrödinger equation with a spherical Woods-Saxon potential

$$V_{\text{WS}}(r)=\{\begin{array}{cc}\frac{V_0}{1+e^{(rR_0)/a_0}},\hfill & r<R_{\text{max}},\hfill \\ \mathrm{},\hfill & rR_{\text{max}},\hfill \end{array}$$
(34)
where $`R_{\text{max}}`$ is introduced for practical reasons to define the box boundary, the eigenfunction can be written as $`\varphi _{nlm_l}(𝒓)=R_{nl}(r)Y_{lm_l}(\theta ,\varphi )`$. Its radial Schrödinger equation is derived as

$`\left[{\displaystyle \frac{1}{2M}}\left({\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{r}}r^2{\displaystyle \frac{}{r}}{\displaystyle \frac{l(l+1)}{r^2}}\right)+V_{\text{WS}}(r)\right]R_{nl}(r)`$
$`=E_{nl}R_{nl}(r).`$ (35)
Equation (35) is solved on a discretized radial mesh with a mesh size $`\mathrm{\Delta }r`$. $`R_{\text{max}}`$ ($`\mathrm{\Delta }r`$) should be chosen larger (smaller) enough to make sure that the final results do not depend on it. The radial wave functions thus obtained form a complete basis,

$$\left\{R_{nl}(r);n=0,1,\mathrm{};l=0,1,\mathrm{},n\right\},$$
(36)
in terms of which the radial part of the upper and the lower components of the Dirac spinor in Eq. (22) are expanded respectively as

$$\{\begin{array}{c}G_\alpha ^\kappa (r)=i\underset{n=0}{\overset{n_{\text{max}}}{}}g_{\alpha n}rR_{nl}(r),\hfill \\ F_\alpha ^\kappa (r)=i\underset{\stackrel{~}{n}=0}{\overset{\stackrel{~}{n}_{\text{max}}}{}}f_{\alpha \stackrel{~}{n}}rR_{\stackrel{~}{n}\stackrel{~}{l}}(r).\hfill \end{array}$$
(37)
The radial Dirac equation, Eq. (22), is transformed into the WS basis as

$$\left(\begin{array}{cc}𝒜_{mn}& _{m\stackrel{~}{n}}\\ 𝒞_{\stackrel{~}{m}n}& 𝒟_{\stackrel{~}{m}\stackrel{~}{n}}\end{array}\right)\left(\begin{array}{c}g_{\alpha n}\\ f_{\alpha \stackrel{~}{n}}\end{array}\right)=ϵ_\alpha \left(\begin{array}{c}g_{\alpha n}\\ f_{\alpha \stackrel{~}{n}}\end{array}\right),$$
(38)
where the matrix elements are calculated as follows

$$\{\begin{array}{c}𝒜_{mn}=_0^{R_{\text{max}}}r^2𝑑rR_{ml}(r)\left(V(r)+S(r)+M\right)R_{nl}(r),\hfill \\ _{m\stackrel{~}{n}}=_0^{R_{\text{max}}}r^2𝑑rR_{ml}(r)\left(+\frac{}{r}\frac{\kappa _\alpha 1}{r}\right)R_{\stackrel{~}{n}\stackrel{~}{l}}(r),\hfill \\ 𝒞_{\stackrel{~}{m}n}=_0^{R_{\text{max}}}r^2𝑑rR_{\stackrel{~}{m}\stackrel{~}{l}}(r)\left(\frac{}{r}\frac{\kappa _\alpha +1}{r}\right)R_{nl}(r),\hfill \\ 𝒟_{\stackrel{~}{m}\stackrel{~}{n}}=_0^{R_{\text{max}}}r^2𝑑rR_{\stackrel{~}{m}\stackrel{~}{l}}(r)\left(V(r)S(r)M\right)R_{\stackrel{~}{n}\stackrel{~}{l}}(r).\hfill \end{array}$$
(39)
In practical calculations, an energy cutoff $`E_{\text{cut}}`$ (relative to the nucleon mass $`M`$) is used to determine the cutoff of the radial quantum number $`n_{\text{max}}`$ for each block. In the expansion of the corresponding lower component, we take $`\stackrel{~}{n}_{\text{max}}=n_{\text{max}}+\mathrm{\Delta }n`$ with $`\mathrm{\Delta }n1`$ in order to avoid spurious states Gambhir90 .

The following Woods-Saxon parameters have been used according to Ref. Heyde99

$$\{\begin{array}{c}V_0=(51\pm 33(NZ)/A)\text{ MeV},\hfill \\ R_0=1.27A^{1/3}\text{ fm},a_0=0.67\text{ fm},\hfill \end{array}$$
(40)
where ‘+’ is for the neutron and ‘$``$’ for the proton. As expected, the dependence of final results on the initial WS potential is almost negligible. For example, a variation of $`V_0`$ by 50% gives differences in total binding energies by less than 0.1% and charge radii by less than 0.5% for 16O, 48Ca and 208Pb. Such situation is also checked to be true for the other two parameters in the WS potential, $`R_0`$ and $`a_0`$.

### III.2 Woods-Saxon basis from solving a Dirac equation (the DWS basis)

The radial Dirac equation, Eq. (22), may be solved in the $`r`$ space Horowitz81  with Woods-Saxon-like potentials for $`V_0(r)\pm S_0(r)`$ Koepf91  within a spherical box of the size $`R_{\text{max}}`$, together with the spherical spinor which gives a complete WS basis

$$\left\{[ϵ_{n\kappa m}^0,\psi _{n\kappa m}^0(𝒓,s,t)];ϵ_{n\kappa m}^00\right\},$$
(41)
with $`n=0,1,\mathrm{}`$, $`\kappa =\pm 1,\pm 2,\mathrm{}`$, and $`m=j_\kappa ,\mathrm{},j_\kappa `$. $`\psi _{n\kappa m}^0(𝒓,s,t)`$ takes the form of Eq. (21). We note that states both in the Fermi sea and in the Dirac sea should be included in the basis for the completeness. The nucleon wave function, Eq. (21), can be expanded in terms of this set of basis as

$$\psi _{\alpha \kappa m}(𝒓,s,t)=\underset{n=0}{\overset{n_{\text{max}}}{}}c_{\alpha n}\psi _{n\kappa m}^0(𝒓,s,t),$$
(42)
where $`n_{\text{max}}=n_{\text{max}}^++n_{\text{max}}^{}+1`$ and the summation is over normal levels in the Fermi sea for $`0nn_{\text{max}}^+`$ and over negative levels in the Dirac sea for $`n_{\text{max}}^++1nn_{\text{max}}`$. The negative states is obtained with the same method as the positive ones Horowitz81 . In this WS basis, the Dirac equation, Eq. (7), turns out to be

$$c_{\alpha m}ϵ_m^0+\underset{n=0}{\overset{n_{\text{max}}}{}}c_{\alpha n}H_{mn}^{}=ϵ_\alpha c_{\alpha m},m=1,\mathrm{},n_{\text{max}},$$
(43)
with

$`H_{mn}^{}`$ $`=`$ $`\psi _m^0(𝒓)\left|\left[\mathrm{\Delta }V(𝒓)+\beta \mathrm{\Delta }S(𝒓)\right]\right|\psi _n^0(𝒓)`$ (44)
$`=`$ $`{\displaystyle _0^{R_{\text{max}}}}𝑑rG_m^0(r)\left[\mathrm{\Delta }V(r)+\mathrm{\Delta }S(r)\right]G_n^0(r)`$
$`+`$ $`{\displaystyle _0^{R_{\text{max}}}}𝑑rF_m^0(r)\left[\mathrm{\Delta }V(r)\mathrm{\Delta }S(r)\right]F_n^0(r),`$
where $`\mathrm{\Delta }V(𝒓)=V(𝒓)V_0(𝒓)`$ and $`\mathrm{\Delta }S(𝒓)=S(𝒓)S_0(𝒓)`$. The angular, spin, and isospin quantum numbers are omitted for brevity.

It should be mentioned that Eq. (22) can be solved directly in the $`r`$ space with the same method of generating the DWS basis. It is our aim to test the validity of an efficient solution not only for the spherical RH model but also for its extension to include the deformation and/or the pairing correlation. In fact, if only the SRH theory is concerned, this procedure is just a replacement of the direct solution in the $`r`$ space by a diagnolization of a matrix with some complication introduced by the fact that contributions from states in the Dirac sea must be included.

An energy cutoff $`E_{\text{cut}}`$ (relative to the nucleon mass $`M`$) and the cutoff of radial quantum numbers $`n_{\text{max}}^+`$ are applied to normal levels alternatively according to practical convenience. For the initial Woods-Saxon potentials $`V_0(r)\pm S_0(r)`$, we follow Ref. Koepf91 .

### III.3 Comparison with the $`r`$-space method

In order to check the validity of solving the Dirac equation in the WS basis and to provide numerical experiences for future applications, we compare results of 16O and 208Pb from solving the Dirac equation in the WS basis and those from solving the same equation in the coordinate space. The latter is the most accurate method of solving the Dirac equation for realistic nuclei up to now thus is used as a standard here. The scalar and vector potentials in the Dirac equation are provided by very accurate SRHR calculations with the parameter set NL3 for the Lagrangian, the mesh size $`\mathrm{\Delta }r`$ = 0.05 fm, the box size $`R_{\text{max}}`$ = 30 fm for 16O and $`R_{\text{max}}`$ = 35 fm for 208Pb. Then with thus obtained $`S(r)`$ and $`V(r)`$, the Dirac equation is solved in the coordinate space and in the WS basis with also the parameter set NL3.

To compare results only from solving the Dirac equation avoids errors from other numerical procedures, e.g., the error from the iteration and that from solving the Laplace equations. For the same reason, what we compare between these two methods is not the binding energy which contains the contribution from mesons but the average single particle energy $`E^{\text{sp}}/A`$. $`E^{\text{sp}}=_iϵ_i`$ where $`ϵ_i`$ is the single particle energy and the summation over all occupied states for both neutrons and protons. We also compare the rms radius $`r^2^{1/2}`$ and $`r^4^{1/4}`$. The radius $`r^4^{1/4}`$ reflects the nucleon densities in the tail more than the rms radius.

Table 1 presents the dependence of results of the Dirac equation in the SWS basis on the mesh size $`\mathrm{\Delta }r`$. With $`\mathrm{\Delta }r`$ decreasing, results in the SWS basis approach the standard results, i.e., those in the $`r`$ space. $`\mathrm{\Delta }r`$ = 0.1 fm gives results accurate enough. The dependence on the box size $`R_{\text{max}}`$ and on the basis size determined by $`E_{\text{cut}}`$ are investigated and shown in Figs. 1 and 2 where the deviations of the average single particle energy $`E^{\text{sp}}/A`$, the rms radius $`r^2^{1/2}`$ and $`r^4^{1/4}`$ from the standards are plotted versus $`E_{\text{cut}}`$ for different $`R_{\text{max}}`$. If $`R_{\text{max}}`$ is not large enough, it is difficult to approach the standard results. For example, when $`R_{\text{max}}=3r_0A^{1/3}=9.4`$ fm and $`E_{\text{cut}}`$ = 300 MeV for 16O (correspondingly, $`N_{\text{max}}20`$), the results seem converge, but the discrepancy of the average single particle energy from the standard one remains 0.1 keV. So one must use a large enough box with the size $`R_{\text{max}}`$ around $`4r_0A^{1/3}`$ for light nuclei and $`3r_0A^{1/3}`$ for heavy ones. It is interesting that the convergence of the results does not depend on $`N_{\text{max}}`$ but only on $`E_{\text{cut}}`$. For 16O (208Pb), the results converge to the standard ones at $``$ 300 (400) MeV. From Figs. 1 and 2, we find the radius $`r^4^{1/4}`$ also converges very well which implies that nucleon densities could be calculated accurately even for large $`r`$.

We have made similar investigations for results in the DWS basis and similar conclusions are made. For instances, the deviations of the average single particle energy $`E^{\text{sp}}/A`$, the rms radius $`r^2^{1/2}`$ and $`r^4^{1/4}`$ from the standards are plotted versus $`E_{\text{cut}}`$ for different $`R_{\text{max}}`$ in Figs. 3 and 4.

In the expansion of the nucleon wave function, Eq. (42), one has to take into account not only the levels in the Fermi sea but also those in the Dirac sea because they form a complete basis together. Now the question arises how many levels in the Dirac sea one has to take into account. In the calculations in Figs. 3 and 4, we have used $`n_{\text{max}}^{}=n_{\text{max}}^+`$ with $`n_{\text{max}}^+`$ determined by $`E_{\text{cut}}`$. In Table 3, the dependence of the average single particle energy, the rms radius $`r^2^{1/2}`$ and $`r^4^{1/4}`$ on $`N_{\text{max}}^{}=2n_{\text{max}}^{}+l`$ — a cutoff on the principal quantum number of levels in the Dirac sea — are given for 16O and 208Pb. From Table 3, we find a merit of solving the Dirac equation in the DWS basis: the number of negative states included in the basis could be much smaller that that of the positive states. Let take 16O as an example, $`R_{\text{max}}=4r_0A^{1/3}`$ and $`E_{\text{cut}}`$ = 300 MeV for positive states correspond to $`N_{\text{max}}^+`$ 28. But for negative states, $`N_{\text{max}}^{}`$ = 10 gives very accurate results, e.g., the discrepancy of $`E^{\text{sp}}/A`$ from the standard is smaller than 0.1 keV. This will significantly simplify the deformed problem by decreasing the matrix dimension compared to solve the Dirac equation in the SWS basis.

The above investigations are somehow academic. In practical applications, it is not necessary to go to the accuracy around keV in the single particle energy or 10-4 fm in the radius. So in the following calculations, we will use $`R_{\text{max}}`$ = 20 fm, $`\mathrm{\Delta }r`$ = 0.1 fm and $`E_{\text{cut}}`$ = 60$``$80 MeV for heavy and light nuclei which give reasonable accuracies both for the binding energy and the radius. This set of cutoff’s corresponds approximately $`N_{\text{max}}=2n_{\text{max}}+l25`$ where $`l`$ is the orbital angular momentum of relevant state.

## IV Results and discussions

In this section we present results of SRHWS. Since our main aim is to show the virtues of SRHWS compared to SRHHO and SRHR, we do not include pairing correlations and restrict our study to doubly magic or magic nuclei only. If not specified, the parameter set NLSH is used for the Lagrangian, $`R_{\text{max}}`$ = 20 fm and $`\mathrm{\Delta }r`$ = 0.1 fm throughout this section. Other parameter sets for the Lagrangian do not change the conclusion here. In SRHDWS, the number of normal levels in the Fermi sea and that of negative ones in the Dirac sea are the same for convenience, i.e., $`n_{\text{max}}^+=n_{\text{max}}^{}`$. For SRHHO, $`\mathrm{}\omega _0=41A^{1/3}`$ has been used and cutoff’s of the expansion for fermions and bosons are the same, i.e., $`N_\text{F}=N_\text{B}N_{\text{max}}`$.

### IV.1 Bulk properties of stable nuclei from different SRH theories

In Table 4, the binding energy per nucleon ($`E/A`$) and neutron, proton and charge radii ($`r_\text{n}`$, $`r_\text{p}`$ and $`r_\text{c}`$) of some typical spherical nuclei are presented which are calculated from the present available codes, including SRHR, SRHSWS, SRHDWS, SRHHO. Available data Audi95 ; DeVries87  are also included for comparison. We use approximately the same $`N_{\text{max}}`$ in the SRHHO as that in the SRHWS which is determined by $`E_{\text{cut}}`$.

Generally speaking, for each studied nucleus, the four approaches give almost the same results with an accuracy within 0.1% with few exceptions where the differences are still less than 0.3%. They are in excellent agreement with available data.

With the same parameters of spatial lattice $`R_{\text{max}}`$ and $`\mathrm{\Delta }r`$, SRHWS should reproduce results of SRHR when $`E_{\text{cut}}`$ (or $`N_{\text{max}}`$) is large enough. This is justified in Table 4. One find exactly coincident results between SRHSWS and SRHR for most of the studied nuclei. The remaining differences and those between SRHDWS and SRHR could be diminished by increasing $`E_{\text{cut}}`$.

In Figures 5, 6 and 7, the neutron density distributions are compared between SRHR, SRHSWS, SRHDWS, and SRHHO, in which 16O, 48Ca and 208Pb are chosen as examples. The calculation details are the same as Table 4. For these stable nuclei, all these SRH methods are valid and all calculations are in excellent agreement with each other from the central to outer region of each nucleus. Tiny differences in the central region do not affect the physical observables like the binding energy or nuclear radius as is seen in Table 4. Furthermore, these differences could also be decreased by increasing $`E_{\text{cut}}`$ or $`N_{\text{max}}`$.

From the above discussions, it is clear that SRHWS is equivalent to SRHR and SRHHO for stable nuclei. Thus we conclude that Woods-Saxon basis provide another possibility to solve (non-)relativistic mean field theory.

### IV.2 Neutron density distributions for 72Ca in different SRH theories

As we already discussed in the introduction, one of the merits of SRHR against SRHHO is its proper description of exotic nuclei. In this subsection, we will demonstrate the equivalence between SRHWS and SRHR when reasonably large $`E_{\text{cut}}`$ is applied in SRHWS.

In order to see the results for the unstable nuclei near the neutron drip line, the neutron density distribution for 72Ca is studied here. The nucleus 72Ca is predicted to be the last bound calcium isotope Im00 ; Hamamoto01 ; Zhang02 ; Meng02 . Since it is not a doubly magic nucleus, there might be some uncertainty in present results due to the lack of inclusion of pairing correlations. However, as we stressed in the beginning of this section, the main aim here is to show the virtue of SRHWS compared to SRHHO, it is very unlikely that the pairing would change our conclusion qualitatively.

For stable nuclei, it has been shown that $`R_{\text{max}}`$ 20 fm is large enough. For drip line nuclei, the dependence of the results on $`R_{\text{max}}`$ for 72Ca is presented in Table 5. For both SRHR and SRHWS, $`\mathrm{\Delta }r`$ = 0.1 fm and $`R_{\text{max}}`$ = 20, 25, 30 and 35 fm have been adopted respectively. The energy cutoff $`E_{\text{cut}}`$ = 75 MeV is used to SRHWS. In SRHR and SRHWS calculations, the neutron rms radius $`r_\text{n}`$ and Fermi energy $`\lambda _\text{n}`$ of 72Ca converge around $`R_{\text{max}}`$ = 35 fm while a independence of the binding energy per nucleon $`E/A`$ and proton rms radius $`r_\text{p}`$ on the box size can be seen. These sets of parameters in SRHWS, $`E_{\text{cut}}`$ = 75 MeV and $`R_{\text{max}}`$ = 20, 25, 30 and 35 fm, correspond to cutoff’s on principal quantum number $`N_{\text{max}}`$ = 25, 31, 37 and 43 which are used in SRHHO calculations in order to make fair comparisons between SRHWS and SRHHO. Similar as those from SRHR and SRHWS, $`E/A`$ and $`r_\text{p}`$ depend little on $`N_{\text{max}}`$ in SRHHO. However, the neutron rms radius $`r_\text{n}`$ increases with $`N_{\text{max}}`$ steadily, which shows a much slower convergence. As it is based on a complete basis, SRHHO can reach convergence of $`r_\text{n}`$ if $`N_{\text{max}}`$ is large enough. From Table 5, one finds that for the same $`N_{\text{max}}`$ (or equivalent $`R_{\text{max}}`$), a difference of $`\mathrm{\Delta }r_\text{n}`$ 0.2 fm between SRHHO and SRHWS (SRHR) can be seen. From the slow convergence of $`r_\text{n}`$ with $`N_{\text{max}}`$ in SRHHO ($`\mathrm{\Delta }N_{\text{max}}`$ = 6 gives $`\mathrm{\Delta }r_\text{n}`$ 0.02 fm), we can estimate the lower limit of the $`N_{\text{max}}`$ as $`N_{\text{max}}`$ 90 in order to give $`r_\text{n}`$ = 4.8 fm in SRHR or SRHWS.

We compare the neutron density distribution of 72Ca from different SRH approaches in Figs. 8 and 9. With the same box size, the density distribution from SRHR are almost identical with those from SRHWS, which indicates the equivalence between SRHWS and SRHR. For brevity, only $`\rho _\text{n}(r)`$ from SRHR with $`R_{\text{max}}`$ = 35 fm is displayed in Figs. 8 or 9 which covers the curve corresponding to $`\rho _n(r)`$ from SRHSWS with $`R_{\text{max}}`$ = 35 fm in Fig. 8. On the other hand, $`\rho _\text{n}(r)`$ from SRHHO even with $`N_{\text{max}}`$ = 43 fails to reproduce the result of SRHR due to the localization property of HO potential. In addition, with the same $`N_{\text{max}}`$, the spatial extension of $`\rho _\text{n}(r)`$ from SRHWS are always larger than that from SRHHO. The variational tendency of the curve $`\rho _\text{n}(r)r`$ also explains different convergence behaviors of $`r_\text{n}`$ in SRHWS and SRHHO as given in Table 5. With increasing $`R_{\text{max}}`$, the $`\rho _\text{n}(r)`$ in SRHWS has the correct asymptotic behavior while that in SRHHO decay too quickly.

This result is very encouraging and tell us that even the long tail (or halo) behavior in neutron density distribution for nuclei near the drip line can be described by SRHSWS as well as SRHR, if pairing correlations is incorporated properly.

### IV.3 Contribution from negative levels in the SRHDWS theory

In the expansion of the nucleon wave function, Eq. (42), one has to take into account not only the levels in the Fermi sea but also those in the Dirac sea because they form a complete basis together. We study the effects of negative levels and give results of an example, 16O, in Table 6. Firstly, without negative levels included, the nucleus is over bound and the nuclear size is smaller as seen from Table 6. Secondly, in contrary to the case with negative levels included, the calculated nuclear properties depend on the initial potentials very much if no negative levels included.

It should be noted that the contribution from negative levels depends on the initial Woods-Saxon potentials for generating the DWS basis. So do the cutoff $`N_{\text{max}}`$ or $`E_{\text{cut}}`$ for convergence. If the initial Woods-Saxon potential is exactly identical to the converged potentials, the matrix in Eq. (43) is diagonal, negative states do not contribute because of the no sea approximation. Positive states can also be chosen as less as possible, e.g., $`1s_{1/2}`$, $`1p_{3/2}`$ and $`1p_{1/2}`$ are enough for 16O. From the third column corresponding to $`V_0`$ = 72 MeV in Table 6, one finds that the initial nuclear potential for the Dirac equation proposed in Ref. Koepf91  is a good choice for SRHDWS as the negative states only contribute $``$ 1.25 % to both $`E/A`$ and $`r_{\text{rms}}`$. If we change the initial potentials, e.g., by changing $`V_0`$ by 25%, much larger contributions from negative states are found in Table 6.

In order to know the contribution of negative levels in the Dirac sea to the wave function, the value of $`_n|c_n^{}|^2`$ in the expansion, Eq. (42) has been given in Table 7 for occupied states of 16O. We note that the nucleon wave function is normalized to one. It can be seen that a small component from negative states in the wave functions, about $`10^{45}`$, contributes to the physical observables such as $`E/A`$ and $`r_{\text{rms}}`$ by the magnitude of 1%$``$10% as seen from Table 6. Again we notice that the initial Woods-Saxon potentials differ more from the converged ones, the larger is the contribution from negative levels.

## V Summary

We have solved the spherical relativistic Hartree theory in the Woods-Saxon basis (SRHWS). The Woods-Saxon basis is obtained by solving either the Schrödinger equation (SRHSWS) or the Dirac equation (SRHDWS). Formalism and numerical details for both cases are presented. The WS basis in the SRHDWS theory could be very smaller than that in the SRHSWS theory. This will largely facilitate solving the deformed problem.

The results from SRHWS are compared with those from solving the spherical relativistic Hartree theory in the harmonic oscillator basis, SRHHO, and those in the coordinate space, SRHR. For stable nuclei, all approaches give identical results for properties such as total binding energies and the neutron, proton and charge rms radii as well as neutron density distributions.

For neutron drip line nuclei, e.g. 72Ca, which has a very small neutron Fermi energy $`\lambda _\text{n}0.2`$ MeV, both SRHR and SRHWS easily approach convergence by increasing box size while SRHHO does not. Furthermore, SRHWS can satisfactorily reproduce the neutron density distribution from SRHR, but SRHHO fails with similar cutoff’s.

In SRHDWS calculations, the negative levels in the Dirac Sea must be included in the basis in terms of which nucleon wave functions are expanded. We studied in detail the effects and contributions of negative states. Without negative levels included, the calculated nuclear properties depend on the initial potentials very much. A small component from negative states in the wave functions, about $`10^{45}`$, contributes to the physical observables such as $`E/A`$ and $`r_{\text{rms}}`$ by the magnitude of 1%$``$10%. When the initial potentials differ more from the converged ones, the contribution from negative levels becomes more important.

We conclude that the Woods-Saxon basis provides a reconciler between the harmonic oscillator basis and the coordinate space which may be used to describe exotic nuclei both properly and efficiently.

The extension of spherical relativistic Hartree theory in the Woods-Saxon basis to deformed cases with pairing included is in progress.

###### Acknowledgements.

S.-G. Z. would like to thank the Max-Planck-Institut für Kernphysik for kind hospitality where part of this work is done. J. M. would like to thank Physikdepartment, Technische Universität München for kind hospitality. This work is partly supported by the Major State Basic Research Development Program Under Contract Number G2000077407 and the National Natural Science Foundation of China under Grant No. 10025522, 10047001 and 19935030.


text:
# The physics of pulses in gamma-ray bursts: emission processes, temporal profiles and time lags.

## 1 Introduction

Cosmic gamma-ray burst (hereafter GRBs) exhibit a great diversity of duration and profiles. The distribution of durations is clearly bimodal with two peaks at about 0.2 and 20 seconds. GRB light curves are highly variable but can often be interpreted in terms of a succession of elementary pulses which possibly overlap (Norris et al. 1996). These pulses appear as the building blocks of the profiles and understanding their physical origin would certainly represent a clue for a better description of the whole GRB phenomenon. The pulse temporal evolution has often been described by a fast rise followed by an exponential decay (the so-called FRED shape; see Fishman et al. (1994)) but other mathematical behaviors such as stretched exponentials, gaussian (Norris et al. 1996) or power-law decays (Ryde & Svensson 2000) have been also proposed. Spectral hardness decreases during pulse decay and two relations between the temporal and spectral properties, the HIC (hardness-intensity correlation; Golenetskii et al. (1983)) and the HFC (hardness-fluence correlation; Liang & Kargatis (1996)) appear to be satisfied by a substantial fraction of GRB pulses during the decay phase. Pulse profiles peak earlier in higher energy bands and the corresponding time lags between different energy channels correlate to pulse hardness and peak luminosity (Norris et al. 2000). These observational results must be reproduced by the models and may help to discriminate among different possibilities.

Two distinct mechanisms have been proposed to explain the origin of pulses in GRBs. In the external shock model they are formed when a relativistic shell ejected by the central engine is decelerated by the circumstellar material (Meszaros & Rees 1993). An homogeneous medium leads to a single pulse but an irregular, clumpy environment can produce a complex profile if a large number of small clouds are present (Dermer & Mitman 1999). In the internal shock model (Rees & Meszaros 1994) the central engine generates a relativistic flow with a highly non uniform distribution of the Lorentz factor and the pulses are made by collisions between rapid and slower parts of the flow. In the two scenarios the variability of the profiles has a very different interpretation. In one case it provides a “tomography” of the burst environment while in the second it reveals the activity of the central engine.

In this paper we consider in some details the mechanism of pulse formation by internal shocks. Three characteristic time scales may be relevant during pulse evolution: the time $`t_{\mathrm{rad}}`$ required to radiate the energy dissipated in shocks; the dynamical time $`t_{\mathrm{dyn}}`$, i.e. the time taken by internal shocks to travel throughout the flow and the angular spreading time $`t_{\mathrm{ang}}`$ corresponding to the delay in arrival time of photons emitted from a spherical shell. A short radiative time $`t_{\mathrm{rad}}t_{\mathrm{dyn}},t_{\mathrm{ang}}`$ appears to be mandatory to avoid adiabatic losses and maintain a sufficient efficiency. This condition is satisfied by the synchrotron process which is the most commonly invoked radiation mechanism in GRBs. If the thickness of colliding shells is small compared to their initial separation, $`t_{\mathrm{dyn}}t_{\mathrm{ang}}`$ and the pulse temporal evolution is fixed by geometry; conversely if the source produces a continuous wind rather than a series of discrete, well separated shells, $`t_{\mathrm{dyn}}t_{\mathrm{ang}}`$ and hydrodynamical effects control the pulse shape.

Pulse evolution has been studied extensively when it is dominated by geometry (see e.g. Fenimore et al. (1996); Kobayashi et al. (1997)) but discrepancies between model predictions and the observations (Soderberg & Fenimore 2001) have cast some doubt about the validity of the internal shock model. Our purpose is to see if the situation can be improved when the hydrodynamical point of view is adopted. We first summarize in Sect. 2 some basic informations regarding pulse temporal and spectral evolution. We then develop in Sect. 3 a simple model where pulses are formed when a fast moving wind is decelerated by a comparatively slower shell. Spectral evolution is considered in Sect. 4 where constraints are obtained on the GRB radiation mechanism. Temporal profiles computed with our model are presented in Sect. 5 and time lags are discussed in Sect. 6. Sect. 7 is the conclusion.

## 2 Temporal and spectral evolution during pulse decay

We consider a pulse characterized by a photon flux $`N(t)`$ in the energy range ($`E_1`$, $`E_2`$), a peak energy $`E_\mathrm{p}(t)`$ of the $`E^2𝒩(E,t)`$ spectrum and a photon fluence defined by

$$\frac{\mathrm{d}\mathrm{\Phi }_\mathrm{N}(t)}{\mathrm{d}t}=N(t)=_{E_1}^{E_2}𝒩(E,t)dE.$$
(1)
The HIC and the HFC are then given by

$$E_\mathrm{p}(t)N(t)^\delta $$
(2)
and

$$E_\mathrm{p}(t)e^{a\mathrm{\Phi }_\mathrm{N}(t)},$$
(3)
where $`a`$ is an exponential decay constant. For pulses satisfying both the HIC and the HFC, Ryde & Svensson (2000) have shown that the photon flux and the peak energy follow simple power laws during the decay phase

$$N(t)=\frac{N_0}{1+t/\tau }$$
(4)
and

$$E_\mathrm{p}(t)=\frac{E_{\mathrm{p},0}}{(1+t/\tau )^\delta },$$
(5)
where $`t=0`$ corresponds to the maximum of $`N(t)`$.

Ryde & Svensson (2002) performed a detailed analysis of the decay behavior of a sample of 25 long and bright pulses to check whether it was indeed described by eqs 4 and 5. They found that to account for the temporal and spectral evolution of all the pulses, eq. 4 had to be replaced by the more general expression

$$N(t)=\frac{N_0}{(1+t/\tau )^n}.$$
(6)
If $`n1`$, the $`\delta `$ indices appearing in eqs 2 and 6 are different and following Ryde & Svensson (2002) we then write $`E_\mathrm{p}(t)=E_{\mathrm{p},0}/(1+t/\tau )^\delta _{}`$ with $`\delta _{}=n\delta `$. Ryde & Svensson (2002) found that the distribution of $`n`$ in their sample was sharply peaked at $`n=1`$ with however a secondary bump at $`n3`$. The values of $`\delta _{}`$ were all smaller than 1.5 for the $`n=1`$ pulses but could reach 3.5 when $`n3`$ . The distribution of $`\delta `$ was narrower with $`0.5\delta 1`$ in most of the sample.

Once the decay behavior of $`N(t)`$ and $`E_\mathrm{p}(t)`$ has been specified, it becomes possible to obtain the evolution of the bolometric energy flux $`F_\mathrm{E}(t)`$ since

$$F_\mathrm{E}(t)=_0^{\mathrm{}}𝒩(E,t)EdE=E_\mathrm{p}^2_0^{\mathrm{}}𝒩(x,t)xdx,$$
(7)
where $`x=E/E_\mathrm{p}`$. We suppose that the temporal and spectral behavior can be separated in $`𝒩(x,t)`$ :

$$𝒩(x,t)=A(t)(x),$$
(8)
$`(x)`$ representing the spectrum shape. The photon flux in the energy range $`(E_1,E_2)`$ is then given by

$$N(t)=_{E_1}^{E_2}𝒩(E,t)dE=\frac{F_\mathrm{E}(t)\phi (E_\mathrm{p})}{E_\mathrm{p}(t)\phi _0},$$
(9)
so that

$$F_\mathrm{E}(t)=N(t)E_\mathrm{p}(t)\frac{\phi _0}{\phi (E_\mathrm{p})}=\frac{N_0E_{\mathrm{p},0}}{(1+t/\tau )^{1+\delta _{}}}\frac{\phi _0}{\phi (E_\mathrm{p})}$$
(10)
with

$$\phi (E_\mathrm{p})=_{E_1/E_\mathrm{p}}^{E_2/E_\mathrm{p}}(x)dx$$
(11)
and

$$\phi _0=_0^{\mathrm{}}x(x)dx.$$
(12)
The derivative of $`\phi (E_\mathrm{p})`$ has been represented in Fig. 1 for the BATSE spectral range (20, 1000 keV), using a standard Band function (Band et al. 1993) with low and high energy indices $`\alpha =2/3`$, $`1`$ or $`1.5`$ and $`\beta =2.5`$. At low (resp. high) $`E_\mathrm{p}`$, $`\phi (E_\mathrm{p})`$ is given by a simple power-law $`E_\mathrm{p}^{(\beta +1)}`$ (resp. $`E_\mathrm{p}^{(\alpha +1)}`$) but for intermediate values ($`E_1<E_\mathrm{p}<E_2`$) which are representative of the decay phase in the Ryde & Svensson (2002) sample, $`\phi (E_\mathrm{p})`$ does not have a simple analytical form. Assuming that it can still be approximated by a power-law, $`\phi (E_\mathrm{p})E_\mathrm{p}^{(\zeta +1)}`$ where $`\zeta `$ is a weighted average of the low and high energy spectral indices ($`2\zeta 1`$) the bolometric energy flux also follows a power-law

$$F_\mathrm{E}(t)\frac{1}{(1+t/\tau )^ϵ}$$
(13)
with

$$ϵ=n+(2+\zeta )\delta _{}.$$
(14)
The slope of the HIC is then given by

$$\delta =\frac{\delta _{}}{n}=\frac{1}{ϵ/\delta _{}(2+\zeta )}.$$
(15)
If the temporal and spectral evolution during pulse decay is due to geometrical effects alone $`ϵ=3`$ and $`\delta _{}=1`$ (Granot et al. 1999) which leads to $`\delta =\frac{1}{1\zeta }`$. With $`2\zeta 1`$ the resulting value $`0.3\delta 0.5`$ lies below what is found in most observed pulses (Soderberg & Fenimore 2001).

Geometrical effects govern pulse evolution if the shell thickness is small compared to their initial separation. But if a continuous outflow emerges from the central engine the hydrodynamical time scale can play a dominant role during pulse decay. We have then developed a simple model to check whether a better agreement can be found with the observations when the hydrodynamical aspect of the flow is taken into account.

## 3 A simple pulse model

We consider a relativistic wind where a slow shell of mass $`M_0`$ and Lorentz factor $`\mathrm{\Gamma }_0`$ decelerates a more rapid part of the flow characterized by a constant mass flux $`\dot{M}`$ (in the source frame) and Lorentz factor $`\mathrm{\Gamma }_1>\mathrm{\Gamma }_0`$. We do not solve the true hydrodynamical problem but rather approximate the flow evolution by considering that fast material is “accreted” by the slow shell. The accretion rate is given by

$$\frac{\mathrm{d}M}{\mathrm{d}t}=\dot{M}(1\gamma ^2),$$
(16)
where $`t`$ is the observer time and $`\gamma =\mathrm{\Gamma }/\mathrm{\Gamma }_1`$ ($`\mathrm{\Gamma }`$ and $`M`$ being the current Lorentz factor and mass of the slow shell). Due to the accretion of fast moving material, the Lorentz factor of the slow shell increases. When a mass element $`\mathrm{d}M`$ is accreted the Lorentz factor becomes

$$\mathrm{\Gamma }+\mathrm{d}\mathrm{\Gamma }=\left(\mathrm{\Gamma }_1\mathrm{\Gamma }\frac{\mathrm{\Gamma }_1\mathrm{d}M+\mathrm{\Gamma }M}{\mathrm{\Gamma }\mathrm{d}M+\mathrm{\Gamma }_1M}\right)^{1/2},$$
(17)
so that

$$\frac{\mathrm{d}\gamma }{\mathrm{d}M}=\frac{1\gamma ^2}{2M},$$
(18)
which can be integrated to give

$$\mu =\left(\frac{1+\gamma }{1\gamma }\right)/\left(\frac{1+\gamma _0}{1\gamma _0}\right),$$
(19)
where $`\mu =M/M_0`$ and $`\gamma _0=\mathrm{\Gamma }_0/\mathrm{\Gamma }_1`$. Introducing $`t_0=M_0/\dot{M}`$ and $`\tau =t/t_0`$, eqs 1619 yield

$$\frac{\mathrm{d}\gamma }{\mathrm{d}\tau }=Q(1\gamma ^2)(1\gamma )^2$$
(20)
with $`Q=\frac{1}{2}\left(\frac{1+\gamma _0}{1\gamma _0}\right)`$. Equation 20 has the analytical solution

$$\tau =\frac{1}{Q}\left[F(\gamma )F(\gamma _0)\right],$$
(21)
where the function $`F(\gamma )`$ is given by

$$F(\gamma )=\frac{1}{8}\mathrm{log}\left(\frac{1+\gamma }{1\gamma }\right)+\frac{1}{4(1\gamma )}+\frac{1}{4(1\gamma )^2}.$$
(22)
The solution $`\gamma (\tau )`$ corresponding to eq. 21 has been represented in Fig. 2 for $`\gamma _0=0.25`$. When $`\tau 2`$, it is well approximated by

$$\gamma (\tau )1\frac{1}{2\sqrt{Q\tau }}.$$
(23)
Once $`\gamma (\tau )`$ is known it is possible to calculate the dissipated power

$$\dot{}(\tau )=\frac{\dot{M}\mathrm{\Gamma }_1c^2}{2}(1\gamma ^2)(1\gamma )^2,$$
(24)
which has been represented in Fig. 3. At large $`\tau `$, it behaves as $`\tau ^{3/2}`$ since

$$\dot{}(\tau )(1\gamma )^3(1+\gamma )\tau ^{3/2}$$
(25)
for $`\tau 2`$. Pulse evolution is essentially completed at $`\tau 10`$ when $`\mathrm{\Gamma }/\mathrm{\Gamma }_1>0.8`$ and $`\dot{}`$ has decreased by more than an order of magnitude.

The dissipated power given by eq. 24 is slightly different from what the observer will see since the energy released at time $`t`$ is spread over an interval $`\mathrm{\Delta }t`$ corresponding to the difference in arrival time for photons emitted by a shell of radius $`r`$ moving at a Lorentz factor $`\mathrm{\Gamma }`$

$$\mathrm{\Delta }t=\frac{r}{2c\mathrm{\Gamma }^2}.$$
(26)
The solution for $`\dot{}`$ including angular spreading has been obtained numerically and is also shown in Fig. 3. It differs from the analytical expression (eq. 24) at early times but preserves the power law decay of slope $`ϵ=3/2`$ at late times.

## 4 Spectral evolution and emission processes

We now use the analytical model to follow the spectral evolution during pulse decay. If the dissipated energy is radiated by the synchrotron process the peak energy $`E_\mathrm{p}`$ is

$$E_\mathrm{p}=E_{\mathrm{syn}}\mathrm{\Gamma }B\mathrm{\Gamma }_\mathrm{e}^2,$$
(27)
where $`B`$ is the magnetic field and $`\mathrm{\Gamma }_\mathrm{e}`$ the characteristic electron Lorentz factor behind the shock. With classical equipartition assumptions $`B`$ and $`\mathrm{\Gamma }_\mathrm{e}`$ can be expressed as

$$B=(8\pi \alpha _\mathrm{B}\rho ϵc^2)^{1/2}$$
(28)
and

$$\mathrm{\Gamma }_\mathrm{e}=\frac{\alpha _\mathrm{e}}{\zeta }\frac{m_\mathrm{p}}{m_\mathrm{e}}ϵ,$$
(29)
where $`\rho `$ is the density and $`ϵc^2`$ the dissipated energy per unit mass (both in the comoving frame); $`\alpha _\mathrm{B}`$ and $`\alpha _\mathrm{e}`$ are the equipartition parameters and $`\zeta `$ is the fraction of electrons which are accelerated. Finally,

$$E_{\mathrm{syn}}\mathrm{\Gamma }\rho ^{1/2}ϵ^{5/2},$$
(30)
where the comoving density $`\rho `$ is proportional to $`r^2`$ ($`r`$ being the shock radius $`r\mathrm{\Gamma }^2ct`$) and $`ϵ`$ is obtained from $`\dot{}=\frac{dM}{dt}\mathrm{\Gamma }ϵc^2`$ and eq. 16

$$ϵ=\frac{(1\gamma )^2}{2\gamma }.$$
(31)
This leads to the following expression for $`E_{\mathrm{syn}}`$

$$E_{\mathrm{syn}}\frac{(1\gamma )^5}{\gamma ^{7/2}t},$$
(32)
which behaves as a power law ($`E_\mathrm{p}t^{7/2}`$) when $`(1\gamma )t^{1/2}`$. This is much steeper than the observed spectral evolution of pulses which satisfy both the HIC and the HFC. Instead of using eq. 30 we therefore parametrize the peak energy with the more general phenomenological expression

$$E_\mathrm{p}\mathrm{\Gamma }\rho ^xϵ^y\frac{(1\gamma )^{2y}}{\gamma ^{4x+y1}t^{2x}},$$
(33)
which becomes

$$E_\mathrm{p}\frac{1}{t^{2x+y}}$$
(34)
at late times. The exponents $`x`$ and $`y`$ can be different from their standard synchrotron values $`1/2`$ and $`5/2`$ if the equipartition parameters $`\alpha _\mathrm{B}`$, $`\alpha _\mathrm{e}`$ or $`\zeta `$ vary with $`\rho `$ or/and $`ϵ`$. For example, Daigne & Mochkovitch (1998) adopted a fraction $`\zeta `$ of accelerated electrons proportional to $`ϵ`$ so that $`\mathrm{\Gamma }_\mathrm{e}`$ remains constant which leads to $`x=y=1/2`$ and $`E_\mathrm{p}t^{3/2}`$. However, since most of the observed values of $`\delta _{}=2x+y`$ are smaller than 1.5, it seems necessary to further reduce the $`x`$ and $`y`$ indices and we have therefore considered below the case $`x=y=1/4`$, i.e. $`\delta _{}=0.75`$.

## 5 Temporal profiles

We obtain the temporal profile of synthetic pulses from eqs 9, 24 and 33 of our model. We have represented in Fig. 4 a pulse formed when a wind of Lorentz factor $`\mathrm{\Gamma }_1=400`$ and power $`\dot{M}\mathrm{\Gamma }_1c^2=10^{52}`$ erg.s-1 is decelerated by a slow shell with $`\mathrm{\Gamma }_0=100`$. We adopt $`t_0=0.4`$ s, $`x=y=1/4`$ and $`z=1`$. The profile is computed in the BATSE range (20 – 1000 keV) and the constant of proportionality in eq. 33 is fixed to get a peak energy $`E_\mathrm{p}=300`$ keV for the whole pulse spectrum. The pulse duration is close to $`10(1+z)t_0`$ as expected from the results obtained in Sect.3.

The evolution after maximum is initially close to a $`1/t`$ decay (i.e. $`n1`$ in eq. 6) as can be seen in Fig. 4 where $`1/N(t)`$ has been also represented. This can be simply understood from eq. 14 which, for the decay slopes of the dissipated power $`ϵ=1.5`$ (eq. 25) and of the peak energy $`\delta _{}=2x+y=0.75`$ gives

$$n=0.75\zeta .$$
(35)
With $`2\zeta 1`$, the central value of $`n`$ is indeed close to unity. Since the decay phase of our synthetic pulse can be decribed by eqs. 4 and 5 it should also satisfy both the HIC and the HFC. This is checked in Fig. 5 where the two relations have been plotted in the time interval delimited by the two vertical lines in Fig. 4. The HIC is not a strict power law but its average slope $`\delta 0.9`$ is in global agreement with the observations. The HFC is satisfied to a better accuracy since the relation between $`\mathrm{log}E_\mathrm{p}`$ and the photon fluence is quasi-linear in the considered interval.

## 6 Time lags

Norris et al. (2000) have shown that time lags between different energy channels correlate with spectral hardness and possibly also with the burst peak luminosity. GRBs are distributed in a triangular domain of the time lag–hardness ratio diagram (the hardest bursts having the shortest time lags; see Fig. 6) and the time lag–luminosity relation obtained from 6 bursts with known redshifts takes the form

$$L_{51}130\left(\frac{\mathrm{\Delta }t_{31}}{0.01\mathrm{s}}\right)^{1.14},$$
(36)
where $`\mathrm{\Delta }t_{31}`$ is the time lag between BATSE channels 3 and 1 and $`L_{51}`$ is the luminosity in units of $`10^{51}`$ erg.s-1. In our model, we estimate time lags by cross-correlating profiles in different energy channels obtained from eq. 9. The first factor in eq. 9 is

$$\psi (t)=\frac{F_\mathrm{E}(t)}{E_\mathrm{p}(t)},$$
(37)
which behaves as $`t^{2x+y1.5}`$ during pulse decay. The sign of $`\mathrm{\Delta }=2x+y1.5`$ is of great importance in determining the time $`t_{\mathrm{max}}`$ of maximum count rate and the related time lags. If $`\mathrm{\Delta }<0`$ the function $`\psi (t)`$ has a maximum at some early time $`t_\mathrm{m}`$ before decreasing as $`t^\mathrm{\Delta }`$ while it steadily increases for $`\mathrm{\Delta }>0`$. The second factor in eq. 9 is $`\phi (E_\mathrm{p})`$ so that $`t_{\mathrm{max}}`$ is solution of the implicit equation

$$\frac{\dot{\psi }(t_{\mathrm{max}})}{\psi (t_{\mathrm{max}})}+\frac{\dot{E}_\mathrm{p}^{\mathrm{max}}}{E_\mathrm{p}^{\mathrm{max}}}\frac{\mathrm{d}\mathrm{log}\phi }{\mathrm{d}\mathrm{log}E_\mathrm{p}}|_{E_\mathrm{p}^{\mathrm{max}}}=0,$$
(38)
where $`E_\mathrm{p}^{\mathrm{max}}`$ and $`\dot{E}_\mathrm{p}^{\mathrm{max}}`$ are the values of $`E_\mathrm{p}`$ and its time derivative at $`t=t_{\mathrm{max}}`$. Since in most pulses the evolution of $`E_\mathrm{p}`$ precedes the count rate, $`\dot{E}_\mathrm{p}^{\mathrm{max}}`$ is negative. If the low energy slope of the Band spectrum $`\alpha 1`$, the derivative $`\frac{\mathrm{d}\mathrm{log}\phi }{\mathrm{d}\mathrm{log}E_\mathrm{p}}|_{E_\mathrm{p}^{\mathrm{max}}}`$ is positive (see Fig. 1) and eq. 38 then shows that $`\dot{\psi }(t_{\mathrm{max}})>0`$ which, for $`\mathrm{\Delta }<0`$, leads to

$$t_{\mathrm{max}}<t_\mathrm{m},$$
(39)
which provides a strict upper limit on the time lags between different energy channels

$$\mathrm{\Delta }t<t_\mathrm{m}.$$
(40)
(If $`\alpha >1`$, $`\frac{\mathrm{d}\mathrm{log}\phi }{\mathrm{d}\mathrm{log}E_\mathrm{p}}|_{E_\mathrm{p}^{\mathrm{max}}}`$ can be weakly negative at large $`E_\mathrm{p}`$ but even in this case $`t_{\mathrm{max}}`$ never greatly exceeds $`t_\mathrm{m}`$.)

When $`\mathrm{\Delta }0`$, no constraint such (40) applies and the time lags can be quite large. Daigne & Mochkovitch (1998) who adopted $`x=y=1/2`$ ($`\mathrm{\Delta }=0`$) obtained lags of several seconds between BATSE channels 1 and 3 (or 4) for a pulse lasting about 10 s, while currently observed values are in the range $`10^2`$ to a few $`10^1`$ s. Moreover, for $`\mathrm{\Delta }>0`$ the time lags increase with pulse hardness, in contradiction with the observations.

Conversely, with $`x=y=1/4`$ the time lags are short ($`\mathrm{\Delta }t_{31}0.5`$ s) even for long pulses and they decrease with increasing hardness. Figure 6 shows the time lag – hardness ratio relation given by our model superimposed to the Norris et al. (2000) results for BATSE bursts. The thick grey lines in Fig. 6 correspond to sequences of pulses of comparable duration ($`t_{90}10`$ s) obtained with similar distributions of the Lorentz factor (a slow shell with $`\mathrm{\Gamma }_0=100`$ decelerating a fast wind with $`\mathrm{\Gamma }_1=300`$, 400 or 600) and a varying value of the (isotropic) injected power (from $`\mathrm{5\; 10}^{51}`$ to $`10^{54}`$ erg.s-1). Even if differences in duration and redshift will contribute to add scatter to their distribution it appears that synthetic pulses populate the same triangular domain as observed ones.

We finally checked if our model was able to reproduce the time lag – luminosity correlation (eq. 36). When $`x=y=1/4`$, we do find that the lags decrease with increasing luminosity. This is a consequence of the HIC and the time lags – hardness ratio relation discussed above. The results are shown in Fig. 7 where the three lines correspond to the wind cases with $`\mathrm{\Gamma }_1=300`$, 400 and 600 already considered in Fig. 6. It can be seen that there is an overall agreement between the model predictions and eq. 36. However, at low luminosities and large time lags we obtain a rather wide strip instead of a single relation as eq. 36. If this is confirmed by the analysis of more GRBs with known redshifts, the time lag – luminosity correlation will still be useful for statistical studies of large burst samples but may be quite inaccurate to estimate the luminosity of a specific event.

## 7 Conclusion

We have developed a simple model where GRB pulses are produced when a rapid part of a relativistic outflow is decelerated by a comparatively slower shell. We do not solve the true hydrodynamical problem but rather assume that the slow shell “accretes” the fast moving material which allows to obtain an analytical solution for the dissipated power $`\dot{}`$. During pulse decay $`\dot{}t^{3/2}`$ as $`\dot{}t^3`$ when the evolution is fixed by shell geometry. To compute the spectral evolution of our synthetic pulses we parametrize the peak energy as $`E_\mathrm{p}\rho ^xϵ^y\mathrm{\Gamma }`$ where $`\rho `$, $`ϵc^2`$ and $`\mathrm{\Gamma }`$ are respectively the post shock values of the density, dissipated energy (per unit mass) and Lorentz factor. At late times, we get $`E_\mathrm{p}t^{(2x+y)}`$ which constraints $`x`$ and $`y`$ since in most observed bursts $`E_\mathrm{p}t^\delta _{}`$ with $`\delta _{}1.5`$. The synchrotron process with standard equipartition assumptions corresponds to $`x=1/2`$ and $`y=5/2`$ (i.e. $`2x+y=3.5`$) and gives a much too steep spectral evolution. One has then to suppose that the equipartion parameters $`\alpha _\mathrm{e}`$, $`\alpha _\mathrm{B}`$ and $`\zeta `$ vary with $`\rho `$ or/and $`ϵ`$ to reduce $`x`$ and $`y`$ (a possible alternative being that energy is radiated by another process – different from synchrotron – but which can still be approximated by eq. 33).

We have considered the case $`x=y=1/4`$ and the resulting pulses then have temporal and spectral properties in excellent agreement with the observations. They follow both the HIC and the HFC during the decay phase and the time lags between energy channels decrease with increasing pulse hardness and peak luminosity. We therefore conclude that if GRB pulses are produced by internal shocks, their temporal and spectral properties are probably governed by the hydrodynamics of the flow rather by the geometry of the emitting shells.