text:
# Probing Anderson localization of light via decay rate statistics

## Abstract

We have studied the distribution of resonant widths $`P\left(\mathrm{\Gamma }\right)`$ in one-, two- and three dimensional multiple light scattering systems. $`P\left(\mathrm{\Gamma }\right)`$ should follow a universal power law $`P\left(\mathrm{\Gamma }\right)\mathrm{\Gamma }^1`$ in the localized regime as confirmed by extensive numerical calculations. This behavior can be interpreted as an unambiguous signature of exponential Anderson localization of light in open systems.

The research on Anderson localization of light has been of great interest cargese  since it was originally proposed as the optical counterpart of electronic localization john . Localization, as proposed by Anderson, is defined as an inhibition of wave diffusion in infinite disordered media due to interference of multiple scattered waves anderson . A much stronger definition is that the eigenfunctions in an infinite disordered medium are characterized by an exponential decay in space, $`\left|\psi \left(𝐫\right)\right|\mathrm{exp}\left(\left|𝐫𝐫^{}\right|/\xi \right)`$, where $`\xi `$ is the localization length. In finite, open media, waves can “leak” through the sample boundaries. Anderson localization must thus relate to manifestations of leakage in observables quantities. For optical systems, they are typically the emerging intensity, the total transmission or the coherent backscattering cone. The observation of an exponential scaling of transmission azgarcia ; nature , as well as the rounding of the backscattering cone cone , may not have definitively established localization since absorption could be responsible for these same effects.

There are several criteria to determine the onset of the localized regime. The Ioffe-Regel criterion states that, in three dimensions (3D), localization occurs for $`k\mathrm{}1`$ (with $`k`$ the light wavenumber inside the medium and $`\mathrm{}`$ the mean free path). Another approach to define electromagnetic localization is based on the variance of fluctuations of transmission, even in the presence of absorption natureaz . In open systems, the “eigenstates” are resonances with a finite energy width $`\mathrm{\Gamma }`$ (or, equivalently, with a finite lifetime $`t1/\mathrm{\Gamma }`$) due to leakage. The Thouless criterion establishes that localization can be set to occur when the typical time that an excitation needs to propagate through the entire system of size $`R`$, $`t_T1/\mathrm{\Gamma }_TR^2/D`$ (Thouless time), exceeds the maximal time scale of the system, $`t_H1/\mathrm{\Gamma }_H1/\mathrm{\Delta }E`$ (Heisenberg time) thouless . Here $`D`$ is the diffusion constant and $`\mathrm{\Delta }E`$ the mean level spacing.

This Thouless criterion applies to the average leakage width. Hence it is reasonable to assume that the statistical properties of resonance widths are strongly affected by localization. The aim of the present paper is to investigate how localization manifests itself in the distribution of resonance widths $`P(\mathrm{\Gamma })`$ in multiple light scattering in open systems. We will show that $`P(\mathrm{\Gamma })`$ exhibits the universal power law $`P(\mathrm{\Gamma })\mathrm{\Gamma }^1`$ in 1D, 2D and 3D optical disordered systems, thereby generalizing recent theoretical titovfyo  and numerical terraguarnieri  studies in 1D models of mesoscopic transport. We assert that the algebraic decay $`P(\mathrm{\Gamma })\mathrm{\Gamma }^1`$ represents a universal property of Anderson localization of light.

Although the statistical properties of resonance widths in open systems have been extensively studied over the last years, in particular for chaotic/ballistic systems distgamma ; fyosomm ; sommfyotit , their behavior for systems exhibiting localization has received considerably less attention. As argued by Casati et al., $`P(\mathrm{\Gamma })`$ should follow a power law $`P(\mathrm{\Gamma })\mathrm{\Gamma }^1`$ in localized, classically chaotic systems dima . $`P(\mathrm{\Gamma })`$ was analytically obtained for 1D disordered systems, showing a slightly different power law $`P(\mathrm{\Gamma })\mathrm{\Gamma }^{1.25}`$ titovfyo . This prediction was corroborated later by numerical calculations in 1D and quasi-1D tight binding models terraguarnieri . The $`P(\mathrm{\Gamma })\mathrm{\Gamma }^1`$ behavior was also reported in 1D and 3D strongly driven atomic Rydberg states in the context of dynamical localization andreas . Exactly at the Anderson transition distcritical  and in the diffusive regime distdifusive , $`P(\mathrm{\Gamma })`$ was shown to follow a power law with a power different from $`1`$. Concerning the study of $`P(\mathrm{\Gamma })`$ for optical systems, the only work on the subject is, to the best of our knowledge, due to Patra patra , who mainly focused on the small $`\mathrm{\Gamma }`$ regime and its application to random lasers. For small $`\mathrm{\Gamma }`$, it is known that $`P(\mathrm{\Gamma })`$ is different from a power law, for both the diffusive fyosomm ; patra  and the localized regimes sommfyotit ; patra . It should be emphasized that the power law decay of $`P(\mathrm{\Gamma })`$ is expected to occur only for intermediate values of $`\mathrm{\Gamma }`$, i.e., for not too large $`\mathrm{\Gamma }`$, where the resonances are strongly coupled to the continuum and where $`P(\mathrm{\Gamma })`$ decays faster than algebraically titovfyo ; terraguarnieri , and for $`\mathrm{\Gamma }\mathrm{\Gamma }`$.

We will present a simple physical argument, inspired by Refs. dima ; titovfyo ; terraguarnieri , to explain the universal $`P(\mathrm{\Gamma })\mathrm{\Gamma }^1`$ behavior for the localized regime, i.e., independent on the dimensionality of the system. Due to the opening of the system, exponentially localized eigenstates of the corresponding closed system (linear size $`R`$) acquire a finite frequency width $`\mathrm{\Gamma }^{}`$, $`\mathrm{\Gamma }^{}e^{r^{}/\xi }`$, with $`r^{}`$ the distance to the boundaries. Near the system edges, the leakage is strong and the resonances are large compared to $`\mathrm{\Gamma }_T`$. On the other hand, far from the boundaries the leakage is small and the typical $`\mathrm{\Gamma }`$ in this region is much smaller than $`\mathrm{\Gamma }_T`$. Assuming that the resonances are uniformly distributed in space, it follows that the (integrated) probability of finding a resonance width $`\mathrm{\Gamma }`$ smaller than $`\mathrm{\Gamma }^{}`$, $`P_{int}\left(\mathrm{\Gamma }<\mathrm{\Gamma }^{}\right)`$, is equal to the probability of finding a resonance situated at a distance $`r`$ from the boundaries larger than $`r^{}`$, $`P\left(r>r^{}\right)`$, i.e., $`P_{int}\left(\mathrm{\Gamma }<\mathrm{\Gamma }^{}\right)=P\left(r>r^{}\right)`$. Since $`P\left(r>r^{}\right)\mu _D\left(Rr^{}\right)/\mu _D\left(R\right)`$ with $`\mu _D`$ the D-dimensional volume, we conclude that the probability density is

$`P\left(\mathrm{\Gamma }^{}\right)`$ $`=`$ $`{\displaystyle \frac{dr^{}}{d\mathrm{\Gamma }^{}}}{\displaystyle \frac{d}{dr^{}}}\left[P(r>r^{})\right]`$ (1)
$``$ $`{\displaystyle \frac{\xi }{\mathrm{\Gamma }^{}}}{\displaystyle \frac{d}{dr^{}}}\left[{\displaystyle \frac{\mu _D\left(Rr^{}\right)}{\mu _D\left(R\right)}}\right].`$
The purely geometrical factor $`d/dr^{}\left[\mu _D\left(Rr^{}\right)/\mu _D\left(R\right)\right]`$ depends on the dimensionality of the system but does not affect the exponent in $`\mathrm{\Gamma }`$.

To test the validity of (1) in the context of Anderson localization of light, we will consider scalar wave propagation in disordered media using the method introduced by Rusek and Orlowski ro ; ro00 . This approach is based on the analysis of the Green matrices spectra, which describe light scattering from randomly distributed pointlike dipoles (i.e., particles much smaller than the wavelength of light). For an incident plane wave $`\psi _0\left(𝐫\right)`$ in a system of $`N`$ identical dipoles with scattering matrix $`t`$, the field acting in the dipole at $`𝐫_i`$ is given by ro ; ro00 :

$$\psi \left(𝐫_i\right)=\psi _0\left(𝐫_i\right)+t\underset{ji}{\overset{N}{}}G\left(𝐫_{ij}\right)\psi \left(𝐫_j\right).$$
(2)
The complex-valued $`N\times N`$ matrix $`G\left(𝐫_{ij}\right)`$ describes light propagation of the wave scattered by the dipole at $`𝐫_i`$ to the dipole at $`𝐫_j`$. Since the eigenvalues $`\lambda _M`$ of $`𝐌𝐈t𝐆`$ and $`\lambda _G`$ of $`𝐆`$ are related by $`\lambda _M=1t\lambda _G`$, and $`t`$ depends on frequency $`\omega `$ via the scattering phase shift $`\delta \left(\omega \right)`$ green , an eigenvalue $`\lambda _G`$ with $`\text{Re}\lambda _G=1`$ will permit an appropriate choice for $`\delta \left(\omega \right)`$ in order to obtain an eigenvalue $`\lambda _M=0`$, which would correspond to a genuinely localized state somewhere inside the random medium ro . Assuming a Breit-Wigner model for the scatterers, (with one resonance of width $`\mathrm{\Gamma }_0`$ at the position $`\omega _0`$) and for which $`\delta \left(\omega \right)`$ has a simple form, it is possible to obtain, in a good approximation, the resonance widths $`\mathrm{\Gamma }`$ via $`\lambda _G`$, $`\mathrm{\Gamma }/\mathrm{\Gamma }_01+\text{Re}\lambda _G`$ ro00 . In the following, we will numerically diagonalize $`𝐆`$ in 1D, 2D and 3D and calculate the distribution of resonance widths $`P\left(\mathrm{\Gamma }\right)`$ using the above approximation.

In Fig. (1), $`P\left(\mathrm{\Gamma }\right)`$ is calculated for 1D systems composed of $`100`$ randomly distributed scatterers in a linear segment for two different values of the uniform optical density $`\rho `$: $`\rho =1`$ and $`\rho =10`$ scatterers per wavelength. In 1D, all eigenstates are known to be exponentially localized even for weak disorder and $`\xi `$ is of the order of the mean free path $`\mathrm{}`$. $`P\left(\mathrm{\Gamma }\right)`$ exhibits a power law with an exponent very close to $`1`$, in good agreement with (1). In addition, the exponent does not change by increasing $`\rho `$, i.e., by decreasing $`\xi `$. This demonstrates that the algebraic decay $`P\left(\mathrm{\Gamma }\right)\mathrm{\Gamma }^1`$ in the localized regime is valid not only for 1D models of mesoscopic transport titovfyo ; terraguarnieri , but also for our model of wave propagation in disordered media. At large $`\mathrm{\Gamma }`$, $`P\left(\mathrm{\Gamma }\right)`$ decays faster than algebraically. This can be explained by the fact that this region is dominated by short living resonances, typically close to the boundaries, for which the prediction (1) breaks down.

Fig. (2) shows $`P\left(\mathrm{\Gamma }\right)`$ for 2D systems of $`2500`$ scatterers randomly distributed in a square for $`\rho =1`$ and $`\rho =10`$ scatterers per wavelength squared. According to the scaling theory of localization, all eigenstates in 2D closed systems should be localized, but $`\xi `$ scales exponentially with $`\mathrm{}`$ costas . The $`\mathrm{\Gamma }^1`$ decay of $`P\left(\mathrm{\Gamma }\right)`$ in Fig. (2) is clearly visible for both values of $`\rho `$ used, with an exponent very close to $`1`$, in excellent agreement with (1). Remark however that the range of the power law is broader when $`\rho `$ is higher. Increasing $`\rho `$ means decreasing $`\mathrm{}`$ and, consequently, decreasing $`\xi `$. The range of the algebraic decay $`P\left(\mathrm{\Gamma }\right)\mathrm{\Gamma }^1`$ is expected to be broader as more and more states become localized. Such a behavior was also reported in numerical calculations within the Anderson model terraguarnieri . For large $`\mathrm{\Gamma }`$, $`P\left(\mathrm{\Gamma }\right)`$ decays again faster than algebraically as in the 1D case.

In Fig. (3) the 3D case is considered, where $`P\left(\mathrm{\Gamma }\right)`$ is calculated for systems composed by $`1000`$ point scatterers randomly distributed in a sphere (radius $`R`$) for $`\rho =1`$, $`\rho =10`$, $`\rho =30`$ and $`\rho =60`$ scatterers per wavelength cubed. In 3D, the system is expected to undergo, upon varying the degree of disorder, a transition from extended states to localized states. It is therefore interesting to investigate if and how this transition manifests itself in $`P\left(\mathrm{\Gamma }\right)`$. As in the 2D case, we notice that, as $`\rho `$ increases, the range of the algebraic decay $`P\left(\mathrm{\Gamma }\right)\mathrm{\Gamma }^\alpha `$ increases. We also remark that, as $`\rho `$ increases, the associated exponents tend more and more to the value $`1`$. The exponents, obtained by a linear fit in the range where the power law is present, are $`\alpha 0.76`$ for $`\rho =1`$, $`\alpha 0.83`$ for $`\rho =10`$, $`\alpha 0.95`$ for $`\rho =30`$ and $`\alpha 1.1`$ for $`\rho =60`$. This suggests, according to (1), the onset of the localized regime for higher $`\rho `$. In fact, the Ioffe-Regel criterion for localization ($`k\mathrm{}<1`$) is estimated to be satisfied for $`\rho >2\pi ^220`$ for scatterers at resonance. This condition is fulfilled for $`\rho =30`$ and $`\rho =60`$, for which $`\alpha `$ is very close to $`1`$, showing that the system with these densities are indeed in the localized regime and confirming that the power law $`P\left(\mathrm{\Gamma }\right)\mathrm{\Gamma }^1`$ can be considered a genuine signature of Anderson localization of light. The fact that the exponents for $`\rho =30`$ and $`\rho =60`$ are not exactly equal to $`1`$ can probably be attributed to finite-size effects. Once again, note that $`P\left(\mathrm{\Gamma }\right)`$ decays faster than a power law for large $`\mathrm{\Gamma }`$.

Let us now compare the mean width of the distribution to the inverse Thouless time $`\mathrm{\Gamma }_T=1/t_T=6D_B/R^2`$, where $`D_B`$ is the Boltzmann diffusion constant. To estimate $`\mathrm{\Gamma }_T`$ one should recall that in 3D, $`D_B=v_E\mathrm{}^{}/3`$, where $`\mathrm{}^{}`$ is the transport mean free path (which is, for point scatterers, equal to $`\mathrm{}`$) and $`v_E`$ the energy transport velocity, $`v_Ec_0/[1+\tau _{dwell}/\tau _{mf}]`$ reports , with $`\tau _{dwell}=1/\mathrm{\Gamma }_0`$ the dwell time in a single scattering and $`\tau _{mf}=\mathrm{}/c_0`$ the mean free time. Applying these considerations, $`\mathrm{\Gamma }_T`$ can be written as $`\mathrm{\Gamma }_T/\mathrm{\Gamma }_02\left(\mathrm{}/R\right)^2`$. Fig. (4) exhibits $`P\left(\mathrm{\Gamma }\right)`$ for the same optical densities of Fig. (3) but now with $`\mathrm{\Gamma }`$ normalized to $`\mathrm{\Gamma }_T`$. For low $`\rho `$ ($`\rho =1`$ and $`\rho =10`$), $`P\left(\mathrm{\Gamma }\right)`$ is peaked near $`\mathrm{\Gamma }_T`$, showing that the system is in the diffusive regime. Notice that there is a non-vanishing probability to find modes that live much longer than $`t_T`$ even in the diffusive regime, the so-called “prelocalized” modes preloc . As $`\rho `$ increases, we observe that $`P\left(\mathrm{\Gamma }\right)`$ is no longer centered at $`\mathrm{\Gamma }_T`$ and that the probability to find a mode with resonance width smaller than $`\mathrm{\Gamma }_T`$ also increases. This means that, on average, the modes live longer than $`t_T`$. At the same time, Fig. (3) shows that localization manifests itself in $`P\left(\mathrm{\Gamma }\right)`$ not only via the broadening of the power law range but also via the fact that the associated exponents approach to $`1`$. We conclude again that the $`P\left(\mathrm{\Gamma }\right)\mathrm{\Gamma }^1`$ behavior is an unambiguous signature of Anderson localization of light in open media. It must be mentioned that the present 3D study may be relevant for recent multiple light scattering experiments in atomic media nice , for which modeling the scatterers by pointlike dipoles constitutes an excellent approximation.

In summary, we have studied the distribution of resonance widths $`P\left(\mathrm{\Gamma }\right)`$ in 1D, 2D and 3D multiple light scattering systems composed of randomly distributed pointlike scalar dipoles. We have developed a simple physical argument, based on the exponential decay of localized eigenfunctions, to show that $`P\left(\mathrm{\Gamma }\right)`$ should follow an universal power law $`P\left(\mathrm{\Gamma }\right)\mathrm{\Gamma }^1`$ decay in the localized regime. This prediction was confirmed by extensive numerical calculations and demonstrates that the $`P\left(\mathrm{\Gamma }\right)\mathrm{\Gamma }^1`$ behavior can be interpreted as an unambiguous signature of Anderson localization of light in open media.

###### Acknowledgements.

Most of the computations presented in this paper were performed on the cluster PHYNUM (CIMENT, Grenoble). This work was supported in part by the Polonium contract 03290 VK. M.R. and A.O. were supported by the Polish State Committee for Scientific Research (KBN) under Grant no. 2 P03B 044 19. F.A.P. wishes to thank CNPq/Brazil for financial support.


format: text

text:
# Naked Singularities in Higher Dimensional Szekeres space-time

## I Introduction

In recent years, there is an exhaustive study of gravitational collapse of inhomogeneous spherical dust \[1-4\]. It is mainly concentrated to central shell focusing singularity. The behaviour of the singularity (black hole or naked singularity) has been analyzed mainly by the studying of outgoing radial null geodesic and the strength is measured using due to Tipler . It has been shown that appearance of naked singularity strongly depends on the initial data \[6-10\].

In contrast, there is very little information about singularity formation and structure in non-spherical collapse. Here horizon will be formed if and only if the hoop conjecture  is satisfied. Most of the studied so far in non-spherical collapse are by numerical analysis and for some definite shape of the gravitating mass \[12-16\].

However, the study of gravitational collapse for quasi-spherical Szekeres space-time  was started by Szekeres himself  long ago. Afterwards, it was studied further by Joshi et al  and extensively by Goncalves . Recently, we have obtained solutions for ($`n+2`$)-dimensional Szekeres space-time with perfect fluid (or dust) as the matter content . Subsequently, we have also analyzed the gravitational collapse for the above dust solution to examine the local behaviour of the singularity so formed .

In this paper, we extend our study for global characteristic of the singularity by studying both null and time like geodesic originated from the singularity using the formalism of Barve et al . The paper is organized as follows: In section II we have written the basic equations and regularity conditions. In section III and IV, we have studied the radial null and time-like geodesics originated from the singularity. We have studied the local visibility of the singularity for radial null and time-like geodesics due to collapse in section V. Finally the paper ends with a short discussion.

## II Basic Equations and regularity conditions

The metric ansatz for ($`n+2`$)D Szekeres space-time is given by

$$ds^2=dt^2e^{2\alpha }dr^2e^{2\beta }\underset{i=1}{\overset{n}{}}dx_i^2$$
(1)
where $`\alpha `$ and $`\beta `$ are functions of all the ($`n+2`$) space-time co-ordinates with expression  (where $`\beta ^{}0`$)

$$e^\alpha =R^{}+R\nu ^{}$$
(2)
and

$$e^\beta =R(t,r)e^{\nu (r,x_1,\mathrm{},x_n)}$$
(3)
Here $`R`$ satisfied the differential equation

$$\dot{R}^2=f(r)+\frac{F(r)}{R^{n1}}$$
(4)
and the expression for $`\nu `$ is

$$e^\nu =A(r)\underset{i=1}{\overset{n}{}}x_i^2+\underset{i=1}{\overset{n}{}}B_i(r)x_i+C(r)$$
(5)
with the restriction

$$\underset{i=1}{\overset{n}{}}B_i^24AC=f(r)1$$
(6)
for the arbitrary functions $`A(r),B_i(r)`$ and $`C(r)`$. Also in the expression (4), $`f(r)`$ and $`F(r)`$ are arbitrary functions of $`r`$ alone. But due to complexity of the problem we shall restrict ourselves to marginally bound case only (i.e., $`f(r)=0`$). In this case equation (4) can easily integrated to give

$$R=\left[r^{\frac{n+1}{2}}\frac{n+1}{2}\sqrt{F(r)}t\right]^{\frac{2}{n+1}}$$
(7)
and the energy density for the matter field

$$\rho (t,r,x_1,\mathrm{},x_n)=\frac{n}{2}\frac{F^{}+(n+1)F\nu ^{}}{R^n(R^{}+R\nu ^{})}$$
(8)
takes the form

$$\rho _i(r,x_1,\mathrm{},x_n)=\rho (0,r,x_1,\mathrm{},x_n)=\frac{n}{2}\frac{F^{}+(n+1)F\nu ^{}}{r^n(1+r\nu ^{})}$$
(9)
initially at $`t=0`$, where we have chosen the scaling of $`R`$ as

$$R(0,r)=r.$$

Now assuming that the collapse starts from a regular initial hypersurface so we take the following series forms for $`F(r)`$, $`\rho _i(r)`$ and $`\nu ^{}(r)`$

$$F(r)=\underset{j=0}{\overset{\mathrm{}}{}}F_jr^{n+j+1},$$
(10)
$$\rho _i(r)=\underset{j=0}{\overset{\mathrm{}}{}}\rho _jr^j$$
(11)
and

$$\nu ^{}(r)=\underset{j=1}{\overset{\mathrm{}}{}}\nu _jr^j,(\nu __11)$$
(12)
where the coefficients are related as

$`\rho _0={\displaystyle \frac{n(n+1)}{2}}F_0,\rho _1={\displaystyle \frac{n}{2}}\left(n+1+{\displaystyle \frac{1}{1+\nu __1}}\right)F_1,`$
$`\rho _2={\displaystyle \frac{n}{2}}\left[\left(n+1+{\displaystyle \frac{2}{1+\nu __1}}\right)F_2{\displaystyle \frac{F_1\nu __0}{(1+\nu __1)^2}}\right],`$
$$\rho _3=\frac{n}{2}\left[\left(n+1+\frac{3}{1+\nu __1}\right)F_3\frac{2F_2\nu __0}{(1+\nu __1)^2}\frac{(1+\nu __1)\nu __1\nu __0^2}{(1+\nu __1)^3}F_1\right],$$
(13)
$$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{},$$

$$OR$$

$`\rho _0={\displaystyle \frac{n}{2}}\left[{\displaystyle \frac{F_1}{\nu _0}}+(n+1)F_0\right],\rho _1={\displaystyle \frac{n}{2}}\left[{\displaystyle \frac{2F_2}{\nu _0}}+\left\{(n+1){\displaystyle \frac{\nu _1}{\nu _0^2}}\right\}F_1\right],`$
$`\rho _2={\displaystyle \frac{n}{2}}\left[{\displaystyle \frac{3F_3}{\nu _0}}+\left\{(n+1){\displaystyle \frac{2\nu _1}{\nu _0^2}}\right\}F_2+\left({\displaystyle \frac{\nu _1^2}{\nu _0^3}}{\displaystyle \frac{\nu _2}{\nu _0^2}}\right)F_1\right],`$
$$\rho _3=\frac{n}{2}\left[\frac{4F_4}{\nu _0}+\left\{(n+1)\frac{3\nu _1}{\nu _0^2}\right\}F_3+2\left(\frac{\nu _1^2}{\nu _0^3}\frac{\nu _2}{\nu _0^2}\right)F_2+\left(\frac{2\nu _1\nu _2}{\nu _0^3}\frac{\nu _3}{\nu _0^2}\frac{\nu _1^3}{\nu _0^4}\right)F_1\right],$$
(14)
$$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{},$$

according as $`\nu __1>1`$ or $`\nu __1>=1`$.

The singularity curve $`t=t_s(r)`$ for the shell focusing singularity is characterized by

$$R(t_s(r),r)=0$$

and we have

$$t_s(r)=\frac{2}{n+1}\frac{r^{\frac{n+1}{2}}}{\sqrt{F(r)}}$$
(15)
where $`t_0=\frac{2}{(n+1)\sqrt{F_0}}`$ is the time for the central shell focusing singularity. However, near $`r=0`$, above singularity curve can be approximately written as

$$t_s(r)=t_0\frac{F_m}{(n+1)F_0^{3/2}}r^m$$
(16)
where $`m(1)`$ is an integer and $`F_m`$ is the first non-vanishing term beyond $`F_0`$.

## III Radial Null Geodesics

The equation of the outgoing radial null geodesic (ORNG) which passes through the central singularity in the past is taken as (near $`r=0`$)

$$t_{ORNG}=t_0+ar^\xi ,$$
(17)
to leading order in $`t`$-$`r`$ plane with $`a>0,\xi >0`$. Now to visualize the singularity the time $`t`$ in the geodesic (17) should be less than $`t_s(r)`$ in equation (16). Hence comparing these times we have the restrictions

$$\xi m\text{and}a<\frac{F_m}{(n+1)F_0^{3/2}}$$
(18)
Further, for the metric (1) an outgoing radial null geodesic should satisfy

$$\frac{dt}{dr}=R^{}+R\nu ^{}$$
(19)
We shall now examine the feasibility of the null geodesic starting from the singularity with the above restrictions for the following two cases namely, (i) $`\xi >m`$ and (ii) $`\xi =m`$.

When $`\xi >m`$ then near $`r=0`$ the solution for $`R`$ in (7) simplifies to

$$R=\left(\frac{F_m}{2F_0}\right)^{\frac{2}{n+1}}r^{\frac{2m}{n+1}+1}$$
(20)
Now combining (17) and (20) in equation (19) we get (upto leading order in $`r`$)

$$a\xi r^{\xi 1}=\left(\nu __1+1+\frac{2m}{n+1}\right)\left(\frac{F_m}{2F_0}\right)^{\frac{2}{n+1}}r^{\frac{2m}{n+1}}$$
(21)
which implies

$$\xi =1+\frac{2m}{n+1}\text{and}a=\frac{(\nu __1+1+\frac{2m}{n+1})}{(1+\frac{2m}{n+1})}\left(\frac{F_m}{2F_0}\right)^{2/(n+1)}$$
(22)
We note that if $`\nu __1=0`$ then the above results are same as TBL model , so we restrict ourselves to $`\nu __10`$. As $`\xi >m`$ so from (22)

$$m<\frac{n+1}{n1}$$

Since $`m`$ is an integer so we have the following possibilities:

$`n=2(4D\text{ space-time}):m=1,2,\xi =\frac{4}{3},\frac{7}{3}`$
$`n3(5D\text{ or  higher dim. space-time}):m=1,\xi =1+\frac{2}{n+1}`$.

Now for $`\nu __1>1`$, if we assume that the initial density gradient falls off rapidly to zero at the centre then we must have $`\rho _1=0,\rho _2<0`$ i.e., $`F_1=0,F_2<0`$. Hence the least value of $`m`$ should be 2. Thus for a self consistent outgoing radial null geodesic starting from the singularity we must have $`n=2`$ only. So the singularity will be naked only in four dimension. However, if we relax the restriction on the initial density and simply assume that the initial density gradually decreases as we go away from the centre (i.e., $`\rho _1<0`$) then $`F_1<0`$ and $`m=1`$ is the possible solution. Therefore, naked singularity is possible in any dimension. However for $`\nu __1=1`$, the assumption about the initial density gradient ($`\rho _1=0,\rho _2<0`$) does not implying $`F_1=0`$. Hence here $`m=1`$ is a possible solution for all dimension with $`F_1<0`$. Therefore naked singularity may occur in any dimension for this choice of the parameter $`\nu _1`$.

Now we shall study the second case namely $`\xi =m`$. Here also for $`\nu __1=0`$ we have identical results as in TBL model . For $`\nu __10`$, using the same procedure as above we get

$$m=\frac{n+1}{n1}$$

and

$$a=\frac{1}{m}\left(\frac{n+1}{2}\sqrt{F_0}a\frac{F_m}{2F_0}\right)^{\frac{1n}{1+n}}\left[(\nu __1+m)\frac{F_m}{2F_0}+\frac{1}{2}(n+1)(\nu __1+1)\sqrt{F_0}a\right]$$
(23)
Hence integral solution for $`m`$ is possible only for $`n=2`$ and 3. So consistent outgoing radial null geodesic through the central shell focusing singularity is possible only upto five dimension i.e., we cannot have naked singularity for higher dimensional space-time greater than five. It is to be noted that this conclusion is not affected by the value of the parameter $`\nu __1`$.

Now we shall examine whether the restriction (18) for $`a`$ is consistent with the expression $`a`$ in equation (23). In fact equation (23) takes the form

$$2^{\frac{2}{n+1}}bm=\left[(n+1)b\zeta \right]^{\frac{1}{m}}\left[(\nu __1+m)\zeta +(n+1)(\nu __1+1)b\right]$$
(24)
with the transformation

$$a=bF_0^{\frac{1}{n1}},F_m=\zeta F_0^{\frac{m}{2}+1}.$$
(25)
Since equation (24) is a real valued equation of $`b`$, so we must have

$$(n+1)b+\zeta <0,$$

which using (25) gives us the restriction on $`a`$ in equation (18). Hence the geodesic (17) will have consistent solution for $`a`$ and $`m`$. So the above conclusion regarding the formation of naked singularity is justified. Further, introducing the variable $`\varphi `$ by the relation

$$\varphi =(n+1)b\zeta $$
(26)
We have seen from equation (24)

$$4\varphi ^{n1}(\zeta +\varphi )^{n+1}=[2\zeta (n1)(\nu __1+1)\varphi ]^{n+1}$$
(27)
with the restriction $`0<\varphi <\zeta `$.

Finally, we investigate whether there is only one null geodesic emanating from the singularity (i.e., when singularity is visible for an infinitesimal time) or it is possible to have an entire family of geodesic through the singularity (this occurs when singularity is visible for an infinite time). So we write the equation for the outgoing radial null geodesics to next order as

$$t=t_0+ar^\xi +hr^{\xi +\sigma }$$
(28)
where as before $`a,\xi ,h,\sigma `$ are all positive. Then proceeding in the same way we have

(i) when $`\xi >m`$, $`a`$ and $`\xi `$ have the same expressions as in equation (22) while

$$\sigma =\frac{m(1n)}{1+n}+1$$
(29)
and

$$h=\frac{a}{\xi +\sigma }(\nu __1\sigma )\sqrt{F_0}\left(\frac{F_m}{2F_0}\right)^{\frac{1n}{1+n}}$$
(30)
We note that since here $`m`$ can have values 1 and 2 for $`n=2`$ and $`m=1`$ for $`n3`$ so $`\sigma `$ is always positive. But for $`h`$ to be positive definite, $`\nu __1`$ is restricted within the range $`1\nu __1\sigma `$ i.e., $`\nu __1`$ is negative definite.

(ii) for $`\xi =m`$, $`a`$ and $`\xi `$ have the same expressions as in equation (23) and $`\sigma `$ is obtained from the equation

$$m+\sigma =2^{\frac{n1}{n+1}}b(n+1)\left[\zeta (n+1)b\right]^{\frac{2n}{n+1}}2^{\frac{n1}{n+1}}\nu __1F_0^{1/2}\left[\zeta (n+1)b\right]^{\frac{1n}{1+n}}$$
(31)
with expressions for $`b`$ and $`\zeta `$ from equation (25), while the other constant $`h`$ is totally arbitrary. Thus it is possible to have an entire family of outgoing null geodesics terminated in the past at the singularity, provided $`\sigma `$ obtained from equation (31) is positive definite.

## IV Radial Time-like Geodesics

We shall now investigate whether it is possible to have any outgoing time-like geodesic originated from the singularity. For simplicity of calculation we shall consider only outgoing radial time-like geodesic (ORTG). Let us denote by $`K^a=\frac{dx^a}{d\tau }`$, a unit tangent vector field to a ORTG with $`\tau `$ an affine parameter along the geodesic. Hence from the geodesic equation we have

$$K^t\left(K^r\right)^.+2\dot{\alpha }K^tK^r+K^r\left(K^r\right)^{}+\alpha ^{}\left(K^r\right)^2=0$$
(32)
with

$$K^t=\pm \sqrt{1+e^{2\alpha }\left(K^r\right)^2}.$$
(33)
In order to satisfy the above two equations the simplest choice for ($`K^r,K^t`$) is

$$K^t=\pm 1,K^r=0.$$
(34)
and so we have the solution

$$tt_0=\pm (\tau \tau _0),r=r_0=\text{constant}.$$
(35)
Here $`+`$ (or $``$) sign corresponds to ORTG (or ingoing RTG) and $`\tau _0`$ is the proper time at which radial time-like geodesic passes through the central singularity.

Now similar to null geodesic let us choose the radial time-like geodesic near the singularity to be of the form (to leading order)

$$t_{ORTG}(r)=t_0+cr^p$$
(36)
with $`c`$ and $`p`$ as positive constants. Further, consistent with equations (34) and (35) we assume

$$K^r(t,r)=A(tt_0)^\lambda r^\delta $$
(37)
where $`A(>0),\lambda `$ and $`\delta `$ are constants. Hence we have

$$K^r(t_{ORTG},r)=Ar^q,q=\lambda p+\delta .$$
(38)
Also from the geodesic equation (36) using (38) we get

$$\frac{dt_{ORTG}}{dr}=cpr^{p1}=\frac{K^t}{K^r}=\sqrt{A^2r^{2q}+(R^{}+R\nu ^{})^2}$$
(39)
Using the solution for $`R`$ near the singularity in equation (39) and equating equal powers of $`r`$ we have

(i) $`p=1q`$, $`c=\frac{1}{A(1q)}`$ if $`\frac{2m}{n+1}<q<1`$.

(ii) $`p=1+\frac{2m}{n+1}`$, $`c=\frac{(\nu __1+1+\frac{2m}{n+1})}{(1+\frac{2m}{n+1})}\left(\frac{F_m}{2F_0}\right)^{\frac{2}{n+1}}`$ if $`q<\frac{2m}{n+1}`$

(iii) $`p=1+\frac{2m}{n+1}`$, $`c=\frac{1}{p}\left[A^2+\left(\frac{F_m}{2F_0}\right)^{\frac{4}{n+1}}\left(\nu __1+1+\frac{2m}{n+1}\right)^2\right]^{1/2}`$ if $`q=\frac{2m}{n+1}`$

Therefore equation (36) has consistent solution for $`c`$ and $`p`$ depending on the parameters involved. So it is possible to have outgoing radial time-like geodesic originated from the singularity.

## V Local Visibility

We have shown the existence of outgoing both null and time-like geodesics which were originated in the past from the singularity. We shall now examine whether these geodesics are visible to any non-space like observer. Now for visibility of the singularity (local nakedness) any future directed geodesics (null or time-like) starting from the singularity should be outside the domain of dependence of any trapped surfaces both before and at the time of formation of apparent horizon (AH) which is the outer boundary of the trapped surface. Now if $`t_{ah}(r)`$ is the time of formation of apparent horizon so we must have

$$\dot{R}(t_{ah}(r),r)=1$$

i.e.,

$$R^{n1}(t_{ah}(r),r)=F(r)$$
(40)
The apparent horizon is given by the curve

$$t_{ah}(r)=t_0\frac{1}{(n+1)F_0^{3/2}}\left(F_mr^m+O(r^{m+1})\right)\frac{2}{n+1}F_0^{\frac{1}{n1}}\left(r^{\frac{n+1}{n1}}+O(r^{\frac{n+1}{n1}+1})\right)$$
(41)
As the apparent horizon and the singularity curve form at the same time at $`r=0`$, so the visibility of the singularity is determined by the relative slopes of the curve $`t_{ah}(r)`$ and the curve for outgoing radial non-space like geodesics (denoted by $`t_{ORG}(r)`$). Hence the necessary and sufficient condition for singularity to be at least locally naked is that

$$\begin{array}{c}lim\\ r0\end{array}\begin{array}{c}\left\{(\frac{dt_{ah}}{dr})/(\frac{dt_{ORG}}{dr})\right\}>1\end{array}$$
(42)
Now we shall examine this condition for both ORNG and ORTG separately.

### V.1 Outgoing Radial Null Geodesic

In this case the ratio of the slopes is

$`{\displaystyle \frac{(\frac{dt_{ah}(r)}{dr})}{(\frac{dt_{ORNG}(r)}{dr})}}={\displaystyle \frac{1}{a\xi (n+1)F_0^{3/2}}}\left(mF_mr^{m\xi }+O(r^{m\xi +1})\right)`$
$$\frac{2}{a\xi (n1)}F_0^{\frac{1}{n1}}\left(r^{\frac{n+1}{n1}\xi }+O(r^{\frac{n+1}{n1}\xi +1})\right)$$
(43)
So we have the following possibilities:

(i) If $`m>\frac{n+1}{n1}`$ then the above ratio of the slope is always negative, hence the condition is violated. So we always get black hole.

(ii) If $`m<\frac{n+1}{n1}`$ then the above ratio of the slopes approaches to $`+\mathrm{}`$ as $`r0^+`$ for $`m<\xi `$ and while the limit will be finite and greater than unity for $`m=\xi `$ provided $`a<\frac{m}{(n+1)\xi }\frac{F_m}{F_0^{3/2}}`$. Thus the singularity will be locally naked for any dimension if $`m=1`$ and $`a<\frac{1}{(n+1)\xi }\frac{F_1}{F_0^{3/2}}`$ while naked singularity appears only for four dimension if $`m=2`$ and $`a<\frac{2}{3\xi }\frac{F_2}{F_0^{3/2}}`$.

(iii) If $`m=\frac{n+1}{n1}`$ then we have shown that naked singularity appears only upto five dimension and the singularity will be locally visible for

$$F_m<2F_0^{\frac{m}{2}+1},m<\xi $$

or
$$F_m<2F_0^{\frac{m}{2}+1}a(n+1)F_0^{3/2},m=\xi .$$

### V.2 Outgoing Radial Time-like Geodesic

In this case the ratio of the slopes is

$`{\displaystyle \frac{(\frac{dt_{ah}(r)}{dr})}{(\frac{dt_{ORTG}(r)}{dr})}}={\displaystyle \frac{1}{cp(n+1)F_0^{3/2}}}\left(mF_mr^{mp}+O(r^{mp+1})\right)`$
$$\frac{2}{cp(n1)}F_0^{\frac{1}{n1}}\left(r^{\frac{n+1}{n1}p}+O(r^{\frac{n+1}{n1}p+1})\right)$$
(44)
We shall now study the following possibilities for local visibility:

(i) If $`m>\frac{n+1}{n1}`$ then as before only black hole appears.

(ii) If $`m<\frac{n+1}{n1}`$ we have the conclusion as in ORNG except that instead of restricting $`a`$, we have here $`c<\frac{m}{(n+1)p}\frac{F_m}{F_0^{3/2}}`$.

(iii) Similar is the situation for $`m=\frac{n+1}{n1}`$.

Hence the conditions for local visibility of naked singularity is consistent for both time-like and null geodesic.

## VI Discussion and Concluding Remarks

In this paper we have studied quasi-spherical gravitational collapse in $`(n+2)`$-D Szekeres space-time by studying both the null and time-like geodesics originated from the central shell focusing singularity. Following the approach of Barve et al  for null geodesic we have shown that if we choose the parameter $`\nu __1`$ to be greater than $`1`$, then it is possible to have a class of outgoing radial null geodesic for four and five dimensions only with the restriction that initial density falls off rapidly to the centre. However for $`\nu __1=1`$, naked singularity is possible in all dimensions irrespective of the assumption on the initial density. On the other hand, for time-like geodesic the results are very similar to the study of null geodesic. Regarding local visibility we have consistent results for both ORNG and ORTG with some restrictions on the parameters involved in the equation of the geodesics. Therefore we conclude that it is possible to have at least local naked singularity for the given higher dimensional non-spherical space-time.

Acknowledgement:

One of the authors (U.D) thanks CSIR (Govt. of India) for the award of a Senior Research Fellowship.

References:
$`[1]`$ A. Ori and T. Piran, Phys. Rev. Lett. 59 2137 (1987); J. P. S. Lemos Phys. Rev. Lett. 68 1447 (1992).
$`[2]`$ P.S. Joshi, Global Aspects in Gravitation and Cosmology(Oxford Univ. Press, Oxford, 1993).
$`[3]`$ P. S. Joshi and I. H. Dwivedi, Commun. Math. Phys. 166 117 (1994).
$`[4]`$ P. S. Joshi and I. H. Dwivedi, Class. Quantum Grav. 16 41 (1999).
$`[5]`$ F. J. Tipler, Phys. Lett. A 64 8 (1987).
$`[6]`$ F. C. Mena, R. Tavakol and P. S. Joshi, Phys. Rev. D 62 044001 (2000).
$`[7]`$ P.S. Joshi, N. Dadhich and R. Maartens, Phys. Rev. D 65 101501(R)(2002).
$`[8]`$ A. Banerjee, U. Debnath and S. Chakraborty, gr-qc/0211099 (2002)(accepted in Int. J. Mod. Phys. D) .
$`[9]`$ U. Debnath and S. Chakraborty, gr-qc/0211102 (2002) .
$`[10]`$ R. Goswami and P.S. Joshi, gr-qc/02112097 (2002) .
$`[11]`$ K. S. Thorne, 1972 in Magic Without Magic: John Archibald Wheeler, Ed. Klauder J (San Francisco: W. H. Freeman and Co.).
$`[12]`$ S. L. Shapiro and S. A. Teukolsky, Phys. Rev. Lett. 66 994 (1991).
$`[13]`$ T. Nakamura, M. Shibata and K.I. Nakao, Prog. Theor. Phys. 89 821 (1993) .
$`[14]`$ T. Harada, H. Iguchi and K.I. Nakao, Phys. Rev. D 58 041502 (1998) .
$`[15]`$ H. Iguchi, T. Harada and K.I. Nakao, Prog. Theor. Phys. 101 1235 (1999); Prog. Theor. Phys. 103 53 (2000).
$`[16]`$ C. Barrabes, W. Israel and P. S. Letelier, Phys. Lett. A 160 41 (1991); M. A. Pelath, K. P. Tod and R. M. Wald, Class. Quantum Grav. 15 3917 (1998).
$`[17]`$ P. Szekeres, Commun. Math. Phys. 41 55 (1975).
$`[18]`$ P. Szekeres, Phys. Rev. D 12 2941 (1975).
$`[19]`$ P. S. Joshi and A. Krolak, Class. Quantum Grav. 13 3069 (1996); S. S. Deshingkar, S. Jhingan and P. S. Joshi, Gen. Rel. Grav. 30 1477 (1998).
$`[20]`$ S. M. C. V. Goncalves, Class. Quantum Grav. 18 4517 (2001).
$`[21]`$ S. Chakraborty and U. Debnath , gr-qc/0304072 (2003).
$`[22]`$ U. Debnath, S. Chakraborty and J. D. Barrow, gr-qc/0305075 (2003).
$`[23]`$ S. Barve, T. P. Singh, C. Vaz and L. Witten, Class. Quantum Grav. 16 1727 (1999).


text:
# Ab initio Simulations of Fe-based Ferric Wheels

## Introduction

There has been a remarkable increase in the number of studies on molecular magnets during the last decade which clearly is driven by the progress in the chemistry of such systems. Many different families of molecular magnetic systems composed of all possible transition-metal atoms in various combinations with stabilizing organic ligands have been synthesized, their crystal structure determined and their magnetic properties analyzed .

While there is a large amount of experimental data available, the theoretical description in terms of quantitative – and often also qualitative – understanding is missing in some cases, in others it stays behind the experimental progress. The traditional approach to describe these systems is based on model Hamiltonians (“exact” diagonalization within, e.g., Heisenberg model). One the other hand, with the advance in computer technology and the development in theoretical methods, most notably in density-functional theory (DFT) , first-principles calculations became feasible. Both approaches simplify real systems and are prone to certain shortcomings. The model Hamiltonians neglect the true chemical environment and bonding. All interactions are reduced to just few model parameters. One is dependent on experimental input in order to estimate these parameters, and the accuracy of quantitative predictions depends on the parameters chosen. Even more problematic, it is not *a priori* clear which interactions are important and should be included in the model and which are negligible. In contrast, first-principles calculations based on quantum mechanics use the atomic positions as the only experimental input and, as chemical bonding is included in a straightforward way, they are able to provide microscopic foundations of the magnetic behavior. Unfortunately, the solution of Schrödinger’s equation is limited to a few electrons so that one has to use approximations. The chemical complexity of molecular magnets (or, their large number of atoms) excluded so far the “straightforward” simulations using approaches of superior accuracy, like quantum Monte Carlo. Density-functional based calculations have become the working horse in today’s first-principles simulations. DFT has been proven to be very accurate and reliable in many cases, although electron correlation effects beyond the mean-field approximation in DFT are neglected, which limits the accuracy and may lead to failures in the description of a given system.

Both approaches, the model-based and DFT-based, are useful and important for our understanding, and are best used together. While DFT can provide the input parameters for models (e.g. the exchange parameters $`J`$ for a Heisenberg model), the model Hamiltonians can deliver information on many-body effects, temperature dependence, and quantum fluctuations, which are beyond the reach of DFT. By describing the relevant physics from different points of view one develops a complete picture where essential features do not remain unexplored.

Concentrating in the following on DFT-based strategies, we note that molecular magnets turn out to be hard cases for computational treatment – for different reasons in different established calculation approaches. The openness of the crystal structure creates problems for traditional muffin-tin- or atomic sphere approximation-based tools, and requires requires a “full-potential” formalism instead. The presence of transition metals and/or other constituents with deep semicore states, which can not be neglected, creates problems for pseudopotential planewave methods. Moreover, the need for a large (albeit largely empty) “simulation box” around a molecule makes planewave-related methods expensive. An abundance of hydrogen and hence strong disbalance in atom “sizes” makes the linearized augmented plane wave scheme numerically demanding. The often unusual chemical composition and coordination reduces the usefulness of *a priori* tuned basis sets in those tight-binding-type methods which heavily depend on pre-adjusted basis functions for their efficiency. Further, the option to treat non-collinear spin-density distribution and to have accurate forces on atoms for studying structure relaxation and dynamics may be advantageous and sometimes necessary. As a consequence, flexible, efficient and accurate methods (and codes) are in demand, of which several have been successfully applied to systems of reasonable complexity.

## Calculations

As an example of a system of practical interest we consider the hexanuclear “ferric wheels” $`M`$Fe6\[N(CH2CH2O)3\]6Cl ($`M`$ = Li, Na, see Fig. 1), synthesized at the Institut für Organische Chemie in Erlangen  and labeled as substances 4 and 3 in the latter publication. These materials belong to a larger family of Fe-coronates, of which 8- and 10-member rings are also known. A common feature of these substances is the bonding of Fe atoms via doubled oxygen bridges, somehow resembling the 90 coupling of magnetic atoms in transition-metal oxides. In the vicinity of each Fe atom one finds a nitrogen atom, connected via C2H4 chains to oxygen atoms. The Li- and Na-centered molecules have a very similar structure, with only slightly different Fe–O–Fe angles (101.1 for Li-wheel and 103.3 for Na-wheel).

The interest in the magnetic properties of the system is the presence of large magnetic moments, traditionally argued to correspond to $`S`$=5/2 on the Fe site, thus implying a highly ionized Fe(III) state. Such an identification of the iron magnetic moment is consistent with magnetization and torque measurements , from which Waldmann *et al.* derived the values of exchange parameters $`J`$ for the spin Hamiltonian of the Heisenberg model:

$$H=J\left(\underset{i=1}{\overset{5}{}}𝐒_i𝐒_{i+1}+𝐒_6𝐒_1\right)+(\text{anisotropy term})+(\text{Zeeman term}).$$
(1)
The obtained values of $`J`$ are $``$18 to $``$20 K for the Li-wheel (depending on sample and method) and $``$22.5 to $``$25 K for the Na-wheel, thus implying an antiferromagnetic (AFM) ground state. While the data from bulk magnetic measurements seem to be well established by now and reliable, virtually no spatially resolved probing of the magnetic structure has been carried out yet. Recent X-ray photoelectron and X-ray emission spectroscopy studies  allow to probe the electronic structure in the valence band and on the Fe site, albeit without resolution in spin. This is exactly where *ab initio* calculations may yield a valuable contribution in providing site- and spin-resolved densities of states, as well as the spatial distribution of charge and magnetization densities. The verification of these predictions demands the comparison with (atom- and/or spin-averaged) spectroscopic data, and eventually comparison with other calculational approaches at a different level of approximations.

In the present study, we compare the results of electronic structure calculations by two different methods – Siesta and NRLMOL . Both methods are based on the DFT; from the various exchange-correlation functionals implemented we used the generalized gradient approximation after Perdew–Burke–Ernzerhof  for the present study. A further similarity is that their basis sets are atom-centered functions – numerical pseudoatomic orbitals in Siesta  and contracted Gaussian-type orbitals in NRLMOL . The main difference is that Siesta uses norm-conserving pseudopotentials whereas NRLMOL is an all-electron code. NRLMOL uses the symmetry, which is relatively high in case of the Fe6-wheels, of a molecule in question in all numerical tasks. Moreover, all computationally intensive parts of the code are massively parallelized , which makes DFT calculations on molecular magnets (with a typical size of 100–200 atoms) feasible. Siesta is competitive due to the compactness and efficiency of its numerical basis sets. Moreover, it can work as an order-$`N`$ method (albeit not yet for magnetic systems), solving the DFT problem without explicit diagonalization. Other features of the two calculational methods, which are potentially important in application to molecular magnets but have not been actually used in the present study, are the availability of accurate forces in both codes, the spin-orbit coupling treated within second order perturbation theory in NRLMOL , the possibility of general-form (i.e., non-collinear) magnetization in Siesta . Being an *ab initio* pseudopotential code, Siesta is usually able to provide very good accuracy for the electronic structure and total energies without fitting to any exterior information, but still benchmark calculations make sense for critical cases and/or new systems. The present study serves in part exactly this task, as the accuracy of all-electron NRLMOL method is expected to be superior to that of Siesta. As the latter is routinely used for dynamical simulations in large systems with virtually all symmetry lost, no special treatment is provided for the use of possible point group symmetries available in the molecule. For the systems in question, due to high degeneracy of energy levels close to the chemical potential, the neglect of symmetry leads to a certain instability in the convergence, which however can be managed by keeping the mixing parameter very low (10-3–10-4) and imposing the fixed spin moment scheme (FSM, Ref. ).

## Results

In the following, we discuss the results obtained by Siesta for the Li-centered molecule, and by NRLMOL – for the Na-centered one, with the coordinates as reported from the crystallographic study of corresponding materials . The NRLMOL treatment was restricted to the ground-state AFM configuration (alternating orientations of Fe magnetic moments over the ring); the Siesta calculation addressed in addition different magnetic configurations, that allowed to make an estimate of exchange parameters.

Fig. 2 displays the partial densities of states (DOS) on Fe and its several neighbors in the AFM configuration, as calculated by both methods. The discrete levels of the energy spectra are weighted (with the charge density integrated over atom-centered spheres in NRLMOL, or according to Mulliken population analysis in Siesta), and broadened for better visibility, with broadening parameter of 0.15 eV in Siesta and 0.14 eV in NRLMOL. The local moments corresponding to integrating such partial DOS over occupied states are given in Table 1.

One notes a remarkable agreement between both calculations up to the finest details in the distribution of state densities at Fe and O sites. This is the more astonishing as the local DOS is a loosely defined property, normally dependent on the choice of region it attributes to, or, in the present context, to the choice of basis functions which are quite differently constructed in both calculational approaches. It also indicates that there is no significant influence of the central atom (Li or Na) on the electronic structure. Using NRLMOL we calculated the total charge within a sphere with a radius of 2.97 Bohr units around the Na atom. There are already 10.98 electrons contained inside the sphere and zero spin polarization which clearly indicates that the Na atom does not play a role for the electronic structure of the ferric wheel.

Notably, both methods find the local magnetic moments on Fe sites are very close to 4 $`\mu _\text{B}`$ and *not* to 5 $`\mu _\text{B}`$ as is generally assumed, based on the above mentioned magnetization data. The maximal magnetization $`S`$=5/2 of the Fe atom corresponds to a Fe(III)-ion with in 3$`d_{}^5d_{}^0`$ configuration. Our first-principles calculations suggest a somewhat different picture: the minority-spin DOS has a non-zero occupation due to the hybridization (chemical bonding) of Fe$`3d`$ with O$`2p`$ states. However, the magnetic polarization in the organic ligand which provides the octahedral coordination for the iron atoms, due to Fe is substantial, the most pronounced effect being on the apical oxygen atom (which is not participating in the bonding to the next Fe neighbor). Taken together with the (smaller) polarization of the bridging oxygen atoms and magnetization at the nitrogen site, the distributed magnetic moment *per* Fe atom yields 5 $`\mu _\text{B}`$, recovering the agreement with the magnetization results.

An important consequence is that the charge state of iron is not Fe(III) but more close to Fe(II), according to our calculations. Moreover, the distributed magnetic moment behaves like a rigid one, in a sense that it can be inverted, following a spin flip on a Fe site. This conclusion has been derived from our analysis of other magnetic configurations, done with Siesta . The local DOS does not change considerably when switching from AFM to FM configuration – only the HOMO/LUMO gap becomes less pronounced, and a slight ferro-magnetic shift between the two spin bands appears. Further, the change of total energies in the sequence of FSM, from the maximal spin moment of 30 $`\mu _\text{B}`$ to one Fe spin inverted (20 $`\mu _\text{B}`$), and then to two (second-neighbors) Fe spins inverted is linear, supporting the conclusion that the system maps well onto the Heisenberg model. The corresponding exchange parameter, as estimated over both 30$``$20 and 20$``$10 $`\mu _\text{B}`$ flips, is around $``$80 K. This yields a correct sign and correct order of magnitude, relatively to experimental estimates of Ref. . However, the magnitude is overestimated by a factor of $``$4. The exchange parameters $`J`$ depend on the spatial overlap of the $`d`$-orbitals on different Fe-sites. It is well known that the $`d`$-orbitals within DFT are not localized enough compared to experiment, consequently the $`J`$ values will be overestimated. There are two main reasons for this shortcoming. First, possible on-site correlations as known from atomic physics are underestimated in case of “conventional” DFT. Second, DFT is not free from spurious self-interactions due to the replacement of the point-like electrons by corresponding densities. Bringing in the atomic physics in the form of LDA+$`U`$ (adding a local orbital dependent atomic Coulomb interaction parameter $`U`$ to DFT, ) or self-interaction corrections (SIC)  will improve the results by lowering the $`d`$-orbitals in energy and therefore localizing them stronger. SIC only affects occupied states, whereas LDA+$`U`$ plunges the occupied $`d`$-states and shifts the unoccupied ones to higher energies. By increasing, on the average, the magnetic excitation energy across the spin majority-minority gap, both mechanisms help to effectively reduce the magnitude of $`J`$. In what regards the LDA+$`U`$ approach, this has been shown by Boukhvalov *et al.* for another molecular magnet, Mn12 .

It is worth noting that the inversion of the third Fe spin that yields the AFM configuration does not follow the same linear trend. The origin of this is not clear at the moment and should be tested in further high-precision calculations. We recently came aware of calculations by E. Ruiz on similar systems , done also with the use of Siesta. According to them, the $`J`$ values extracted just from the FM–AFM total energy difference along the procedure described in Ref.  yielded the $`J`$ values of $``$ $``$24 K (Li system) and $``$ $``$45 K (Na-system). A similar estimate from FM to AFM transition with our total energy data gives $``$38 K, i.e., such estimates seem to be in better agreement with the reported experimental data. It remains to be tested, however, whether the FM to AFM transition can be reasonably well fitted to the Heisenberg model with a single (nearest-neighbors) exchange parameter, in view of above mentioned nonlinear dependence of the total energy on the total spin. In V15, another molecular magnet, a complex system of six exchange parameters was found to be necessary  in an otherwise conceptually similar fit to the Heisenberg model.

A clear visualization of the above discussed delocalized (or, rather, distributed) magnetic moment associated with the Fe atom comes from the map of spin density, obtained from the NRLMOL calculation (Fig. 3). One should take into account that the volume enclosed by the iso-surfaces is not directly correlated to the total moment at the site. One sees moreover an absence of magnetization on carbon and hydrogen sites. The fact that the magnetization is noticeable and changes its sign when passing through bridge oxygen atoms emphasizes the failure of methods depending the spherical averaging of atom-centered potentials. From the side of experiment, it would be interesting to probe the spin polarization on oxygen and nitrogen atoms.

## Conclusion

We presented a study of the electronic structure of a six-membered ferric wheel with Li or Na as a central atom. The local densities of states on the iron and neighboring sites clearly indicate that the magnetic moment is distributed over several sites. The local moment at the iron site is still the largest one with about 4 $`\mu _\text{B}`$, although indicating that the iron in the ferric wheels is closer to Fe(II) than Fe(III) as generally expected. The local spin of $`S`$=5/2 per iron site as deduced from magnetization measurements is reproduced if one takes the moments on the neighbor atoms into account. The largest moment is found on the apical oxygen atom, followed by smaller moments on nitrogen and the bridging oxygen atoms. This picture has been confirmed by a plot of the spin density map. A further goal of the present study has been a comparative test of two different DFT implementations. In most cases the local (energy-resolved) DOS and (spatially resolved) magnetization density show a very good quantitative agreement. The most pronounced deviations are related to nitrogen, indicating less localized magnetic moment there and thus higher sensitivity of numbers to the exact normalization of DOS.

In order to explore the range of systems where the presented first-principles methods give reliable results further studies on more systems are required. In particular studies of experimentally well characterized classes of compounds where only the ligands are changed in a controlled way but not the magnetic core would reveal if our calculational approach can reproduce the observed changes in the magnetic behavior. Our current results make us very confident in the predictive power of the presented methods. This should allow for a microscopic understanding based on the electronic structure of single molecule magnets.

## Acknowledgments

The authors thank the Deutsche Forschungsgemeinschaft for financial support (Priority Program “Molecular Magnetism”) and appreciate useful discussions with Jürgen Schnack, Oliver Waldmann, Sorin Chiuzbăian, Manfred Neumann, Eliseo Ruiz and Mikhail Katsnelson. The crystal structure data have been kindly provided by the group of Prof. Dr. R. W. Saalfrank from the Institute for Organic Chemistry of the University Erlangen-Nürnberg.


text:
# ANTIBOUND AND RESONANT STATES IN HALO NUCLEI

## Abstract

An unified shell model scheme to evaluate simultaneously the contribution of antibound states and Gamow resonances to the nuclear halos is presented. The calculations, performed in the complex energy plane, are applied to the case of 11Li. It is found that 11Li may develop a resonant state excitation near the breakup threshold.

It is by now well established that halos in nuclei are produced by particles moving in single-particle states which extend far in space. The setting for this to occur requires neutron configurations consisting mainly of unbound $`s`$\- and $`p`$-waves since in this case no barrier large enough to trap the nucleons inside the nuclear core will be present. It is the pairing interaction that holds the escaping neutrons together, granting them their Borromean character .

Stated in this simplified fashion the explanation of halos seems perhaps obvious. However it took many years of intense work, experimental as well as theoretical, to reach the understanding of halos. It is not our intention here to quote the vast amount of papers written on this subject. A recent review, including abundant references, can be found in Ref. . For our purpose it is enough to point out that for the description of 11Li, which is the typical halo nucleus, the two relevant single-particle states, as specified by the experimental spectrum of 10Li , consist of a low- lying $`s_{1/2}`$ state, which probably is an antibound (or virtual) state at about -50 keV , and a $`p_{1/2}`$-resonance at about 540 keV. A second $`p`$ resonance at around 0.250 MeV may also exist . The two-body correlations induce a bound ground state in 11Li at about -0.295 MeV. Notice that we are here using the shell-model language, where the core (9Li) is considered as inert, the single-particle states are given by 10Li and the two-body nucleus is 11Li .

The existence of a very low-lying virtual $`s`$-state in 10Li has important consequences for the correlations developed in 11Li . As discussed below, an antibound state close to the continuum threshold enhances the localisation of the low-lying scattering states. Therefore the $`s`$-wave content of the ground state of 11Li is also increased, reaching the corresponding (large) experimental value. Moreover, the antibound state in 10Li can affect the excited spectrum of 11Li as well as the ground state. Since the dineutron system has also a virtual state close to threshold, in each of the two-body subsystems of 11Li there is an antibound state near zero energy. Thus 11Li fulfils the Efimov conditions  and consequently it may form low-lying excited states close to the breakup threshold. The posibility of such excitations was suggested in Ref.  in connection to the momentum distribution of the fragments of 11Li.

In all the studies done until now the role of the antibound $`s`$-state in 11Li was discussed only indirectly, i. e. through the associated scattering length. One of the aims of this letter is to present for the first time a formalism in which antibound states are treated as ordinary single-particle states. That is, the antibound states will form part of the single-particle representation in the same fashion as bound states do in standard shell model calculations. For this we will use an extended version of the shell model on the complex energy plane (CXSM) . Within this formalism we will investigate the existence of low-lying excited resonant states of Efimov type, for which there are no accurate calculations yet.

The CXSM is a powerful method to analyse the influence of the single-particle resonances upon the two-particle correlated states as well as to understand the developing of bound states by the two-body force acting upon states embedded in the continuum. In its initial form  the representation corresponding to the CXSM consists of Gamow resonances, bound states and complex scattering states. These, together with the antibound states, are poles of the one-particle Green function and, therefore, we will refer to them as ”poles”.

The space spanned by this basis (Berggren space)  is determined by a metric (Berggren metric) which is non-Hermitian. Since on the real energy axis the Berggren space coincides with the Hilbert space, any physical quantity evaluated within the CXSM coincides with the corresponding one evaluated within the standard shell model. In this sense one can assert that the CXSM is a generalization of the shell-model and, as in this model, the many-body energies and corresponding wavefunctions can readily be evaluated and analysed.

The CXSM can be easily extended to take into account the antibound states. The corresponding single-particle representation contains the antibound states on the same footing as the other discrete elements, i. e. bound states and Gamow resonances. We will show how to construct this representation for the case of 11Li.

In the first step of the calculation one evaluates the single-particle states of the unbound nucleus 10Li. As in Ref. , for the central field we choose a Woods-Saxon potential with different depths for even and odd orbital angular momenta $`l`$. The single-particle states as well as the scattering waves will be evaluated by using the high precision piecewise perturbation method . Since for the resonant state $`0p_{1/2}`$ the experimental data are uncertain we will perform two calculations corresponding to the energies 200 keV and 500 keV. At the same time, this will allow us to assess the influence of the $`p_{1/2}`$ resonance upon the halo formation. For clarity of presentation we will start studying the 200 keV case and afterwards the 500 keV case will be presented. At the end we will compare the results and discuss the similarities and differences between the two cases.

For the 200 keV case we use the Woods-Saxon potential given by $`a`$ = 0.67 fm, $`r_0`$ = 1.27 fm, $`V_0`$ = 50 (36.9) MeV and $`V_{so}`$ = 16.5 (12.624) MeV for $`l`$ even (odd). With these parameters we found the single-particle bound states $`0s_{1/2}`$ at -23.278 MeV and $`0p_{3/2}`$ at -2.589 MeV forming the 9Li core. The valence poles are the low lying resonances $`0p_{1/2}`$ at (0.195,-0.047) MeV and $`0d_{5/2}`$ at (2.731,-0.545) MeV and the wide resonance $`0d_{3/2}`$ at (6.458,-5.003) MeV. Besides, the state $`1s_{1/2}`$ appears as an antiboud state*** The principal quantum number $`n`$ labelling the single-particle states indicates that the corresponding wave functions are localized in a region inside the nucleus and that its real part has in that region $`n`$ nodes, excluding the origin. at -0.050 MeV. We thus reproduce the experimental single-particle energies giving from the very beginning unequivocal endorsement to the low lying $`s`$-state as due to an antibound state.

We also found other resonances at high energies. However we include in the basis single-particle states lying up to 10 MeV of excitation energy only. We found that expanding the basis from this limit does not produce any effect upon the calculation up to the six digits of precision that we require.

The non-resonant continuum in the CXSM representation is given by the points in the contour embracing the resonances . One can use different contours for different partial waves, but it is important to keep in mind that the corresponding Berggren space and therefore the calculated quantities on the complex energy plane will depend on the contours. Only on the real energy axis the calculated quantities do not depend upon the shapes of the contours. This property will allow us to check the precision of our calculation. That is, the calculation of physical quantities by using any contour should coincide with those evaluated on the real energy axis, i. e. by means of the continuum shell-model . Let us stress again that differences between the CXSM and the standard shell-model appear only when evaluating resonances .

In the calculations presented below we will use the contours shown in Fig. 1. One sees in this figure that for the partial waves containing a Gamow resonance we take the conventional CXSM contour defined by the vertices B  (see also Ref. ). Instead, for the antibound state we introduce the new type of countour defined by the vertices A. The general requirements which such contours should fulfil are discussed in Ref. . With proper values for the energies corresponding to the vertices of these contours one can readily include in the basis the three resonances and the antibound state mentioned above. Notice that in our basis no single-particle bound state is present.

The next step in the calculation is to adopt a residual interaction. We will use a separable force and, therefore, the Hamiltonian matrix reduces to a dispersion relation (for details see Ref. ). To evaluate the ground state of 11Li we adjust the strength $`G`$ of the separable interaction to reproduce the corresponding energy, i. e. -295 keV. We thus obtained G = 0.00194 MeV.

With the mean field and the two-body interaction thus established we evaluated the ground state wave function. First we performed the calculations by choosing the real energy as a contour. In this case the wave function is spread over many components. The largest of these components correspond to configurations $`p_{1/2}p_{1/2}`$ lying close to 400 keV (i. e. about twice the energy of the $`0p_{1/2}`$ resonance) and $`s_{1/2}s_{1/2}`$ lying close to threshold (i. e. close to twice the energy of the antibound state). The wave function consists of 44 % $`s`$-states, 48 % $`p`$-states and 8 % $`d`$-states, as expected .

A remarkable feature of the calculation is that the antibound state exerts such a strong effect upon the two-particle wavefunction. From a CXSM point of view this is because the energy corresponding to the configuration $`(1s_{1/2})^2`$ is very close (in the complex energy plane) to the two-particle energy. But it can also be understood from a continuum shell-model point of view. That is, the radial scattering wave function with energies $`E=\mathrm{}^2k^2/2\mu `$ ($`k`$ real and positive) close to a bound or antibound state with energy $`E_0`$ ($`k_0=\pm i|k_0|`$, respectively) lying near threshold can be written as ,

$$_l(kr)\sqrt{\frac{2k|k_0|}{k^2+|k_0|^2}}_l(|k_0|r)$$
(1)
which shows that close to threshold the radial dependence of the scattering wave functions is independent of the energy, except the square root factor. This factor is maximum for $`k=|k_0|`$. Therefore, the matrix elements of the two-body interaction are large within a relative large interval close to threshold. This induces large components of the two-body wave function in that energy region. Notice that it does not matter whether the S-matrix pole $`E_0`$ corresponds to a bound or to an antibound state. The effect is exactly the same.

This feature can be seen in Fig. 2, where we give the localization of the scattering states, $`L(E)`$, defined as

$$L(E)=_0^{1.2R_N}_{l=0}^2(kr)r^2𝑑r$$
(2)
where $`R_N`$ = 2.6 fm is the core nuclear radius.

One sees in Fig. 2 a strong increase of the localization as the energies of the poles approaches threshold. This propery is also responsible for a similar increase in the elastic cross section, as seen in Figs. 12 and 13 of Ref. . Therefore one expects large elastic cross sections at low energies in halo nuclei, particularly in 11Li.

So far we have shown the advantages of the CXSM to evaluate the effect of the antibound state on already known properties of the ground state of the halo nucleus 11Li. However, the transparency of the method becomes essential in the search for other physically meaningful two-particle states in the continuum. In what follows we discuss the problem of low-lying $`0^+`$ excitations and their influence upon the neutron halo.

The possibility of low-lying excitations in 11Li, with energies below the first experimentaly measured excitation at 1.2 MeV, was considered in Ref.. In this reference it is argued that a low-lying narrow resonance at around $`0.20.4`$ MeV would be consistent with the momentum distribution of the 11Li fragments. To our knowleadge the existence of such a low-lying resonance in 11Li was not investigated in microscopical calculations so far.

Within the CXSM two-particle resonances are easy to calculate since they appears as a result of the diagonalization of the Hamiltonian (which in our case reduces to the solution of the dispersion relation) in the complex energy plane. Thus we found that the first excited state (i. e. the state $`0_2^+`$) appears at the complex energy (0.202,-0.137) MeV. The corresponding wave function consists of nearly 100 % $`p`$-states, with a small admixture of $`s`$-states.

It is interesting to analyse how this state is built up by the two-body interaction starting from the zeroth-order configuration $`(0p_{1/2})^2`$. For this we increased the interaction gradually starting from G=0, as seen in Fig. 3. As the attractive interaction increases the resonance becomes narrower and approaches threshold, as expected from perturbation theory. However, a point is reached where continuum configurations become important and the resonance widens. This happens at G=0.0005 MeV in the figure. Up to this point the resonance is a purely $`(0p_{1/2})^2`$ state and, therefore, it is localized inside the nucleus. That is, it is a physically meaningful resonance. But from here on other configurations become important. These configurations are overwhelmingly those where one of the particles moves in the continuum and the other in the resonance $`0p_{1/2}`$. More specifically, in the two-particle resonace that we are studying the most important states in the continuum are those corresponding to $`p_{1/2}`$ waves with energies corresponding to the points on the segment $`B_1B_2`$ in Fig. 1. At $`G=G_0`$=0.00194 MeV, corresponding to the G-value fitting the ground state, there is a strong mixing with the continuum configurations. Increasing $`G`$ farther the configuration $`(0p_{1/2})^2`$ looses its importance, the resonance is split in a number of pieces and eventually dissolves into the continuum. However, since these configurations virtually include only $`p`$-waves the wavefunctions still consist of only $`p`$-states.

The analysis that we have done so far is based upon the assumption that the $`p`$ resonance in 10Li is located at 200 keV. To see the influence of this resonance on the structure of the halo we will now analyse the ground and the excited states of 11Li by using the 500 keV case. For this we adopted the Woods-Saxon depth $`V_0`$ = 35.366 MeV for odd $`l`$-values, keeping all other parameters as before. We thus obtained the energy (0.470,-0.197) MeV for the state $`0p_{1/2}`$. The state $`0p_{3/2}`$, belonging to the core, is found now at -2.016 MeV. The other odd l-value poles lie beyond the range of energies included here. But, nevertheless, we have checked that they do not affect the results.

We kept the two-particle interaction used in the previous case, except that the strength necessary to adjust the energy of 11Li(gs) is now G = 0.00694 MeV.

As before, we found that on the real energy axis the ground state wave function is spread in many components. The largest of them lie close to threshold for the configurations $`s_{1/2}s_{1/2}`$ and around 1 MeV for the $`p_{1/2}p_{1/2}`$ configurations. The wave function consists of 49 % $`s`$-states, 39 % $`p`$-states and 12 % $`d`$-states, which is also within the range of accepted values .

Since the position of the $`p_{1/2}`$ pole seems likely to correspond to the present 500 keV case , we will analyse here the effects of the antibound and the Gamow poles upon the ground state of 11Li by using the contours of Fig. 1. We therefore present in Table I the contribution of different configurations to that ground state. The corresponding complex amplitudes depend on the chosen contours and have no direct physical meaning. But the total content of a given partial wave in the bound ground state wave function, which is a physical quantity, does not depend upon the chosen contour. From Table I we can see that for the $`p`$ and $`d`$ waves the configurations are built mainly on the corresponding Gamow resonances. The situation is different for the $`s`$-wave since apart from the configurations built upon the antibound state there is also an important contribution coming from the complex scattering states. This contribution is given mainly by those $`s`$ scattering states located on the segments $`(0,0)A_1`$ and $`A_1A_2`$ of Fig. 1, which are the closest to the antibound state.

Up to this point there is not much difference between the 200 and the 500 keV cases, which may explain why various studies of the halo structure of 11Li(gs) with the common feature of having low-lying s- and p-states, provide similar results . This is because the wave function of 11Li(gs) is mainly controlled by low-spin single-particle states lying close to the continuum threshold. The exact positions of the resonances do not influence the wave function very much. Even the two-particle interaction (if reasonable) is not so important in the evaluation of the properties of 11Li(gs) since this state is a Cooper pair  and as such is mainly induced by the Pauli principle acting upon the valence particles and those in the core. However, the two-body interaction as well as the position of the single-particle poles may have a fundamental importance to determine the physically meaningful excited states. The states arising from the particles moving in the continuum are not localized inside the nucleus and, therefore, will be weakly affected or not affected at all by the interaction. This can be seen in Fig. 4, where we present the evolution of the state $`0_2^+`$ as a function of the strength G. The $`0p_{1/2}`$ resonance is now wider and higher in energy than before. As a result, the point corresponding to G=0 in the figure is closer to points coming from the continuum contour. Yet, these continuum states do not seem to affect the resonance as G increases. That is, the behaviour of the resonance as G is varied is very similar to that in Fig. 3 as well as to resonances in non-halo nuclei . The reason for this is that those points, label ”s-states” in the figure, correspond to configurations of the type $`cs_{1/2}1s_{1/2}`$, where ”c” labels points in the segments $`(0,0)A_1`$ and $`A_1A_2`$ of Fig. 1. The overlap between these configurations and the mainly $`(0p_{1/2})^2`$ configuration of the resonance is small. Only at large values of G (above G=0.004 MeV in the figure), the resonance starts to feel the presence of the continuum states. The remarkable feature of the figure is the sudden turning down of the curve corresponding to the physical resonance at G=0.005 MeV. As G increases in this region the continuum plays a mounting role. As in the 200 keV case above, the most important of the continuum configurations are those in which one particle moves in the continuum and the other in the $`0p_{1/2}`$ resonance. There could be many comparatively large configurations of these type and it would not be useful to give all of them. More instructive is to show their contribution to the normalization of the wave function in this case where, in contrast to the ground state case of Table I, the zeroth order energies are not very close to the energy of the state $`0_2^+`$. We thus define $`S(l)`$ as the sum of the squares of the amplitudes corresponding to configurations where at least one of the two particles moves in continuum states. On the real energy axis the sum of $`S(l=1)`$ and $`X^2((0p_{1/2})^2)`$, where $`X`$ is the wave function amplitude, is the probability of the p-content of the wave function. This is a quantity that we evaluated above for the ground state. The dependence of these quantities upon G corresponding to the state $`0_2^+`$ is shown in Table II. Since we are studying states with complex energies the numbers $`S`$ as well as $`X^2`$ are complex in this Table. Moreover, their absolute values could be larger than 1 although the sum of all possible l-contributions is normalized to (1,0). This is a good example of the non-Hermitian character of the Berggren metric.

One sees in this Table that the two-particle resonance starts to mix with the continuum at G= 0.005 MeV and at G=G0 it is composed mainly of continuum configurations. Therefore, at this point it has already lost its localization features. It has become a part of the continuum background.

We are now in a position to recognize other systems where halos may be present. We thus looked for nuclei which may be considered shell-model cores lying on the neutron drip line with low-lying single-particle resonances carrying low-spin. Following the trend of single-particle states in the relativistic mean field calculations we found that Z=20, N=50 may be such a core. In order to simulate the order of the single particle states given by the relativistic calculations we used a Wood-Saxon potential defined by $`a`$ = 0.67 fm, $`r_0`$ = 1.27 fm, $`V_0`$ = 39 MeV and $`V_{so}`$ = 22 MeV. With this potential the antibound $`2s_{1/2}`$ state (note that $`n`$=2) appears again at -0.050 MeV. But now the next valence shells are $`1d_{5/2}`$ at (0.469,-0.048) MeV, $`1d_{3/2}`$ at (2.080,-1.525) MeV, $`0g_{7/2}`$ at (6.739,-0.738) MeV and $`0h_{11/2}`$ at (5.344,-0.102) MeV. The states in the core are ordered as usual. As expected, the highest of these is the state $`0g_{9/2}`$, lying at -2.276 MeV. Using the same separable interaction as before and assuming again that the ground state of 72Ca lies at -295 keV, we obtained for the strength of the interaction the value $`G_0`$=0.00174 MeV. Close to this $`G_0`$ value we found also a low-lying two-particle resonance with the energy of about (0.550, -0.350) MeV. The behaviour of this $`0_2^+`$ resonance as a function of $`G`$ is very similar to the 500 keV case of Fig. 4. The discussion performed there is also valid here and, therefore, we will not analysed this rather academic case farther. But it is important to point out that in this and the other cases presented here, we have been careful to choose contours that leave the region around the two-particle resonances in the complex energy plane free of continuum configurations. We thus established an ”allowed” region . Otherwise the two-particle resonance would be embedded in a see of continuum states, making the calculations difficult and the evaluated quantities unreliable.

In conclusion, we have presented in this letter a formalism that treats the antibound states exactly in the same fashion as bound states and Gamow resonances in the framework of a shell model scheme in the complex energy plane. We found, as expected, that antibound states lying close to the continuum threshold are of a fundamental importance to build up the halo. But we also found that an excited low-lying two-particle resonance may exist in halo nuclei. For the case of 11Li this low-lying resonant appears in the energy range of $`0.20.5`$ MeV. This energy range is the same as the one suggested in Ref.  in connection to the momentum distribution of 11Li fragments. However, the calculations exibit a very drastic change in the structure of the resonant excited state when the strength of the force is approaching the value used for the determination of the ground state. Do to this it is rather difficult to conclude whether this state is a physical resonance or not. Further efforts directed to the measuring of the eventual low-lying resonance excitations in 11Li would be essential in order to settle this issue.

###### Acknowledgements.

This work has been supported by FOMEC and Fundación Antorcha (Argentina), by the Hungarian OTKA fund Nos. T37991 and T29003 and by the Swedish Foundation for International Cooperation in Research and Higher Education (STINT).