text:
# Many-body approach to proton emission and the role of spectroscopic factors

## I Introduction

One of the classic problems in quantum mechanics is that of tunneling through a classically forbidden region or, more specifically, the decay of a quasi-stationary state to the continuum. In nuclear physics, this manifests itself in the processes of $`\alpha `$-decay in heavy nuclei and proton emission by proton drip-line nuclei. Of particular current interest are the lifetimes of proton emitters, especially in the lighter region of the nuclear chart, and the implications of this in nuclear astrophysics.

Over the years, a number of different theoretical approaches have been used to describe the decay process in nuclear physics, either by means of perturbation theory of decaying states or by time reverse study of resonance states via scattering theory Jackson and Rhoades-Brown (1977); Gurvitz et al. ; Mang (1964); Thomas (1954); Arima and Yoshida (1974); Vogt (1996). Some authors solve the time-dependent problem while others use a stationary picture and make use of approximation methods such as the distorted-wave Born approximation or the semi-classical Wentzel-Kramers-Brillouin approach to evaluate the width Åberg et al. (1997). Other more accurate methods, such as R-matrix theory, are sometimes very sensitive to the channel radius giving dramatic variation in the calculated widths Mang (1964); Thomas (1954); Arima and Yoshida (1974); Vogt (1996). The method of Gurvitz and Kalbermann Gurvitz and Kalbermann (1987); Gurvitz (1988); Gurvitz et al. , also known as the two-potential approach (TPA) Jackson and Rhoades-Brown (1977); Åberg et al. (1997), is based on splitting the barrier potential into an interior and an exterior components. The inner potential binds the particle, which can then be described by a bound eigenstate of the relative Hamiltonian, while the outer potential acts as a perturbation that converts it into a quasi-stationary state (a wave packet), which can decay.

An important shortcoming of all the above approaches, however, and in common with the descriptions of so many nuclear processes, is the approximate treatment of the many-body structure effects. In most descriptions of the proton-emission process the initial ($`A+1`$)-body wave function is written as a product of an $`A`$-body wave function, describing the daughter nucleus, and the proton’s single-particle wave function. The decay width is then written in the form of a single-particle width multiplied by a spectroscopic factor, which contains the many-body information of the system. This procedure, however, makes various assumptions about the relationship between the many-body problem and the effective one-body problem that have to be tested. In this work we consider the TPA of Gurvitz and Kalbermann and extend it to properly account for the many-body correlations.

The standard reduction from a many-body problem to an effective one-body picture has been revisited in a recent study of radiative proton capture Escher et al. (2001). The work focused on one-body overlap functions and their associated equations of motion. The one-body overlap functions are obtained by integrating the product of the wave functions for an ($`A+1`$)-body system and its $`A`$-body subsystem over the coordinates of the latter. While the overlap functions are unambiguously defined, it was demonstrated in Ref. Escher and Jennings (2002) that useful ‘auxiliary’ one-body functions can be defined in several different ways. Naturally, the associated equations of motion differ for the three approaches considered in Ref. Escher and Jennings (2002). In the current work, we derive expressions for the proton decay width using two of the three approaches mentioned. The resulting decay widths have formally the same structure in both approaches, but the overall normalization factors differ. Only one of these normalization factors can be interpreted as the square root of a spectroscopic factor. This has consequences for the interpretation of experimental results and in particular for the determination of spectroscopic factors from decay widths.

We start, in Sec. II, with the time-dependent Schrödinger equation and follow the standard theory of decaying states Goldberger and Watson (1964). This method is briefly compared with the scattering approach to decay problems. In Secs. III and IV we use projection operator techniques and perturbation theory to derive an expression for the decay width in terms of the imaginary part of the pole in the Green’s function matrix element. In Sec. V, we describe the two alternative routes for reducing the many-body problem to an effective one-body problem. We compare the resulting two expressions for the width in Sec. VI. An alternative expression for the decay width is given in Sec. VII which shows more clearly the relation of the width to the spectroscopic factor.

## II Proton emission formalism

In this section we begin to recast the formalism of Gurvitz and Kalbermann in a form which is more convenient for our purposes. The approach of Refs. Gurvitz and Kalbermann (1987); Gurvitz (1988) starts with a square-integrable wave function, $`|\psi _0`$, which corresponds to the quasi-bound nucleus whose decay we are interested in. The initial wave function is close to the resonance state in the nuclear interior but decays rather than oscillates in the exterior region. This wave function cannot be an eigenstate of the full Hamiltonian or it would have a trivial time dependence and no decay would take place. Taking $`|\psi _0`$ as the wave function at $`t=0`$, we follow its time evolution using the time-dependent Schrödinger equation:

$`i\mathrm{}{\displaystyle \frac{}{t}}|\psi (t)=H|\psi (t).`$ (1)
This initial value problem is solved by using the one-sided Fourier transform (sometimes called a Laplace transform). We obtain

$`i\mathrm{}|\psi _0+E|\stackrel{~}{\psi }(E)`$ $`=`$ $`H|\stackrel{~}{\psi }(E),`$ (2)
where $`|\stackrel{~}{\psi }(E)=_0^{\mathrm{}}𝑑te^{i(E+iϵ)t/\mathrm{}}|\psi (t)`$ and $`ϵ`$ is a positive infinitesimal real number. Solving for $`|\stackrel{~}{\psi }(E)`$ we have

$`|\stackrel{~}{\psi }(E)`$ $`=`$ $`{\displaystyle \frac{i\mathrm{}}{EH+iϵ}}|\psi _0.`$ (3)
The probability amplitude for the nucleus remaining in the initial state after a time $`t`$ is given by

$`\psi _0|\psi (t)`$ $`=`$ $`{\displaystyle \frac{i}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑Ee^{iEt/\mathrm{}}\psi _0|{\displaystyle \frac{1}{EH+iϵ}}|\psi _0,`$ (4)
which is obtained by taking the matrix element of the previous equation with $`\psi _0|`$ and carrying out the inverse Fourier transform. This is the Fourier transform of one particular matrix element of the many-body Green’s function. In general the Lehmann representation of the latter contains contributions from many poles. However, the overlap $`\psi _0|\psi (t)`$ takes a simple form if the right-hand side of Eq. (4) is dominated by the contribution of only one pole. In this case the decay rate can be extracted from the imaginary part of the pole location and we obtain a simple exponential function which describes the decay of the initial state. Thus the initial state $`|\psi _0`$ should be chosen to minimize the contributions from other poles.

The same result can be found by considering the related scattering problem. A quasi-bound state can also be thought of as a resonance in the scattering amplitude. We consider scattering of a proton off the $`(A1)`$-body system and we define a t-matrix, $`T`$, through the equation

$`Ta^{}(k)|\psi _{A1}=V|\psi _A,`$ (5)
where $`a^{}(k)`$ is the creation operator for a particle with momentum $`k=\sqrt{2mE/\mathrm{}^2}`$ and $`m`$ is the reduced mass of the system. By standard techniques it can be shown that

$`T`$ $`=`$ $`V+VG_0T`$ (6)
$`=`$ $`V+VGV,`$ (7)
where $`G_0=1/(EH_0+iϵ)`$, $`G=1/(EH+iϵ)`$, and $`H_0`$ describes the free motion of the ejected particle with respect to the final state of the $`(A1)`$-body system, i.e., $`H_0a^{}(k)|\psi _{A1}=Ea^{}(k)|\psi _{A1}`$ and $`H=H_0+V`$. We observe that the poles of the t-matrix are given by the poles of $`V`$ and $`G`$. Assuming $`V`$ has no nearby pole we see that the poles of $`T`$ are just those of the many-body Green’s function. The width of the state is then given by the imaginary part of the pole location as in the previous case.

In principle, the potential $`V`$ is the sum of the interactions of the $`A`$-th particle with each of the particles in the $`(A1)`$-body system. In nuclear physics this is generally approximated by a nucleon-nucleus optical potential. This approach, however, has the disadvantage of loosing track of the Pauli exchange correlations and other many-body effects. For bound states one of the principal effects is included through the use of the spectroscopic factor. However for scattering states there is, strictly speaking, no spectroscopic factor. Using Eq. (4) as the starting point we can formally take the many-body effects into account while deriving an effective one-body equation.

## III Projection operator formalism

The main ingredient of Eq. (4) is the matrix element of the Green’s function, $`M=\psi _0|\frac{1}{EH+iϵ}|\psi _0`$. Its expression can be simplified as done in Ref. Escher and Jennings (2002) by using the projection operator formalism. We define a projection operator $`P=|\psi _0\psi _0|`$ and the complementary operator $`Q=1P`$. The equation for the Green’s function can be written as

$`(EH){\displaystyle \frac{1}{EH+iϵ}}=(EH)G=1.`$ (8)
Acting on the left by $`P`$ or $`Q`$ and on the right by $`P`$ we get the two equations

$`P(EH)(P+Q)GP`$ $`=`$ $`P,`$
$`Q(EH)(P+Q)GP`$ $`=`$ $`0.`$ (9)
Solving the second equation for $`QGP`$ and substituting this into the first equation we have

$`\left(EPHPPHQ{\displaystyle \frac{1}{EQHQ+iϵ}}QHP\right)PGP=P.`$ (10)
Taking the matrix element of this equation with $`|\psi _0`$ and using the explicit form for $`P`$ we find

$`\left[E\psi _0|\left(H+HQ{\displaystyle \frac{1}{EQHQ+iϵ}}QH\right)|\psi _0\right]\psi _0|G|\psi _0=1.`$ (11)
Since the left-hand side of this equation is a product of two terms, the poles of $`\psi _0|G|\psi _0`$ must coincide with the zeros of the multiplying factor. Thus the poles of the Green’s function are given by

$`E=\psi _0|\left(H+HQ{\displaystyle \frac{1}{EQHQ+iϵ}}QH\right)|\psi _0,`$ (12)
or, if we use a spectral representation of $`\frac{1}{EQHQ+iϵ}`$, by

$`E`$ $`=`$ $`\psi _0|H|\psi _0+{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑E^{}{\displaystyle \frac{|\psi _0|HQ|\zeta _E^{}|^2}{EE^{}+iϵ}},`$ (13)
where $`|\zeta _E^{}`$ is the solution of the equation

$`E|\zeta _E=QHQ|\zeta _E,`$ (14)
normalized according to $`\zeta _E^{}|\zeta _E=\delta (E^{}E)`$. Note that $`|\psi _0`$ is also an eigenstate of $`QHQ`$, but with energy $`E=0`$. This state is excluded from the sum (integral) in Eq. (13) by the projection operator. In fact the only role of $`Q`$ in this equation is to exclude the discrete state $`|\psi _0`$. The residue, $`R`$, of the pole of $`\psi _0|G|\psi _0`$ is given by

$`R=\left[1{\displaystyle \frac{d}{dE}}\psi _0|HQ{\displaystyle \frac{1}{EQHQ+iϵ}}QH|\psi _0\right]^1.`$ (15)
This follows from a Taylor series expansion of Eq. (11) in $`E`$. For the present problem $`R`$ is close to one since we will choose $`|\psi _0`$ to be an eigenvalue of $`H`$ inside the nucleus and different from an eigenstate only in a region where the wave function is exponentially damped. Where $`|\psi _0`$ is an eigenstate $`\psi _0|HQ`$ is zero so the contribution from the second term will be exponentially small.

More insight into Eq. (12) can be obtained by an alternative derivation. Consider the equation $`E=\psi _0|H|\psi _E/\psi _0|\psi _E`$ where $`|\psi _E`$ is an eigenstate of the full Hamiltonian, $`H`$. The wave function can be written as $`|\psi _E=(P+Q)|\psi _E`$. Following the standard Feshbach projection operator technique we can write $`Q|\psi _E=Q\frac{1}{EQHQ+iϵ}QHP|\psi _E`$. Using the explicit form of the projection operator $`P=|\psi _0\psi _0|`$ yields Eq. (12) immediately.

## IV Perturbative approximation

In general, Eq. (12) is highly non-linear and has many solutions. Every solution of this equation will give a pole of the Green’s function, however, not every pole of the Green’s function will necessarily be found by using this equation. For example, if $`|\psi _0`$ has a definite angular momentum only poles with that angular momentum can be found with that particular choice of $`|\psi _0`$.

For the problem at hand, namely proton emission, we are not interested in the complete complexity of the Green’s function. To obtain the decay width, only the imaginary part of the pole location of the nearby pole is relevant. We define an Hermitian Hamiltonian $`H_0`$ such that $`H_0|\psi _0=E_0|\psi _0`$ with $`H=H_0+\delta H`$. Expanding the right hand side of Eq. (12) about $`E_0`$ and neglecting terms of order $`(EE_0)^2`$ we obtain

$`EE_0`$ $``$ $`\psi _0|\delta H|\psi _0+\psi _0|HQ{\displaystyle \frac{1}{E_0QHQ+iϵ}}QH|\psi _0`$ (16)
$`+(EE_0)\left[{\displaystyle \frac{d}{dE}}\psi _0|HQ{\displaystyle \frac{1}{EQHQ+iϵ}}QH|\psi _0\right]_{E=E_0}`$
$``$ $`R\left(\psi _0|\delta H|\psi _0+\psi _0|HQ{\displaystyle \frac{1}{E_0QHQ+iϵ}}QH|\psi _0\right),`$ (17)
where the last line has been obtained by solving for $`(EE_0)`$ in the first equation. This expression for the energy, Eq. (17), contains a factor of $`R`$, the residue of pole of $`\psi _0|G|\psi _0`$. Note that $`\psi _0|HQ=\psi _0|\delta HQ`$ implies that $`R1`$ is of the order of $`\delta H^2`$. Since it multiplies a factor of order $`\delta H`$, it will introduce terms of order higher than $`(EE_0)^2\delta H^2`$ and therefore it can be neglected. In the following, we thus take $`R=1`$.

Let us emphasize that $`E_0`$ is real since it is an the eigenvalue of the Hermitian operator $`H_0`$. The pole location, on the other hand, occurs at a complex energy. If the width of the state is large, for any reason not just proton emission, then $`EE_0`$ does not hold (since $`E_0`$ is real and the pole location has a large imaginary part) and the approximation fails. However, when the width is narrow then it will be possible to chose $`|\psi _0`$ or equivalently $`H_0`$ such that perturbation theory is applicable.

The equation for the energy is now linear in $`E`$ and has the form typical of perturbation theory. The explicit connection with perturbation theory can be made when the unperturbed Hamiltonian is taken to be $`H_{00}=H_0+Q\delta HQ`$, and not $`H_0`$ as one might have expected. This definition is not only necessary in order to cast Eq. (13) into the form of a perturbation expansion but at the same time eliminates the problem of non-compactness Gurvitz and Kalbermann (1987); Gurvitz (1988). The non-compactness arises since $`H_0`$ is chosen in such a manner that it does not go to zero at infinity. However both $`H`$ and $`H_{00}`$ go to zero at infinity.

The states $`|\zeta _E`$ occurring in the integral in Eq. (13) are eigenstates of $`H_{00}`$, as is $`|\psi _0`$. In contrast to Eq. (14) the state $`|\psi _0`$ now occurs with energy $`E_0`$ rather than $`0`$. We can also rewrite $`\psi _0|HQ|\zeta _E^{}`$ as $`\psi _0|(HH_0)Q|\zeta _E^{}`$ or $`\psi _0|(HH_{00})Q|\zeta _E^{}`$ since $`\psi _0|`$ is an eigenstate of both $`H_0`$ and $`H_{00}`$ or, equivalently, since $`Q`$ commutes with both $`H_0`$ and $`H_{00}`$. Using the previously defined $`\delta H=HH_0`$, Eq. (13) can be rewritten as

$`E`$ $``$ $`E_0+\psi _0|\delta H|\psi _0+{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑E^{}{\displaystyle \frac{|\psi _0|\delta HQ|\zeta _E^{}|^2}{E_0E^{}+iϵ}}`$ (18)
$``$ $`E_0+\psi _0|\delta H^{}|\psi _0+{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑E^{}{\displaystyle \frac{|\psi _0|\delta H^{}Q|\zeta _E^{}|^2}{E_0E^{}+iϵ}},`$ (19)
where

$`\delta H^{}=HH_{00}=P\delta H+\delta HPP\delta HP.`$ (20)
Equation (18) for $`E`$ is not in the form of perturbation theory since $`|\zeta _E`$ is an eigenfunction of $`H_{00}`$ and not $`H_0`$. However, Eq. (19) is. The only approximation made in deriving Eqs. (18) and (19) was the replacement of $`E`$ in the denominator on the right hand side by $`E_0`$, i.e., we assumed that perturbation theory is valid.

We emphasize that neither $`|\psi _0`$ nor $`|\zeta _E`$ is an eigenstate of the full Hamiltonian $`H`$ and our derivation depends crucially on this point. Instead, $`|\psi _0`$ is an eigenstate of $`H_0`$ and $`H_{00}`$ while $`|\zeta _E`$ is an eigenstate of $`QHQ`$ and $`H_{00}`$. The fact that both wave functions are eigenstates of $`H_{00}`$ allows us to interpret the right-hand side of Eq. (19) as the first few terms of a perturbative expansion.

## V Choice of $`H_0`$

The proper choice of $`H_0`$ was discussed in a very convincing and direct way by Gurvitz and Kalbermann Gurvitz and Kalbermann (1987); Gurvitz (1988) for the problem of a single particle in a local potential well. The reader is referred to those papers for specific details. The basic idea is to take an $`H_0`$ such that at infinity the potential goes to a finite value, larger than $`E_0`$, rather than to zero. This causes the state under consideration to be bound, but if the decay width is small this new state should be very close to the scattering state. We now generalize this concept to the case of $`A`$ interacting particles.

We consider the case where there is one open channel. Asymptotically the wave function describes a free proton and the bound $`(A1)`$-body system. Thus we take

$`\delta H={\displaystyle 𝑑𝐫a^{}(𝐫)|\mathrm{\Phi }_{A1}V(𝐫)\mathrm{\Phi }_{A1}|a(𝐫)},`$ (21)
where $`|\mathrm{\Phi }_{A1}`$ is the ground state of the $`(A1)`$-body system, $`𝐫`$ is the relative coordinate of the proton and the $`(A1)`$-body system. The integral covers the whole space. If $`V(𝐫)`$ is taken to be greater than $`E_0`$ outside the range of the nuclear potential, $`r_0`$, it will prevent the initial state $`|\psi _0`$ from decaying. At the same time, we want the perturbing potential $`V(𝐫)`$ to be zero inside the nucleus so that it does not modify the wave function in the interior. With this ansatz expression (18) for $`E`$ becomes

$`E`$ $``$ $`E_0{\displaystyle 𝑑𝐫\varphi _0^{}(𝐫)V(𝐫)\varphi _0(𝐫)}+{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑E^{}{\displaystyle \frac{|𝑑𝐫\varphi _0^{}(𝐫)V(𝐫)\varphi _E^{}(𝐫)|^2}{E_0E^{}+iϵ}},`$ (22)
where $`\varphi _0(𝐫)=\mathrm{\Phi }_{A1}|a(𝐫)|\psi _0`$ is a spectroscopic amplitude (i.e. a one-body overlap function involving bound many-body states) and $`\varphi _E^{}(𝐫)=\mathrm{\Phi }_{A1}|a(𝐫)|\zeta _E^{}`$ is an optical model wave function (cf. the discussion in Refs. Escher et al. (2001); Escher and Jennings (2002)). The projection operator $`Q`$ does not have to be included explicitly since the integral does not include the discrete state $`\varphi _0(𝐫)`$.

A miracle has occurred here. Due to the choice of $`\delta H`$, which is physically motivated, the expression has reduced to an effective one-body problem. All the many-body aspects of the problem are contained in the one-body overlap functions, $`\varphi _0(𝐫)`$ and $`\varphi _E^{}(𝐫)`$. We stress that the only approximation made so far is that second order perturbation theory was used to justify replacing $`E`$ with $`E_0`$ in Eq. (13). For states with narrow widths this should be acceptable. We have not assumed that $`|\psi _0`$ is a product state or made any other assumptions regarding its structure.

Equation (14) for $`|\zeta _E`$ can be written as $`E|\zeta _E=H_{00}|\zeta _E`$. In the neighborhood of $`E=E_0`$ the Hamiltonian in this equation for $`|\zeta _E^{}`$ can be approximated as

$`\delta H^{}`$ $``$ $`{\displaystyle 𝑑𝐫a^{}(𝐫)|\mathrm{\Phi }_{A1}V_{00}(𝐫)\mathrm{\Phi }_{A1}|a(𝐫)}.`$ (23)
For an appropriately chosen $`V_{00}(𝐫)`$, the approximation made here is the same as that of Eq. (2.15) in Ref. Gurvitz and Kalbermann (1987). The basic argument given there proceeds as follows: the eigenfunction of $`H`$ at resonance will be large when the $`A`$th particle is in the nuclear interior. However the state $`|\zeta _E^{}`$ will be small due to the projection operators in $`QHQ`$. It will look like the real state in exterior region but be suppressed in the interior. The form of $`V_{00}(𝐫)`$ will ensure this if we take it to be $`0`$ in the exterior region and repulsive (greater than $`E_0`$) in the interior region (see Fig. 3 in Ref. Gurvitz and Kalbermann (1987)). This is the opposite of what we did for $`V(𝐫)`$ which was large in the exterior and $`0`$ in the interior.

Despite its nice form, Eq. (22) is not useful until we specify how to calculate the functions $`\varphi _0(𝐫)`$ and $`\varphi _E^{}(𝐫)`$. Following Ref. Escher and Jennings (2002) we set up one-body equations for these functions using the Feshbach projection operator technique. We start with the exact scattering state $`|\mathrm{\Phi }_A`$ and write the equation-of-motion for the corresponding one-body overlap function (which in this case is the Feshbach optical model wave function) and see how it is modified when the perturbing Hamiltonian $`\delta H`$ is added. The one-body overlap function for the exact scattering state is given by Feshbach (1958, 1962); Block and Feshbach (1963); Escher and Jennings (2002)

$`E\varphi (𝐫)`$ $`=`$ $`{\displaystyle 𝑑𝐫^{}𝑑𝐫^{\prime \prime }\mathrm{\Phi }_{A1}|a(𝐫)\left(H+HQ_F\frac{1}{EQ_FHQ_F}Q_FH\right)a^{}(𝐫^{\prime \prime })|\mathrm{\Phi }_{A1}}`$ (24)
$`𝒩(𝐫^{\prime \prime },𝐫^{})^1\varphi (𝐫^{})`$
$`=`$ $`{\displaystyle 𝑑𝐫^{}(𝐫,𝐫^{})\varphi (𝐫^{})},`$ (25)
where $`P_F=𝑑𝐫𝑑𝐫^{}a^{}(𝐫)|\mathrm{\Phi }_{A1}𝒩(𝐫,𝐫^{})^1\mathrm{\Phi }_{A1}|a(𝐫^{})`$, $`Q_F=1P_F`$ and the norm operator is $`𝒩(𝐫,𝐫^{})=\mathrm{\Phi }_{A1}|a(𝐫)a^{}(𝐫^{})|\mathrm{\Phi }_{A1}`$. Since $`\delta HQ_F=Q_F\delta H=0`$, it follows that $`\varphi _0(𝐫)`$ satisfies the equation

$`E\varphi _0(𝐫)={\displaystyle 𝑑𝐫^{}\left[(𝐫,𝐫^{})+𝒩(𝐫,𝐫^{})V(𝐫^{})\right]\varphi _0(𝐫^{})}.`$ (26)
Note the explicit asymmetry of this Hamiltonian. Even in the absence of inelastic processes it is not Hermitian. An analogous potential in alpha emission has been strongly criticized Varga and Lovas (1991); Lovas et al. (1998) in the study of alpha emission. Two methods of addressing the problem were discussed in Ref. Escher and Jennings (2002). The first is to define a new amplitude $`\overline{\varphi }(𝐫)=𝑑𝐫^{}𝒩(𝐫,𝐫^{})^{1/2}\varphi (𝐫^{})`$ and take

$`\delta H={\displaystyle 𝑑𝐫𝑑𝐫^{}𝑑𝐫^{\prime \prime }a^{}(𝐫)|\mathrm{\Phi }_{A1}𝒩(𝐫,𝐫^{})^{1/2}\overline{V}(𝐫^{})𝒩(𝐫^{},𝐫^{\prime \prime })^{1/2}\mathrm{\Phi }_{A1}|a(𝐫^{\prime \prime })}.`$ (27)
The equation for $`\overline{\varphi }(𝐫)`$ is then

$`E\overline{\varphi }_0(𝐫)={\displaystyle 𝑑𝐫^{}\left[\overline{}(𝐫,𝐫^{})+\delta (𝐫,𝐫^{})\overline{V}(𝐫^{})\right]\overline{\varphi }_0(𝐫^{})},`$ (28)
where $`\overline{}(𝐫,𝐫^{})`$ is given in Ref. Escher and Jennings (2002) as

$`\overline{}(𝐫,𝐫^{})={\displaystyle 𝑑𝐫^{\prime \prime }𝑑𝐫^{\prime \prime \prime }𝒩(𝐫,𝐫^{\prime \prime })^{1/2}(𝐫^{\prime \prime },𝐫^{\prime \prime \prime })𝒩(𝐫^{\prime \prime \prime },𝐫^{})^{1/2}}.`$ (29)
A similar development holds for $`\varphi _E(𝐫)`$ and we obtain the corresponding equation

$`E\overline{\varphi }_E(𝐫)={\displaystyle 𝑑𝐫^{}\left[\overline{}(𝐫,𝐫^{})+\delta (𝐫,𝐫^{})\overline{V}_{00}(𝐫^{})\right]\overline{\varphi }_E(𝐫^{})}.`$ (30)
The equations we find for $`\overline{\varphi }_0(𝐫)`$ and $`\overline{\varphi }_E(𝐫)`$ are quite remarkable. As previously noted, we have reduced the many-body problem to an effective one-body equation in which all many-body effects are contained in $`\overline{}(𝐫,𝐫^{})`$. The perturbing potential is taken to be local in both cases. This is not really an additional approximation since the justifications given previously still hold.

Another way of getting a symmetric form for the additional potential is by using a projection operator given in Ref. Escher and Jennings (2002)

$`P_F={\displaystyle 𝑑𝐫[a(𝐫)+a^{}(𝐫)]|\mathrm{\Phi }_{A1}\mathrm{\Phi }_{A1}|[a(𝐫)+a^{}(𝐫)]},`$ (31)
for which Eq. (26) reduces to the Dyson equation of many-body theory. For time reversal invariant states the associated norm operator reduces to the unity operator. Repeating the derivation with this new projection operator we obtain

$`E\varphi _0(𝐫)`$ $`=`$ $`{\displaystyle 𝑑𝐫^{}\left[_{}(𝐫,𝐫^{})+\delta (𝐫,𝐫^{})V(𝐫^{})\right]\varphi _0(𝐫^{})},`$ (32)
$`E\varphi _E(𝐫)`$ $`=`$ $`{\displaystyle 𝑑𝐫^{}\left[_{}(𝐫,𝐫^{})+\delta (𝐫,𝐫^{})V_{00}(𝐫^{})\right]\varphi _E(𝐫^{})},`$ (33)
where $`_{}(𝐫,𝐫^{})`$ is the mass operator that occurs in the particle-hole Green’s function Escher and Jennings (2002). In general there is no simple relation between $`_{}(𝐫,𝐫^{})`$ and $`\overline{}(𝐫,𝐫^{})`$. The equations for $`\varphi (𝐫)`$ and $`\overline{\varphi }(𝐫)`$ are now formally the same. At this point we would expect the widths calculated with these functions to be numerically similar. The only difference in the two treatments is the form of $`V(𝐫)`$ and $`V_{00}(𝐫)`$. When the argument leading to the justification of perturbation theory is correct we expect both approaches to work equally well.

The argument for Eq. (23) can be restated using the one-body functions defined above. The resonant state, that is the solution of the exact Schrödinger-like equation $`E_0\varphi _0^{\mathrm{true}}(𝐫)=𝑑𝐫^{}(𝐫,𝐫^{})\varphi _0^{\mathrm{true}}(𝐫^{})`$, grows exponentially as we go into the classical forbidden region under the Coulomb barrier from the outside and the wave function is exponentially large in the interior. On the other hand, $`\varphi _E(𝐫)`$ coming from $`|\zeta _E`$, the solution of the many body-state defined with the Hamiltonian $`QHQ`$, is exponentially suppressed due to the orthogonality of $`|\psi _0`$ and $`|\zeta _E`$. If the suppression is large (i.e. the state under consideration is narrow) the precise form of $`|\zeta _E`$ in the interior is not important as it is essentially zero. The approximate form of $`H_{00}`$ is chosen so the approximate $`|\zeta _E`$ is also zero in the interior. As we will see later we do not need $`|\zeta _E`$ in the interior but only outside at some radius, $`r_0`$.

## VI Two expressions for the width

The decay width is given in terms of the imaginary part of the pole location in Eq. (22). We find

$`\mathrm{\Gamma }_0`$ $``$ $`2\pi \left|{\displaystyle 𝑑𝐫\varphi _0^{}(𝐫)V(𝐫)\varphi _{E_0}(𝐫)}\right|^2`$ (34)
or

$`\mathrm{\Gamma }_0`$ $``$ $`2\pi \left|{\displaystyle 𝑑𝐫\overline{\varphi }_0^{}(𝐫)\overline{V}(𝐫)\overline{\varphi }_{E_0}(𝐫)}\right|^2,`$ (35)
depending on whether we use the bar-ed amplitudes or the standard overlap functions. As discussed below, we expect the approximations leading to the two equations to be valid simultaneously, so that the resulting two expressions for the width will agree. We stress again that the reduction to an effective one-body problem is not an approximation but emerges rather naturally from the formalism. In a pure one-body problem $`\varphi _0(𝐫)`$ and $`\overline{\varphi }_0(𝐫)`$ would be normalized to unity whereas here they are normalized to the factors $`\sqrt{S_0}`$ and $`\sqrt{\overline{S}_0}`$, respectively. The normalization factors contain the many-body aspects of the problem. For the standard overlap function $`\varphi _0(𝐫)`$, this normalization is the well-known spectroscopic factor (see also Ref. Escher and Jennings (2002)). It is useful to define normalized one-body functions $`\widehat{\varphi }_0(𝐫)`$ and $`\widehat{\overline{\varphi }}_0(𝐫)`$ via $`\varphi _0(𝐫)=\sqrt{S_0}\widehat{\varphi }_0(𝐫)`$ and $`\overline{\varphi }_0(𝐫)=\sqrt{\overline{S}_0}\widehat{\overline{\varphi }}_0(𝐫)`$ (i.e. $`𝑑𝐫|\widehat{\varphi }_0(𝐫)|^2=𝑑𝐫|\widehat{\overline{\varphi }}(𝐫)_0|^2=1`$). The wave-functions $`\varphi _{E_0}(𝐫)`$ and $`\overline{\varphi }_{E_0}(𝐫)`$ describe scattering states and are normalized asymptotically at large $`r`$. If we assume the one-body Hamiltonians $``$ and $`\overline{}`$ to be local (more will be said about this later) we can take over verbatim the development given in Refs. Gurvitz and Kalbermann (1987); Gurvitz (1988); Gurvitz et al. . The only difference between the treatment presented there and the approach shown here is the presence of the factors $`S_0`$ (for the standard overlap functions) and $`\overline{S}_0`$ (for the bar-ed function). To make the dependence on these normalization factors explicit we rewrite the last equations as

$`\mathrm{\Gamma }_0`$ $``$ $`2\pi S_0\left|{\displaystyle 𝑑𝐫\widehat{\varphi }_0^{}(𝐫)V(𝐫)\varphi _{E_0}(𝐫)}\right|^2`$ (36)
and

$`\mathrm{\Gamma }_0`$ $``$ $`2\pi \overline{S}_0\left|{\displaystyle 𝑑𝐫\widehat{\overline{\varphi }}_0^{}(𝐫)\overline{V}(𝐫)\overline{\varphi }_{E_0}(𝐫)}\right|^2.`$ (37)
To further understand the results we peruse one particular line of development. As previously noted, we take $`V(𝐫)`$ to be zero inside the nucleus, $`V(𝐫)=0`$ for $`rr_0`$. For $`rr_0`$, we have $`𝑑𝐫^{}_{}(𝐫,𝐫^{})\varphi _E(𝐫^{})=E\varphi _E(𝐫)`$ since $`V_{00}(𝐫)`$ vanishes for large $`r`$. For definiteness we follow Refs. Gurvitz and Kalbermann (1987); Gurvitz (1988) and take $`r_0`$ to the maximum of the potential. For the example of fluorine discussed below, this corresponds to about 5 fm. Thus we write the integral in Eq. (34) as

$`I`$ $`=`$ $`{\displaystyle 𝑑𝐫\theta (rr_0)\varphi _0^{}(𝐫)V(𝐫)\varphi _{E_0}(𝐫)}`$ (38)
$`=`$ $`{\displaystyle 𝑑𝐫𝑑𝐫^{}\theta (rr_0)\varphi _0^{}(𝐫)\left\{_{}(𝐫,𝐫^{})+[E_0+V(𝐫)]\delta (𝐫,𝐫^{})\right\}\varphi _{E_0}(𝐫)},`$ (39)
with a similar equation holding for the bar-ed quantity. Naively one might expect this to vanish and indeed it would if the step function $`\theta `$ were not there. In the following we shall consider the case in which the non-locality in $`_{}(𝐫,𝐫^{})`$ is small enough to be properly described by a local potential and an effective mass $`m_k(𝐫)`$. In this case one has Ma and Wambach (1983)

$$_{}(𝐫,𝐫^{})=\left[\frac{\mathrm{}^2}{2m_k(𝐫)}+V_{}(𝐫)\right]\delta (𝐫,𝐫^{}).$$
(40)
In this equation $`_{}`$ is correctly evaluated at the energy $`E_0`$ of the scattered state, therefore $`m_k(𝐫)`$ accounts for the sole spatial non-locality and differs from the usual definition of effective mass $`m^{}(𝐫)`$ Negele and Yazaki (1981); Jaminon and Mahaux (1989). When Eq. (40) holds, the integral, Eq. (39), can be simplified through integration by parts and, for spherically symmetric potentials $`V(r)`$, $`V_{}(r)`$ and effective mass $`m_k(r)`$, one obtains

$`I`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m_k(r_0)}}\left[\varphi _{\mathrm{r0}}^{}(r_0)\varphi _{\mathrm{r}E_0}^{}(r_0)\varphi _{\mathrm{r}E_0}(r_0)\varphi _{\mathrm{r0}}^{}(r_0)\right]`$ (41)
$`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2\overline{m}_k(r_0)}}\left[\overline{\varphi }_{\mathrm{r0}}^{}(r_0)\overline{\varphi }_{\mathrm{r}E_0}^{}(r_0)\overline{\varphi }_{\mathrm{r}E_0}(r_0)\overline{\varphi }_{\mathrm{r0}}^{}(r_0)\right],`$ (42)
where the wave functions are assumed to factorize as $`\varphi (𝐫)=\frac{\varphi _\mathrm{r}(r)}{r}Y_{\mathrm{}}^m(\theta ,\phi )`$ and $`Y_{\mathrm{}}^m(\theta ,\phi )`$ are spherical harmonic functions of the angular variables $`\theta ,\phi `$. The prime denotes the derivation with respect to the radial coordinate $`r`$.

In deriving these last two equations we made an additional approximation, namely that the potentials are local in the vicinity of $`r_0`$ \[or at least the non-locality is restricted to an effective mass $`m_k(r)`$\]. The two equations will be simultaneously valid only if the norm operator is local, i.e. $`𝒩(𝐫,𝐫^{})=𝒩(𝐫)\delta (𝐫𝐫^{})`$, in the vicinity of $`r_0`$. In this case $`\overline{m}_k(r_0)/m_k(r_0)=𝒩(r_0)`$ and the two equations are identical. Since $`r_0`$ is in the tail of the density distribution, one expects the norm operator to be unity in its vicinity, which means that the two expressions for the width should give the same result for realistic models. This is verified numerically in Fig. 1, where the function $`\varphi _\mathrm{r}(r)`$ (solid lines) and $`\overline{\varphi }_\mathrm{r}(r)`$ (dashed lines) are displayed for the $`3/2^{}`$ and $`3/2^+`$ decaying states of <sup>17</sup>F at respective excitation energies of 4.64 and 5.00 MeV (the proton threshold is at 0.60 MeV). The functions are calculated with the self-consistent Green’s function method of Refs. Barbieri and Dickhoff (2001, 2002); Barbieri (2002). The $`3/2^+`$ state is well explained by the nuclear mean field approach and provides a typical example of a strong state, for which $`S_0\overline{S}_0`$. The $`3/2^{}`$ state, on the other hand, is a typical example of a weak state for which Ref. Barbieri (2002) gives $`S_00.04`$ and $`\overline{S}_00.14`$. It is worth to emphasize the difference between these two cases. For strong states the quenching of the spectroscopic factor is due both to short-range correlations and to the coupling to other excitations of the system Barbieri and Dickhoff (2002). Nevertheless they maintain a strong single-particle structure and the orbital occupancy is of the order of unity. Weak states, instead, have a more complicated structure and can be seen as collective excitations on their own rather than having a single-particle character. As a consequence the one-body spectroscopic factor can be one order of magnitude smaller or less. In this case the functions $`\varphi (𝐫)`$ and $`\overline{\varphi }(𝐫)`$ sample the one-body substructure in a different way, whence the difference in their normalizations with $`\overline{S}_0`$ larger than $`S_0`$ Escher and Jennings (2002). For <sup>16</sup>O + p, the nuclear interaction becomes negligible beyond 5 fm typically and Fig. 1 shows that $`\varphi (𝐫)`$ and $`\overline{\varphi }(𝐫)`$ are equal beyond this radius, which means that the norm operator is unity. In principle, there is no problem calculating the width from a microscopic model: one may define two different functions $`\varphi (𝐫)`$ and $`\overline{\varphi }(𝐫)`$, which have two different normalization factors $`\sqrt{S}`$ and $`\sqrt{\overline{S}}`$, respectively, but the estimate for the width will be the same with both of them. In fact, a similar formula would be valid for any amplitude of the form $`\varphi ^\nu (𝐫)=𝑑r^{}𝒩(𝐫,𝐫^{})^\nu \varphi (𝐫^{})`$ for arbitrary $`\nu `$, since $`𝒩(𝐫,𝐫^{})`$ is unity around $`r_0`$.

Let us emphasize that for the actual calculation of the width by the above formulas, a microscopic model based on the harmonic oscillator basis, as in Refs. Barbieri and Dickhoff (2002); Barbieri (2002) may not be the best choice since the wave function is not expected to be very precise in the vicinity of $`r_0`$. This is demonstrated in Fig. 1, where the $`3/2^+`$ resonant wave function of the phenomenological mean field Woods-Saxon potential of Ref. Sparenberg et al. (2000) is shown for comparison. The agreement with the microscopic wave function is good in the interior but deteriorates above 4 fm (where the approximation of the potential by a harmonic oscillator breaks down, see Fig. 2). The difference between both models is particularly large in this case since the $`3/2^+`$ state is wide ($`\mathrm{\Gamma }_0=1.5`$ MeV); for narrow states, the harmonic oscillator approximation could be sufficient, a possibility which will be explored elsewhere. Let us finally stress that the phenomenological mean field potential of Ref. Sparenberg et al. (2000) does not reproduce the $`3/2`$ state because its structure is not well approximated by an <sup>16</sup>O core plus a proton.

Equations (36) and (37) indicate that there may be a serious problem in extracting spectroscopic factors from measured decay widths. The standard method for determining a spectroscopic factor involves dividing an experimental width by a single-particle width calculated with a phenomenological model. However, it is not clear a priori whether the phenomenological wave functions are good approximations to $`\widehat{\varphi }_0(𝐫)`$ and $`\varphi _E(𝐫)`$, or to $`\widehat{\overline{\varphi }}_0(𝐫)`$ and $`\overline{\varphi }_E(𝐫)`$; hence, it is not obvious whether dividing the experimental width by the result of a single-particle calculation provides the spectroscopic factor $`S_0`$ or the normalization $`\overline{S}_0`$ of the auxiliary function $`\overline{\varphi }(𝐫)`$ (or the norm of yet another one-body function). Normally, one assumes in proton emission studies that $`\widehat{\varphi }_0(𝐫)`$ and $`\varphi _E(𝐫)`$ can be equated with the wave functions obtained from phenomenological potentials (see for example Ref. Åberg et al. (1997)). On the other hand, in the context of some alpha emission studies it has been argued very strongly that $`\widehat{\overline{\varphi }}_0(𝐫)`$ and $`\overline{\varphi }_E(𝐫)`$ correspond to phenomenological wave functions Varga and Lovas (1991); Lovas et al. (1998) . If the latter is true then the experiments would be sensitive to $`\overline{S}_0`$ rather than to $`S_0`$. For strong states, with a clear core-plus-particle structure, this is mainly a philosophical issue, since $`S_0\overline{S}_0`$ holds. For weak states, however, $`S_0`$ and $`\overline{S}_0`$ can be significantly different from each other.

We have attempted to resolve the ambiguity outlined above by calculating the effective local potentials corresponding to $`\varphi (𝐫)`$ and $`\overline{\varphi }(𝐫)`$ and comparing them with typical phenomenological potentials. This is done by inversion of the local (radial) Schrödinger equation:

$`V_{\mathrm{eff}}(r)=E+{\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{\varphi _\mathrm{r}^{\prime \prime }(r)}{\varphi _\mathrm{r}(r)}},`$ (43)
where the effective potential $`V_{\mathrm{eff}}(r)`$ is the sum of the interaction potential (nuclear + Coulomb) and the centrifugal term. The effective potentials corresponding to the radial wave functions of Fig. 1 are shown in Fig. 2. For the strong $`3/2^+`$ state ($`\mathrm{}=2`$ centrifugal term) the three potentials are in reasonable agreement, except above 4 fm where the potential extracted from the microscopic functions asymptotically approaches the harmonic oscillator potential that generated them. For the weak $`3/2^{}`$ state ($`\mathrm{}=1`$ centrifugal term) the potentials deduced from $`\varphi (𝐫)`$ and $`\overline{\varphi }(𝐫)`$ display a singularity and are very different from traditional phenomenological potentials. This singularity occurs because the zeroes of $`\varphi _\mathrm{r}(r)`$ and $`\varphi _\mathrm{r}^{\prime \prime }(r)`$ occur at different radii. This is probably not an artifact of the model but a real effect and arises from the relative sign of the 0 and $`2\mathrm{}\omega `$ contributions to the spectroscopic amplitude. Let us remark that Eq. (43) assumes a constant effective mass $`m_k(r)=m`$. We have checked that introducing a realistic effective mass does not lead to a simpler potential for the weak state; more will be said about this elsewhere. Let us finally emphasize that, though their energies are close to one another, the strong and weak states correspond to very different potentials, which suggests that a strong energy dependence is a necessary feature for such potentials (see the discussion in the next section). In conclusion, this example shows that constructing reliable phenomenological local potentials for extracting spectroscopic factors from experimental cross sections is non-trivial. Moreover, since the characteristics of the potential depend strongly on the state it is difficult to determine which normalization ($`S_0`$ or $`\overline{S}_0`$) would be extracted from a comparison with the experimental data.

## VII Alternative expressions for the width

Let us now return to the expression for the decay width and establish a link with known results. In Refs. Gurvitz and Kalbermann (1987); Gurvitz (1988); Gurvitz et al.  the spatial derivatives were evaluated using a special form for $`V(r)`$ and $`V_{00}(r)`$. In those references these auxiliary potentials were chosen such that the total potential for $`\varphi _{\mathrm{r0}}(r)`$ was constant outside a given radius and the one for $`\varphi _{\mathrm{r}E_0}(r)`$ was constant inside. In Refs. Gurvitz and Kalbermann (1987); Gurvitz (1988) the separation radius is at the maximum of the potential barrier, whereas in Ref. Gurvitz et al.  the advantage of using a radius of the order of the nuclear range (our $`r_0`$) is pointed out. For this last case, the width is then given by

$`\mathrm{\Gamma }_02\pi \left[{\displaystyle \frac{\mathrm{}^2\alpha }{2m_k(r_0)}}\right]^2|\varphi _{\mathrm{r0}}(r_0)\varphi _{\mathrm{r}E_0}(r_0)|^2,`$ (44)
where $`\alpha =\sqrt{2m_k(r_0)[V(r_0)E_0]/\mathrm{}^2}`$.

An alternative approach is to exploit the Wronskian form of Eqs. (41) or (42). If $`V(r)`$ is zero for radii less than the outer turning point and $`V_{00}(r)`$ zero for radii greater than the inner turning point then there is a region where $`\varphi _{\mathrm{r0}}(r)`$ and $`\varphi _{\mathrm{r}E_0}(r)`$ satisfy the same differential equation. If the potential is local the Wronskian is a constant. Assuming a thick barrier, there will be a region where $`\varphi _{\mathrm{r0}}(r)`$ is an irregular Coulomb function, the regular Coulomb function having decayed away, and where $`\varphi _{\mathrm{r}E_0}(r)`$ is a regular Coulomb function, the irregular Coulomb function having decayed away. All that is required is to determine the proportionality constants. These have a simple expression for a constant effective mass $`m_k(r)=m`$: for $`\varphi _{\mathrm{r}E_0}(r)`$ we have $`\varphi _{\mathrm{r}E_0}(r)=\sqrt{2m/\mathrm{}^2k_0\pi }F(k_0r)`$, with $`k_0=\sqrt{2mE_0/\mathrm{}^2}`$, while for $`\varphi _{\mathrm{r0}}(r)`$ the proportionality constant can be written as $`\varphi _{\mathrm{r0}}(r_0)/G(k_0r_0)`$. Here $`F(k_0r)`$ and $`G(k_0r)`$ are the regular and irregular Coulomb functions respectively, the Wronskian of which equals $`k_0`$. The resulting width is then, as in Ref. Gurvitz et al. ,

$`\mathrm{\Gamma }_0{\displaystyle \frac{\mathrm{}^2k_0}{m}}|\varphi _{\mathrm{r0}}(r_0)/G(k_0r_0)|^2.`$ (45)
This can be simplified further. The wave function $`\varphi _{\mathrm{r0}}(r)`$ is normalized to the appropriate spectroscopic factor, $`S_0`$. We can also consider the true scattering wave function at resonance, $`\varphi _{\mathrm{r0}}^{\mathrm{true}}(r)`$. In the interior it will behave like $`\varphi _{\mathrm{r0}}(r)`$ while in the exterior region it will behave like $`G(k_0r)`$. Normalizing it to $`G(k_0r)`$ in the exterior region we obtain for the width:

$`\mathrm{\Gamma }_0{\displaystyle \frac{S_0\mathrm{}v_0}{_0^{r_\mathrm{t}}𝑑r|\varphi _{\mathrm{r0}}^{\mathrm{true}}(r)|^2}},`$ (46)
where $`v_0=\mathrm{}k_0/m`$ is the asymptotic velocity. The exact value of the upper limit on the integral is not crucial and we take it to be the outer turning point (see numerical justification below).

Equation (46) can also be derived in a more transparent way (see also Refs. Breit (1959); Iliadis (1997)). Consider the Schrödinger-like equation $`[H(E)E]\varphi (𝐫,E)=0`$ for the overlap function. Since the effective one-body equation we have been considering can depend on the energy we keep an explicit energy dependence in the Hamiltonian. Differentiating this equation with respect to $`E`$ we find

$`[H(E)E]{\displaystyle \frac{\varphi (𝐫,E)}{E}}=\left[1{\displaystyle \frac{H(E)}{E}}\right]\varphi (𝐫,E).`$ (47)
Next we multiply by $`\varphi ^{}(𝐫,E)`$ and integrate up to some radius $`r_l`$. If $`H(E)`$ is not hermitian then $`\varphi ^{}(𝐫,E)`$ should be replaced by the time reversed state. If the potential is local in the vicinity of $`r_l`$ we can integrate by parts on the left hand side. Assuming spherical symmetry, this gives us

$`{\displaystyle \frac{\mathrm{}^2}{2m_k(r_l)}}\left[\varphi _\mathrm{r}^{}(r_l,E){\displaystyle \frac{\varphi _\mathrm{r}^{}(r_l,E)}{E}}\varphi _\mathrm{r}^{}(r_l,E){\displaystyle \frac{\varphi _\mathrm{r}(r_l,E)}{E}}\right]`$
$`={\displaystyle _0^{r_l}}𝑑r\varphi _\mathrm{r}^{}(r,E)\left[1{\displaystyle \frac{H_\mathrm{r}(E)}{E}}\right]\varphi _\mathrm{r}(r,E),`$ (48)
where the prime denotes a partial derivative with respect to $`r`$ and $`H_\mathrm{r}(E)`$ is the radial Hamiltonian. If we take $`r_l`$ to be outside the range of nuclear force, $`\varphi _\mathrm{r}(r,E)`$ can be written as

$`\varphi _\mathrm{r}(r,E)\underset{r>r_0}{}\mathrm{cos}\delta (E)F(kr)+\mathrm{sin}\delta (E)G(kr).`$ (49)
At a narrow resonance the phase shift $`\delta (E)`$ will be rapidly varying so we expect that the largest part of the energy dependence will come from the phase shift and not from the Coulomb functions. The energy dependence of the Coulomb functions will be minimized if the radius is chosen to be near the outside turning point. For example, for large $`r`$ the regular Coulomb function will have a $`\mathrm{sin}(kr)`$ dependence. Differentiating with respect to $`E`$ will give $`r(dk/dE)\mathrm{cos}(kr)`$, which diverges for large $`r`$. As $`r`$ decreases to the turning point this asymptotic form for the wave function breaks down. However a similar argument using exponentials holds inside the turning points. Thus near the outside turning point we have

$`{\displaystyle \frac{\varphi _\mathrm{r}(r,E)}{E}}{\displaystyle \frac{d\delta (E)}{dE}}[\mathrm{sin}\delta (E)F(kr)+\mathrm{cos}\delta (E)G(kr)].`$ (50)
We have checked this relation numerically for Woods-Saxon plus Coulomb potentials and verified that for widths less then 15 KeV the error does not exceed 3%. That the radius should be chosen near the outer turning point was also confirmed numerically. The Wronskian relation can now be written as

$`{\displaystyle \frac{\mathrm{}^2k}{2m_k(r_\mathrm{t})}}{\displaystyle \frac{d\delta (E)}{dE}}{\displaystyle _0^{r_\mathrm{t}}}𝑑r\varphi _\mathrm{r}^{}(r,E)\left[1{\displaystyle \frac{H_\mathrm{r}(E)}{E}}\right]\varphi _\mathrm{r}(r,E),`$ (51)
where we have taken $`r_l=r_\mathrm{t}`$ to be the outside turning point radius.

At the resonant energy we expect the energy variation of the phase shift to be a maximum so the resonance energy occurs when $`_0^{r_\mathrm{t}}𝑑r\varphi _\mathrm{r}^{}(r,E)\left[1\frac{H_\mathrm{r}(E)}{E}\right]\varphi _\mathrm{r}(r,E)`$ is a maximum. At the resonance energy $`d\delta (E_0)/dE=2/\mathrm{\Gamma }_0`$, and we have

$`\mathrm{\Gamma }_0{\displaystyle \frac{\mathrm{}v_0}{_0^{r_\mathrm{t}}𝑑r\varphi _\mathrm{r}^{}(r,E_0)\left[1\frac{H_\mathrm{r}(E_0)}{E}\right]\varphi _\mathrm{r}(r,E_0)}},`$ (52)
where now $`v_0=\mathrm{}k_0/m`$ is the asymptotic value of the velocity. For bound states the spectroscopic factor can be written as $`S_0=\left[_0^{\mathrm{}}𝑑r\varphi _\mathrm{r}^{}(r)\varphi _\mathrm{r}(r)\right]/\left\{_0^{\mathrm{}}𝑑r\varphi _\mathrm{r}^{}(r)\left[1\frac{H_\mathrm{r}(E)}{E}\right]\varphi _\mathrm{r}(r)\right\}`$ (see Ref. Wegmann (1969)), which does not depend on the specific normalization of the overlap function $`\varphi (𝐫)`$. By extending this relation to define the spectroscopic factor for resonant states, we recover Eq. (46).

We finally note that the expression (52) of the width is independent of the choice of $`\varphi _\mathrm{r}(r,E)`$ or $`\overline{\varphi }_\mathrm{r}(r,E)`$. Its denominator can be rewritten as

$`D={\displaystyle _0^{r_\mathrm{t}}}𝑑r\varphi _\mathrm{r}^{}(r,E)\left\{{\displaystyle \frac{}{E}}\left[EH_\mathrm{r}(E)\right]\right\}\varphi _\mathrm{r}(r,E).`$ (53)
The bar-ed amplitudes are defined as $`\overline{\varphi }(𝐫,E)=𝑑𝐫𝒩(𝐫,𝐫^{})^{1/2}\varphi (𝐫,E)`$ and a corresponding expression for the Hamiltonian is $`\overline{H}(E)=𝒩^{1/2}\left[H(E)E\right]𝒩^{1/2}+E`$. Thus the denominator is invariant under this transformation, as well as under a general transformation with an arbitrary power of $`𝒩(𝐫,𝐫^{})`$. All that matters is that $`H_\mathrm{r}(E)`$ and $`\varphi _\mathrm{r}(r,E)`$ are consistent with one another. These wave functions are phase equivalent and any of them can be used.

## VIII Discussion

We have embedded the elegant Gurvitz-Kalbermann approach of proton emission Gurvitz and Kalbermann (1987); Gurvitz (1988); Gurvitz et al.  into a full many-body picture. We have reduced the formalism to an effective one-body problem and demonstrated that the decay width can be expressed in terms of a one-body matrix element multiplied by a normalization factor. At first sight, this result agrees with the standard procedure for extracting spectroscopic factors from measurements via dividing an experimental width by a calculated single-particle width (see, for example, Ref. Åberg et al. (1997)). The present work, however, clearly demonstrates that this procedure for determining spectroscopic factors is only valid if the phenomenological potential used to generate the single-particle width corresponds to the potential in $`_{}`$ \[see Eqs. (32) and (33)\]. It is not a priori clear that this is actually the case. In fact, the authors of Refs. Varga and Lovas (1991); Lovas et al. (1998) (and prior to that the authors of Ref. Fliessbach and Mang (1976)) have argued strongly that the phenomenological potential approximates the potential in $`\overline{}`$ \[see Eq. (29)\]. While the studies of Refs. Varga and Lovas (1991); Lovas et al. (1998) were carried out for alpha decay, the arguments given there can be carried over to a description of the proton emission process. Furthermore, Eq. (52) suggests that $`_0^{r_t}𝑑r\varphi _\mathrm{r}^{}(r)\left[1\frac{H_\mathrm{r}(E)}{E}\right]\varphi _\mathrm{r}(r)`$ is the appropriate observable that can be extracted from proton emission experiments.

Besides the ambiguity regarding whether the standard spectroscopic factor or an auxiliary normalization is extracted from the experimental procedure, it has been demonstrated that constructing a reliable phenomenological potential is non-trivial. The situation is quite complicated since the interaction with the nuclear medium strongly depends on the initial state of the ejected proton. For states with a clear core-plus-particle structure (i.e., with a large spectroscopic factor), traditional phenomenological potentials seem to provide good approximations to the nuclear mean field and spectroscopic factors can be determined from proton emission studies. In this case, the spectroscopic factor extracted from the experiment can be safely compared to the results of nuclear many-body calculations (note also that for large $`S_0`$ values, $`S_0`$ and $`\overline{S}_0`$ are approximately equal Escher and Jennings (2002) and the distinction between the two approaches discussed here becomes irrelevant).

For weak states, which have a more complicated many-body structure, standard phenomenological potentials do not give a proper approximation to the nuclear medium, as shown by the radial shape of the 3/2<sup>-</sup> states in Fig. 2. Also, as discussed in Ref. Escher and Jennings (2002), the dependence of the spectroscopic factor on the energy derivative of the effective one-body Hamiltonian implies that the nuclear medium must be strongly energy dependent. This feature, which is missing in most phenomenological optical potentials, is also confirmed by the numerical results displayed in Fig. 2. Thus for weak states simple potential models are probably not valid for either $`\varphi (𝐫)`$ or $`\overline{\varphi }(𝐫)`$. States that are neither weak nor strong will also have to be dealt with on a case by case basis.

Since $`\varphi (𝐫)`$ and $`\overline{\varphi }(𝐫)`$ are identical for large radii they are phase equivalent and elastic scattering experiments cannot distinguish between them. We conclude that additional experimental input, together with an accurate derivation of the optical potential based on first principles, is required in order to resolve the question regarding which one-body Hamiltonian is most appropriately approximated by a phenomenological model.

###### Acknowledgements.

The authors thank C. Chandler for valuable discussions and comments on a preliminary draft of the paper. Financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) is appreciated. This work was performed in part under the auspices of the U. S. Department of Energy by the University of California, Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48.


format: text

text:
# Electronic properties of NiCl2 tubular nanostructures

## Abstract

Atomic models of zigzag (n, 0)- and armchair (n,n)-like NiCl2 nanotubes (n = 4 - 29) formed by rolling (100) single layers of the bulk NiCl2 which crystallizes in the CdCl2-type structure, are constructed and their electronic properties and bond indices are investigated using the tight-binding band theory. The calculations performed show that all the considered nanotubes in non-magnetic state are uniformly metallic-like. The density of states at the Fermi level contains a considerable contribution from Ni3d states, its value depends on the atomic configuration and diameter of the tubes. The Ni-Cl covalent bonds were found in NiCl2 tubes, whereas Ni-Ni and Cl-Cl covalent interactions are almost absent. According our estimations, the zigzag (n, 0)-like NiCl2 nanotubes are more energetically favorable.

Quasi-one-dimensional (1D) inorganic nanotubes (NTs) attract much interest due to wide prospects of technological application of materials produced on their basis dresselhaus ; saito ; moriarty ; knupfer ; ajayan ; tenne ; ivanovskii ; pokropivnyi . Besides carbon (see dresselhaus ; saito ; moriarty ), the first inorganic compounds, which were the basis for nanotube production, were layered hexagonal boron nitride copra ; louseau ; rubio  and d-metal dichalcogenides MX2 (M = Mo, W, Ta; X = S, Se) tenne ; nath ; rosentsveig ; tenne2 ; tenne3 , having structures analogous to that of graphite. Recently, NTs of other sulfides were successfully synthesized: InS hollingsworth , ZnS wang , Bi2S3 ye , ReS2 bronson , WxMoyCzS2 hsu , ZrS2, HfS2 nath2 ; nath3 .

Along with the above-mentioned NTs, a large number of other inorganic nanotubes have been produced or predicted. Among them there are new 1D tubular structures of metals (Na, Cu opitz , Ni bao , Ag, Au lay ; ahn , Bi li ; su ), boron boustani , and silicon zhang . Some efforts are reported in synthesis and modeling of nanotubes based on inorganic compounds with planar or pseudoplanar structures: BC3, BC2N, CNx (see reviews tenne ; ivanovskii ; pokropivnyi ), metal diborides MB2 (M = Mg, Al, Sc, Ti, Zr) chernozatonskii ; quandt ; ivanovskaya , AlB2-like ternary silicides Sr(GaxSi1-x)2, Ca(AlxSi1-x)2 shein  etc. A large number of 1D nanomaterials (nanowires, nanorods, nanotubes) based on metal oxides such as TiO2 imai ; shimizu ; kasuga ; seo ; miao ; yao , V2O5 satishkumar ; krumeich ; muhr ; niederberger ; dobley ; pan ; chandrappa , Co3O4 shi ; shi2 , MnO2, WO3 lakshmi , MoO3, ReO2 satishkumar2  and also on ZnO hu , SiOx zhang2 ; wang2  etc. have been obtained. Important progress was achieved in producing nanotubes of semiconducting materials: SiGe, InAs/GaAs, InGaAs/GaAs, SiGe/Si InGeAs/GaAs schmidt ; prinz ; deneke .

Recently, the synthesis of closed-cage nanoclusters and nanotubes NiCl2 which is the first reprsentative of layered metal-halogen compounds, was reported hacohen ; hacohen2 . High-resolution transmission electron microscopy analysis (HRTEM) and electron diffraction measurements showed that NiCl2 layers forming the walls of those multi-walled NTs were completely crystalline. Taking into account that bulk NiCl2 is of ferro- and weak antiferromagnetic Ni-Ni coupling type in and across the layers, respectively, it is suggested hacohen  that unique magnetic properties may be expected for NiCl2 NTs.

In this communication, the atomic simulation of the nanotubes constructed from NiCl2 single layers is performed and their electronic band structure and bond indices are calculated and analyzed as a function of the tubes diameters (D) in armchair\- and zigzag-like forms.

NiCl2 crystallizes in the CdCl2 structural type (space group P3m) with the lattice parameters a = 0.3483, c =1.7400 nm. The bulk NiCl2 consists of layers of a Ni sheet sandwiched between two chlorine sheets (Cl-Ni-Cl). Within the layers there exist ionic-covalent Ni-Cl bonds; the neighboring (Cl-Ni-Cl) layers interact with each other due to weak Van der Waals forces. The atomic models of tubular NiCl2 structures were constructed analogously to the nanotubes of layered Mo, W, Nb dichalcogenides seifert ; seifert2 ; seifert3 ; seifert4 ; seifert5  by mapping a ”triple” (Cl-Ni-Cl) layer onto the tube surface. Such ”single-walled” tubes consist of three coaxial Cl-Ni-Cl cylinders built up of wrapped hexagonal Ni, Cl monoatomic sheets, Fig. 1. As for single-walled carbon NTs, see dresselhaus ; saito ; moriarty ; knupfer ; ajayan ; tenne ; ivanovskii ; pokropivnyi , three groups of NiCl2 NTs can be classified as non-chiral (armchair (n,n)-, zigzag (n,0)-like) and chiral (n,m)-like nanotubes. The majority of the synthesized hacohen2  NiCl2 NTs were non-chiral.

We calculated the electronic structure of non-chiral (n,0) and (n,n) NiCl2 NTs as a function of n in ranges from (8,0) to (29,0) and from (4,4) to (29,29), which correspond to the intervals of the diameters of ”outer” chlorine cylinders (Dout) 2.272 - 7.096 and 1.984 - 11.902 nm, respectively, Table 1.

The tight binding band structure method within the extended Huckel theory (EHT) approximation hoffmann  was employed. Besides the electronic spectra, the densities of states (DOS), crystal orbital overlap populations (COOPs), and total band energies of the nanotubes (Etot) were obtained.

Figure 2 shows the calculated values of Etot (per NiCl2 unit) versus NT diameters. The Etot dependence follows a 1/D2 behavior indicative of a decrease in NTs stability with decreasing D. Analogous dependences of strain energy (the difference between the energies of the plane atomic layer and the corresponding NT characterizes the chemical stability of tubular structures) are known for carbon and the majority of non-carbon NTs dresselhaus ; saito ; moriarty ; knupfer ; ajayan ; tenne ; ivanovskii ; pokropivnyi . It is worth noting that according to our results for all the considered tubes, zigzag-like configurations of the NiCl2 nanotubes are more stable.

Within non-magnetic approach, all NTs are metallic-like. The electronic DOSs for the (15, 15) and (15, 0) NiCl2 nanotubes are showed in Fig. 3. The pronounced very sharp DOS peak near the Fermi level consists mainly of Ni3d states. Hence the magnetic instability can occur and it may be assumed that the near-Fermi Ni3d states will play a large role in the magnetic properties of NiCl2 NTs. However this aspect remains to be considered.

It is significant that the shape of the electronic spectrum is essentially modified by the tube geometry, see Fig. 3. For example, for armchair (15, 15) NTs, the near-Fermi bands are separated from the filled valence bands by a forbidden gap, which is absent for zigzag (15, 0) NTs.

To analyze the covalence of intra-atomic interactions, the COOPs of paired bonds were examined. It was established that Ni-Cl covalent ineractions are the strongest bonds (Table 2), whereas Ni-Ni and Cl-Cl covalent-type interactions are absent (COOPs less 0). There is also sharp anisotropy of separate Ni-Cl bonds depending both on their orientation relative to the NT axis and the coordination of chlorine atoms belonging to ”inner” or ”outer” cylinders, as well as on the tube geometry. It is seen that the ratio of COOPs values for Ni-Clin and Ni-Clout bonds differs drastically for armchair (n,n) - and zigzag (n,0)-like tubes. In particular, for zigzag-like tubes the COOPs values of Ni-Clin are much greater than those for Ni-Clout. This suggests higher reactivity of the ”outer” chlorine atoms and a possibility of their rearrangement. The above features of interatomic bonds depending on the tube geometry may become an important factor in developing microscopic models of NiCl2 tubes growth both in zigzag\- and armchair-like (and also chiral hacohen2 ) forms.

In conclusion we have constructed atomic models of single-walled NiCl2 nanotubes in zigzag\- and armchair-like forms and investigated their electronic properties and bond indices by the tight-binding band method. We show that both zigzag\- and armchair-like nanotubes in non-magnetic state are metallic-like, with a very sharp DOS peak near the Fermi level predominantly of Ni3d character. This points to magnetic instability and a possibility of magnetic phenomena in NiCl2 NTs. Zigzag-like NiCl2 nanotubes were found to be more stable. It was established that Ni-Cl covalent bonds are the strongest interactions in nickel dichloride NTs.


text:
# Does the Sun shrink with increasing magnetic activity?

## 1 INTRODUCTION

Temporal variations in the solar radius has been a controversial topic, as direct measurements of the solar radius have given conflicting results (Delache, Laclare & Sadsaoud 1985; Wittmann, Alge & Bianda 1993; Fiala, Dunham & Sofia 1994; Laclare et al. 1996; Noël 1997; Emilio et al. 2000). The reported temporal variations in the solar radius ranges from 0 to 700 km. It is important to estimate the radius variations with solar cycle as these can provide useful constraint for models to explain the luminosity variation with solar cycle (Gough 2001). In particular, the ratio of the radius variation to the luminosity variation, $`W=(\mathrm{\Delta }R/R)/(\mathrm{\Delta }L/L)`$ depends on the theoretical model of luminosity variations. The luminosity variation is known to be about 0.001 (e.g., Mecherikunnel 1994) between the maximum and minimum of solar activity. Thus, it is important to obtain a reliable estimate of radius variation over the solar cycle so that we can distinguish between these models.

Recently, Schou et al. (1997) and Antia (1998) have demonstrated that the frequencies of f-modes can be used to estimate the solar radius. Since these frequencies have been measured with a relative accuracy of $`10^5`$ we may expect to determine the solar radius to similar accuracy. However, there are systematic errors of order of 100 km in calibration of photospheric radius from the measured frequencies (Tripathy & Antia 1999). If these systematic errors are independent of time, then it would be possible to determine the temporal variation in the solar radius using f-mode frequencies. The attempts so far (Dziembowski et al. 1998, 2000, 2001; Antia et al. 2000, 2001) give conflicting results. Using the first few data sets from the Michelson Doppler Imager (MDI), Dziembowski et al. (1998) found that the solar radius is increasing with solar activity. They found an increase by about 4 km in 6 months just after the solar minimum. If this variation was indeed correlated to solar activity we would expect a much larger variation in radius during the solar cycle. Subsequently, using more data Dziembowski et al. (2000) found no systematic variation in the solar radius. This work used all data sets from MDI that were obtained before the contact with SOHO satellite was lost. Using a few data sets from the Global Oscillation network Group (GONG), Antia et al. (2000) found the solar radius to be decreasing with activity, but subsequently using more extensive data sets from GONG and MDI Antia et al. (2001) found no evidence for any variation in the solar radius. However, using essentially same data sets from MDI, Dziembowski et al. (2001) (hereinafter DGS) have found a decrease in the solar radius. Unfortunately, the claimed variation, if any, in all these works is of the order of a few km and even a small change in systematic errors can give rise to spurious variations of this order. Clearly, a more careful analysis of f-mode frequencies is required before drawing any conclusions about variation of the solar radius.

Antia et al. (2001) have shown that the variation in f-mode frequencies is more complex than what is assumed in other studies. These variations can be decomposed into at least two components. One of these components is oscillatory with a period of 1 yr, while the second component is correlated with solar activity. The amplitudes of both these components increase with frequency and hence are not likely to arise from radius variations. Variation in the solar radius will cause frequency shifts that are proportional to frequency, but the observed variations have much steeper dependence on frequency. The oscillatory component is most likely to be an artifact introduced by orbital period of the Earth. Antia et al. (2001) have also shown that, most of the discrepancy between different results about radius variation using f-mode frequencies can be explained if these two components are invoked in the temporal variations. In particular, Antia et al. (2000) failed to detect the oscillatory component as they used only 5 data sets covering a period of 3 years. Further, after accounting for these two components in temporal variations there is no evidence for any variation in the solar radius. DGS have claimed that the solar radius decreases at a rate of 1.5 km yr-1 during 1996–2000. However, they have not removed the oscillatory component in f-mode frequency variation and hence their claim needs to be examined carefully.

## 2 RADIUS VARIATION FROM F-MODE FREQUENCIES

It can be easily shown that if the solar radius varies by even 1 km over the solar cycle, the rate of resulting release or absorption of gravitational energy would be larger than the solar luminosity. Hence, we can rule out such radius variations. Thus, any possible variation in the solar radius must be confined to the outermost layers of the Sun. DGS have argued that since observed f-modes are trapped in a layer beneath the visible surface, they would measure the radius variation at this depth. In particular, the fractional variation in radius, could be a function of radial distance. Such a variation is, of course, realistic, but the problem is, it may not be possible to analyze such variations easily. For example, DGS have split the f-mode frequency variations into two parts, one arising from radius variation and the other from some variations in the outermost layers, which scales inversely as the mode inertia. Thus, they express the change in f-mode frequencies as

$$\mathrm{\Delta }\nu _{\mathrm{}}=\frac{3}{2}\frac{\mathrm{\Delta }R}{R}\nu _{\mathrm{}}+\frac{\mathrm{\Delta }\gamma }{I_{\mathrm{}}},$$
(1)
where $`\nu _{\mathrm{}}`$ is the frequency of the f-mode of degree $`\mathrm{}`$, $`\mathrm{\Delta }R`$ is the change in radius while $`\mathrm{\Delta }\gamma `$ measures the contribution from surface term and $`I_{\mathrm{}}`$ is the mode inertia. While fitting the expression to the observed data DGS assume $`\mathrm{\Delta }R/R`$ to be constant, which implies that the radius variations are homologous, at least, in the region where the observed modes are trapped. Thus as far as f-modes are concerned they have assumed that $`\mathrm{\Delta }R/R`$ is constant, presumably because otherwise it is difficult to proceed with the analysis. Subsequently, they claim that these radius variations arise from magnetic field variation in a layer below the outermost surface layers. This is certainly conceivable, but if that is the case then there should be an additional term in Eq. (1) that arises from the direct effects of the magnetic field. The effect of the magnetic field on f-mode frequencies cannot be assumed to be solely due to those arising from radius variation. Frequency shifts due to magnetic fields (e.g., Campbell & Roberts 1989) are not in general, proportional to frequency as is implied by Eq. (1) (the surface term cannot arise from such fields in the interior). Thus, a more complex model will be required to fit the frequency differences arising from magnetic field. The same applies to frequency variations due to density perturbations (Chitre et al. 1998). Basically, if f-mode frequency variations are due to magnetic field or density perturbation, then we need to calculate these shifts explicitly, rather than modeling them via radius variation, which cannot account for the entire effect.

To estimate the depth at which f-modes are trapped, we can consider the kinetic energy density from the eigenfunctions of f-modes in the relevant range of $`\mathrm{}=140`$–300 which are mainly used in this study. If we assume that the trapping region of each mode covers the layers where the kinetic energy density is greater than $`1/e`$ of its peak value, then the upper limit of $`\mathrm{}=300`$ f-mode is at a depth of about 1 Mm, while the lower limit of $`\mathrm{}=140`$ f-mode extends to a depth of about 12 Mm. Hence, a depth range of 1–12 Mm is expected to be covered by this study. The f-modes are not trapped between a pair of rigid boundaries and the extent of region covered by them will depend on the definition of boundary as well as on the precise mechanism responsible for frequency variations.

## 3 RESULTS USING MDI DATA

The inconsistencies pointed out in the previous section arise in attempting to find a physical model which gives rise to the radius variation inferred from f-mode frequencies. This model is relevant only if the data actually show any evidence for change in radius from f-mode frequencies. Thus, in this section, I ignore the cause of radius variation and just address the question whether the changes in f-mode frequencies imply any change in the solar radius (as defined by DGS). The f-mode frequency variation is expressed using two terms, one arising from variation in the solar radius and another from unspecified variations in the outermost surface layers. There are good reasons to look for such a term since it is known that p-mode frequency variations largely arise from variations in outer layers (e.g., Basu & Antia 2000). Thus, following DGS, I assume the frequency variation to be given by Eq. (1), where $`\mathrm{\Delta }R/R`$ and $`\mathrm{\Delta }\gamma `$ are constants. In order to test whether the observed data fit this form, I have used the same f-mode frequency data (Schou 1999) that DGS have used except for some additional data sets that are now available. These data consist of 30 sets each covering a period of 72 days starting from May 1, 1996 and ending on August 21, 2002. Note that there is a gap in data sets between 1998.5 and 1999.2, when the contact with SOHO satellite was lost. The only data set during this period has significantly worse fit and may be ignored. For each data set, I take the difference in frequency with respect to a standard solar model (in the sense observed $``$ model) with radius $`R_{}=695.78`$ Mm. The model radius is chosen to ensure that the frequency differences are small. These frequency differences are then fitted to Eq. (1). All modes with $`\mathrm{}>140`$ are used in these fits. Even if I ignore the data set taken immediately after recovery of SOHO satellite ($`\chi ^2=3.54`$), the $`\chi ^2`$ per degree of freedom in these fits varies from 1.2 to 2.5. One such fit for data obtained around 1997.0 is shown in Fig. 1. It is clear that the fit is not good and that the variation in frequency differences is more complicated than what is modeled by Eq. (1). This should be expected from the results of Antia et al. (2001), since depending on the phase of the oscillatory component in the frequency variations, it can have a different sign as compared to the other component of the variation. Furthermore, the oscillatory component has roughly the same frequency dependence as the second term in Eq. (1), while the non-oscillatory component has a significantly less steep frequency dependence and hence, may not be adequately represented by Eq. (1). Thus, at times when both components are in phase the net frequency variation may be approximated by Eq. (1), but six months later when the oscillatory component changes sign the two will be opposite, as is the case around 1997.0, data for which are shown in Fig. 1. Since the oscillatory component has a steeper frequency dependence as compared to the non-oscillatory component, at high frequencies the trend appears to reverse. Such a behavior cannot be modeled by Eq. (1).

Fig. 2 shows the results from fits to all data sets from MDI. The upper panel shows the inferred radius variation, which is similar to Fig. 2 of DGS. The middle panel shows the fitted variation in surface term $`\mathrm{\Delta }\gamma `$. This figure also looks similar to Fig. 2 of DGS, though the y-axis is different. The cause of this difference is not clear as the surface term cannot be positive as shown in Fig. 2 of DGS. The lowest panel shows the $`\chi ^2`$ per degree of freedom for each of the fits. The oscillatory trend is quite clear in all these panels. Further, comparing different panels it is clear that the best fits are obtained when the oscillatory and non-oscillatory components are of the same sign and the magnitude of the oscillatory component is close to maximum, which is the case when $`\mathrm{\Delta }\gamma `$ is the lowest. The non-oscillatory component appears to be reducing with increasing solar activity and as a result, fits improve during high activity period. For example, the best fits during pre-recovery period have $`\chi ^21.5`$ per degree of freedom, while during 2000, 2001 this reduces to 1.2. In their Fig. 1, DGS have shown fits to some data sets at interval of 1 yr. It can be seen from Fig. 2 in this paper that these correspond to times when $`\chi ^2`$ is close to a local minima. The fits shown by DGS correspond to filled squares in Fig. 2. Some of the intermediate data sets give bad fits as can be seen from Fig. 1. About 20% of the fits have $`\chi ^22`$ and most of these show clear deviation from the assumed form.

These oscillations in $`\mathrm{\Delta }R`$ arise because the expression is inadequate to fit the data and do not represent real variations in $`\mathrm{\Delta }R`$. Basically there is additional contribution to $`\delta \nu _{\mathrm{}}`$ which can not be represented by either of the terms in Eq. (1) and this gets projected on to the two terms giving spurious results. In particular, the oscillatory trend is also projected on to both terms in Eq. (1) and further, the division between the two terms is also a function of time. As a result, the oscillations get modulated and it is not straightforward to isolate the oscillatory part in the fitted results for $`\mathrm{\Delta }R`$. If all these oscillations are ignored and a linear function in time is fitted to the inferred radius variation as DGS have done, then the $`\chi ^2`$ per degree of freedom for this fit is 8.7. This fit is shown by the continuous line in the upper panel of Fig. 2 and corresponds to a radius variation of $`1.2`$ km yr-1, slightly less than that inferred by DGS. This difference is because of the additional data that have become available. A slightly smaller $`\chi ^2`$ of 5.5 per degree of freedom is obtained if instead, these points are fitted by a step function with a discontinuity around 1999.2. This fit is shown by heavy line in upper panel of Fig. 2. Looking at the top panel of Fig. 2, it appears that inferred radius has suddenly changed around 1999 and the step function fit appears to support this hypothesis. The large $`\chi ^2`$ is to be expected as Eq. (1) that is used to calculate $`\mathrm{\Delta }R`$, does not really fit the observed data at all times.

Ideally, one should remove the oscillatory component in frequency variation before considering longer period variations, but for simplicity, I consider fits at interval of 1 year which will be at the same phase of oscillatory component and further, select the phase such that the fits are the best in some sense. These points are marked by filled squares in Fig. 2. Table 1 gives the results obtained for these sets, which includes the $`\chi ^2`$ per degree of freedom as well as the average 10.7 cm radio flux during the time interval covered by the data set, which is a measure of solar activity. Looking at this table it is clear that the radius is not changing continuously. In fact, most of the radius variation has occurred between 1998.4 and 1999.4. Possible radius variation during 1996.4–1998.4 and 1999.4–2002.4 is less than 1 km. The solar activity did increase significantly during the period 1998.4–1999.4, but there has been comparable change in activity during other periods too. Hence that cannot explain the variation seen in Table 1. This happens to be the period during which contact with SOHO satellite was lost and it is most likely that this variation reflects systematic errors arising from changes in the MDI instrument that may have occurred during recovery of the satellite. Even if we assume that this variation is real, the rate of shrinking is not 1.5 km yr-1 as claimed by DGS, but something like 3 km yr-1 during 1998.4–1999.4 and essentially no variation at other times. Thus, any model to explain this frequency change by a radius variation must explain why there is little radius variation during most of the time and why all variation is confined to less than 1 yr at some intermediate phase of solar cycle.

In order to study the robustness of the inferred radius variation, I attempt the fits by restricting the mode set or the data sets and the results are summarized in Table 2. If high degree modes are neglected, then the fits to data using Eq. (1) improve to some extent, which is mainly because the total variation in frequencies reduces with degree. Nevertheless, the fit to linear variation in $`\mathrm{\Delta }R`$ is still bad and its slope keeps reducing as the upper limit on $`\mathrm{}`$ is reduced. Thus if only modes with $`140<\mathrm{}<250`$ are used the radius variation comes out to be $`0.57\pm 0.08`$ km yr-1 (with a $`\chi ^2=2.6`$), while if the upper limit on $`\mathrm{}`$ is reduced to 200, it becomes $`0.47\pm 0.17`$ km yr-1 ($`\chi ^2=2`$). In these cases if a step function is fitted the $`\chi ^2`$ comes out to be 2.1. Figure 3 shows the fits in these cases. It can be seen that the magnitude of possible discontinuity around 1999 reduces as the upper limit on $`\mathrm{}`$ is reduced and is hardly visible when the upper limit is reduced to $`\mathrm{}=200`$. In this case the errors in inferred radius are rather large and the one year oscillations are essentially wiped out by statistical fluctuations. The reduction in $`\chi ^2`$ is mainly due to increase in estimated errors in $`\mathrm{\Delta }R`$. Antia et al. (2001) have shown that the amplitude of oscillatory term reduces with decreasing $`\mathrm{}`$ and that also contributes to improvement in fits. If the $`\mathrm{}`$ range is reduced still further, the errors in fitted quantities are too large to make any meaningful deductions and the results are not shown in Table 2. This mainly arises because the frequency variation are very small and it is not easy to distinguish between contributions of the two terms in Eq. (1) over a small range of $`\mathrm{}`$.

If only the first term in Eq. (1) is used for fitting, then it is equivalent to taking an average of relative frequency variations over all modes. In that case, any variation in frequency will imply a variation in radius. This may not be realistic and the resulting fits (to Eq. (1)) are always bad. Since both components in the frequency variations increase steeply with $`\mathrm{}`$, as the upper limit on $`\mathrm{}`$ is reduced the inferred radius variation should reduce and the limiting value at the lowest frequency range would give an upper limit to any possible radius variation. In this case, since only one parameter is fitted, it is possible to get some fits with only a few low degree modes and hence it is possible to reduce the upper limit on $`\mathrm{}`$. The estimated rate of reduction in radius decreases from 2.2 km yr-1 when all modes are used to 0.74 km yr-1 when the upper limit on $`\mathrm{}`$ is reduced to 160. If the upper limit on $`\mathrm{}`$ is reduced still further, there are very few modes in some data sets and it is not possible to obtain any meaningful fits. However, recently the MDI data sets have been updated and the new data has more f-modes. With these revised data sets it is possible to reduce the upper limit to $`\mathrm{}=120`$ and the inferred radius variation comes out to be $`+0.08\pm 0.11`$ km yr-1 (Fig. 4). As demonstrated by Antia et al. (2001), at these low frequencies the oscillatory component in frequency variation is also not observed. This fit appears to be consistent with results of Basu & Antia (2002) who found that systematic error in MDI data is restricted to modes with $`\mathrm{}>120`$. When these modes are eliminated no radius variation is found. In order to enable a direct comparison with DGS, these revised data sets are not generally used in this work, but Table 2 lists the results obtained using these sets also over the full time interval. It is clear that the results are not significantly different from earlier results. If the fits are restricted to include only the post gap data sets (i.e., after 1999) and still use only the first term in Eq. (1), then the resulting rate of radius decrease is 1.1 km yr-1 with all modes and comes down to $`0.03\pm 0.22`$ km yr-1 when modes with $`\mathrm{}<160`$ are considered. This is to be expected as the frequency variation increases rapidly with degree. Thus the inferred radius variation is maximum when high degree modes are used and is negligible when only relatively low degree modes are used. If there is any component in frequency differences which corresponds to radius variation, the limiting radius variation will tend to this value when the modes in low $`\mathrm{}`$ range are used. Since this limiting value happens to be consistent with zero, we can conclude that there is no radius variation during 1999–2002 (or during the entire period when the revised MDI data are used). On the other hand, if the second term in Eq. (1) is also included and only data sets after 1999.0 are used then the resulting fit does not show any significant variation in radius irrespective of the upper limit on $`\mathrm{}`$. Similar results are obtained when only data sets before the gap (i.e., before 1998.6) are used. Thus it is clear that most of the inferred variation in solar radius has taken place during the gap in MDI data.

From the top panel in Fig. 2 it can be seen that in the pre-gap data (before 1999.0) the filled squares are close to the minimum in $`\mathrm{\Delta }R`$, while after the gap the filled squares are close to the maximum in $`\mathrm{\Delta }R`$. The reason for this flip is not clear. It could be due to instrumental variations during recovery or alternately it may be because before the gap the non-oscillatory component had larger amplitude as compared to the oscillatory component, while after the gap the amplitude of non-oscillatory component is less. There is also some variation in the number of modes and the set of modes between different data sets. If the expression fits the data well this variation will not matter, but unfortunately that is not true. It is difficult to assign much significance to these results as during the time of these data sets oscillatory component has maximum magnitude and this can give rise to spurious results in the fits. Comparing Figs. 2 and 4 it can be seen that the inferred radius using $`\mathrm{}<120`$ modes is in between the pre-gap and post-gap values using all modes. Thus it is clear that averaging over the oscillations in Fig. 2 does not give the correct estimate of radius. Ideally, one should subtract out the oscillatory component in the frequency variation itself before analyzing the data, as has been done by Antia et al. (2001).

## 4 DISCUSSION AND CONCLUSIONS

There have been a few claims in the recent times about variation in the solar radius from f-mode frequencies (Dziembowski et al. 1998, 2001; Antia et al. 2000). Unfortunately, the variation in f-mode frequencies is more complex than what is modeled in these studies. From the results shown in Table 2, it is clear that there is no evidence to support a decrease in radius since the result depends on the set of modes included in the study and the linear fit is generally bad. From Table 2 it can be seen that the estimated rate of radius change varies between $`2.2`$ and $`+0.1`$ km yr-1 depending on the range of $`\mathrm{}`$ and the number of terms in Eq. (1) used for fitting. This large variation is simply because Eq. (1) does not fit the observed data properly. More terms will be required to obtain a proper estimate of radius variation. In particular, when only $`\mathrm{}<120`$ f-modes are used the inferred radius variation always comes out to be consistent with zero. Similarly, when only data sets before (or after) the gap are used the inferred radius variation comes out to be negligible. Basically, the observed data sets do not appear to have any component of relative frequency variation which is independent of degree as would be required for radius variation (see Eq. (1)).

A large part of the inferred variations in solar radius is most probably due to instrumental effects. For example, oscillatory trend with a period of 1 year is probably due to orbital motion of the Earth and SOHO satellite around the Sun. Similarly, the sharp variation seen between 1998.4 and 1999.4, (depending on degree $`\mathrm{}`$) is most likely a result of changes in instrumental characteristics during the recovery of the SOHO satellite. These two effects can account for all claims of radius variation made earlier. All these instrumental errors need to be eliminated before any claim can be made about the cause of frequency variation. From the results presented above it is clear that after eliminating these instrumental effects there is no significant variation in the solar radius as determined by f-mode frequencies. Similar conclusion was obtained by Antia et al. (2001) using a more detailed analysis of both GONG and MDI data.

The systematic error between MDI data sets before and after recovery also manifests in other studies (Antia 2002; Basu & Antia 2002; Antia, Chitre & Thompson 2003). In particular, it is found that this systematic error is mostly confined to modes with $`\mathrm{}>120`$, which is consistent with the results in this work. When these modes are neglected no radius variation is found, while if these are included then we find varying amount of radius variation around 1999. If the radius variation is real, it cannot depend on $`\mathrm{}`$. It is quite likely that systematic errors are present in all MDI data sets, but their magnitude has changed during recovery.

If we assume that the inferred radius variation between 1998.4 and 1999.4 is of instrumental origin, then we can put some limits on radius variation. From Table 1 it can be seen that the maximum variation between the 3 points before the data gap is 1.3 km, while that in the 4 points after the gap is 1.1 km. Considering an error of about 0.6 km in each data point, this variation is consistent with no radius variation. This would suggest an upper limit comparable to error bars in each point on radius variation during half of the solar cycle. Similar conclusion can be obtained from Fig. 4, which shows the results using only $`\mathrm{}<120`$ modes which are not expected to be affected by the systematic error in MDI data. From these results we can put a conservative upper limit of 2 km on radius variations during the last 6 years. This would yield $`\mathrm{\Delta }R/R<3\times 10^6`$ and $`W<0.003`$ as the ratio of radius to luminosity variation. Such a small value should favor models involving changes in the outer layers to explain the observed luminosity variations (Gough 1981; Däppen 1983; Balmforth, Gough & Merryfield 1996; Gough 2001). Of course, the value of $`W`$ in these models is determined by radius variation at the photosphere, while f-modes are sensitive to variations at depths of about 1–12 Mm. But if the photospheric radius variations are much larger than those inferred by the f-modes, then the cause is most likely to be near the surface. Emilio et al. (2000) find a much larger increase in photospheric radius by about $`15\pm 2`$ km during the solar cycle. These direct measurements from MDI are also affected by a number of systematic errors and as they have pointed out this value should be regarded as an upper limit to radius variation. Thus our results from f-mode frequencies which effectively measure the solar radius in the subsurface layers, are probably not inconsistent with these measurements.

## ACKNOWLEDGMENTS

This work utilizes data obtained by the Solar Oscillations Investigation / Michelson Doppler Imager on the Solar and Heliospheric Observatory (SOHO). SOHO is a project of international cooperation between ESA and NASA.