text:
# Recombination of low energy electrons with U28+

## 1 Introduction

Atomic processes are of great importance from astrophysical and other various point of view. Especially the low energy recombination cross sections and rate coefficients are required for understanding the fusion and astrophysical plasmas. In particular dielectronic recombination (DR) plays a significant role in astrophysics because it is the dominant electron-ion recombination process for most ions in low density, photo-ionized and electron-ionized cosmic plasmas . Accurate rate coefficients are needed to calculate reliably the ionization balance, thermal structure and line emission of cosmic plasmas. Most importantly the low energy positron antiproton recombination provides a challenging scheme for the production of antihydrogen.

It is well-known that due to the presence of additional recombination channels such as DR, the rate coefficients are found to be larger than radiative recombination (RR) rates for many electron complex ions. In DR process the incident electron is captured in a doubly-excited state of the compound ion, which is then stabilized by photoemission. This process originally suggested by J. Sayers and was first considered by Massey and Bates in the study of ionospheric oxygen. Electron-ion recombination has been measured directly in the laboratory since early 1980’s . More recently the use of heavy-ion accelerators and electron coolers of ion storage rings has greatly advanced the experiment . Recombination rates for various ions have been measured at electron energies from threshold to hundreds of electron volts (eV) with a fraction-of-an-eV resolution . For few-electron ions the measured rates were found to be in good agreement with theoretical predictions which included the contribution of DR resonances on top of the RR background, e.g., in He+ , Li-like C4+  and Ar15+ , and B-like Ar13+ . However, more complicated ions, e.g., Au50+ , U28+ and Au25+, showed complicated resonance spectra and strongly-enhanced recombination rates at low electron energies. The Au25+ ion has been studied in detail by Gribakin et al. and Flambaum et al. using statistical methods. They suggested that the strongly enhanced low energy electron recombination observed in this ion is mediated by complex multiply-excited states rather than simple dielectronic resonances and the dielectronic resonances play the role of doorways to the electron capture process. The statistical method developed by Flambaum et al. is based on the assumption of strong (chaotic) configuration mixing. This assumption has been verified by Gribakin and Sahoo  in a recent study. However, for U28+, no theory so far has described the low energy DR process successfully.

Historically the recombination rate enhancement was observed first in U28+. This measurement has been performed in a merged beam experiments at UNILAC accelerator in Darmstadt and at heavy-ion storage ring TSR in Heidelberg. The experiment found rate enhancement in the U28+ spectrum exceeds the theoretical calculation by at least a factor 20 in the energy range below 10 eV. Later in 1998, they extended the experiment to high energies up to 420 eV  and a comparison has been made with the theoretical calculations which is based on distorted wave approximations. The DR cross sections are calculated in this method are able to explain the main resonant features in the range 80–180 eV, but failed to identify the resonances and reproduce the rate at smaller energies. The cross sections involving the excitations from the 5s25p2 ground state configuration calculated in this method are well described by using either semirelativistic wave functions as found in AUTOSTRUCTURE codes or fully relativistic wave functions in HULLAC codes. However, the resonance structure observed at low energies i.e. below 80 eV largely remains unexplained. They finally concluded that in complex ions, particularly in U28+ ion, what are the resonances just above the threshold and how they contribute to the low energy recombination remain a ’mystery’. In this paper we performed an extensive study of the excited spectra and eigenstates of U28+ near its ionization threshold and calculated recombination rate coefficients for electron recombination with U28+. We identified some of the resonances near the threshold those contribute significantly to the low energy recombination. The present results are found to be in good agreement with the experiment. This work develops further a statistical theory towards the full understanding of the mechanism of low energy electron recombination with U28+ and other similar complex ions.

## 2 Many-electron excitations

We consider the recombination of an electron with U28+. Due to electron correlation the slow electron can be captured in one of the excited states of the compound ion U27+. It has 65 electrons and its ground state configuration belongs to 1$`s^2`$……..5$`s^2`$5$`p^3`$. Figure 1 shows the spectrum of relativistic $`nlj`$ orbitals obtained by relativistic Hartree-Fock calculations. Atomic units (a.u.) are used unless otherwise stated. All the occupied orbitals below fermi level ($`31.19a.u.`$) are obtained by a self consistent calculation of U28+ground state. Each of excited orbitals above the fermi level are calculated by placing one electron on to it in the frozen 1$`s^2`$……..5$`s^2`$5$`p^2`$ (U28+) core. Our configuration interaction (CI) calculation shows that the ground state of U28+ and U27+ ions are characterized by their total angular momentum J $`=`$ 0.0 and 1.5 respectively. The difference between their total energies $``$27741.40 a.u. and $``$27771.37 a.u., gives us an estimation of ionization threshold (I.T.) $`=`$ 29.97 a.u. (815.184 eV).

The excited states of the ion are generated by transferring one, two, three, etc., electrons from the ground state into the empty orbitals above fermi level. Since we are interested in the excited spectrum of the ion near its ionization threshold (29.97 a.u.), we consider it as a system of having only 29 electrons above the frozen Kr-like 1$`s^2`$….4$`p^6`$ core. The number of many-electron states obtained by distributing 29 electrons over 40 relativistic orbitals, 4$`d_{3/2}`$ through to 8$`g_{9/2}`$ are huge in number. It is practically impossible to perform a CI calculation for all of them. We construct the excite spectrum by using mean field approach by calculating their mean energies $`E_i`$, and number of many electron states $`N_i`$ associated with each of them:

$$E_i=E_{core}+\underset{a}{}ϵ_an_a+\underset{a<b}{}\frac{n_a(N_b\delta _{ab})}{1+\delta _{ab}}U_{ab},$$
(1)
$$N_i=\underset{a}{}\frac{g_a!}{n_a!(g_an_a)!},$$
(2)
where $`n_a`$ are the orbital occupation numbers of the relativistic orbitals in a given configuration and $`_an_a=n`$. $`ϵ_a=<a|H_{core}|a>`$ is the s

ingle-particle energy of the orbital $`a`$ in the field of the core, $`g_a=2j_a+1`$, and $`U_{ab}`$ are the average Coulomb matrix elements for the electrons in orbitals $`a`$ and $`b`$ ( direct minus exchange):

$$U_{ab}=\frac{g_a}{g_a\delta _{ab}}\left[R_{abab}^{(0)}\underset{\lambda }{}\delta _pR_{abba}^{(\lambda )}\left\{\genfrac{}{}{0pt}{}{j_a}{\frac{1}{2}}\genfrac{}{}{0pt}{}{ȷ_b}{\frac{1}{2}}\genfrac{}{}{0pt}{}{\lambda }{0}\right\}^2\right]$$
(3)
$`R_{abba}^{(\lambda )}`$ is the two-body radial Coulomb integral of $`\lambda `$ multipole, and $`\delta _p=1`$ when $`l_a+l_b+\lambda `$ is even and 0 otherwise. Using Eqs. (1)-(3) we obtained about 353 configurations within $`\pm `$ 1 a.u. of ionization threshold. They comprise a total of $`1.9\times 10^5`$ many electron states. The single-particle spectrum of U28+ does not show large gaps. Owing to the “gapless”single-particle spectrum, the density increases rapidly as a function of energy, as described by the Fermi-gas-model ansatz

$$\rho (E)=AE^\nu \mathrm{exp}(a\sqrt{E}),$$
(4)
with $`A`$ = 0.0885, $`\nu `$ = 2.33, and $`a`$ = 3.99 a.u., where $`E`$ is the energy above the ground state in atomic units. Figure 2 shows the level density calculated by averaging with a Gaussian with 1 a.u. variance. This figure depicts the level densities for both odd and even parity configurations. Also included in the figure are only the doubly excited configurations of even and odd parity found within $`\pm `$ 1 a.u. of ionization threshold. It is found that the level density of the odd configurations near its ionization threshold is of the order is of $`1\times 10^5`$ and

those of even configurations it is of the order of $`1.3\times 10^4`$. This provides an evidence that this system is characterized by huge level density. This system may be compared with Au25+, in which the level density near its ionization threshold is about $`1\times 10^7`$. This shows that U28+ ion is less complicated than Au25+ from the view point of dense excited spectrum. Apart from this one can estimate the mean level spacing from the distribution of total angular momentum J as shown in Figure 3.

Ionic eigenstates are characterized by their total angular momentum and parity $`J^\pi `$, and are $`2J+1`$ times degenerate. Therefore the total level density can be broken into a sum of partial level densities: $`\rho (E)=_{J^\pi }(2J+1)\rho _{J^\pi }(E)`$. The excitation spectrum of U28+ near the ionization threshold, $`E=I29.97`$ a.u., contains many $`J`$ ranging from $`\frac{1}{2}`$ to $`\frac{25}{2}`$. Their distribution is in agreement with statistical theory , which predicts that at a given energy $`\rho _{J^\pi }`$ are proportional to the function

$$f(J)=\frac{2(2J+1)}{(2J_m+1)^2}\mathrm{exp}\left[\frac{(2J+1)^2}{2(2J_m+1)^2}\right],$$
(5)
where $`J_m`$ is the most abundant $`J`$ value. Numerically for U27+ we find $`J_m\frac{7}{2}`$.

Using Eq. (5) one can estimate the partial densities by $`\rho _J=\rho _{J^+}+\rho _J^{}f(J)\rho /2J+1`$, where $`2J+1`$ is an average over $`f(J)`$. For the most abundant angular momenta $`JJ_m`$, and assuming $`\rho _{J^+}\rho _J^{}`$, we have $`\rho _{J^\pi }(E)=A_{J^\pi }E^\nu \mathrm{exp}(a\sqrt{E})`$. Near the ionization threshold this gives $`\rho _{J^\pi }1.8\times 10^3`$ a.u. (Fig. 3), which means that the spacing between the multiply-excited states with a given $`J^\pi `$ is small: $`D=1/\rho _{J^\pi }15`$ meV. Whereas in Au25+ it is about 1 meV, i.e. 15 times larger than that found in Au25+. This situation explains why individual resonances appear in recombination rates of U28+ but the same is not observed in Au25+ experimentally even at an energy resolution of 0.1 eV. However, the large density of multiply-excited states is only a “kinematic” reason behind the experimental observation. To explain the fact that the electron can actually be captured into these states, we need to analyze the dynamics of electron capture and show that the residual Coulomb interaction between the electrons (i.e. that beyond the mean field) makes for an efficient capture and accounts for the observed enhanced recombination rate.

## 3 Configuration mixing

Taking into account the fact that the rasidual Coulomb interaction is the key problem in many-electron processes, we construct the basis of many-electron states $`\mathrm{\Phi }_k`$ from single-particle (e.g., Dirac-Fock) orbitals, and solve the eigenvalue problem for the Hamiltonian matrix $`H_{ik}=\mathrm{\Phi }_i|\widehat{H}|\mathrm{\Phi }_k`$, which yields the eigenvalues $`E_\nu `$ and eigenstates $`|\mathrm{\Psi }_\nu =_kC_k^{(\nu )}|\mathrm{\Phi }_k`$ of the system (configuration interaction method). We performed two sets of model calculations. One includs all the configurations within $`\pm 1`$ a.u. of ionization threshold of U28+, which produces 2516 states for $`J^\pi `$=$`\frac{7}{2}^{}`$ sequence. Similarly the second calculation includes all the dielectronic (doorway) states within $`\pm 1`$ a.u. of ionization threshold and produces 108 states for $`J^\pi `$=$`\frac{7}{2}^{}`$ sequence. This shows that the number of states associated with the dielectronic configurations for a given $`J^\pi `$ are not large in number. As a result we performed the full CI calculations for Hamiltonian matrix of size 108 and 2516 respectively and obtained the eigen values and eigenstate components. To study the mixing between the doubly excited states we analyze the eigenstate components by calculating the weight of a given doorway configuration shown in Figure 4(a)-(e). The weight ($`w`$) of a doorway state can be calculated as $`\mathrm{\Sigma }_{k=1,Nc}|C_k^\nu |^2`$, where $`N_c`$ is the number of states in each configuration. It has been shown in Ref that when there is strong and uniform configuration mixing the weights significantly reduce from 1, but in the present case the weights go down from 1 but not significantly. So one can say that these doubly excited (doorways) states weakly mix with each other. When they are included in the large calculation (2516 states), their weights significantly go down from 1 as shown in Figure 4(f)-(j). This gives a signature of strong mixing i.e., the doorways show a significant mixing with multiply excited states. It has also been observed that the multiply excited configurations mix with each other quite comfortably. This mixing is mainly responsible for an enhancement of recombination rates over RR. In the recombination cross sections, the doorways which do not mix completely with either multiply excited states or with other doorways (in a sense they remain ’isolated’) appear in the form of narrow peaks (resonances) which explains the experiment as well as the theory.

On the other hand, when the level density is high and the two-body interaction is sufficiently strong the system is driven into a regime of many-body quantum chaos, where the effect of configuration mixing can be described statistically in the case of Au25+. This regime is characterized by the following . (i) Each eigenstate contains a large number $`N`$ of principal components $`C_k^{(\nu )}1/\sqrt{N}`$, corresponding to the basis states $`\mathrm{\Phi }_k`$ which are strongly mixed together. (ii) Owing to the strong mixing, the only good quantum numbers that can be used to classify the eigenstates, are the exactly conserved total angular momentum and parity $`J^\pi `$ and the energy. (iii) The degree of mixing in this regime is in some sense complete, i.e. all basis states that can be mixed (within a certain energy rage, see below) are mixed. The notion of configurations based on the single-particle orbitals becomes largely irrelevant for the purpose of classifying the eigenstates. Each eigenstate contains substantial contributions of a few nearby configurations. As mentioned above, in the present case we found that there is a weak configuration mixing between the dielectronic doorway configurations and complicated multiply excited states, and the multiply excited states show a substantially strong configuration mixing with each other. But the degree of mixing is not sufficiently strong to drive the system into chaotic regime.

This feature can be studied from the inverse participation ratio (IPR): $`\mathrm{\Sigma }_{j=1,Nc}|C_k^\nu |^4`$.

Figure 5 shows the inverse participation ratio for the two sets of model calculations. The top figure depicts the inverse participation ratio of the doorways only and the middle one for the large (2516 states) calculation that includes all the configurations (both dielectronic and multiply excited states). We found most of the doorways are associated with number of states either 4 or 2 or even less. The top figure shows that the inverse participation ratio lies flat between 2 and 4 except at the energy range where the eigen energies are close to the ionization threshold. This indicates that a few number of the doorways having energies close to the threshold participate in mixing. However the middle figure shows a picture of strong but non uniform mixing which involves a lots of multiply excited states. When we included the IPR of doorways with the IPR obtained from large calculation as shown in the bottom figure, it lies well below as indicated by solid circles. It also provides us with an information that some of the doorways take part in mixing with multiply excited configurations and other doorways are weakly mixed or even remain isolated.

This strong mixing takes place in a certain energy range $`|E_kE_\nu |\mathrm{\Gamma }_{\mathrm{spr}}`$, where $`E_kH_{kk}`$ is the mean energy of the basis state and $`\mathrm{\Gamma }_{\mathrm{spr}}`$ is the so-called spreading width. More precisely, the mean-squared value of $`C_k^{(\nu )}`$ as a function of $`E_kE_\nu `$, can be described by a Breit-Wigner (BW) formula

$$\overline{\left|C_k^{(\nu )}\right|^2}=N^1\frac{\mathrm{\Gamma }_{\mathrm{spr}}^2/4}{(E_kE_\nu )^2+\mathrm{\Gamma }_{\mathrm{spr}}^2/4},$$
(6)
with $`N=\pi \mathrm{\Gamma }_{\mathrm{spr}}/2D`$ fixed by normalization $`_k\left|C_k^{(\nu )}\right|^2\overline{\left|C_k^{(\nu )}\right|^2}𝑑E_k/D=1`$.

The mean-squared components are obtained by averaging over the basis states associated with each of the doorway configurations and are plotted as a function of eigen energies. We obtained $`\mathrm{\Gamma }_{\mathrm{spr}}`$ of doorways (close to the ionization threshold) from the BW fit which is shown in Figure 6(a)-(e). From the BW fit we observed that the $`\mathrm{\Gamma }_{\mathrm{spr}}`$ is not constant and it varies from 0.1 a.u. to 0.2 a.u. in these doorways. There may be doorways which even show much less spreading. Roughly one can say the mixing takes place within 0.1 a.u.. So we calculated the recombination rates with two different values of spreading width which can be found in section IV. It may be pointed out that though the single BW fit is not as accurate as it should be, still it gives an estimation of the important quantity i.e. $`\mathrm{\Gamma }_{\mathrm{spr}}`$. It is worth mentioning that the spreading width found in Au25+ is about 0.5 a.u. and these value does not change if one performs calculations by adding more configurations.

## 4 Recombination

### 4.1 Theory

For low-energy electrons the contribution of the autoionising states (resonances) to the recombination cross section is given by (see, e.g., Ref. )

$$\sigma _r=\frac{\pi }{k^2}\underset{\nu }{}\frac{2J+1}{2(2J_i+1)}\frac{\mathrm{\Gamma }_\nu ^{(r)}\mathrm{\Gamma }_\nu ^{(a)}}{(\epsilon \epsilon _\nu )^2+\mathrm{\Gamma }_\nu ^2/4},$$
(7)
where $`\epsilon =k^2/2`$ is the electron energy, $`J_i`$ is the angular momentum of the initial (ground) target state, $`J`$ are the angular momenta of the resonances, $`\epsilon _\nu =E_\nu I`$ is the position of the $`\nu `$th resonance relative to the ionization threshold of the compound (final-state) ion, and $`\mathrm{\Gamma }_\nu ^{(a)}`$, $`\mathrm{\Gamma }_\nu ^{(r)}`$, and $`\mathrm{\Gamma }_\nu =\mathrm{\Gamma }_\nu ^{(r)}+\mathrm{\Gamma }_\nu ^{(a)}`$ are its autoionisation, radiative, and total widths, respectively . When the resonance spectrum is dense, $`\sigma _r`$ can be averaged over an energy interval $`\mathrm{\Delta }\epsilon `$ which contains many resonances, $`D\mathrm{\Delta }\epsilon \epsilon `$, yielding

$$\overline{\sigma }_r=\frac{2\pi ^2}{k^2}\underset{J^\pi }{}\frac{2J+1}{2(2J_i+1)D}\frac{\mathrm{\Gamma }_\nu ^{(r)}\mathrm{\Gamma }_\nu ^{(a)}}{\mathrm{\Gamma }_\nu ^{(r)}+\mathrm{\Gamma }_\nu ^{(a)}},$$
(8)
where $`\mathrm{}`$ means averaging. If the fluorescence yield, $`\omega _f\mathrm{\Gamma }_\nu ^{(r)}/(\mathrm{\Gamma }_\nu ^{(r)}+\mathrm{\Gamma }_\nu ^{(a)})`$, fluctuates weakly from resonance to resonance (see below), one can write $`\overline{\sigma }_r=\overline{\sigma }_c\omega _f`$, where

$$\overline{\sigma }_c=\frac{\pi ^2}{k^2}\underset{J^\pi }{}\frac{(2J+1)\mathrm{\Gamma }^{(a)}}{(2J_i+1)D}$$
(9)
is the energy-averaged capture cross section, and $`\mathrm{\Gamma }^{(a)}`$ is the average autoionisation width.

In a situation when there is a strong configuration mixing between the dielectronic doorway states and multiply excited states, the capture cross sections can be obtained as a sum over the single-electron excited states $`\alpha ,\beta `$ and hole states $`\gamma `$, as well as the partial waves $`lj`$ of the continuous-spectrum electron $`\epsilon `$. As a result, we have

$`\overline{\sigma }_c`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{k^2}}{\displaystyle \underset{\alpha \beta \gamma ,lj}{}}{\displaystyle \frac{\mathrm{\Gamma }_{\mathrm{spr}}}{(\epsilon \epsilon _\alpha \epsilon _\beta +\epsilon _\gamma )^2+\mathrm{\Gamma }_{\mathrm{spr}}^2/4}}{\displaystyle \underset{\lambda }{}}{\displaystyle \frac{\alpha ,\beta V_\lambda \gamma ,\epsilon lj}{2\lambda +1}}`$ (10)
$`\times `$ $`\left[\alpha ,\beta V_\lambda \gamma ,\epsilon lj(2\lambda +1){\displaystyle \underset{\lambda ^{}}{}}(1)^{\lambda +\lambda ^{}+1}\left\{{\displaystyle \genfrac{}{}{0pt}{}{\lambda }{\lambda ^{}}}{\displaystyle \genfrac{}{}{0pt}{}{j_\alpha }{j_\beta }}{\displaystyle \genfrac{}{}{0pt}{}{j}{j_\gamma }}\right\}\alpha ,\beta V_\lambda ^{}\epsilon lj,\gamma \right],`$
where $`\epsilon _\alpha `$, $`\epsilon _\beta `$ and $`\epsilon _\gamma `$ are the orbital energies, the two terms in square brackets represent the direct and exchange contributions, and $`\alpha ,\beta V_\lambda \gamma ,\epsilon lj`$ is the reduced Coulomb matrix element (see Ref.).

It is assumed that the energies of dielectronic doorway states relative to the threshold is given by $`\epsilon _\alpha `$ \+ $`\epsilon _\beta `$ \- $`\epsilon _\gamma `$. A more accurate value can be obtained by using mean field energies (configuration energies) of doorway configurations in Eq.(10). The effect of using these two different energies can be found in next subsection. We have also shown the the qualitative difference between the results obtained by using two different values of $`\mathrm{\Gamma }_{spr}`$. Because $`\mathrm{\Gamma }_{spr}`$ is well defined in the case of a strong and chaotic configuration mixing. However, it can not be properly defined if the mixing of the configurations is weak and non-uniform.

Equation (10) is directly applicable to targets with closed-shell ground states. If the target ground state contains partially occupied orbitals, a factor

$$\frac{n_\gamma }{2j_\gamma +1}\left(1\frac{n_\alpha }{2j_\alpha +1}\right)\left(1\frac{n_\beta }{2j_\beta +1}\right),$$
(11)
where $`n_\alpha `$, $`n_\beta `$, and $`n_\gamma `$ are the orbital occupation numbers in the ground state $`\mathrm{\Phi }_i`$, must be introduced on the right-hand side of Eq. (10). Steps similar to those that lead to Eq. (10) were used to obtain mean-squared matrix elements of operators between chaotic many-body states .

The chaotic nature of the multiply-excited states $`\mathrm{\Psi }_\nu `$ can also be employed to estimate their radiative widths $`\mathrm{\Gamma }_\nu ^{(r)}`$. Electron-photon interaction is described by a single-particle dipole operator $`\widehat{d}`$. Any excited electron in $`\mathrm{\Psi }_\nu `$ may emit a photon, thus leading to radiative stabilization of this state. The total photo-emission rate $`\mathrm{\Gamma }_\nu ^{(r)}`$ can be estimated as a weighted sum of the single-particle rates,

$$\mathrm{\Gamma }_\nu ^{(r)}\underset{\alpha ,\beta }{}\frac{4\omega _{\beta \alpha }^3}{3c^3}|\alpha \widehat{d}\beta |^2\frac{n_\beta }{2j_\beta +1}\left(1\frac{n_\alpha }{2j_\alpha +1}\right)_\nu ,$$
(12)
where $`\omega _{\beta \alpha }=\epsilon _\beta \epsilon _\alpha >0`$, $`\alpha \widehat{d}\beta `$ is the reduced dipole operator between the orbitals $`\alpha `$ and $`\beta `$, and $`\mathrm{}_\nu `$ is the mean occupation number factor. Since $`\mathrm{\Psi }_\nu `$ have large numbers of principal components $`N`$, their radiative widths display small $`1/\sqrt{N}`$ fluctuations. This can also be seen if one recalls that a chaotic multiply-excited state is coupled by photo-emission to many lower-lying states, and the total radiative width is the sum of a large number of (strongly fluctuating) partial widths. A similar effect is known in compound nucleus resonances in low-energy neutron scattering .

There is a certain similarity between Eqs. (10) and (12) and those for autoionisation and radiative rates obtained in a so-called configuration-average approximation . In both cases the answers involve squares or products of two-body Coulomb matrix elements \[see the direct and exchange terms in Eq. (10)\], or single-particle dipole amplitudes \[Eq. (12)\]. However, there are a number of important differences between the present results and the configuration-average approximation. The latter considers dielectronic recombination and introduces averaging over configurations as a means of simplifying the calculation. The DR cross section is averaged over an arbitrary energy interval $`\mathrm{\Delta }\epsilon `$, and only the configurations within this energy range contribute to the average. Effects of configuration mixing as well as level mixing within a configuration are neglected.

It is important to compare the radiative and autoionisation widths of chaotic multiply-excited states. Equation (12) shows that $`\mathrm{\Gamma }^{(r)}`$ is comparable to the single-particle radiative widths. On the other hand, the autoionisation width $`\mathrm{\Gamma }^{(a)}`$, is suppressed by a factor $`\left|C_k^{(\nu )}\right|^2N^1`$ relative to that of a typical dielectronic resonance. Therefore, in systems with dense spectra of chaotic multiply-excited states the autoionisation widths are small. Physically this happens because the coupling strength of a two-electron doorway state to the continuum is shared between many complicated multiply-excited eigenstates. As a result, the radiative width may dominate the total width of the resonances, $`\mathrm{\Gamma }^{(r)}\mathrm{\Gamma }^{(a)}`$, making their fluorescence yield close to unity. However, in the case of U28+, $`\mathrm{\Gamma }^{(r)}\mathrm{\Gamma }^{(a)}`$. Our numerical results for the recombination of U28+ presented in Sec. 4.2, confirm this scenario.

The resonance recombination cross section should be compared with the direct radiative recombination cross section

$$\sigma _d=\frac{32\pi }{3\sqrt{3}c^3}\frac{Z_i^2}{k^2}\mathrm{ln}\left(\frac{Z_i}{n_0k}\right),$$
(13)
obtained from the Kramers formula by summing over the principal quantum number of the final state . Here $`Z_i`$ is the ionic charge ($`Z_i=28`$ for U28+), and $`n_0`$ is the principal quantum number of the lowest unoccupied ionic orbital ($`n_0=5`$). Note that the direct and energy-averaged resonance recombination cross sections of Eqs. (13) and (8) have similar energy dependences.

### 4.2 Numerical results

Numerical calculations of the cross section from Eqs. (10) and (11) involve summation over the orbitals shown in Fig. 1 with electron partial waves up to $`i_{11/2}`$. The results of the calculations for the recombination rates are displayed in Figure 7. We calculated the recombination rates using two different energies, i. e., configuration energies and relativistic energies. It is found that both the calculations are in good agreement. However, the magnitude of recombination rates obtained from both the present calculations lie above the experimental rates. This is due to the fact that we used $`\omega _f`$ = 1 in the present calculations similar to that in Au25+. It is clear that in the present case $`\omega _f`$ is less than 1 and remains constant throughout the energy range considered. The exchange contributions are found to be about 200 times smaller than the direct one so we did not include these when calculated the final rates. The radiative rates are found to be smaller in magnitude in comparison to DR rates and have almost negligible interference with the resonances. It may be pointed out that the present calculation with spreading width $`\mathrm{\Gamma }_{spr}`$ $`0.15`$ a.u, does not show enough resonance structures as have been observed in the experiment. Hence, we performed another calculation with $`\mathrm{\Gamma }_{spr}`$ $`0.05`$ a.u. and compared with the results obtained using spreading with 0.15 a.u. shown in Figure 8. It is clear that the results do not show any change in magnitude. However, the calculation with $`\mathrm{\Gamma }_{spr}`$ $`0.05`$ interprets the resonance peaks very well which are in reasonable agreement with the experiment, though the position of peaks are different from the experiment. This is because the present energies of the doorway states are approximate and a few eV relative error is expected in this approach. The presence of narrow peaks in rate coefficients can be interpreted as: the dielectronic states which play the role doorway to the electron capture process weakly mix with each other as has been discussed in previous section and appear as single peaks. Because in statistical calculations one peak corresponds to one doorway state. It may be recalled that the distorted wave calculation, predicts the recombination rates in agreement with the experiment above 80 eV and below this energy the results are smaller in magnitude.

Table 1 lists the most important dielectronic doorway contributions to the dimensionless sum in Eq. (10), which also determines the ratio of the autoionisation width to the spacing between the resonances. In total they account for about two thirds of the total cross section. Although these transitions have been selected according to the size of their contribution, their energies from two different calculations are close to the threshold, as seen in the last two columns in Table 1. Indeed, the spreading of configurations discussed in Sec. 2 allows configurations near the threshold, $`|\mathrm{\Delta }E|\mathrm{\Gamma }_{\mathrm{spr}}`$, to contribute. On the other hand, the contribution of configurations lying far away from threshold, $`|\mathrm{\Delta }E|\mathrm{\Gamma }_{\mathrm{spr}}`$, is suppressed. It may be pointed out that the theory of Mitnik et. al. considers the excitations from 5$`s`$ orbital only. However, the present calculations show that excitations involving 4$`f`$ orbitals are more important especially at threshold since they contribute significantly to the recombination cross sections. As expected, the energy dependences of the resonance and direct recombination rates are very close, although the latter is about 200 times smaller. There is a good overall agreement between the resonance rate and experimental data at electron energies between 1 eV and 100 eV.

To compare with Au25+, we calculated the rate coefficients for the electron energy ranging from 1 eV to 100 eV.

Because the previous calculation was restricted to the energy range below 1 eV. Figure 9 shows results three different calculations and the experimental data. It may be seen that the present (we call our best calculation) calculation uses configuration energies shows a very good qualitative agreement with the experiment in comparison to the other two calculations. The RR is found be of smaller in magnitude. The most striking feature in this graph is that the experimental data quickly departs from the theoretical values as one goes to higher energies. This can be interpreted as that in this energy range the fluorescence yield certainly goes down from 1 which can be well understood from the Figure 10. In this figure we show the configuration energies as a function of total angular momentum J. It shows that due to the presence of a lots levels, many inelastic channels open up and they participate strongly. This situation makes the fluorescence yield less than 1 and hence the theory predicts recombination rates higher in magnitude.

Which in turn just opposite to the case of U28+ as shown in Figure 11.

Finally, we compare in table II, the present energies of the lowest configurations in U28+ with the energies calculated by using HULLAC code. As mentioned earlier that in Ref, energies have been calculated by using different codes such as multiconfiguration Hartree-Fock (MCHF), the AUTOSTRUCTURE code in the perturbative-relativistic \[AS(PR)\] and in the semirelativistic \[AS(SR)\] mode and the HULLAC code. On comparison we found the present energies are close to those obtained by Mitnik et. al.  using HULLAC code since both the calculations are fully relativistic.

## 5 Summary and outlook

We have shown that the dielectronic states weakly mix with each other and show a substantial mixing with complicated multiply excited states. This explains the mechanism of low energy recombination of electron with U28+. The present results are found to be in good agreement with the experimental data. We found that one must not ignore the excitations from 4$`f`$ orbitals while considering the low energy recombination of U28+ ion. On detailed study, we found the configuration mixing between the doubly excited states and multiply excited states is not complete or uniform. Although this theory is valid in a system having an extreme degree of configuration mixing , still it predicts quite good results for other systems where there is not enough strong configuration mixing, for example U28+. This work develops further a statistical theory to study the low energy recombination process.

## 6 acknowledgements

We thank Prof. A. Müller for providing experimental data of Au25+ and U28+ in numerical form and for his stimulating discussions. We also thank to Dr. C. Harabati for allowing us to make use of his codes as well as for his useful comments. Financial help from EPSRC is highly acknowledged.

## 7 References


format: text

text:
# Strongly Non-Arrhenius Self-Interstitial Diffusion in Vanadium

## Abstract

We study diffusion of self-interstitial atoms (SIAs) in vanadium via molecular dynamics simulations. The $`111`$-split interstitials are observed to diffuse one-dimensionally at low temperature, but rotate into other $`111`$ directions as the temperature is increased. The SIA diffusion is highly non-Arrhenius. At $`T<600`$ K, this behavior arises from temperature-dependent correlations. At $`T>600`$ K, the Arrhenius expression for thermally activated diffusion breaks down when the migration barriers become small compared to the thermal energy. This leads to Arrhenius diffusion kinetics at low $`T`$ and diffusivity proportional to temperature at high $`T`$.

The creation and migration of self-interstitial atoms (SIAs) are critical for microstructural evolution of materials in a variety of situations, such as in the high energy radiation environment of nuclear reactors Young et al. (1978) and in ion implantation Nordlund et al. (1999). Although SIA formation energies are much larger than typical thermal energies, they form in abundance during collision cascades induced by impinging energetic particles. SIAs in metals are typically very mobile (i.e., their migration barriers are relatively small) and hence play an important role in controlling the rates of many microstructural processes in such applications, in particular the phenomenon of void swelling.

Since SIA properties and mobilities are very difficult to determine experimentally, one often employs computer simulations Díaz de la Rubia and Guinan (1991); Wirth et al. (1997); Soneda and T. Díaz de la Rubia (1998); Wirth et al. (1999). For example, simulations of body centered cubic (bcc) iron (and several other bcc metals), have shown that SIAs preferentially lie along $`110`$ orientations but rotate into $`111`$-directions, where they can migrate easily using the crowdion configuration as transition state. Other simulations have suggested that SIA migration in vanadium is very similar to that in Fe Minashin and Ryabov (1996); Morishita et al. (1997, 2000, 2001). However, these empirical interatomic potential-based simulations are not consistent with recent first principles calculations that clearly show that the lowest energy SIA configuration in V is a $`111`$-split interstitial, rather than the $`110`$-split configuration found in Fe Han et al. (2002a). Interestingly, the first principles calculations also revealed that the $`111`$-oriented SIA migration energy is extraordinarily small ($`0.01`$eV), which explains the experimental observation of diffusion down to 4 K Coltman Jr. et al. (1975). SIA migration in Fe and V must therefore differ in their microscopic mechanisms.

We perform a series of molecular dynamics (MD) simulation of SIA migration in V using an improved interatomic potential for V Han et al. (2002b) (refit to experimental and first principles data Han et al. (2002a) to reproduce the stable interstitial configuration) to address this discrepancy. In particular, we examine SIA diffusion as a function of temperature to determine the SIA migration mechanisms. We find that while SIA migration in V is similar to that in bcc Fe in many respects, its temperature dependence is highly unusual, exhibiting strongly non-Arrhenius behavior and correlation effects. It is this non-Arrhenius behavior that is the focus of the present Letter.

The SIA is introduced in the form of a stable $`111`$-split interstitial and equilibrated for 10 ps at fixed temperature using a Langevin thermostat. The simulation was then switched to a microcanonical (NVE) ensemble in order to study SIA migration dynamics. Simulations were run at temperatures between 100 and 2000 K in a cubic simulation cell of edge length $`10a_0`$, where $`a_0`$ is the temperature dependent-lattice parameter (the linear thermal expansion coefficient was $`8.4\times 10^6`$ K<sup>-1</sup>). The diffusivity was measured by averaging over several 1 ns simulations. The interstitial position was identified by dividing the space into Wigner-Seitz (W-S) cells centered around each perfect crystal lattice site. Interstitials are located in W-S cells containing more than one atom.

Representative trajectories of the $`111`$-split interstitial center of mass, collected over the whole 1 ns simulations, are shown in Fig. 1. For each temperature, more than 1000 jumps were observed, where an SIA jump is the exchange of an atom between neighboring W-S cells. As seen in Fig. 1, the interstitial migration mechanism is strongly temperature dependent. For low and intermediate temperatures (100-600 K) the $`111`$-split interstitial executes a fully one-dimensional (1D) random walk along a $`111`$-direction during the $`1`$ ns simulation, as shown in Fig. 1(a). At $`T`$ 700 K, the $`111`$-split interstitial begins to make infrequent rotations from one $`111`$\- to another $`111`$-direction. This results in a three-dimensional (3D) trajectory that consists of long 1D random walk segments with abrupt reorientations, as seen in Fig. 1(b). As the temperature increases, the frequency of the rotation events increases and the lengths of the 1D trajectory segments decrease. At high temperatures, the rotation events become very frequent, leading to nearly isotropic diffusion (Fig. 1(d)).

Although these trajectories appear to be qualitatitively similar to those reported for other bcc metals (i.e. Fe and Mo) Wirth et al. (1997); Yu. N. Osetsky et al. (1997); Pasianot et al. (2000); Marian et al. (2001) (cf. Fig. 1(d) and Fig. 5 in Pasianot et al. (2000)), they differ in the elementary migration mechanism. The stable interstitial is the $`111`$-split configuration in V, but the $`110`$-split configuration in Fe and Mo. In the Fe and Mo cases, the split interstitial sits in the $`110`$-orientation until it is thermally activated into one of the $`111`$-directions where it can migrate easily before returning to a $`110`$-orientation Pasianot et al. (2000). There are no relaxation events of this type in interstitial migration in V. Here, the stable $`111`$-split interstitial migrates long distances and only requires significant thermal activation to reorient or rotate.

The temperature dependence of the rate of rotation of the split interstitials from one $`111`$ direction to another, $`\omega _r`$, in V is shown in Fig. 2. The data is well described by a conventional Arrhenius fit of the form $`\omega _r=\nu _0\mathrm{exp}[\mathrm{\Delta }E_r/k_BT]`$, suggesting that rotation is a simple thermally activated process. The activation energy, $`\mathrm{\Delta }E_r=0.44`$ eV, is consistent with the energy difference between the $`111`$ and $`110`$ configurations computed in first-principles (0.35 eV) and static calculations using the new interatomic potential (0.4 eV) Han et al. (2002b).

The diffusivity $`D`$ (solid symbols) of the $`111`$-split interstitial is shown in Fig. 3 for a temperature range between 100 K and 2000 K. $`D`$ was determined from $`D=R^2(t)/2t`$ for 1D diffusion, where the mean squared displacement $`R^2(t)`$ was calculated from the trajectory using standard averaging procedures Guinan et al. (1977). Note that although the rotations at higher $`T`$ change the dimensionality of the diffusion path, they do not contribute to transport and hence the diffusion mechanism is always one-dimensional. If the diffusivity were Arrhenius ($`D/a_0^2=\nu _0\mathrm{exp}[\mathrm{\Delta }E_d/k_BT]`$), the data in Fig. 3 would lie along a straight line. This is clearly not the case; Fig. 3 shows pronounced curvature - especially at high temperature. Although Arrhenius behavior is widely expected for diffusion in the solid state, it is clearly inapplicable here.

Non-Arrhenius behavior can have several different origins. The energy barrier could simply change with temperature as a result of thermal expansion, as has been argued for the self-diffusion in bcc metals via a vacancy mechanism Eftaxias and Hadjicontis (1989); Neumann and Beke (1990). The magnitude of the observed deviations from Arrhenius behavior is too large to attribute to such a mechanism. The existence of multiple reaction pathways with different energy barriers can also lead to curvature in Fig. 3. However, detailed examination of the atomic configuration during diffusion shows that there is no change in diffusion mechanism over the entire temperature range. Although rotations are first observed at $`700`$ K within the 1 ns duration of the simulations, strong deviations from Arrhenius behavior are observed already at lower temperatures. A third possibility is that the degree of correlations in the diffusion process (i.e., particle jumps retain memory and the random walk is non-Markovian) is temperature dependent. Indeed, examination of the SIA trajectories show that the $`111`$-interstitial has a higher probability of jumping back in the direction from whence it came, rather than forward along the same direction. At high temperatures, by contrast, this effect appears to be reversed.

We quantified this observation by measuring a correlation factor for split interstitial diffusion $`f`$, defined as $`f=2D/D_b`$, where $`D_b`$ is the “bare” diffusion constant defined as $`D_b=l_0^2n`$. Here, $`n`$ is the mean number of jumps/second and $`l_0=\sqrt{3}a_0/2`$ is the jump length (nearest neighbor distance in the bcc lattice). If the interstitial trajectory is described by a sequence of jump vectors $`\stackrel{}{l}_i`$, the correlation factor $`f`$ is alternatively given by $`f=1+2_i^{n1}\stackrel{}{l}_i\stackrel{}{l}_0/l_0^2`$, i.e. $`f=1`$ for an uncorrelated random walk. The inset to Fig. 3 shows the variation of $`f`$ with 1/$`k_BT`$. At low $`T`$ ($`T600`$ K), SIA motion is anticorrelated $`(f<1)`$. The effect of the correlations on the diffusivity can be isolated by plotting $`D/f`$ rather than $`D`$ in Fig. 3 (open symbols). The low temperature $`D/f`$ data $`(T500K)`$ lie along a nearly straight line with slope $`\mathrm{\Delta }E_d=0.018`$ eV. Hence, the temperature dependent correlation factor explains the relatively weak deviations from Arrhenius behavior at low $`T`$. The unusually small migration energy is consistent with experimental observations Coltman Jr. et al. (1975) and the first principle estimate Han et al. (2002a). One possible origin of the anticorrelations is that the finite relaxation time of the local environment around the SIA becomes less important as the thermal energy of the system increases.

As the temperature is increased beyond 300 K, the correlation factor rises quickly to a value greater than unity. $`f>1`$ is very unusual and may be thought of as correlated interstitial hopping over several barriers without completely thermalizing in between. This interpretation is consistent with the fact that the correlation corrected diffusivity, $`D/f`$, only yields a straight line at low $`T`$ but not at high $`T`$. We note that where $`D/f`$ is rising quickly, the energy barrier $`\mathrm{\Delta }E_d`$ obtained in the low temperature regime (where the behavior is Arrhenius) is smaller than the thermal energy. This is the source of the multiple interstitial hops at high $`T`$. Conventional derivations of activated escape over barriers Haenggi et al. (1990) usually assume $`\mathrm{\Delta }E_dk_BT`$ .

Rather than trying to apply the Arrhenius description to interstitial self-diffusion in this system, a more general model is the motion of a particle in a periodic potential at all temperatures (from $`k_BT\mathrm{\Delta }E_D`$ to $`k_BT\mathrm{\Delta }E_D`$). In the limit that the barrier height is completely negligible relative to the energy of the heat bath (i.e., a free particle), standard arguments predict $`D=k_BT/m\gamma ,`$ where $`m`$ is the mass of the particle and $`\gamma `$ a relaxation time scale associated with a velocity-dependent friction force. Since this free particle diffusion coefficient is expected to be a linear function of $`T`$, we replot the data from Fig. 3 on a linear temperature scale (Fig. 4). In this representation, the data is nearly linear, albeit with weak curvature at low temperature. This suggests that SIA diffusion in V is free particle like at high temperature ($`DT`$), but follows the normal Arrhenius, hopping dynamics at low $`T`$ ($`De^{\mathrm{\Delta }E_D/k_BT}`$). The deviation from linearity at low temperature and the deviation from Arrhenius behavior at high temperature suggests that a cross-over is occuring between the free and hopping particle limits.

In order to better understand this transition, we explicitly consider a simple one-dimensional particle of mass $`m`$ diffusing in a sinusoidal potential by numerically solving the Langevin equation

$$m\ddot{x}\gamma \dot{x}=ϵ/2\mathrm{cos}[x/\sigma ]+\eta ,$$
(1)
where $`\eta `$ is a Brownian white noise. Inserting values of $`m`$ and $`\sigma `$ for V and $`ϵ=\mathrm{\Delta }E_d`$, this model yields nearly linear diffusivity for $`1<k_BT/ϵ<10`$ for reasonable values of $`\gamma `$, followed by a crossover into Arrhenius behavior (see inset of Fig. 4) for $`k_BT/ϵ1`$. Changing the value of $`\gamma `$ shifts the temperature at which the transition from the particle hopping to free particle behavior is observed. This excellent agreement between the MD results and model predictions demonstrates that the observed strongly non-Arrhenius interstitial diffusion in V is a direct consequence of the relative magnitudes of the activation energy and the thermal energy.

It is interesting to compare the situation described here for V to that for bcc Fe. The crowdion mechanism enabling the easy interstitial migration along $`111`$ directions in V is also available in Fe, and estimates for the migration barrier $`\mathrm{\Delta }E_d=0.04`$ eV Wirth et al. (1997) are similar to the V case. However, measurements of the apparent activation energy for diffusivity, analogous to the one presented here (albeit over a smaller temperature range between 700 K and 1200 K) yield much larger values of $`\mathrm{\Delta }E_d=0.12`$ eV Marian et al. (2001) or $`\mathrm{\Delta }E_d=0.17`$ eV Soneda and T. Díaz de la Rubia (2001). This larger effective barrier is due to the fact that the $`111`$-split interstitial must be thermally excited from the $`110`$-state, i.e. the easy-diffusion configuration is not populated at all times as in V. The fraction of time during which the $`111`$-split interstitial is available for transport is given by $`F=P_{111}/(P_{111}+P_{110})`$, where $`P_{111}\mathrm{exp}[\mathrm{\Delta }E_f/k_BT]`$ and $`P_{110}\mathrm{exp}[\mathrm{\Delta }E_b/k_BT]`$ are probabilities for the interstitial to transform from the $`110`$ state to the $`111`$ state and back, respectively. Therefore, the interstitial diffusivity in Fe is the product of the diffusivity that the interstitial would have if it was always in the $`111`$ orientation (like in V) and $`F`$. For $`\mathrm{\Delta }E_f\mathrm{\Delta }E_b`$, the effective activation energy for interstitial diffusivity in Fe is $`(\mathrm{\Delta }E_D+\mathrm{\Delta }E_f)`$, but for $`\mathrm{\Delta }E_b\mathrm{\Delta }E_f`$, it is $`\mathrm{\Delta }E_D`$. Clearly, this implies for Fe that $`\mathrm{\Delta }E_f\mathrm{\Delta }E_b`$.

Molecular dynamics simulations of self-interstitial diffusion in bcc V were performed over an unusually wide temperature range (100-2000 K). Interstitial atoms in the $`111`$-split configuration migrate very fast one-dimensionally along $`111`$ directions during the 1 ns simulations. As $`T`$ is increased above $`600`$ K, rotations of the split-interstitial from one $`111`$ orientation to another occur with increasing regularity. The rotations can be described by Arrhenius kinetics with activation energy $`\mathrm{\Delta }E_r=0.44eV`$. At temperatures $`T<600`$ K, the diffusion exhibits significant anticorrelations. An Arrhenius analysis of the data (corrected for these anticorrelations) yields a very small migration energy barrier $`\mathrm{\Delta }E_d=0.018`$ eV. For $`T>600K`$, $`\mathrm{\Delta }E_d`$ is much smaller than the thermal energy and the Arrhenius expression is no longer applicable. The diffusivity then crosses over from Arrhenius to free particle type diffusion with increasing $`T`$. The fact that this type of diffusion is a linear function of $`T`$ rather than Arrhenius, as usually assumed, could have important implications for predicting the lifetimes of reactor components in vanadium and other bcc metals.

###### Acknowledgements.

We thank A. F. Voter, G. H. Gilmer, and J. A. Caro for useful discussions. This work was performed under the auspices of the U. S. Department of Energy, Office of Fusion Enery Sciences (DE-FG02-01ER54628) and Lawrence Livermore National Laboratory under Contract No. W-7405-Eng-48.


text:
# 1 Introduction

## 1 Introduction

The discovery of synchronized oscillations in the visual cortex and other brain regions has triggered significant research in artificial spike-based computation . While neurons in the central nervous system are about six orders of magnitude ”slower” than silicon-based elements, in both elementary computation time and signal transmission speed, their performance in networks often compares very favorably with their artificial counterparts even when reaction speed is concerned. In a sense, evolution may have been forced to develop extremely efficient computational schemes given available hardware limitations.

In a recent paper , we proposed new models for two common instances of such neural computation, winner-take-all and coincidence detection, featuring fast convergence and $`O(n)`$ network complexity. We saw that both computations could be achieved using a similar architecture, using global feedback inhibition in the first case, and global excitation in the second. In this paper, we further extend this computational architecture to $`k`$-winners-take-all and soft-winner-take-all.

Fast Winner-take-all (WTA) computation in  is based on the FitzHugh-Nagumo model, a well-known simplified version of the classical Hodgkin-Huxley model. Compared to previous WTA networks , it has significant computational advantages. The network’s initial states can be set arbitrarily, and convergence is guaranteed in at most two spiking periods, with a high computation resolution. The network’s complexity is linear in the number of inputs and its size can be adjusted at any time during the computation. As this paper shows, by modifying the starting point of the global inhibitory neuron’s charging mode, $`k`$-Winners-Take-All ($`k`$-WTA) can be computed instead. By running the charging mode independently, soft-Winner-Take-All (soft-WTA) can be computed. Both extensions inherit the advantages of the original WTA network.

After a brief review of the basic WTA network in section 2, $`k`$-WTA computation and soft-WTA computation are studied in Sections 3 and 4. Brief concluding remarks are offered in Section 5.

## 2 Winner-Take-All Network

The WTA network in  is based on the FitzHugh-Nagumo (FN) model :

$$\{\begin{array}{cc}\dot{v}=v(\alpha v)(v1)w+I\hfill & \\ \dot{w}=\beta v\gamma w\hfill & \end{array}$$

For appropriate parameter choices, there exists a unique equilibrium point for any given value of $`I`$, which is stable except for a finite range $`I_lII_h`$ where the system tends to a limit cycle. The steady-state value of $`v`$ at the stable equilibrium point increases with $`I`$.

The network structure is illustrated in Figure 1 where $`n`$ FN neurons receive external stimulating inputs $`I_i`$ and a global inhibition $`z`$. The dynamics of the FN neurons ($`i=1,\mathrm{},n`$) are

$$\{\begin{array}{cc}\dot{v}_i=v_i(\alpha v_i)(v_i1)w_i+I_iz\hfill & \\ \dot{w}_i=\beta v_i\gamma w_i\hfill & \end{array}$$

The dynamics of the global inhibitory neuron is

$$\dot{z}=\{\begin{array}{cc}k_c(zz_0)\mathrm{charging}\mathrm{mode}\hfill & \\ k_dz\mathrm{discharging}\mathrm{mode}\hfill & \end{array}$$
(1)
which starts charging if there is any FN neuron spiking in the network (i.e., if $`i`$, $`v_iv_0`$ for some given threshold $`v_0`$) and switches to discharging if the state $`z`$ is saturated. With a fast charging rate $`k_c`$ and a slow discharging rate $`k_d`$, the network computes the largest input (corresponding to the only spiking FN neuron) in at most two periods. Initial conditions can be set arbitrarily and the computation resolution is very high. Detailed analysis and discussions can be found in .

## 3 $`K`$-Winners-Take-All Network

$`K`$-WTA, a common variation of WTA computation where the output indicates for each neuron whether its input is among the $`k`$ largest, has been studied in such fields as competitive learning, pattern recognition and pattern classification . As Maass argued in , in principle a $`k`$-WTA network can replace a two-layer threshold circuit to perform most standard nonlinear computational operations.

Most $`k`$-WTA studies are based on steady-state stability analysis. Many models define the winners as the neurons with the largest initial states  or require initial conditions to be set precisely , making the networks not well suited to time-varying inputs. Others adopt particular design methodologies  but the network size or the number of winners is limited. $`K`$-WTA is also implemented in analog VLSI circuits , which extend the elegant WTA model in  but inherit its low resolution limit as well.

The neural network described in Section 2 can be easily extended to $`k`$-WTA computation, where an FN neuron spikes if and only if its input is among the $`k`$ largest. Indeed, as the global inhibition force decreases, the FN neurons enter the oscillation region rank-ordered by their inputs. Thus, while for WTA computation, the global inhibition neuron is charged after the first arrival, for $`k`$-WTA computation the charging moment is simply modified to capture the $`k^{\mathrm{th}}`$ arrival instead.

To this effect, we augment the dynamics of the FN neuron with an additional state variable $`u_i`$ (for simplicity, we shall still call FN neuron such a generalized element)

$$\{\begin{array}{cc}\dot{v}_i=v_i(\alpha v_i)(v_i1)w_i+I_iu_iz\hfill & \\ \dot{w}_i=\beta v_i\gamma w_i\hfill & \\ \dot{u}_i=k_u(\zeta _iu_0u_i)\hfill & \end{array}$$

where $`u_0`$ is a constant saturation value and $`k_u`$ the charging/discharging rate. The variable $`\zeta _i`$ takes two values, namely it switches to $`0`$ whenever $`z`$ approaches a saturation value $`z_0`$, else it switches to $`1`$ if $`v_i`$ exceeds a given threshold $`v_0`$. This make the dynamics of $`u_i`$ a local self-inhibition, which starts charging if the basic FN neuron spikes and discharges whenever the global inhibitory neuron spikes. Note that the specific form of the dynamics of $`u_i`$ can be more general, as long as the value of $`u_i`$ varies between $`0`$ and $`u_0`$, and the transition periods are very fast (which is satisfied here by choosing a large $`k_u`$).

The dynamics of the global inhibitory neuron is the same as (1), except that we start its charging mode if any $`k`$ FN neurons in the network spike. Such a moment can be captured by determining that $`_{i=1}^nu_i`$ approaches $`ku_0`$. Thus, if any FN neuron spikes, it excites only the corresponding local inhibitory portion but has no effect on the rest of the network. If there are $`k`$ local inhibitions turned on, the global inhibitory neuron is charged, which then releases all the local inhibitions and starts a new period.

Compared to the WTA network in Section 2, the basic principle underlying the $`k`$-WTA network described above is the same, exploiting the simple properties of the FN model. Thus, most of the computational advantages of the WTA network  are inherited by the $`k`$-WTA extension. In particular

* The initial conditions of the network can be set arbitrarily.
* With appropriate parameters, the computation can be completed at most in two periods, where the first period is affected by the initial conditions but the $`k`$ spiking neurons during the following periods are guaranteed to be those with the largest inputs. If the initial inhibitions are strong, the computation is completed in one period.
* Since initial conditions are immaterial and the computation speed is very fast, the $`k`$-WTA network is able to track time-varying inputs. Moreover, since the network complexity is $`O(n)`$, individual FN neurons can be added or removed at any time during the computation.
* The inputs $`I_i`$ should be lower-bounded by $`I_l`$, the lower threshold of the FN oscillation region. It should also be upper-bounded to set inhibition saturations, although the upper bound value is not restricted.
* The computation resolution also follows that of the WTA network. It can be improved by decreasing the discharging rate $`k_d`$, as well as the relaxation time of the FN neurons.
* FN neurons receiving equal inputs behave identically, which means that the $`k`$-WTA computation may generate more than $`k`$ winners in this particular case.

* ###### Example 3.1

:  The result is illustrated in simulation in Figure 2, with $`n=10`$ and $`k=3`$. The parameters of the FN neurons are set as $`\alpha =5.32,\beta =3,\gamma =0.1`$, with spiking threshold $`v_0=5`$. The parameters of the local inhibition are $`u_0=160,k_u=100`$. The inputs $`I_i`$ are chosen randomly from $`20`$ to $`125`$. The parameters of the global neuron are $`z_0=240,k_c=100,k_d=1/40`$. All initial conditions are chosen arbitrarily. The three spiking neurons after the first charging of the global neuron are those with the three largest inputs.

Note that the output frequency is determined mainly by the global neuron’s dynamics and the value of the $`k^{\mathrm{th}}`$ largest input. It can be increased by increasing the global neuron’s discharging rate after the first winner spikes so as to facilitate the other winners’ spiking. $`\mathbf{}`$

* ###### Example 3.2

:  Figure 3 illustrates a simulation result with $`n=3`$ and $`k=2`$. The parameters are the same as those in Example 3.1. The inputs keep varying and switch winning positions several times. The spiking neurons always track the two largest inputs. $`\mathbf{}`$

## 4 Soft-Winner-Take-All

Soft-WTA  (or softmax) is another variation of WTA computation, where the outputs reflect the rank of all inputs according to their size. Although soft-WTA is a very powerful primitive  in that it can be used to compute any continuous function, its “neural” implementation is complex. Recently,  studied soft-WTA as an optimization problem not based on a biologically plausible mechanism;  presented a hardware model of selective visual attention which lets the attention switch between the selected inputs, but whose switching order does not completely reflect the input ranks. In this section, we develop a simple neural network which computes soft-WTA very fast and generates spiking outputs rank-ordered by their inputs.

Let $`k=n`$ in the $`k`$-WTA network described in Section 3. Then we get a pre-ordered spiking sequence in each stable period, since all the FN neurons enter the oscillation region and then spike rank-ordered by their inputs. However, such an $`n^{\mathrm{th}}`$ arrival moment may not be measurable if the number of inputs $`n`$ is unknown or time-varying. To avoid this problem and make the solution more general, we let the charging mode of the global inhibitory neuron start only if the inhibition $`z`$ is lower than a given bound $`z_{low}`$. The spikings of all the FN neurons in the network are guaranteed by the condition that

$$z_{low}<I_l+I_{min}$$

where $`I_l`$ is the lower bound of the oscillation region of the FN model and $`I_{min}`$ is the minimum input value.

* ###### Example 4.1

:  Figure 4 illustrates the result in simulation with $`n=10`$. The parameters are the same as those in Example 3.1. The inputs $`I_i`$ are distributed uniformly between $`80`$ and $`120`$. The inhibition lower bound is $`z_{low}=60`$. Initial conditions are chosen arbitrarily. The computation is completed in the second period, during and after which the spiking times of the FN neurons are ranked by their inputs. $`\mathbf{}`$

The simple soft-WTA network presented above inherits the main computational advantages of our WTA and $`k`$-WTA networks. Initial conditions can be arbitrary, the computation is completed in at most two periods, network complexity is linear and neurons can be added or removed at any time. We expect it to be effective in many applications such as selective attention, associative memory and competitive learning, and also to provide an efficient desynchronization mechanism for perceptional binding .

## 5 Concluding Remarks

Basic neural computations such as winner-take-all, $`k`$-winners-take-all, soft-winner-take-all, and coincidence detection can all be implemented using a common architecture and biological plausible neuron models. Fast and robust convergence is guaranteed, and time-varying inputs can be tracked. Further research will study models of higher level brain functions, such as perception, based on the WTA networks and on general nonlinear synchronization mechanisms derived in .

Acknowledgments: This work was supported in part by a grant from the National Institutes of Health. The authors benefited from stimulating discussions with Matthew Tresch.


text:
# DIAGNOSTICS OF NEUTRON STAR AND BLACK HOLE X-RAY BINARIES WITH X-RAY SHOT WIDTHS

## 1 Introduction

Complex X-ray emission variability at wide timescales has been observed in accreting stellar mass black hole candidate (BHC) and weakly magnetized neutron stars (NS) systems (for reviews see Tanaka & Lewin, 1995; van der Klis, 1995). On the subsecond timescale, aperiodic fluctuations are common in these systems, as first discovered by Oda et al. (1971) in the prototype BHC Cyg X-1 in its hard state with a large amplitude (up to $`40\%`$ rms of the mean flux). Many properties of X-ray rapid variability from Cyg X-1, e.g., the Fourier power spectrum, the mean and variance, and the autocorrelation function of the lightcurve, could be described as the superposition of individual shots with the characteristic duration of a few tenths of a second (Terrel, 1972; Oda, 1977; Sutherland, Weisskopf & Kahn, 1978; Nolan et al., 1981; Meekings et al., 1984; Miyamoto et al., 1988, 1992; Lochner, Swank, & Szymkowiak, 1991). The average shot profiles of Cyg X-1 have been obtained by superposing many shots by aligning their peaks in the lightcurves obtained with Ginga (Negoro, Miyamoto, & Kitamoto, 1994; Negoro, Kitamoto, & Mineshige, 2001) and RXTE (Feng, Li, & Chen, 1999).

Studying the short timescale variability of the X-ray emission of accreting BHC or NS systems is important for understanding the emitting region and emission process of high energy photons, as well as the nature of the compact object. Timing studies carried out in the time domain can make direct measurements of, e.g., power, spectral lag, of a stochastic process at high frequencies (short timescales). Li (2001) has developed a technique of timescale analysis, with which timescale spectra can be derived directly in the time domain. The timescale analysis of the correlation coefficient between intensity and hardness ratio of Cyg X-1 (Li, Feng, & Chen, 1999) is consistent with the results from superposed shots (Feng, Li, & Chen, 1999; Negoro, Kitamoto, & Mineshige, 2001). The time domain power density at a certain timescale can be defined as the differential of the variation power in a lightcurve with the corresponding time step (Li, 2001). The timescale spectra of power density of a sample of NS and BHC binaries were analyzed and the results indicate that the characteristic timescale of hard X-ray shots is $`0.1`$ s for BHCs, and much shorter ($`1`$ ms) for NSs (Li & Muraki, 2002).

In the present work, we propose a statistic $`w`$, the differential coefficient of the mean absolute difference of an observed lightcurve, to make timescale analysis of shot width. The timescale spectrum of $`w`$ can reflect the width distribution of a stochastic shot process with noise. We introduce the method of $`w`$ analysis in §2. We studied the $`w`$ spectra of a sample of NS and BHC binaries and made an energy dependent analysis of the shot widths for three BHC binaries with RXTE/PCA data; the results are presented in §3. In §4 we discuss the possible implications of our results.

## 2 Method

### 2.1 Algorithm

Let $`\{x(j;\delta t)\}`$ denote the originally observed lightcurve, where $`x(j;\delta t)`$ (cts/s) is the counting rate during the time interval $`[j\delta t,(j+1)\delta t)`$ ($`j=0,1,2,\mathrm{}`$), and $`\delta t`$ is the time resolution. To study the variability on a larger timescale $`\mathrm{\Delta }t=M\delta t`$, we need to construct a lightcurve with the time step $`\mathrm{\Delta }t`$ from the original time series by combining its M successive bins as

$$x_m(i;\mathrm{\Delta }t)=\frac{1}{M}\underset{j=iM+m}{\overset{(i+1)M+m1}{}}x(j;\delta t)\text{cts/s},$$
(1)
where $`i=0,1,\mathrm{},N1`$ and $`m`$ ($`[0,M1]`$) is the discrete phase of the time series. As the lightcurve $`\{x_m(i;\mathrm{\Delta }t)\}`$ does not include any information about the variation on timescales shorter than $`\mathrm{\Delta }t`$, it can be used to study the variability over the region of timescales $`\mathrm{\Delta }t`$. A mean absolute difference of the observed lightcurve at the timescale $`\mathrm{\Delta }t`$ can be defined as

$$W(\mathrm{\Delta }t)=\frac{1}{M}\underset{m=0}{\overset{M1}{}}(\frac{1}{N1}\underset{i=0}{\overset{N2}{}}\frac{|x_m(i+1;\mathrm{\Delta }t)x_m(i;\mathrm{\Delta }t)|}{\mathrm{\Delta }t}).$$
(2)
The quantity $`W(\mathrm{\Delta }t)`$ includes information of variability with timescales larger than or equal to $`\mathrm{\Delta }t`$. To extract the information at the timescale $`\mathrm{\Delta }t`$, we can calculate the rate of change of $`W(\mathrm{\Delta }t)`$

$$w(\mathrm{\Delta }t)=\frac{dW(\mathrm{\Delta }t)}{d\mathrm{\Delta }t}.$$
(3)
To detect the signal in a spectrum $`w(\mathrm{\Delta }t)`$ against a noise background, we need to know the timescale spectrum $`w_{noise}(\mathrm{\Delta }t)`$ of a time series consisting of only noise. We assume two independent random variables $`X_1`$ and $`X_2`$ follow the Poisson distribution with the same mean $`\overline{X}`$ (cts), then the expectation of their absolute difference value is given by

$$E(|X_1X_2|)=2e^{2\overline{X}}\underset{n=0}{\overset{\mathrm{}}{}}nI_n(2\overline{X}),$$
(4)
where $`I_n`$ is the modified Bessel function of the first kind with the order number $`n`$. Therefore the mean absolute difference for the background can be derived as

$$W_{noise}(\mathrm{\Delta }t)=\frac{E(|X_1X_2|)/\mathrm{\Delta }t}{\mathrm{\Delta }t}=\frac{2e^{2\overline{x}\mathrm{\Delta }t}_{n=0}^{\mathrm{}}nI_n(2\overline{x}\mathrm{\Delta }t)}{\mathrm{\Delta }t^2}$$
(5)
with $`\overline{x}`$ being the average counting rate of the lightcurve, then

$$w_{noise}(\mathrm{\Delta }t)=\frac{dW_{noise}(\mathrm{\Delta }t)}{d\mathrm{\Delta }t}.$$
(6)
Finally the signal’s value of the statistic $`w`$ can be derived as

$$w_s(\mathrm{\Delta }t)=w(\mathrm{\Delta }t)w_{noise}(\mathrm{\Delta }t)$$
(7)
Calculating $`w_s(\mathrm{\Delta }t)`$ by Eq. (7) at different step $`\mathrm{\Delta }t`$ for an observed lightcurve we can get a timescale spectrum of $`w_s`$. The differential coefficients in Eqs. (3) and (6) can be calculated numerically. In practice, a lightcurve may be divided into $`L`$ segments, we can calculate $`w_{s,l}`$ for each segment and then obtain the average value

$$\overline{w}_s=\underset{l=0}{\overset{L1}{}}w_{s,l}/L$$
(8)
and its standard error

$$\sigma (\overline{w}_s)=\sqrt{\underset{l=0}{\overset{L1}{}}(w_{s,l}\overline{w}__s)^2/L(L1)}.$$
(9)
### 2.2 Simulation

We make simulation study in order to understand the timescale spectra of $`w`$. The $`w`$ spectra of two time series consisting of square shots are calculated. The shot width is 0.1 s. One time series (top panel of 1) consists of periodic shots with a period of 1 s, and in the other time series (middle panel of 1) shots are randomly distributed with shot separations following an exponential distribution with a mean of 1 s. The time resolution is 5 ms and the total time length is 1250s with each segment of 250 s. A peak around the timescale 0.1 s appears clearly in each spectrum (see the bottom panel of Fig. 1). The other characteristic timescale of the periodic series, 0.5 s, is also shown in its $`w`$ spectrum (the dashed line in Fig. 1) with a smaller amplitude. At small timescales, the random shot series has stronger variability than the periodic one, indicated by the larger amplitude in the $`w`$ spectrum.

The $`w_s`$ spectra of a set of simulated lightcurves consisting of both white noise and signal shots randomly occurring with different widths and amplitudes are also calculated. The shot width is randomly sampled from a uniform distribution between $`[\tau _1,\tau _2]`$ and its profile has either square or Gaussian form. The peak amplitude of a shot is drawn from a uniform distribution between 0 and 1000 cts/s. The shot occurring time is randomly distributed following an exponential distribution with a mean of 1 s. The time resolution is 10 ms and the total duration is 5000 s with each segment of 500 s. A background of 1000 cts/s is added. Poisson counting statistics are also included in the simulations. The results, shown in Fig. 2 (panel (a) for square shots and (b) for Gaussian shots), indicate that we can use the position of short timescale cut-off in the $`w_s`$ spectrum to measure the shortest width of random shots in an observed lightcurve with noise. For comparison, we also calculate the Fourier power spectrum density (PSD) for each shot series. The results are plotted in panel (c) and (d) of Fig. 2 for square and Gaussian shots, respectively. It is clear by comparing the PSD with the $`w_s`$ spectra in Fig. 2 that the $`w`$ spectrum can be used to detect the shot width more conveniently and accurately.

To measure the lower cut-off in a $`w_s`$ spectrum, we fit a polynomial to the front edge and obtain the root value. The error of the cut-off timescale is evaluated from the fitting error with 90% confident level. In this way, the lower cut-off timescales for the above simulations are $`0.0614\pm 0.0006`$, $`0.0743\pm 0.0008`$ and $`0.0893\pm 0.0009`$ (s) respectively for different lightcurves in the left panel of fig. 2, and $`0.0794\pm 0.0007`$, $`0.0923\pm 0.0009`$ and $`0.1092\pm 0.0010`$ (s) for the right panel of fig. 2. The measured lower cut-off timescale in the $`w`$ spectrum is very close to the input lower scale of shot widths, which are 0.05, 0.075 and 0.1 (s) respectively for both square and gaussian shots. There are larger differences for Gaussian shots because of the Gaussian FWHM may not equal to $`w`$ spectrum’s shot width. We therefore conclude that we can detect the lower cut-off timescale of X-ray shots from lightcurves with the $`w`$ spectrum.

To address the question whether $`w`$ spectrum may also detect the frequency breaks often seen in the PSD of the lightcurves many astrophysical sources, we simulate two lightcurves with different PSDs (using the inverse FFT method of Timmer & König, 1995). The PSD consists of power-law forms at different segments: power-law index $`\tau =0`$ when $`f<3`$ Hz, $`\tau =1`$ when $`f3`$ and $`<15`$ Hz, and $`\tau =1.5`$ or $`2`$ when $`f15`$ Hz (left panel of Fig. 3). The $`w`$ spectra of these two lightcurves are presented in the right panel of Fig. 3. We can see that the frequency breaks in the PSD do not correspond to timescale cut-offs, but manifest as different slopes in the $`w`$ spectra.

In sum, the $`w`$ spectrum is only sensitive to the width distribution of a shot series, but insensitive to other shot parameters, such as the shot profile, shot separation distribution and shot amplitude. In particular, the $`w`$ spectrum is especially powerful in detecting the smallest shot width of a shot series.

## 3 Data Analysis

We apply the timescale analysis method to observations of a number of X-ray binaries with the Proportional Counter Array (PCA) on board Rossi X-ray Timing Explorer (RXTE). The data are screened when five PCUs are all on and the original time bin size $`\delta t`$ of the lightcurve is binned to 2 ms or 4 ms according to the data mode, with each segment of 50000 bins (100 s or 200 s). Because this analysis method can tolerate high level of noise, we did not avoid PCA channels 0-7, which have higher level of instrumental background. The observation IDs and other details of used data are listed in Table 1.

The top panel of Figure 4 shows the normalized $`w_s`$ spectra of five accretion BHCs, GRS 1915+105, GRO J1655-40, GX 339-4, XTE J1550-564, and Cyg X-1. All spectra show a cut-off at a timescale about 0.05 s. But for five selected NS binaries with weak magnetic field, i.e., 4U 1705-440, GS 1826-24, 4U 0614+091, 4U 1608-522, and Cyg X-2, the $`w_s`$ spectra have no obvious cut-off above $`10^2`$ s (see the bottom panel of Fig. 4).

The energy dependent analysis is carried out for three BHC binaries: Cyg X-1, XTE J1550-564, and GRO J1655-40. All sources are at their hard state. The choice of energy bands for each observation reflects the different data mode of each observation. From Figure 5, one can see a common trend: below about 10-20 keV the cut-off timescale $`\mathrm{\Delta }t_c`$ decreases with increasing energy, and above about 10-20 keV $`\mathrm{\Delta }t_c`$ increases with increasing energy. An emission feature at iron K$`\alpha `$ line appears clearly in the energy spectra of Cyg X-1 and XTE J1550-564, and a broad peak around $``$ 6 keV can also be seen in their $`\mathrm{\Delta }t_c`$ spectra. For GRO J1655-40, the Fe reflection component is very weak in the energy spectrum of this observation, and this broad peak in its $`\mathrm{\Delta }t_c`$ spectra is also not seen.

## 4 Discussion

Extracting physical information from lightcurves for multi-timescale and aperiodic processes is difficult. However, it is important to study the underlying physics in high-energy processes around compact objects. The timescale analysis technique has the freedom of selecting or designing a proper statistic for a particular purpose. Our results of simulation and data analysis for accreting X-ray sources show that the proposed statistic $`w`$ is useful in timescale analysis of shot widths. From their $`w`$ spectra (the top panel of Fig. 4), we find that the BHCs have shortest shot width $``$ 0.05 s, consistent with that shots in BHCs have characteristic timescale $``$ 0.1 s, which was previously derived by a timescale analysis of power density (Li & Muraki, 2002). The fact that no short timescale cut-off appears in the $`w`$ spectra of NS binaries (the bottom panel of Fig. 4) is also consistent with what indicated by their timescale spectra of power density (Li & Muraki, 2002), and by Fourier spectral analysis (Sunyaev & Revnivtsev, 2000) of these systems in the hard state. The timescale analysis technique of shot width proposed in this work, like the time domain power spectrum analysis, can help us distinguish BHCs from NSs. Compared with the timescale analysis of power density spectrum, the statistic $`w`$ can measure more directly and sensitively the lower cut-off of shot timescale, with which we can study the energy dependence of the cut-off timescale quantitatively.

The variation timescale is usually taken as an indicator of the spatial dimensions of the underlying physical process. The remarkable result from our timescale analysis with the statistic $`w`$ is the energy dependent timescale cut-off of BHCs (Fig. 5) is hard to explain by a simple correspondence between the variation timescale and the radius where shots take place. The monotonic increase of the timescale cut-off with energy in the region of $`E10`$-20 keV is consistent with previous time domain analysis of variability in Cyg X-1 and other BHCs (e.g., Feng, Li, & Chen, 1999; Li, Feng, & Chen, 1999; Maccarone, Coppi, & Poutanen, 2000; Negoro, Kitamoto, & Mineshige, 2001), which is expected from the common understanding that hard X-rays come from inverse Compton scattering of low-energy seed photons by high-energy electrons in a hot corona (e.g., Eardley, Lightman, & Shapiro, 1975; Poutanen & Svensson, 1996; Dove et al., 1997) or hot advection-dominated accretion flow (ADAF) (e.g., Esin et al., 1998), where a shot of low energy seed photons are hardened and broadened by the Comptonization process.

A rough estimate from Figure 5 can be made to evaluate the corona size $`R_c`$ of a BHC. We assume a typical temperature of $``$100 keV for the hot corona at the hard spectral state. The mean free path of X-rays in corona can be considered as the corona’s spatial dimension $`R_c`$, because its optical depth is $`<`$ 1 from the energy spectral fitting. We assume the turning point energy $`E_0`$ at $`w`$ spectrum as the energy of input seed shot photons, and after $`k`$ times of inverse Compton scattering, the mean photon energy is increased to $`E_k`$, the timescale $`\mathrm{\Delta }t_{ck}`$ of shots at energy $`E_k`$ can then be located in the $`\mathrm{\Delta }t_c`$ spectrum. By a simple relationship of $`\mathrm{\Delta }t_{ck}=\mathrm{\Delta }t_{c0}+kR_c/c`$, where $`\mathrm{\Delta }t_{c0}`$ is the turning point timescale in the $`\mathrm{\Delta }t_c`$ spectrum and $`c`$ is the light speed, $`R_c`$ is then estimated. From Figure 5, the approximately derived $`R_c`$ are in order of 30, 140 and 70 $`R_s`$ (Schwarzschild radius) for Cyg X-1, XTE J1550-564 and GRO J1655-40, respectively, assuming these BHC masses are 10, 10 and 7 $`M_{}`$ (Quirrenbach et al., 2002).

It is obvious that the general trend of timescale cut-off vs. energy in the energy region of $`E10`$-20 keV is opposite to that in higher energy region. We first ignore the broad peak around 6 keV in the $`\mathrm{\Delta }t_c`$ spectrum and focus on the decreasing trend of $`\mathrm{\Delta }t_c`$ with increasing energy at $`E10`$-20 keV. One probability is that shots of different energies are produced in different regions of the disk: the high energy shots are located at inner disk and low energy ones are from outer region. The turning point at $``$ 10-20 keV may indicate that the maximum energy of the initial shot photons, which are most probably produced at the turbulent region of the cold disk joining the hot corona, can be as high as up to $``$ 10-20 keV (Manmoto et al., 1996; Li, Feng, & Chen, 1999).

Another possibility to explain the $`\mathrm{\Delta }t_c`$ spectrum at energy below the turning point is to consider the shots with energies below the turning point are produced by down-Comptonization, i.e., original shots with energy at the turning point are scattered by cold electron cloud covering the accretion disk (Zhang et al., 2000). From the top panel of figure 5 (Cyg X-1), we can see that from 11.6 keV (the turning point energy) to 2.4 keV (the lowest energy) the timescale is increased by 0.037 s. If the mean kinetic energy of cold electrons is much less than the turning point energy, about 170 times of scattering are needed for the X-ray energy increased from the turning point energy to 2.4 keV. Assuming the relationship between the increase of timescale $`\mathrm{\Delta }t`$, the number of scattering $`N`$ and mean free path $`l`$ under the optically thick condition, is described by $`\mathrm{\Delta }t=lN/c`$, the derived mean free path $`l`$ is $`6.5\times 10^6`$ cm (2.2 $`R_s`$ of Cyg X-1).

Under the above assumption that an optically thick cold electron cloud covers the initial shots (with energy at the turning point of the $`\mathrm{\Delta }t_c`$ spectrum) located near the innermost region of the disk, we can evaluate the geometry of the cold electron cloud, from the shape of the the $`\mathrm{\Delta }t_c`$ spectrum (top panel of fig. 5 for Cyg X-1). The distance between the location where the down-scattered photons escape to the seed shot photon region can be described as $`D=\tau l`$, where $`\tau `$ is the optical depth and equal to $`\sqrt{N}`$ in an optically thick system. We find that $`\tau `$ increases and $`l`$ decreases when the photon energy decreases, and consequently $`D30R_s`$ of $`D`$ is obtained, almost independent of the escaped photon energy. Similarly we also applied the above procedure to the data on XTE J1550-564 and GRO J1655-40. In Figure 6 we show the relationship between $`D`$ and the escaped photon energy.

Our calculation for Cyg X-1 also indicates that if the mean temperature of electrons is larger than 0.03 keV, the above down-Comptonization process cannot decrease the average photon energy from the turning point to below 2.4 keV. For example, if the electron temperature is 0.03 keV, after about 1000 Compton scattering the average photon energy saturates at 2.4 keV, down from the seed shot photon energy of 11.6 keV. So if the above assumed down-Comptonization scenario is true for the $`\mathrm{\Delta }t_c`$ spectrum below the turning point, we predict that a sharp increase of timescale in the $`\mathrm{\Delta }t_c`$ spectrum below the saturation energy, which is uniquely related to the temperature of the cold electron cloud.

Finally we examine the broad peak in the $`\mathrm{\Delta }t_c`$ spectrum around $``$ 6 keV for Cyg X-1 and XTE J1550-564. It is most probably produced by the disk reflection component: reprocessed iron K$`\alpha `$ fluorescent photons in the cold disk irradiated by hard X-rays with energy above the $``$7keV iron K-edge will further broaden the shots. This interpretation is also consistent with the explanation of time delays and the behavior of the iron line (Gilfanov, Churazov & Revnivtsev, 2000; Poutanen, 2002) with frequency-resolved energy spectra (Revnivtsev, Gilfanov & Churazov, 1999). Besides producing the fluorescence iron K$`\alpha `$ emission, the refection also has a strong component at energies above 10 keV, which may contribute to the increase of shot width above the tuning point in addition to the up-Comptonization process we discussed above.

We are very grateful to the expert and insightful comments by the anonymous referees; panels (c) and (d) in figure 2 and the whole figure 3 are added in responding to the comments, as well as several clarifications in the manuscript. This study is supported in part by the Special Funds for Major State Basic Research Projects and by the National Natural Science Foundation of China under project No. 10233030. SNZ also acknowledges supports by NASA’s Marshall Space Flight Center and through NASA’s Long Term Space Astrophysics Program.


text:
# The Effect of aberration on polarization position angle of pulsars

## 1 Introduction

Pulsars exhibit broad-band coherent radio emission ranging from tens of MHz to tens of GHz and is highly polarized. It is usually observed that pulse widths at higher frequencies are smaller than that at lower frequencies (see Mitra & Rankin 2002 for a recent review). This phenomenon is termed as radius to frequency mapping (RFM) according to which emission at progressively higher frequencies arise closer to the stellar surface. The polarization position angle (PPA) of the linear polarization for several pulsars is seen to execute a smooth ‘S-shaped’ curve which has been interpreted as evidence for highly beamed emission arising from dipolar magnetic field lines as per the rotating vector model (RVM) proposed by Radhakrishnan & Cooke (1969). Fitting the RVM model to the PPA traverse can be used to find the geometrical angles of the neutron star namely $`\alpha `$ the angle between the rotation axis and the dipolar magnetic axis and $`\beta `$ the angle between the rotation axis and the observers line-of-sight. RVM and RFM together has been extensively used to find radio emission heights in pulsars and the shape of the pulsar beam. Canonical values for emission heights are found to be 50 times the stellar radius and the emission beam is thought to be organized in the shape of nested cones of emission, with a central core emission.

Whilst RFM is more commonly observed in pulsars (barring a few exceptions) there are significant deviations from the RVM observed in PPA traverses in pulsars. In some cases across the pulse the PPA traverse is seen to have 90 jumps and/or two orthogonal PPA tracks. Careful study has shown that each of the PPA track are in agreement with the RVM (Gil & Lyne 1994). The other kind of deviation is in the form of wiggles and kinks (non-orthogonal departures) observed in the PPA traverse. In this paper we will concentrate on the latter kind of deviation observed and put forward a possible scenario which can explain such deviations.

## 2 The BCW model and the RVM curve

The RVM model has been further advanced by Blaskiewicz et al. (1991, hereafter BCW) where they included first order special relativistic effects into account. The model predicts that due to aberration and retardation (A/R) effects the PPA traverse depends on the emission height $`r_{em}`$ and the pulsar period $`P`$ (in sec) as,

$$\psi =\psi _0+\mathrm{arctan}\left(\frac{\mathrm{sin}\alpha \mathrm{sin}(\varphi \varphi _0)(6\pi /P)(r_{em}/v_c)\mathrm{sin}\xi }{\mathrm{sin}\xi \mathrm{cos}\alpha +\mathrm{sin}\alpha \mathrm{cos}\xi \mathrm{cos}(\varphi \varphi _0)}\right).$$
(1)
Here $`\psi _0`$ and $`\varphi _0`$ are arbitrary position angle and longitude phase offsets and $`\varphi `$ is the longitude of the pulse. The angle $`\xi =\alpha +\beta `$ and $`v_c`$ is the velocity of light. Note that the above equation reduces to the RVM for $`r_{em}=0`$. It is immediately obvious from the above equation that changes in $`r_{em}`$ across the pulse will result in distortion of the PPA traverse. Also if $`r_{em}`$ is fixed, then the distortions will be larger for faster pulsars, as expected from the A/R effects.

However $`r_{em}`$ across the pulse is an unknown. It is thought that the central core emission arises much below the conal emission. In Figure 1 we show an example of how the PPA curve will be distorted if emission from different parts of the pulse originated from two distinctly different heights. To refine the above concept of emission arising from different heights we put forward the following model,

(a) We note a pulse profile can be decomposed as a set of Gaussians (Kramer 1994). Further each Gaussian can be attributed to either cone or core emission. These core and cone can arise from different heights.

(b) For every ith longitude in the pulse, the resultant emission height $`r_{em}(i)`$ is an intensity (I) weighted quantity given as,

$$r_{em}(i)=\frac{_{c=1}^N(I(i,c)\times r(c))}{_{c=1}^NI(i,c)}$$
(2)
Here the index $`c`$ refers to the Gaussian component and $`r(c)`$ corresponds to the emission height attributed to that component (which correspond to height of the cone or core component). In Figure 1 we show an example on how the PPA traverse behaves when equation 2 is applied. It should be noted that the PPA traverse strongly depends on the parameters $`\alpha `$, $`\beta `$ and $`r_{em}`$ which are highly correlated to each other. Distortions of the PPA traverse due to height dependent emission has also been investigated by Xu. et. al. (1997).

## 3 Application of BCW model to pulsars

Most of the time the distortions in the PPA traverse seem to occur below the central core component in pulsars. While it is seen in a few normal pulsars (NP), it is more commonly seen in millisecond pulsars (MSP). Also in MSP’s the PPA traverse is seen to be flat.

PSR J1022+1001 a millisecond pulsar with a period of 16 msec at 408 MHz is a classical example where the PPA traverse shows distortions below the central core component, as shown in Figure 2 (left plot). The kinky nature of the PPA traverse seen below the central core component. This pulsar can be adequately described by fitting three Gaussians, two outer ones comprising of one cone and the inner core. On the right hand plot of Figure 2 we show simulation of the PPA traverse using equation 2 assuming that the central core emission arises lower than the outer cone emission. This way we are able to generate the kinky feature in the PPA traverse. However we emphasize here that given the simple model that is being used we are only able to match the gross properties of the PPA traverse and do not claim to find the exact solution for $`\alpha `$, $`\beta `$ and $`r_{em}`$. For this a more sophisticated model needs to be considered.

## 4 Conclusion

Using a simple model we have been able to produce the kinkiness in the PPA traverse towards the central region of the pulse profiles which is quite often observed. As a consequence we are able to conclude the following:

(a) The central (or core) emission arises lower than the outer conal emission confirming a suggestion put forward by several people (see for example Rankin 1993).

(b) The PPA traverse results from emission arising over a broad range of heights as has been recently suggested by Mitra & Li (2003).

Recently the central kinky feature seen in PSR J1022+1001 was interpreted as evidence for ‘return currents’ in pulsars by Ramachandran & Kramer (2003). Here we show that by considering height dependent emission regions similar kinky features can be obtained and it is not necessary to invoke return current effects to explain this feature.

## References

Blaskiewicz, M., Cordes, J. M., & Wasserman, I. 1991, ApJ, 370, 643

Gil, J. A. & Lyne, A. G. 1995, MNRAS, 276, L55

Gould, D. M. & Lyne, A. G. 1998, MNRAS, 301, 235

Mitra, D. & Li, X. H., 2003, A&A in press

Mitra, D. & Rankin, J. M. 2002, ApJ, 577, 322

Kramer, M. 1994, A&AS, 107, 527

Radhakrishnan, V. & Cooke, D. J. 1969, ApJL, 3, 225

Ramachandran, R. & Kramer, M., 2003, A&A, 407, 1085

Rankin, J. M. 1993, ApJ, 405, 285

Xu, R. X., Qiao, G. J. & Han, J. L., 1997, A&A, 323, 395


text:
# Conditionally exactly solvable potential and dual transformation in quantum mechanics

## Abstract

We comment that the conditionally exactly solvable potential of Dutt et al (1995 J. Phys. A: Math. Gen. 28 L107) and the exactly solvable potential from which it is derived form a dual system.

PACS numbers: 03.65.Ca, 03.65.Ge, 02.90.+p

Keywords: Quantum mechanics; Conditionally exactly solvable potentials; Dual systems

Conditionally exactly solvable (CES) potentials have received considerable attention in recent times . The main feature of such potentials is that one or more coupling constants in them are fixed to a specific value. There have been instances of CES potentials running into inconsistencies with threshold boundary conditions , but there do exist some which possess valid asymptotic behaviour. One such acceptable class is the one proposed by Dutt et al  sometime ago and which reads 111There is also another class of CES potentials proposed in  but by a redefinition of parameters  it can be made equivalent to (1) and so is not considered here.

$$V(y)=\frac{A}{1+e^{2y}}\frac{B}{(1+e^{2y})^{1/2}}\frac{3}{4(1+e^{2y})^2}$$
(1)
where $`y(\mathrm{},\mathrm{})`$ and $`A`$, $`B`$ are some real parameters defining the shape of the potential. Note the presence of the fixed numerical value $`\frac{3}{4}`$ for one of the coupling constants in $`V(y)`$ that provides its identification as a CES. Potential (1) has for its associated eigenfunctions

$$\psi _n(y)=z^{\frac{1}{4}}(z1)^{\left(\frac{c}{2}\frac{B}{4c}\right)}(z+1)^{\left(\frac{c}{2}\frac{B}{4c}\right)}P_n^{(\frac{B}{2c}c,\frac{B}{2c}c)}(z),$$
(2)
where $`z=1+e^{2y}`$ and $`c`$ is related to the energy eigenvalue $`ϵ_n`$ as $`c=n+\frac{1}{2}+\sqrt{ϵ_n}`$. Actually $`\sqrt{ϵ_n}`$ satisfies a complicated cubic equation but it has been observed  that only one of its roots is compatible with the normalizability condition.

Some remarks are in order concerning the derivation of the eigenfunctions (2) and the energy eigenvalue equation. The Schrödinger equation is subjected to a coordinate transformation and the transformation function is chosen in such a manner that corresponding to an exactly solvable (ES) potential one has a new analytically solvable one. It turns out that for a half-line-full-line mapping function, one can generate the CES potential for some known shape-invariant potential as an input. The available energy eigenvalues and eigenfunctions of the latter then furnish the corresponding ones for the former.

The purpose of this comment is to establish that the CES potential (1) and its accompanying shape-invariant ES potential are actually dual partners in the sense that the corresponding time-independent Schrödinger equations are mapped to each other under appropriate space transformations called the dual transformations . The latter are known to relate some apparently unconnected problems both in classical and quantum mechanics. The one-dimensional harmonic oscillator and the Coulomb problem , the latter and the isotropic oscillator , the Pöschl-Teller and infinite potential well problems  are some examples of dual systems (DS).

In the present context, let us write down the following set of DS given by the pair of Schrödinger equations (with $`\mathrm{}=2m=1`$)

$`\left[{\displaystyle \frac{d^2}{dx^2}}+\lambda \left({\displaystyle \frac{dy}{dx}}\right)^2+\nu {\displaystyle \frac{dy}{dx}}\right]\psi =\mu \psi `$ (3)
$`\left[{\displaystyle \frac{d^2}{dy^2}}{\displaystyle \frac{1}{2}}\{x,y\}\mu \left({\displaystyle \frac{dx}{dy}}\right)^2+\nu {\displaystyle \frac{dx}{dy}}\right]\varphi =\lambda \varphi `$ (4)
where $`\mu `$ and $`\lambda `$ are the energies, $`\nu `$ is a constant, $`\{x,y\}`$ is the Schwarzian derivative, which can be written as

$$\{x,y\}=\frac{1}{y^2}\left[\frac{d}{dx}\left(\frac{y^{\prime \prime }}{y^{}}\right)\frac{1}{2}\left(\frac{y^{\prime \prime }}{y^{}}\right)^2\right]$$
(5)
the primes denoting derivatives with respect to the variable $`x`$. The wave functions $`\psi `$ and $`\varphi `$ are related in the manner

$$\psi =\left(\frac{dx}{dy}\right)^{1/2}\varphi .$$
(6)
In DS as above, the energy and the coupling constant $`\mu `$ and $`\lambda `$ exchange roles. Further DS are meaningful when the potential and its partner are expressible as

$`W(x)=\lambda \left({\displaystyle \frac{dy}{dx}}\right)^2+\nu {\displaystyle \frac{dy}{dx}}`$ (7)
$`U(y)=\mu \left({\displaystyle \frac{dx}{dy}}\right)^2+\nu {\displaystyle \frac{dx}{dy}}{\displaystyle \frac{1}{2}}\{x,y\}.`$ (8)
Generalizations to include integral or even fractional powers in the derivatives in (7) and (8) are straightforward.

We now set $`dy/dx=\mathrm{coth}x`$, i.e., $`y=\mathrm{log}\mathrm{sinh}x`$: in other words, for $`y(\mathrm{},\mathrm{})`$ the variable $`x(0,\mathrm{})`$. It implies from (7) that the Schrödinger equation (3) has the potential $`W(x)=\stackrel{~}{W}(x)+\alpha (\alpha 1)`$, where

$$\stackrel{~}{W}(x)=2\beta \mathrm{coth}x+\alpha (\alpha 1)cosech^2xx(0,\mathrm{})$$
(9)
which is ES (for $`\beta >\alpha ^2`$, $`\alpha >0`$) and shape invariant as well. The energy eigenvalues are given by

$$E_n=\left(\frac{\beta }{\alpha +n}\right)^2(\alpha +n)^2n=0,1,2,\mathrm{}.$$
(10)
Indeed the correspondence with (3) is provided by the following identifications

$`\lambda `$ $`=`$ $`\alpha (\alpha 1)`$
$`\nu `$ $`=`$ $`2\beta `$ (11)
$`\mu `$ $`=`$ $`E_n+\alpha (\alpha 1).`$
We next enquire into the dual potential $`U(y)`$. It is easy to work out $`\{x,y\}`$ from (5) as

$$\{x,y\}=\left(sech^2x\mathrm{tanh}^2x+sech^2x\frac{1}{2}sech^4x\right).$$
(12)
Hence we find from (8) and (4) that $`U(y)=\stackrel{~}{U}(y)+\frac{1}{4}`$, where

$`\stackrel{~}{U}(y)=\left[\left({\displaystyle \frac{1}{2}}\mu \right)\mathrm{tanh}^2x+\nu \mathrm{tanh}x{\displaystyle \frac{3}{4}}\mathrm{tanh}^4x\right]_{x=\mathrm{sinh}^1(e^y)}`$ (13)
$`ϵ_n=\alpha (\alpha 1)+{\displaystyle \frac{1}{4}}.`$ (14)
The relation between the energies $`ϵ_n`$ and $`E_n`$ turns out to be

$$ϵ_n+E_n=\mu +\frac{1}{4}.$$
(15)
On elimination of the parameter $`\alpha `$, from (10), (11) and (14) we then get a cubic equation in $`\sqrt{ϵ_n}`$, similar to that given in . Further, in terms of the variable $`y(\mathrm{},\mathrm{})`$, $`\stackrel{~}{U}(y)`$ in (13) translates to

$$\stackrel{~}{U}(y)=\frac{\frac{1}{2}\mu }{1+e^{2y}}+\frac{\nu }{(1+e^{2y})^{1/2}}\frac{3}{4(1+e^{2y})^2}y(\mathrm{},\mathrm{})$$
(16)
which is identical to the potential in (1) for $`\mu =\frac{1}{2}A`$ and $`\nu =B`$. Finally, the wave functions for $`\stackrel{~}{U}(y)`$ can be obtained from (6) and yield the same form as in (2).

To conclude we have demonstrated that the CES potential of Dutt et al and its ES partner form a DS.

One of us (BB) gratefully acknowledges the support of the National Fund for Scientific Research (FNRS), Belgium, and the warm hospitality at PNTPM, Université Libre de Bruxelles, where this work was carried out. CQ is a Research Director of the National Fund for Scientific Research (FNRS), Belgium.