text:
# References

SINP-TNP/03-03

Hidden Degeneracy in the Brick Wall Model of Black Holes

Kumar S. Gupta111Email: gupta@theory.saha.ernet.in

Theory Division
Saha Institute of Nuclear Physics
1/AF Bidhannagar
Calcutta - 700064, India
Siddhartha Sen222Email: sen@maths.tcd.ie

School of Mathematics
Trinity College
Dublin , Ireland
and

Department of Theoretical Physics
Indian Association for the Cultivation of Science
Calcutta - 700032, India

## Abstract

Quantum field theory in the near-horizon region of a black hole predicts the existence of an infinite number of degenerate modes. Such a degeneracy is regulated in the brick wall model by the introduction of a short distance cutoff. In this Letter we show that states of the brick wall model with non zero energy admit a further degeneracy for any given finite value of the cutoff. The black hole entropy is calculated within the brick wall model taking this degeneracy into account. Modes with complex frequencies however do not exhibit such a degeneracy.

February 2003

PACS : 04.70.Dy

1. Introduction

The event horizon plays a central role in various aspects of black hole physics. A simple illustration of this role of the horizon is provided by the brick wall model proposed by ’t Hooft . This model describes the low energy dynamics of a scalar field the in the background of a Schwarzschild black hole. The density of states near the black hole horizon is however divergent, due to the presence of an infinite number of degenerate modes in the near-horizon region . Such a divergence is controlled in the brick wall model through the introduction of a short distance cutoff (i.e. a brick wall). For every finite value of the cutoff, this model provides a well defined expression for the density of states and entropy. The entropy of the black hole as determined from this model arises essentially from the near-horizon features of the system.

The above discussion suggests that the ambiguity in determining the density of states in the brick wall model is completely removed by the short distance cutoff. In this Letter we show that even with such a finite cutoff, there is still an additional degeneracy present in the system. This degeneracy is revealed by a study of the brick wall model in the near-horizon region of the black hole. For a given value of the cutoff, the allowed quantum states of this model are labelled by energy and other quantum numbers. We shall show that for any given finite value of the cutoff, every such allowed quantum state with non-zero energy has infinite degeneracy.

The infinite degeneracy found is expected to contribute to the black hole entropy. In order to analyze this contribution, we show that a dimensionless parameter can be introduced to regulate the degeneracy. The resulting expression for the black hole entropy is found to depend on the ratio of this parameter to the brick wall cutoff. The assumption of finiteness of this parameter is found to be consistent within the framework of the brick wall model.

Our analysis can be extended to the case of complex energy values. Such states are analogous to the quasinormal modes of the Schwarzschild black hole . Analysis of the equation satisfied by these modes shows that the degeneracy disappears due to the presence of complex frequencies.

This Letter is organized as follows. In Section 2 we show how the degeneracy appears in the brick wall model. Similar analysis for modes with complex frequencies is performed in Section 3 and it is shown that such modes admit no degeneracy. In Section 4 we discuss how the density of states and black hole entropy are modified in presence of the new degeneracy. Section 5 concludes this Letter with some discussions.

2. Hidden Degeneracy in the Brick Wall Model

The brick wall model  describes the low energy quantum dynamics of a scalar field $`\varphi `$ in the background of a massive Schwarzschild black hole of finite mass $`M`$. For simplicity we shall take the scalar field to be massless. In the spherical polar coordinates coordinates $`(r,\theta ,\varphi )`$, the field theory description is assumed to be valid in the region $`r>R+h`$ where $`R=2M`$ is the radius of the horizon and the brick wall cutoff $`h(>0)`$ is small compared to $`R`$. The field $`\varphi `$ is assumed to satisfy the boundary condition

$$\varphi (r=R+h,\theta ,\varphi ,t)=0.$$
(2.1)
In addition, the whole system is assumed to be in a box of radius $`L`$ which provides an infrared cutoff with the associated boundary condition

$$\varphi (r,\theta ,\varphi ,t)=0\mathrm{for}rL.$$
(2.2)
In the background of a Schwarzschild black hole, the field equation for modes with angular momentum $`l`$ and energy $`E`$ is given by

$$\left(1\frac{2M}{r}\right)^1E^2\varphi +\frac{1}{r^2}\frac{}{r}\left(r(r2M)\frac{\varphi }{r}\right)\frac{l(l+1)}{r^2}\varphi =0.$$
(2.3)
In order to study the near horizon properties of the brick wall model, it is useful to introduce a new coordinate $`xrR`$. In the near-horizon region, i.e. when $`xR`$, Eqn. (3) can be written as

$$\frac{^2\varphi }{x^2}+\frac{1}{x}\frac{\varphi }{x}+\frac{R^2E^2}{x^2}\varphi \frac{l(l+1)}{Rx}\varphi =0.$$
(2.4)
In terms of a new field $`\chi `$ defined by

$$\chi =\sqrt{x}\varphi ,$$
(2.5)
Eqn. (2.4) can be written as

$$\frac{^2\chi }{x^2}+\frac{1}{x^2}\left[\frac{1}{4}+R^2E^2\right]\chi \frac{l(l+1)}{Rx}\chi =0.$$
(2.6)
For small values of $`x`$, the angular term in the field Eqn. (2.6) can be ignored as has been noted elsewhere as well . Thus the KG equation in the near-horizon region can be written as

$$\frac{^2\chi }{x^2}+\frac{1}{x^2}\left[\frac{1}{4}+R^2E^2\right]\chi =0.$$
(2.7)
The field equation (2.3) has been analyzed by ’t Hooft using the WKB method . Such an analysis with the brick wall boundary conditions is assumed to provide a discrete set of eigenvalues. Let us consider Eqn. (2.7) with such a generic eigenvalue $`E`$, the latter being now treated as a parameter already determined by the WKB analysis. We can think of as Eqn. (2.7) as a Schrodinger equation in presence of a singular potential  with $`\chi `$ now being interpreted as the zero eigenvalue solution of such an equation. The central question that we address here can now be stated as follows: for the allowed values of the parameter $`E`$ and in the presence of the brick wall cutoff $`h`$, does the singular quantum mechanical problem described by Eqn. (2.7) possess a unique solution or is there any additional degeneracy which is perhaps not manifest? In order to address this issue, it is useful to start by considering a different equation given by

$$\frac{^2\chi }{x^2}+\frac{1}{x^2}\left[\frac{1}{4}+R^2E^2\right]\chi =\eta ^2\chi .$$
(2.8)
In what follows, we shall interprete Eqn. (2.7) as the $`\eta 0`$ limit of Eqn. (2.8). We shall show below that in the limit when $`\eta 0`$, Eqn. (2.8) admits an infinite number of solutions even in the presence of the brick wall cutoff. More precisely, we will see that $`\eta =0`$ is an accumulation point for the solutions of Eqn. (2.8).

We now proceed to solve Eqn. (2.8) with the brick wall boundary conditions. Eqn. (2.8) can be related to the equation for modified Bessel functions with imaginary order. The two linearly independent solutions $`u_1`$ and $`u_2`$ of Eqn. (2.8) can be taken as

$`u_1`$ $`=`$ $`\sqrt{x}I_{i\lambda }(\eta x),`$ (2.9)
$`u_2`$ $`=`$ $`\sqrt{x}I_{i\lambda }(\eta x),`$ (2.10)
where $`\lambda =RE`$ and $`I`$ denotes the modified Bessel function . The general solution of Eqn. (2.8) which vanishes at $`x=h`$ can thus be written as

$$\chi (x)=A\sqrt{x}\left[I_{i\lambda }(\eta x)I_{i\lambda }(\eta h)I_{i\lambda }(\eta x)I_{i\lambda }(\eta h)\right],$$
(2.11)
where $`A`$ is a constant. When $`rL`$, the coordinate $`x(LR)`$. In this limit and under the condition that $`(LR)h`$, the leading divergent behaviour of $`\chi (x)`$ is given by

$$\chi (x)=\frac{A\mathrm{e}^{\eta x}}{\sqrt{2\pi \eta }}\left[I_{i\lambda }(\eta h)I_{i\lambda }(\eta h)\right].$$
(2.12)
For $`x(LR)`$, Eqn. (2.2) requires that $`\chi (x)`$ must vanish. From Eqn. (2.12) we see that this is possible only if

$$I_{i\lambda }(\eta h)=I_{i\lambda }(\eta h).$$
(2.13)
Let us now recall that we are interested in the situation where $`h0`$, i.e. where the brick wall is very close to the horizon. When $`h0`$, $`I_{i\lambda }(\eta h)`$ behaves as

$$I_{i\lambda }(\eta h)\frac{(\eta h)^{i\lambda }}{2^{i\lambda }}\frac{1}{\mathrm{\Gamma }(1+i\lambda )}.$$
(2.14)
Using Eqns. (2.13) and (2.14) we see that the eigenvalues of Eqn. (2.8) are given by

$$\eta _p^2=\frac{1}{h^2}\mathrm{e}xp\left(2\mathrm{l}n2+\frac{2\theta (\lambda )}{\lambda }+2p\frac{\pi }{\lambda }\right),$$
(2.15)
where $`pZ`$ is an arbitrary integer and the function $`\theta (\lambda )`$ is the argument of the Gamma function $`\mathrm{\Gamma }(1+i\lambda )`$. As mentioned before, we are interested only in the zero eigenvalue solutions of Eqn. (2.8). From Eqn. (2.15) we see that $`\eta _p0`$ when $`p\mathrm{}`$, or equivalently when $`p=n`$ where $`n`$ is an arbitrarily large positive integer. In that case, Eqn. (2.15) can be written as

$$\eta _n^2=\frac{1}{h^2}\mathrm{e}xp\left(2n\frac{\pi }{\lambda }\right).$$
(2.16)
From Eqn. (2.16) we see that as $`n\mathrm{}`$, there are an infinite number of zero modes solutions of Eqn. (2.8). In other words, the zero is an accumulation point for the eigenvalues of Eqn. (2.8).

We shall now point out several implications of the above analysis : 1) The energy eigenvalue $`E`$ obtained by the WKB analysis of Eqn. (2.3) serves as a parameter in the above analysis. The degeneracy obtained above exists for every real value of $`E`$ except for $`E=0`$. In the latter case we have $`\lambda =0`$ and Eqn. (2.13) reduces trivially to an identity.
2) The infinite number of zero eigenvalue states of Eqn. (2.8) exist for every finite value of the brick wall cutoff parameter $`h`$. This degeneracy is therefore not removed by the brick wall cutoff.

The above analysis thus shows that for any given real value of $`E0`$ and for any finite value of the brick wall cutoff parameter $`h`$, there exist an infinite number of zero eigenvalue solutions of Eqn (2.8), which can be interpreted as solutions of Eqn. (2.7). It is in this sense that we claim about the existence of the additional degeneracy in the near-horizon region as described by the brick wall model. Finally it may be noted that the above result would continue to hold for massive fields as well as the field equations become independent of the mass term in the hear-horizon region of the black hole .

3. Analysis of Modes with Complex Frequencies

In this Section we shall analyze Eqn. (2.7) for the case of complex frequencies. For the purpose of this Section, we use the notation $`\mathrm{\Omega }`$ instead of $`E`$. Eqn. (2.7) is therefore replaced by

$$\frac{^2\chi }{x^2}+\frac{1}{x^2}\left[\frac{1}{4}+R^2\mathrm{\Omega }^2\right]\chi =0,$$
(3.1)
where $`\mathrm{\Omega }\omega _1+i\omega _2`$ is taken to be complex. It is tempting to identify $`\chi `$ in Eqn. (3.1) with the quasinormal mode amplitude for the Schwarzschild black hole. Indeed, the equation for the quasinormal modes for the Schwarzschild case  reduces to Eqn. (3.1) in the near-horizon region. There are however certain differences in the boundary conditions obeyed by the quasinormal modes  to the ones that are used in the brick wall model. In what follows we shall analyze Eqn. (3.1) with the brick wall boundary conditions.

As described in Section 2, we interprete Eqn. (3.1) as the $`\gamma 0`$ limit of the equation

$$\frac{^2\chi }{x^2}+\frac{1}{x^2}\left[\frac{1}{4}+R^2\mathrm{\Omega }^2\right]\chi =\gamma ^2\chi .$$
(3.2)
As before, we are interested only in solutions of Eqn. (3.2) with $`\gamma 0`$. The general solution of Eqn. (3.2) which vanishes at $`x=h`$ is given by

$$\chi (x)=B\sqrt{x}\left[I_{i\mu }(\gamma x)I_{i\mu }(\gamma h)I_{i\mu }(\gamma x)I_{i\mu }(\gamma h)\right],$$
(3.3)
where $`\mu =R\mathrm{\Omega }`$. When $`x(LR)`$ and $`(LR)h`$, the wavefunction is Eqn. (3.3) can vanish only if

$$I_{i\mu }(\gamma h)=I_{i\mu }(\gamma h).$$
(3.4)
In the limit when $`h0`$, Eqn (3.4) gives

$$(\gamma h)^{2i\mu }2^{2i\mu }=\frac{\mathrm{\Gamma }(1+i\mu )}{\mathrm{\Gamma }(1i\mu )}.$$
(3.5)
Since $`\mathrm{\Omega }=\omega _1+i\omega _2`$, the gamma functions in Eqn. (3.5) can be expressed in the polar form as

$`\mathrm{\Gamma }(1+i\mu )=\mathrm{\Gamma }(1R\omega _1+iR\omega _1)`$ $``$ $`\xi _1\mathrm{e}^{i\theta _1}`$ (3.6)
$`\mathrm{\Gamma }(1\mu )=\mathrm{\Gamma }(1R\omega _2iR\omega _1)`$ $``$ $`\xi _2\mathrm{e}^{i\theta _2}.`$ (3.7)
Using Eqns. (3.6) and (3.7), the imaginary and real parts of Eqn. (3.5) can be written as

$`2R\omega _1\left(\mathrm{l}n\gamma h+\mathrm{l}n2\right)`$ $`=`$ $`(\theta _1\theta _2)+2p\pi `$ (3.8)
$`2R\omega _2\left(\mathrm{l}n\gamma h+\mathrm{l}n2\right)`$ $`=`$ $`\mathrm{l}n{\displaystyle \frac{\xi _1}{\xi _2}},`$ (3.9)
where $`pZ`$ is an arbitrary integer. However, the zero eigenvalue solutions of Eqn. (3.2) are obtained only when $`p=n`$ where $`n`$ is an arbitrarily large positive integer. In that case Eqn. (3.8) can be written as

$$\mathrm{l}n\gamma _nh=\frac{n\pi }{R\omega _1}.$$
(3.10)
Using Eqns. (3.9) and (3.10) we get

$$\left(\frac{\omega _1}{2\pi \omega _2}\right)\mathrm{l}n\frac{\xi _1}{\xi _2}=n.$$
(3.11)
Since $`n`$ is an arbitrarily large, for fixed values of $`\omega _1`$ and $`\omega _2`$ Eqn. (3.11) has no solution. We can therefore conclude that there is no additional degeneracy for the case of the modes with complex frequencies.

4. Black Hole Entropy

As discussed in Section 2, the solutions of Eqn. (2.8) for any non-zero value of $`E`$ admit an infinite degeneracy due to the near-horizon effects. Let us also recall that the black hole entropy found in the brick wall model arises essentially from the near-horizon features of the model. It is therefore expected that the degeneracy found in Section 2 would contribute to the black hole entropy, as they are both related to the physics in the near-horizon region. Below we shall discuss the effect of this degeneracy on the density of states and black hole entropy.

We shall first discus how the degeneracy found in Section 2 affects the density of states. For this purpose, consider two eigenvalues $`\eta _{n_1}`$ and $`\eta _{n_2}`$ of the Eqn. (2.8) where the two positive integers $`n_1,n_2\mathrm{}`$. For definiteness let us assume that $`n_2>n_1`$. From Eqn. (2.16) we get

$$\mathrm{\Delta }nn_2n_1=\frac{\lambda }{\pi }\mathrm{l}n\left(\frac{\eta _{n_1}}{\eta _{n_2}}\right)(\lambda =RE).$$
(4.1)
The quantity $`\mathrm{\Delta }n`$ provides an expression for the degeneracy of the quantum state with energy $`E`$. The logarithm on the r.h.s. of Eqn. (4.1) contains the ratio of two numbers both of which are tending towards zero. We shall assume that the logarithm of this ratio tends to a positive constant $`c`$ as $`n_1,n_2\mathrm{}`$. This is an ad hoc assumption which is nevertheless consistent within the framework of this model. Eqn. (4.1) can now be written as

$$\mathrm{\Delta }n=\frac{c\lambda }{\pi }=\frac{cRE}{\pi }.$$
(4.2)
The above expression vanishes for $`E=0`$, which is consistent with the observation made at the end of Section 2. The density of states calculated in the brick wall model would now have to be multiplied with this degeneracy factor. The corresponding expression of the free energy $`F`$ obtained from the brick wall model  is now modified to read as

$$\pi \beta F=\frac{2}{3}\frac{(2M)^4}{h}_0^{\mathrm{}}\frac{(\mathrm{\Delta }n)E^3dE}{e^{\beta E}1}.$$
(4.3)
The corresponding expression of the black hole entropy is now obtained as

$$S=\frac{80}{\pi ^3}\frac{c}{h}\frac{(2M)^4}{\beta ^3}.$$
(4.4)
On the other hand, the application of quantum mechanical scattering theory to this system gives the black hole entropy as

$$S=\frac{4\pi M^2}{G},$$
(4.5)
where $`G`$ is the Newton’s constant in units chosen in Ref. . Comparing Eqns. (4.4) with (4.5) we get

$$h=\frac{320Gc}{\pi ^7M}.$$
(4.6)
We therefore see within the framework presented above, the main effect of the additional degeneracy is to modify the expression of the brick wall cutoff parameter $`h`$ (compare with Eqn. (8.19) of Ref.  with $`N=1`$). Finally, the requirement that physical distance of the brick wall from the horizon be of the order of Planck length , leads to a finite value of the constant $`c`$ which is consistent with assumption made above.

5. Conclusion

In this Letter we have shown that there is a hidden degeneracy in the brick wall model in the near-horizon region of the black hole. This degeneracy has been found by the analysis of the zero eigenvalue solutions of Eqn. (2.8) which can be identified with the solutions of Eqn. (2.7). Usually, any reasonable quantum field theory is expected to predict an infinite number of modes in the near-horizon region of a black hole . The divergence from these modes is handled in the brick wall model by the introduction of a short distance cutoff. The degeneracy that we have discussed is different from the usual one in the sense that it is not removed by the brick wall cutoff. We have in fact shown that for any finite value of the brick wall cutoff, every quantum state of the system with non-zero energy still admits an infinite degeneracy.

The additional degeneracy discussed above is essentially a near-horizon feature and is thus expected to affect the black hole entropy. We have estimated the effect of this degeneracy in terms of a parameter $`c`$, which has been assumed to be finite. The contribution of this degeneracy to the density of states is finite and non-zero except for $`E=0`$. The corresponding expression for the free energy and entropy as calculated within the brick wall depends on the ratio $`\frac{c}{h}`$. The requirement that the physical distance of the brick wall from the horizon is of the order of Planck length leads to a finite value of $`c`$. The assumption of finiteness of $`c`$ is thus consistent within the framework of the brick wall model.

We have also analyzed the system with complex frequencies. Such a situation arises in the description of quasinormal modes of Schwarzschild black holes. It has been shown that for finite values of the brick wall cutoff there is no further degeneracy in this case.

It has been argued in Ref.  that for a wide variety of quantum fields including massive ones, the field equation in the near horizon region looks essentially the same. The analysis presented in this paper is thus expected to hold for all such fields and also for a variety of black hole metrics with similar near-horizon structures. It is plausible that the degeneracy discussed in this Letter is related to the conformal symmetry which is present in the near-horizon region of a wide class of black holes .


format: text

text:
# New treatment of breakup continuum in the method of continuum discretized coupled channels

## I Introduction

The method of continuum discretized coupled channels (CDCC) has been successful in describing nuclear reactions including weakly bound projectiles CDCC-review1 ; CDCC-review2 ; Yahiro1 ; yahiro2 ; sakuragi1 ; sakuragi2 ; sakuragi3 ; sakuragi4 ; iseri1 ; Sakuragi ; Surrey ; Ogata . CDCC has been attracting much attention since the advent of experiments with radioactive beams, because projectile breakup processes are essential to many of such reactions. CDCC plays an important role in the spectroscopic studies of radioactive nuclei through the nuclear reactions involving such nuclei.

In CDCC for reactions with a projectile consisting of two fragments, the states of the projectile are classified by the linear and the angular momenta, $`k`$ and $`\mathrm{}`$, of relative motion of the two fragments of the projectile, which are truncated by $`kk_{\mathrm{max}}`$ and $`\mathrm{}\mathrm{}_{\mathrm{max}}`$. The truncation is the most basic assumption in CDCC, and it is confirmed that calculated $`S`$-matrix elements converge for sufficiently large $`k_{\mathrm{max}}`$ and $`\mathrm{}_{\mathrm{max}}`$ CDCC-review1 ; Yahiro1 ; Piya . It has been shown that CDCC is the first-order approximation to the distorted Faddeev equations, and corrections to the converged CDCC solution are negligible within the region of space in which the reaction takes place CDCC-foundation .

As a consequence of the truncation, the integral equation form of the equations of coupled channels, derived from the three-body Schrödinger equation, has a compact kernel, satisfying a necessary condition for iterative solutions CDCC-foundation . In practice, however, the coupled channels equations thus obtained are impossible to be solved because of the continuously infinite number of coupled breakup channels. The problem is solved by discretizing the $`k`$ continuum. The discretization leads the coupled equations to a set of differential equations with a finite number of channels.

As for the discretization, three methods have been proposed so far: the average (Av) YK1 ; Yahiro1 ; CDCC-review1 ; CDCC-review2 , the midpoint (Mid) CDCC-review2 ; Piya , and the pseudo-state (PS) CDCC-review1 ; PS1 ; PS2  methods. In the Av and Mid methods, the $`k`$-continuum is divided into a finite number of bins. In the former the continuum channels within each bin are averaged into a single channel, while in the latter they are represented by the channel at a midpoint of the bin. It has been confirmed that calculated $`S`$-matrix elements converge as the width $`\mathrm{\Delta }`$ of the bins is decreased, and also that the two methods yield the same converged $`S`$-matrix elements CDCC-review1 ; Yahiro1 ; Piya . From practical point of view, the Av method is more convenient than the Mid one, since the former requires less numerical works than the latter. The Av method therefore is most widely used in practice.

In the PS method, the breakup states are described by superpositions of $`L^2`$-type basis functions. The wave functions of such pseudo breakup state have wrong asymptotic forms. For this reason, the PS method was mainly used in the past to describe virtual breakup processes in the intermediate stage of elastic scattering and ($`d,p`$) reactions. Very recently, the applicability of the PS method for real breakup reactions for a specific set of basis functions. It was shown that the PS method well simulates the angular distribution of the total breakup cross section calculated with the Av method PS1 . In the analysis, discrete breakup cross sections calculated with the PS method are transformed into continuous one, by assuming a specific form for the $`k`$ distribution of the continuous cross section, that is, a histogram with widths estimated in a reasonable way.

The purpose of the present paper is to propose a new PS method for projectile breakup reactions. If the basis functions form in the good approximation a complete set in the configuration space which is important for breakup processes, we can formulate a method of interpolation for generating continuous breakup $`S`$-matrix elements from the discrete ones provided by the PS method, following the formulation of Tostevin  for the Av method in which the completeness of discretized states is well satisfied if the discretization is enough precise. The method does not assume any form a priori for $`k`$ distributions of continuous breakup $`S`$-matrix elements. The method is independent of the type of the basis functions used for the PS method. Only condition is that the basis functions constitute an approximate complete set in the wide range of $`k`$ and its conjugate coordinate $`r`$ which are important for the breakup processes. As basis functions satisfying this condition, we propose two types of them. One is ordinary Gaussian functions K-G , which we refer to as real-range Gaussian function in the present paper. The other is a natural extension of that, complex-range Gaussian function H-Ka-Ki , i.e, Gaussian functions with the complex range parameters the precise definition of which is given in a later section.

In Section II, we recapitulate CDCC based on both the Av and the PS methods of discretization. In Section III, we present a method of interpolation to get continuous breakup $`S`$-matrix elements from discrete ones calculated with the PS method, and introduce the real and the complex-range Gaussian basis. In Section IV, the validity of the present PS method is tested and justified for two realistic cases, $`d`$+58Ni scattering at 80 MeV and 6Li+40Ca scattering at 156 MeV. In Section V, the potentiality of the present PS method for four-body breakup reactions is discussed. Section VI gives a summary.

## II The Method of Continuum Discretized Coupled Channels

We consider a reaction of a weakly bound projectile (B) impinging on a target nucleus (A). We treat a simple system shown in Fig. 1 in which the projectile, B, is composed of two particles (b and c) and the target is inert. The three-body system is described by a model Hamiltonian $`H=H_{\mathrm{bc}}+K_R+U,`$ where $`H_{\mathrm{bc}}=K_r+V_{\mathrm{bc}}(𝐫)`$ and $`U=U_{\mathrm{bA}}(𝐫_{\mathrm{bA}})+U_{\mathrm{cA}}(𝐫_{\mathrm{cA}})`$. Vector $`𝐫`$ is the relative coordinate between b and c, $`𝐑`$ the one between the center-of-mass of the b-c pair and the center-of-mass of A, and $`𝐫_{\mathrm{XY}}`$ denotes the relative coordinate between two particles X and Y. Operators $`K_r`$ and $`K_R`$ are kinetic energies associated with $`𝐫`$ and $`𝐑`$, respectively. $`V_{\mathrm{bc}}(𝐫)`$ is the interaction between b and c. The interaction $`U_{\mathrm{bA}}`$ ($`U_{\mathrm{cA}}`$) between b (c) and A is taken to be the optical potential for b+A (c+A) scattering. For simplicity, the spin dependence of the interactions is neglected. Furthermore, the Coulomb part of $`U`$ is treated approximately as a function only of $`𝐑`$, i.e., we neglect Coulomb breakup processes, and focus our attention on nuclear breakup.

In CDCC, the three-body wave function $`\mathrm{\Psi }_{JM}`$, with the total angular momentum $`J`$ and its projection $`M`$ on $`z`$-axis, is expanded in terms of the orthonormal set of eigenstates $`\mathrm{\Phi }`$ of $`H_{\mathrm{bc}}`$:

$`\mathrm{\Psi }_{JM}(𝐫,𝐑)={\displaystyle \underset{L}{}}𝒴_{JM}^{\mathrm{}_0,L}\mathrm{\Phi }_0(r)\chi _{\mathrm{}_0LJ}(P_0,R)/R`$
$`+{\displaystyle \underset{\mathrm{},L}{}}𝒴_{JM}^{\mathrm{},L}{\displaystyle _0^{\mathrm{}}}\mathrm{\Phi }_{\mathrm{}}(k,r)\chi _{\mathrm{}LJ}(P,R)/R𝑑k,`$ (1)
where

$$𝒴_{JM}^{\mathrm{},L}=[i^{\mathrm{}}Y_{\mathrm{}}(\mathrm{\Omega }_r)i^LY_L(\mathrm{\Omega }_R)]_{JM}.$$
(2)
For simplicity, we assume that the b+c system has one bound state $`\mathrm{\Phi }_0(r)`$ with angular momentum $`\mathrm{}_0`$ and continuum states $`\mathrm{\Phi }_{\mathrm{}}(k,r)`$ with linear momentum $`k`$ and angular momentum $`\mathrm{}`$, both ranging from zero to infinity. The $`\mathrm{\Phi }_{\mathrm{}}`$ are real functions normalized to the $`\delta `$-function in $`k`$ Yahiro1 . The projectile B is initially in the bound state. The coefficient $`\chi _{\mathrm{}LJ}`$ ($`\chi _{\mathrm{}_0LJ}`$) of the expansion describes a center-of-mass motion of the b-c pair relative to A in the state $`\mathrm{\Phi }_{\mathrm{}}`$ with the linear and orbital angular momenta $`P`$ $`(P_0)`$ and $`L`$, respectively.

In CDCC, the sum over $`\mathrm{}`$ is truncated by $`\mathrm{}\mathrm{}_{\mathrm{max}}`$ and the $`k`$ integral by $`kk_{\mathrm{max}}`$. For each $`\mathrm{}`$, furthermore, the continuum states from $`k=0`$ to $`k_{\mathrm{max}}`$ are discretized into a finite number of states, with the wave function $`\widehat{\mathrm{\Phi }}_i\mathrm{}(r)`$ of the i-th state corresponding to the momentum $`\widehat{k}_i`$. Details of the discretization methods are described in the next section.

After the truncation and the discretization, $`\mathrm{\Psi }_{JM}`$ is reduced to an approximate one,

$`\mathrm{\Psi }_{JM}^{\mathrm{CDCC}}={\displaystyle \underset{L}{}}𝒴_{JM}^{\mathrm{}_0,L}\mathrm{\Phi }_0(r)\widehat{\chi }_{\gamma _0}(P_0,R)/R`$
$`+{\displaystyle \underset{l=0}{\overset{l_{\mathrm{max}}}{}}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{L}{}}𝒴_{JM}^{\mathrm{},L}\widehat{\mathrm{\Phi }}_i\mathrm{}(r)\widehat{\chi }_\gamma (\widehat{P}_i,R)/R,`$ (3)
where

$`\widehat{\chi }_{\gamma _0}(P_0,R)`$ $`=`$ $`\chi _{\gamma _0}(P_0,R),\gamma _0=(0,\mathrm{}_0,L,J),`$
$`\widehat{\chi }_\gamma (\widehat{P}_i,R)`$ $`=`$ $`W_\gamma \chi _\gamma (\widehat{P}_i,R),\gamma =(i,\mathrm{},L,J).`$
On the right hand side of Eq. (3), the first term represents the elastic channel denoted by $`\gamma _0`$ and the second one corresponds to the discretized breakup channels, each denoted by $`\gamma `$. The weight factor $`W_\gamma `$ depends on the discretization method used. Each pair of momenta, $`(\widehat{k}_i,\widehat{P}_i)`$, satisfies the total energy conservation:

$$E=\mathrm{}^2P_0^2/2\mu _{\mathrm{AB}}+ϵ_0=\mathrm{}^2\widehat{P}_i^2/2\mu _{\mathrm{AB}}+ϵ_i,$$
(4)
where $`ϵ_0`$ and $`ϵ_i=\mathrm{}^2\widehat{k}_i^2/2\mu _{\mathrm{bc}}`$ are the energies of the ground and the continuum states, respectively. Inserting Eq. (3) into the three-body Schrödinger equation, $`(HE)\mathrm{\Psi }_{JM}=0`$, leads to a set of coupled differential equations for $`\widehat{\chi }_{\gamma _0}(P_0,R)`$ and $`\widehat{\chi }_\gamma (\widehat{P}_i,R)`$:

$`\left[{\displaystyle \frac{d^2}{dR^2}}+\widehat{P}_i^2{\displaystyle \frac{L(L+1)}{R^2}}{\displaystyle \frac{2\mu _{\mathrm{AB}}}{\mathrm{}^2}}V_{\gamma \gamma }(R)\right]\widehat{\chi }_\gamma (\widehat{P}_i,R)`$
$`={\displaystyle \underset{\gamma ^{^{}}\gamma }{}}{\displaystyle \frac{2\mu _{\mathrm{AB}}}{\mathrm{}^2}}V_{\gamma \gamma ^{^{}}}(R)\widehat{\chi }_\gamma ^{^{}}(\widehat{P}_i^{^{}},R)`$ (5)
for all $`\gamma `$ including $`\gamma _0`$, where $`\widehat{k}_0=\sqrt{2\mu _{\mathrm{bc}}ϵ_0}/\mathrm{}`$ and $`\widehat{P}_0=P_0`$. The coupling potentials $`V_{\gamma \gamma ^{^{}}}(R)`$ are obtained as

$$V_{\gamma \gamma ^{^{}}}(R)=𝒴_{JM}^{\mathrm{},L}\widehat{\mathrm{\Phi }}_i\mathrm{}(r)|U|𝒴_{JM}^{\mathrm{}^{^{}},L^{^{}}}\widehat{\mathrm{\Phi }}_i^{^{}}\mathrm{}^{^{}}(r)_{𝐫,\mathrm{\Omega }_R}.$$
(6)
The coupled equations are solved, under the asymptotic boundary condition

$$\widehat{\chi }_\gamma (\widehat{P}_i,R)u_L^{()}(\widehat{P}_i,R)\delta _{\gamma ,\gamma _0}\sqrt{\frac{\widehat{P}_i}{\widehat{P}_0}}\widehat{S}_{\gamma ,\gamma _0}u_L^{(+)}(\widehat{P}_i,R).$$
(7)
Here $`u_L^{()}(\widehat{P}_i,R)`$ and $`u_L^{(+)}(\widehat{P}_i,R)`$ are incoming and outgoing Coulomb wave functions with the momentum $`\widehat{P}_i`$, and $`\widehat{S}_{\gamma ,\gamma _0}`$ is the $`S`$-matrix element for the transition from the initial channel $`\gamma _0`$ to $`\gamma `$.

## III Discretization of $`k`$ continuum

Among the three methods of discretization of the $`k`$ continuum, the relation between the average (Av) and the midpoint (Mid) methods has already been clarified Piya . The present discussion therefore is focused on the relation between the Av and the pseudo-state (PS) methods.

### III.1 The average method

In the Av method, the $`k`$-continuum \[0, $`k_{\mathrm{max}}`$\], for each $`\mathrm{}`$, is divided into a finite number of bins, each with a width $`\mathrm{\Delta }_i\mathrm{}=k_ik_{i1}`$, and the continuum breakup states in the $`i`$-th bin are averaged with a weight function $`f_i\mathrm{}(k)`$ CDCC-review1 ; CDCC-review2 . The resultant orthonormal state is described as

$$\widehat{\mathrm{\Phi }}_i\mathrm{}(r)=\frac{1}{W_\gamma }_{k_{i1}}^{k_i}\mathrm{\Phi }_{\mathrm{}}(k,r)f_i\mathrm{}(k)𝑑k\text{(for Av) ,}$$
(8)
then the weight factor $`W_\gamma `$ is given by

$$W_\gamma ^2=_{k_{i1}}^{k_i}[f_i\mathrm{}(k)]^2𝑑k.$$

For a bin far from a resonance, it is natural to set $`f_i\mathrm{}(k)=1`$, so that $`W_\gamma =\sqrt{\mathrm{\Delta }_i\mathrm{}}`$ since $`\mathrm{\Phi }_{\mathrm{}}(k,r)`$ changes smoothly with $`k`$. On the other hand, $`\mathrm{\Phi }_{\mathrm{}}(k,r)`$ changes rapidly across the resonance. One way of coping with this situation is to take $`\mathrm{\Delta }_i\mathrm{}`$ much smaller than the width of the resonance so that $`\mathrm{\Phi }_{\mathrm{}}(k,r)`$ does not change much within individual bins. This, however, make the number of bins large. Alternatively, one can take a single bin which contains the whole resonance peak and use a weight function of Breit-Wigner type CDCC-review1 ; sakuragi1 ; sakuragi2 ; sakuragi4 ; Sakuragi ,

$$f_i\mathrm{}(k)=\left|\frac{i\mathrm{\Gamma }/2}{ϵ(k)ϵ_{\mathrm{res}}+i\mathrm{\Gamma }/2}\right|,$$
(9)
where $`ϵ(k)`$ is a continuous intrinsic energy of the b+c system. The discretized intrinsic energy, $`ϵ_i=\mathrm{}^2\widehat{k}_i^2/2\mu _{\mathrm{bc}}`$, corresponding to each bin is obtained as $`\widehat{k}_i^2=(k_i+k_{i1})^2/4+\mathrm{\Delta }_i\mathrm{}^2/12`$ for a non-resonance bin and $`ϵ_i=ϵ_{\mathrm{res}}`$ for a resonance one. Comparing the approximate form (3) with the exact one (1) in the asymptotic region $`R\mathrm{}`$, it is natural to assume

$$S_{\mathrm{},L}^{(J)}(k)\frac{\widehat{S}_{\gamma ,\gamma _0}}{W_\gamma }f_i\mathrm{}(k)\text{(for Av) ,}$$
(10)
for $`k`$ belonging to the $`i`$-th bin, i.e., $`k_{i1}<kk_i`$.

### III.2 The pseudo-state method

In the PS method, $`H_{\mathrm{bc}}`$ is diagonalized in a space spanned by a finite number of $`L^2`$ type basis functions. The resultant eigenstates can well reproduce both bound and continuous states within a finite region of $`k`$ and $`r`$ YK1 . The $`k`$ continuum is automatically discretized by identifying the eigenstates of positive energies with $`\widehat{\mathrm{\Phi }}_i\mathrm{}(r)`$. The weight factor $`W_\gamma `$ is unity if the resultant discretized states $`\widehat{\mathrm{\Phi }}_i\mathrm{}(r)`$ are orthonormalized. Among the eigenstates, only low-lying states belonging to the region $`0<ϵ<\mathrm{}^2k_{\mathrm{max}}^2/2\mu _{\mathrm{bc}}`$ are taken as breakup channels in CDCC equation (5), where $`k_{\mathrm{max}}`$ is the maximum $`k`$ in the Av method.

The CDCC equations (5) thus obtained yield discrete breakup $`S`$-matrix elements. If the basis functions form a complete set with good accuracy in the region of $`r`$ and $`k`$ which is important for breakup processes, an accurate transformation from the approximate breakup $`S`$-matrix elements to the continuous (“exact”) ones is possible, as shown in Tostevin  for the Av method. The exact breakup $`T`$-matrix element is given by

$`T_\mathrm{}L^{(J)}(k)=\mathrm{\Phi }_{\mathrm{}}(k,r)j_L(PR)𝒴_{JM}^{\mathrm{},L}|U|\mathrm{\Psi }_{JM}.`$ (11)
Inserting the approximate complete set $`\{\widehat{\mathrm{\Phi }}_i\mathrm{}(r)\}`$ between the bra vector and the operator $`U`$ in Eq. (11), and replacing the ket vector by the CDCC wave function (3), one obtains the following approximate relation,

$`T_\mathrm{}L^{(J)}(k)`$ $``$ $`{\displaystyle \underset{i}{}}f_i\mathrm{}^{\mathrm{PS}}(k)\widehat{\mathrm{\Phi }}_i\mathrm{}(r)j_L(PR)𝒴_{JM}^{\mathrm{},L}|U|\mathrm{\Psi }_{JM}^{\mathrm{CDCC}}`$ (12)
$``$ $`{\displaystyle \underset{i}{}}f_i\mathrm{}^{\mathrm{PS}}(k)\widehat{T}_{\gamma ,\gamma _0},`$
where

$`f_i\mathrm{}^{\mathrm{PS}}(k)`$ $`=`$ $`\mathrm{\Phi }_{\mathrm{}}(k,r)|\widehat{\mathrm{\Phi }}_i\mathrm{}(r)`$ (13)
and

$`\widehat{T}_{\gamma ,\gamma _0}`$ $`=`$ $`\widehat{\mathrm{\Phi }}_i\mathrm{}j_L(\widehat{P}_iR)𝒴_{JM}^{\mathrm{},L}|U|\mathrm{\Psi }_{JM}^{\mathrm{CDCC}}.`$ (14)
The last form of Eq. (12) has been derived by replacing $`P`$ by $`\widehat{P}_i`$ in the spherical Bessel function $`j_L(PR)`$, which is valid since the $`k`$ distribution of $`f_i\mathrm{}^{\mathrm{PS}}(k)`$ is sharply localized at $`k=\widehat{k}_i`$. $`\widehat{T}_{\gamma ,\gamma _0}`$ is a CDCC breakup $`T`$-matrix element calculated with CDCC. Since the $`\widehat{T}_{\gamma ,\gamma _0}`$ are proportional to the corresponding $`S`$-matrix elements $`\widehat{S}_{\gamma ,\gamma _0}`$,

$$S_\mathrm{}L^{(J)}(k)\underset{i}{}f_i\mathrm{}^{\mathrm{PS}}(k)\widehat{S}_{\gamma ,\gamma _0}.$$
(15)
The “interpolation formula” (15) for the PS method agrees with the corresponding one in Tostevin  for the Av method. Thus, the interpolation formula (15) can be used for any method of discretization, if the discretized wave functions constitute an approximate complete set. The interpolation formula (15) is also independent of the type of the basis function taken, as obvious from the derivation, but it is necessary that the basis functions form a complete set in good approximation in the region of $`r`$ and $`k`$ which is important for the breakup $`T`$-matrix elements. As such basis functions, we here propose two types; one is the conventional real-range Gaussian functions

$$r^{\mathrm{}}\mathrm{exp}[(r/a_j)^2],$$
(16)
where $`a_j`$ ($`j=1`$$`n`$) are assumed to increase in a geometric progression K-G :

$$a_j=a_1(a_n/a_1)^{(j1)/(n1)}.$$
(17)
The other is an extension of Eq. (16H-Ka-Ki : pairs of functions

$$r^{\mathrm{}}\mathrm{exp}\left[\left(r/a_j\right)^2\right]\mathrm{cos}\left(b\left(r/a_j\right)^2\right),$$
(18a)
$$r^{\mathrm{}}\mathrm{exp}\left[\left(r/a_j\right)^2\right]\mathrm{sin}\left(b\left(r/a_j\right)^2\right),$$
(18b)
which can be also expressed as $`(\varphi _i\mathrm{}+\varphi _i\mathrm{}^{})/2`$ and $`(\varphi _i\mathrm{}^{}\varphi _i\mathrm{})/(2i)`$, respectively, with

$$\varphi _i\mathrm{}(r)=r^{\mathrm{}}\mathrm{exp}[\eta _jr^2],\eta _j=(1+ib)/a_j^2,$$
(19)
i.e., a Gaussian function with a complex range parameter. We refer to the new basis (18) as the complex-range Gaussian basis. In Eqs. (18) and (19) $`a_j`$ are the same as in Eq. (17); $`b`$ is a free parameter, in principle, but numerical test show that $`b=\pi /2`$ is recommendable. It should be noted that the total number of basis functions is $`2n`$. Note also that the complex-range Gaussian basis agrees with the real-range Gaussian basis when $`b=0`$.

The complex-range Gaussian basis functions are oscillating with $`r`$. They are therefore expected to simulate the oscillating pattern of the continuous breakup state wave functions better than the real-range Gaussian basis functions do. Actually, the latter reproduces the continuous state $`\mathrm{\Phi }_{\mathrm{}}(k,r)`$ in the region $`0kr20`$ YK1 , while the former does in the even larger region, $`0kr35`$.

If necessary, one can calculate, with the complex-range Gaussian basis as well as the real-range Gaussian basis, all the coupling potentials analytically by expanding the potential $`U`$ into a finite number of Gaussian functions. As far as three-body breakup reactions are concerned, however the analytic forms of the coupling potentials are not necessarily, since they can be easily obtained with numerical integration. Direct numerical calculation is difficult in the case of four-body breakup, and the analytic forms of the potentials are very useful in that case, as shown in section V.

## IV Numerical test of the pseudo-state method

In previous CDCC analyses Yahiro1 ; Piya , calculated elastic and breakup $`S`$-matrix elements were found to converge, for sufficiently large model space. In this section, we test the PS method by comparing the calculated $`S`$-matrix elements with those obtained with the Av method that converged within the error of 1% and hence forth called “exact” $`S`$-matrix elements. The test is made for two systems, $`d`$+58Ni scattering at 80 MeV and 6Li+40Ca scattering at 156 MeV.

### IV.1 $`d`$+58Ni scattering at 80 MeV

The model space taken in the present CDCC calculations is $`\mathrm{}=0,2`$ and $`k_{\mathrm{max}}=1.3\mathrm{fm}^1`$. It should be noted that the p-wave ($`\mathrm{}=1`$) breakup is negligible, because couplings $`V_{\gamma \gamma ^{^{}}}(R)`$ between odd and even parity breakup states contain a small factor, $`U=U_{\mathrm{pA}}(𝐫_{\mathrm{pA}})U_{\mathrm{nA}}(𝐫_{\mathrm{pA}})`$, in Eq. (6). The choice of the $`k_{\mathrm{max}}`$ is discussed below. Table I shows the parameters of the potentials used; the interaction between a nucleon and the target is the nucleon-nucleus optical potential of Becchetti and Greenlees BGP  at half the deuteron incident energy. The interaction between proton and neutron is a one-range Gaussian potential, $`v_{np}=v_0\mathrm{exp}[(r/r_0)^2]`$ with $`v_0=72.15`$ MeV and $`r_0=1.484`$ fm, which reproduces the radius and the binding energy of deuteron.

In the Av method the weight function is taken as $`f_i\mathrm{}(k)=1`$, since the projectile (deuteron) has no resonance state. The model space that gives convergence within error of 1% turns out to be $`\mathrm{\Delta }_i\mathrm{}=1.3/30`$ fm-1, as mentioned above, and $`k_{\mathrm{max}}=1.3`$ fm-1. The resulting values of $`k_{\mathrm{max}}`$ and $`\mathrm{\Delta }_i\mathrm{}`$ are different from those used in the previous analysis Piya ; the main purpose of Ref. Piya  was to show that the convergence of the CDCC solution was obtained within a model space of practical use and that the converged solution satisfied an appropriate boundary condition. The model space taken there, $`k_{\mathrm{max}}=1.0`$ fm-1 and $`\mathrm{\Delta }_i\mathrm{}=1/8`$ fm-1, are indeed enough for the elastic $`S`$-matrix elements and the dominant part of the breakup ones with the smaller $`k`$, therefore the elastic cross sections and the total breakup cross sections are well reproduced. However, the model space is found to be insufficient to obtain the “exact” $`S`$-matrix elements in the high $`k`$ region around $`1.0\mathrm{fm}^1`$, hence we take here $`k_{\mathrm{max}}=1.3`$ fm-1 and $`\mathrm{\Delta }_i\mathrm{}=1.3/30`$ fm-1 mentioned above.

In the real-range Gaussian PS method, a similar convergence is found, when the number of breakup channels, $`N_{\mathrm{PS}}`$, is 18 for both s- and d-waves. The number is even smaller when the complex-range Gaussian basis is taken: $`N_{\mathrm{PS}}`$ is 16 for s-wave and 17 for d-wave. The basis functions finally obtained have parameter sets $`(a_1=1.0,a_n=30.0,n=30)`$ for real-range Gaussian basis and $`(a_1=1.0,a_n=20.0,2n=40,b=\pi /2)`$ for complex-range one. For both of them $`N_{\mathrm{PS}}`$ is smaller than the number of basis functions. High-lying states with $`k>k_{\mathrm{max}}`$, which are obtained by diagonalizing $`H_{\mathrm{bc}}`$, do not affect the breakup $`S`$-matrix elements with $`k<k_{\mathrm{max}}`$, because the coupling potentials between the two $`k`$ regions are weak.

Figure 2 shows the discrete momenta $`\widehat{k}_i\mathrm{}`$ translated from the eigenenergies $`ϵ_i\mathrm{}`$ for the real- and complex-range Gaussian bases. One sees that for the real-range Gaussian basis, the discrete momenta are dense in the smaller $`k`$ region and sparse in the larger $`k`$ one. This distribution is not so effective in simulating the $`k`$ continuum, in the higher $`k`$ region in particular. For the complex-range Gaussian basis, on the other hand, the discrete momenta are distributed almost evenly. A similar sequence of the $`\widehat{k}_i`$ is also seen for the case of the transformed harmonic oscillator basis of Ref. PS1 ; Perez . Such a sequence of $`\widehat{k}_i`$ is close to that in the Av method. Thus, the complex-range Gaussian basis, as well as the transformed harmonic oscillator, is well suited for simulating the $`k`$ continuum in the entire region $`0<k<k_{\mathrm{max}}`$.

For the elastic $`S`$-matrix elements, both the real- and complex-range Gaussian PS methods well reproduce the “exact” one calculated with the Av method, as confirmed in Fig. 3 for the differential cross section. The three types of calculations, the real-range Gaussian PS, the complex-range Gaussian PS and the Av methods, yield an identical cross section at all scattering angles. Thus, both of the PS methods proposed here are useful for treating the breakup effects on the elastic scattering.

Figure 4 shows the result for breakup $`S`$-matrix elements at the grazing total angular momentum $`J=17`$, as a function of $`k`$. The real-range Gaussian PS method (dashed line) well simulates the exact solution calculated with the Av method (step line) in the lower $`k`$ region that corresponds to the main components of the breakup $`S`$-matrix elements, but inaccurate in the higher $`k`$ region around $`k=0.8`$ fm-1. The deviation at higher $`k`$ stems from the fact that real-range Gaussian basis poorly reproduces the continuum breakup state $`\mathrm{\Phi }_{\mathrm{}}(k,r)`$ at the higher $`k`$. Figure 4 shows that this problem can be solved by using the complex-range Gaussian basis (solid line) instead.

### IV.2 6Li+40Ca scattering at 156 MeV

Characteristic to this scattering, the projectile (6Li) has d-wave triplet resonance states ($`3^+,2^+,1^+`$). For simplicity, we neglect the intrinsic spin of 6Li, following Refs. CDCC-review1 ; sakuragi1 ; sakuragi2 . Then the projectile has only one d-wave resonance state with $`ϵ_{\mathrm{res}}=2.96`$ MeV and $`\mathrm{\Gamma }=0.62`$ MeV. Obviously the energy and the width do not reproduce experimental data, but at least the elastic cross section of 6Li is not affected much by the neglection of the spin Tanifuji .

In this scattering, the three-body system consists of deuteron, $`\alpha `$ and 40Ca. The interactions between each pair of the constituents are the optical potential of $`\alpha +^{40}`$Ca scattering at 104 MeV HAU72 , that of $`d+^{40}`$Ca scattering at 56 MeV dCa , and $`v_{\alpha \mathrm{d}}=v_0\mathrm{exp}\left[(r/r_0)^2\right]`$ with $`v_0=74.19`$ MeV and $`r_0=2.236`$ fm. Table II shows the parameters of the optical potentials.

The model space sufficient for describing breakup processes in this scattering is $`k_{\mathrm{max}}=2.0`$ fm-1 and $`\mathrm{}_{\mathrm{max}}=2`$; the model space is composed of two $`k`$-continua for $`\mathrm{}=0`$ and 2. Since there exists a resonance in $`\mathrm{}=2`$, the d-wave $`k`$-continuum is further divided in the Av method into the resonant part $`[0<k<0.55]`$ and the non-resonant part $`[0.55<k<2.0]`$. The $`k`$ continuum of $`\mathrm{\Phi }_{i,\mathrm{}=2}(k,r)`$ in the resonant part varies rapidly with $`k`$. The Av method can simulate the rapid change by taking $`f_{i,\mathrm{}=2}(k)=1`$ with bins of an extremely small width. In fact clear convergence is found for both the elastic and the breakup $`S`$-matrix elements, when the resonance part is described by 30 bins and the non-resonance part of the d-wave and the s-wave $`k`$-continua by 20 bins. Another Av discretization, which has been widely used as a convenient prescription CDCC-review1 ; sakuragi1 ; sakuragi2 ; sakuragi4 ; Sakuragi , is also made for comparison, in which the resonance region is represented by a single state with the weight factor of Breit-Wigner type (9). The two sorts of Av discretization are compared with the real- and complex-range Gaussian PS methods. With the PS methods, convergence of the $`S`$-matrix elements is found with 21 s-wave breakup channels and 22 d-wave ones. The level sequence of the resulting discrete eigenstates is shown in Fig. 5 for both the basis functions. The level sequences have the same properties as in Fig. 2. The parameter sets of the basis functions, finally taken in the PS methods, are $`(a_1=1.0,a_n=20.0,2n=40,b=\pi /2)`$ for the complex-range Gaussian basis and $`(a_1=1.0,a_n=30.0,n=30)`$ for the real-range one.

Figure 6 shows the differential cross section of the elastic scattering. The result with the precise Av discretization based on dense bins, considered to be the “exact” solution, is denoted by the solid line. The dotted line represents the result of the Watanabe model, i.e., with no breakup channels. The conventional Av discretization, based on the weight factor of Breit-Wigner type (dash-dotted line), well describes the breakup effect, particularly at very forward angles ($`\theta <20^{}`$), but deviates considerably from the “exact” solution at larger angles ($`\theta >30^{}`$). The complex-range Gaussian PS discretization (dashed line) well reproduces the “exact” solution with a number of channels being suitable for practical use. The real-range Gaussian PS method gives just the same result as the complex-range one.

Figure 7 represents breakup $`S`$-matrix elements at grazing total angular momentum $`J=43`$. The real- and complex-range Gaussian PS discretization well reproduce the “exact” solution calculated by the Av discretization with dense bins. The results of the two PS methods turn out to coincide within the thickness of the line. The resonance peak can be expressed by only 8 (12) breakup channels in the complex-range (real-range) Gaussian PS method, while the corresponding number of breakup channels is 30 in the Av method, as mentioned above. Thus, one can conclude that the real- and complex-range Gaussian PS methods are very useful for describing not only non-resonant states but also resonant ones.

## V Discussions on Four-body breakup reaction

In the past CDCC calculations the projectile was assumed to be a two-body system, dealing only with three-body breakup reactions. In this section, we investigate the applicability of CDCC to four-body breakup reactions of the projectile consisting of three particles, b+c+x (Fig. 8). The Av method needs the exact three-body wave functions being impossible to obtain. We can circumvent this problem with the present PS method; one can prepare an approximate complete set $`\{\widehat{\mathrm{\Phi }}_i\mathrm{}\}`$ by diagonalizing the Hamiltonian of the projectile in a space spanned by a set of basis functions of $`L^2`$ type. With $`\{\widehat{\mathrm{\Phi }}_i\mathrm{}\}`$ as the wave functions of the breakup channels, one can obtain an approximate total wave function $`\mathrm{\Psi }^{\mathrm{CDCC}}`$ by solving CDCC equations (5). Inserting $`\mathrm{\Psi }^{\mathrm{CDCC}}`$ into the exact form of breakup $`T`$-matrix elements in place of the exact total wave function, one reaches an approximate form:

$`T_4=e^{i(𝐏𝐑+𝐤𝐫+𝐪𝐲)}|U_4|\mathrm{\Psi }^{\mathrm{CDCC}}_{𝐑,𝐫,𝐲},`$ (20)
where $`U_4`$ is the sum of all interactions in the four-body system (A+b+c+x), $`𝐫`$ and $`𝐲`$ are and $`𝐤`$ and $`𝐪`$ are, respectively, two Jacobi coordinates of the three-body (b+c+x) system and momenta conjugate to them. The accuracy of Eq. (20) depends on how complete the set $`\widehat{\mathrm{\Phi }}_i\mathrm{}`$ is within the region $`(0rr_{\mathrm{max}},0yy_{\mathrm{max}})`$ which is important for the breakup process considered. An important advantage of the use of the real- and complex-range Gaussian bases is, as mentioned before, that analytic integrations over $`𝐫`$ and $`𝐲`$ can be done in Eq. (20), by expanding $`U_4`$ in terms of a finite number of Gaussian basis functions. This makes the derivation of $`T_4`$ feasible. Analyses based on this formulation are of much interest as a future work.

## VI Summary

The method of continuum discretized coupled channels (CDCC) is an accurate method of treating three-body breakup processes, in which the discretization of the $`k`$ continuum is essential. In this paper, we propose the new method of pseudo-state (PS) discretization which can be used not only for virtual breakup processes in elastic scattering but also for breakup reactions. First we show that an accurate transformation from the discrete breakup S-matrix elements calculated with the PS method to continuous ones is possible, since the PS basis functions can form in the good approximation a complete set in the finite region of $`𝐫`$ and $`𝐤`$ which is important for the breakup processes. As bases satisfying the completeness, we propose the real- and complex-range Gaussian bases. Both bases can treat virtual breakup processes in the elastic scattering with high accuracy, i.e., with the error of calculated cross sections less than 1%. For breakup processes, the complex-range Gaussian basis is accurate throughout the entire region of the $`k`$ continuum. The real-range Gaussian basis also keeps a good accuracy for the dominant part of breakup $`S`$-matrix elements with the lower $`k`$, although it is partially inaccurate for the higher $`k`$ region. Thus, both bases can be used for realistic analyses of elastic scattering and projectile breakup reactions.

The present new PS method has at least two advantages over the widely used momentum bin average method. One is that it does not need the exact wave function of the projectile over the entire region of $`r`$. This is important from a theoretical point of view. The other advantage is that with the real- and complex-range Gaussian bases one can calculate all the coupling potentials semi-analytically, which is very useful in actual calculations. Furthermore, if the projectile has resonances in its excitation spectrum, the new method discretizes the complicated spectrum with a reasonable number of the basis functions, without distinguishing the resonance states from non-resonant continuous states. These advantages of the new method are extremely helpful, sometimes even essential, in applying CDCC to four-body breakup effects of unstable nuclei such as 6He and 11Li. Actually, a CDCC analysis of four-body breakup effect on the 6He elastic scattering is in progress, and the result of the analysis will appear in a forthcoming paper.

## Acknowledgments

The authors would like to thank M. Kawai for helpful discussions. This work has been supported in part by the Grants-in-Aid for Scientific Research (12047233, 14540271) of Monbukagakusyou of Japan.


text:
# Stellar abundance analyses in the light of 3D hydrodynamical model atmospheres

## 1. Introduction

Determining stellar element abundances play an integral role in most endeavours to improve our understanding of stellar, galactic and cosmic evolution. The term observed abundances is somewhat of a misnomer however, since the chemical composition of course can not be inferred directly from an observed spectrum. The obtained abundances are therefore never more trustworthy than the models employed to analyse the observations. The derivation of accurate abundances require realistic models of both the spectrum formation region and the spectrum formation process, which is the topic of the present contribution.

For a long time, understanding stellar convection has been one of the greatest challenges in stellar astrophysics. Convection can provide an efficient means of transporting energy and thereby greatly influence the stellar structure and evolution, while also mixing nuclear-processed material and preventing radiative diffusion and gravitational settling. For late-type stars, the surface convection zone reaches the stellar atmosphere, which directly affects the emergent spectrum. The solar granulation is the observational manifestation of convection: concentrated, rapid, cold downdrafts in the midst of broad, slow, warm upflows. Qualitatively similar granulation properties are expected in other solar-type stars, as indeed confirmed by 3D numerical simulations (e.g. Nordlund & Dravins 1990; Asplund et al. 1999; Asplund & García Pérez 2001; Allende Prieto et al. 2002) and indicated by observed line asymmetries. The up- and downflows have radically different temperature structures (Stein & Nordlund 1998), which can not be approximated by normal theoretical 1D hydrostatic model atmospheres with different effective temperatures $`T_{\mathrm{eff}}`$ (Fig. 1). In fact, neither the steep temperature structures of the ascending material nor the shallow gradients of the descending gas are well described by such 1D models. Because of the photospheric inhomogeneities and the highly non-linear and non-local nature of spectrum formation, it is clear that no single 1D model can be expected to properly describe all aspects of what is inherently a 3D phenomenon.

Here I will describe recent progress in developing 3D hydrodynamical model atmospheres of late-type stars and some applications to stellar abundance analyses. In spite of their still inferior treatment of radiative transfer and line blanketing compared with existing 1D model atmospheres based on the mixing length theory, detailed comparisons with a range of observational diagnostics suggest that these simulations are indeed highly realistic.

## 2. 3D Hydrodynamical Model Atmospheres

The 3D model atmospheres which form the basis of the abundance analyses presented here have been computed with a 3D, time-dependent, compressible, explicit, radiative-hydrodynamics code developed to study solar and stellar surface convection (Stein & Nordlund 1998). The hydrodynamical equations for conservation of mass, momentum and energy are solved on a non-staggered Eulerian mesh with gridsizes of $`100^3`$. The physical dimensions of the grids are sufficiently large to cover many ($`>10`$) granules simultaneously in the horizontal direction and about 13 pressure scale-heights in the vertical. In terms of continuum optical depth the simulations extend at least up to log $`\tau _{\mathrm{Ross}}5`$ which for most purposes are sufficient to avoid numerical artifacts of the open upper boundary on spectral line formation. The lower boundary is located at large depths to ensure that the inflowing gas is isentropic and featureless, while periodic horizontal boundary conditions are employed. The temporal evolution of the simulations cover several convective turn-over time-scales to allow thermal relaxation to be established and to obtain statistically significant average atmospheric structures and spectral line profiles. No magnetic fields have been included in the simulations described here although the code is capable of simulating magneto-convection.

In order to obtain a realistic atmospheric structure, it is crucial to have the best possible input physics, and properly account for the energy exchange between the radiation field and the gas. The adopted equation-of-state is that of Mihalas et al. (1988), which includes the effects of ionization, excitation and dissociation. The continuous opacities come from the Uppsala package (Gustafsson et al. 1975 and subsequent updates) while the line opacities are from Kurucz (1998, private communication). The 3D radiative transfer is solved at each time-step using 9 inclined rays (2 $`\mu `$-angles and 4 $`\phi `$-angles, plus the vertical) under the assumptions of local thermodynamic equilibrium (LTE, $`S_\lambda =B_\lambda `$) and opacity binning (Nordlund 1982). The four opacity bins are designed to correspond to continuum and weak, intermediate and strong lines. The assignment of the original 2748 wavelength points into the different opacity bins follows from detailed monochromatic radiative transfer calculations of the 1D averaged atmospheric structure. The opacity binning thus includes the effects of line-blanketing in a manner reminiscent of opacity distribution functions. In order to improve the numerical accuracy, the radiative transfer is solved on a finer grid and the results subsequently interpolated back to the original grid for the hydrodynamical variables at each timestep.

It is important to realise that the simulations contain no free parameters which are tuned to improve the agreement with observations. The adoption of the numerical and physical dimensions of the simulation box is determined by practical computional time considerations, the need to resolve the most important spatial scales and the wish to place the artificial boundaries as far as possible from the region of interest. The code is stabilized using a hyper-viscosity diffusion algorithm with parameters specified from standard hydrodynamical test cases. These parameters are not changed in the actual stellar convection simulation and are therefore not freely adjustable parameters. It has been verified that the resulting atmospheric structures are insensitive to the adopted effective viscosity at the current highest affordable numerical resolution (Asplund et al. 2000a). The input parameters discriminating different models are the surface gravity log $`g`$, metallicity \[Fe/H\] and the entropy of the inflowing material at the bottom boundary. The effective temperature of the simulation is therefore a property which depends on the entropy structure and evolves with time around its mean value following changes in the granulation pattern.

Further details of the 3D hydrodynamical model atmospheres are available in Stein & Nordlund (1998), Asplund et al. (1999, 2000a,b) and Asplund & García Pérez (2001).

## 3. Spectral Line Formation

### 3.1. General Description

The 3D hydrodynamical model atmospheres described in the previous section form the basis for the 3D spectral line formation calculations. The Doppler shifts introduced by the convective and oscillatory motions cause line asymmetries and shifts as well as broaden the lines significantly. Spatially resolved profiles thus come in an amazing assortment as a result of the atmospheric inhomogeneities and velocities together with their correlations (Fig. 2). The concepts of micro- and macroturbulence, which are introduced in 1D analyses in order to account for the missing line broadening, are not necessary in 3D calculations, as evident from the excellent agreement between predicted and observed line profiles (Fig. 3, see also Asplund et al. 2000b) and the lack of significant trends in derived abundances with line strength (Asplund et al. 2000c; Asplund 2000). The main component of microturbulence does not arise from microscopic turbulent motions but from gradients in convective motions that are resolved with the present simulations, which is also the main physical explanation for macroturbulence.

The original 3D hydrodynamical simulations which includes deep, extremely optically thick layers are interpolated to a finer depth scale which extends only down to log $`\tau _{\mathrm{Ross}}2.5`$ prior to the spectral line calculations for improved numerical accuracy. At the same time, the simulations are interpolated to a coarser horizontal grid to ease the computational burden. Extensive tests have ensured that this practice is acceptable even when considering fine details such as line asymmetries; for abundance analysis purposes the procedure introduces non-noticable differences ($`<0.01`$ dex for individual snapshots and much less for temporal averages). Depending on whether the 3D line formation is computed in LTE or in non-LTE, different assumptions and approximations are necessary, which are described next.

### 3.2. 3D LTE Line Formation

With the assumption of LTE ($`S_\lambda =B_\lambda `$), the level population is uniquely determined by the local gas temperature from the Boltzmann and Saha distributions. With the source function known, it is then straightforward to solve the radiative transfer equation, which typically is performed for about 100 wavelength points across the line profile for 10-20 inclined outgoing rays for more than 50 different snapshots. This allows asymmetries and shifts to be determined to an accuracy on the m s<sup>-1</sup>-level. For abundance determinations even individual snapshots give results accurate to within about 0.03 dex, since at any instant many granules are covered by the simulation box. All in all, a single temporally and spatially averaged flux profile in 3D correspond most of the time to $`N_\mathrm{t}N_\mathrm{x}N_\mathrm{y}N_{\mathrm{angles}}N_\lambda >10^8`$ 1D radiative transfer calculations. In addition, each 3D profile is normally computed for at least three different abundances to enable interpolation to the requested line strength. Even then, such 3D LTE line calculations are achievable on current workstations thanks to efficient numerical algorithms and implementations.

No micro- and macroturbulence enter the calculations but for flux profiles the projected stellar rotational velocity $`v_{\mathrm{rot}}\mathrm{sin}i`$ causes line broadening in addition to the natural, pressure and convective broadening. Flux profiles are obtained after the actual radiative transfer calculations from a disk-integration of the stored profiles of the individual inclined rays for a given $`v_{\mathrm{rot}}\mathrm{sin}i`$ (Dravins & Nordlund 1990). The detailed comparison with observations can then proceed through a $`\chi ^2`$-analysis by interpolating the computed profiles with varying abundances and adopting different $`v_{\mathrm{rot}}\mathrm{sin}i`$.

### 3.3. 3D Non-LTE Line Formation

Recently, even detailed 3D non-LTE line formation calculations have become feasible (Kiselman 1997; Uitenbroek 1998; Botnen & Carlsson 1999; Asplund et al. 2002). Without the much simplifying assumptions of Saha and Boltzmann distribution for the level population, the rate equations assuming statistical equilibrium ($`\mathrm{d}n_\mathrm{i}/\mathrm{d}t=0`$) must be solved simultaneously with the radiative transfer equation. Due to the much increased computational effort in achieving this, sacrifies must be made in other areas, most notably in terms of number of snapshots the calculations can be made for. With individual snapshots it is still feasible to carry out calculations for atoms consisting of hundreds of radiative transitions, each with 50-100 wavelength points. Fortunately, with 3D model atmospheres the non-LTE effects are expected to be quite similar between different snapshots due the relatively large area coverage always containing many up- and downflows, which is also confirmed by detailed calculations. The by far dominating part of 3D non-LTE computations of reasonably sized model atoms is the actual radiative transfer solution from the current estimate of the source function. This formal solution of the radiative transfer equation is performed for all wavelengths of all transitions for typically 24-48 inclined rays to compute mean intensities and flux profiles. The final converged solution is then obtained by a iterative procedure, normally employing an accelerated lambda iteration technique or similar, allowing convergence typically within 10 iterations. For small model atoms, the calculations are possible to perform on modern workstations but for more extended atoms the use of supercomputers are better suited, not the least due to vast memory requirements with 3D non-LTE codes such as multi3d (Botnen & Carlsson 1999; Asplund et al. 2002).

## 4. Stellar Abundance Analyses with 3D Model Atmospheres

The realistic 3D hydrodynamical model atmospheres of late-type stars together with corresponding LTE and non-LTE line formation tools described above have a natural application to stellar abundance analyses. For abundance studies, the different mean photospheric structure between 3D and 1D models, the existence of atmospheric inhomogenities and the broadening of spectral lines due to convective motions are all important ingredients for the line formation. To illustrate the possible impact of the new generation of 3D model atmospheres on abundance analyses, I will describe the situation for the most abundant cosmic metal, oxygen, in two important types of stars: the Sun and metal-poor halo stars. In both cases it seems like the adoption of 3D model atmospheres drastically affects the conclusions.

### 4.1. Oxygen in the Sun

The solar oxygen abundance ought to be well determined but astonishingly the commonly adopted value may be in error by almost a factor of two. For the Sun, five different types of diagnostics are potentially available: permitted high-excitation O i lines, the forbidden \[O i\] lines at 630 nm, the OH vibrational-rotational and pure rotational lines in the IR, and finally the OH electronic transitions in the UV (which have not been included in the present analysis). All have quite different sensitivities to the atmospheric structure and the line broadening. In addition, all have their pros and cons: the \[O i\] lines are very weak and likely blended (Allende Prieto et al. 2001); the O i lines are sensitive to departures from LTE (Kiselman 1993); the OH IR lines are, like all molecular features, very temperature sensitive (Asplund & García Pérez 2001); the OH UV lines may be affected by missing opacity (Balachandran & Bell 1998) and possibly non-LTE effects (Asplund & García Pérez 2001). Unfortunately the different lines imply wildly discortant abundances, the difference amounting to about 0.3 dex both for the Holweger-Müller (1974) and marcs (Asplund et al. 1997) 1D model atmospheres. Normally the conclusions from the OH IR have been embraced (log $`ϵ_\mathrm{O}8.90`$), largely due to the remarkably small line-to-line scatter (Sauval et al. 1984). An application of our new 3D model atmospheres to the problem, however, suggests that the solar oxygen abundance may be significantly lower (log $`ϵ_\mathrm{O}8.65`$) than previously thought (Allende Prieto et al. 2001; Asplund et al., in preparation).

\[O i\]: Due to the exceptionally good agreement between predicted and observed line profiles in general with 3D model atmospheres (Fig. 3), disturbing blends are readily identified. This statement naturally presumes that detailed comparisons have already been performed for a range of lines with similar properties and typical line formation regions, which fortunately is the case for \[O i\]. The very poor agreement between the theoretical and the observed profile, immediately signals that the \[O i\] 630.0 nm line is indeed blended, which have occasionally been suspected but mostly thought to be relatively unimportant. The effect of the blending Ni i line reduces the derived abundance significantly ($`0.13`$ dex), while the transition from 1D to 3D model atmospheres and a revised $`gf`$-value are of less importance ($`0.08`$ and $`0.03`$ dex, respectively) but compounds the overall difference. The forbidden line implies a solar oxygen abundance of log $`ϵ_\mathrm{O}=8.69\pm 0.05`$ (Allende Prieto et al. 2001).

O i: The O i triplet around 777 nm is well-known to be susceptible to significant departures from LTE. We have therefore investigated its line formation in detail with the 3D non-LTE code multi3d (Botnen & Carlsson 1998; Asplund et al. 2002), employing two independent snapshots from the 3D solar simulation and a model atom consisting of 23 levels and 65 radiative transitions (Kiselman 1993). The adopted model atom is sufficiently comprehensive for this particular problem since even a two-level approach catches the essentials of the line formation. Collisional excitation and ionization due to electrons are accounted for but not the corresponding case of hydrogen collisions, since the available classical recipes (Drawin 1968) likely over-estimates their importance by several orders of magnitude. Our 3D non-LTE calculations indeed reveal pronounced departures from LTE, in terms of abundance about 0.2 dex. The dominant non-LTE mechanism is the same as with 1D models, namely photon losses in the line itself. Preliminary the triplet implies a solar oxygen abundance of log $`ϵ_\mathrm{O}=8.67\pm 0.03`$ (Asplund et al., in preparation).

OH vibrational and rotational lines: The slightly cooler temperature structure in our 3D solar model compared with the standard Holweger-Müller (1974), which is normally the 1D model of choice for solar abundance determinations, directly translates to a lower derived oxygen abundance from the OH IR lines. Equally important, however, is the presence of temperature inhomogeneities due to the non-linear line formation. It therefore comes as no surprise that 1D analyses typically will significantly over-estimate the abundances derived from molecular lines. We have reanalysed 70 and 127 OH vibration-rotation and pure rotation lines, respectively, which have previously been studied by Grevesse et al. (1984) and Sauval et al. (1984) using up-to-date molecular line data. The 3D line formation calculations are based on the assumption of LTE, which is likely a reasonable approximation in view of the efficient collisional thermalization of the relevant levels. A potential source of error comes from the assumed chemical equilibrium in the OH molecule formation, which would tend to artificially decrease the oxygen abundance slightly. The magnitude of this effect is not known at this stage but will be investigated further in the near future. Keeping this in mind, preliminary unweighted means yield log $`ϵ_\mathrm{O}=8.61\pm 0.03`$ for the vibrational lines and log $`ϵ_\mathrm{O}=8.67\pm 0.03`$ for the rotational lines.

Summary: Although the exact value for the solar oxygen abundance depends on the relative weights the different diagnostics are given, it is clear that our detailed 3D investigation unanimously points to a significantly lower solar oxygen abundance than previously thought based on classical 1D model atmospheres. Preliminary we estimate log $`ϵ_\mathrm{O}=8.65\pm 0.05`$ which is lower by almost a factor of two. Such a radical re-evaluation for the most abundant metal in the star we know the best is certainly a sobering realisation when discussing derived stellar abundances based on inferior spectra and where the modelling uncertainties can be much larger. Our revision directly impacts also on stellar evolution calculations and helioseismology among other topics.

### 4.2. Oxygen in Metal-Poor Stars

The oxygen abundances in stars born at different epochs and locations carry vital information about the formation and evolution of galaxies. In particular the temporal (with metallicity as a practical although not perfect proxy) evolution of the oxygen-to-iron ratio can yield estimates of the time-scales involved (O mainly from SNe II while Fe largely from SNe Ia) as well as give important clues to the nature of the elusive Pop III stars. It has long been known that halo stars show an over-abundance of oxygen (\[O/Fe\] $`>0`$) but an exact estimate has been notoriously difficult to achieve due to the large discrepancies between the different abundance indicators. Typically the \[O i\] line implies a more or less plateau at \[O/Fe\] $`+0.4`$ (Barbuy 1988; Nissen et al. 2002) but due to weakness of the line its use for \[Fe/H\] $`<2.0`$ is restricted to analysis of giants and subgiants. Recent analyses of the OH UV lines in halo dwarfs on the other hand seem to suggest an almost linear trend with metallicity, reaching \[O/Fe\]$`+1.0`$ at \[Fe/H\]$`=3.0`$ (Israelian et al. 1998; Boesgaard et al. 1999) while the OH IR vibrational lines in giants again more support the conventional plateau (Melendez & Barbuy 2002). The confusing situation is completed by the fact that the oxygen abundances derived from the high-excitation O i triplet at 777 nm tend to lay in between the OH UV and the \[O i\] results.

We have recently applied the new generation of 3D hydrodynamical model atmospheres of metal-poor stars (Asplund et al. 1999) in an attempt to investigate whether a possible resolution to the oxygen conundrum can be found. Due to the lack of radiative heating from spectral lines in halo stars, the adiabatic cooling from the rising gas in the upflows dominate the energy balance. Consequently, the temperature is much below the radiative equilibrium value in the optically thin layers, which is enforced in 1D models.

OH UV: As noted before, molecule formation is extremely temperature sensitive. The lower temperatures in 3D models therefore directly translates to much larger number densities of molecules and thus to stronger lines. The application of 3D model atmospheres appears to bring down the high advocated linear trend to something which better resembles the canonical plateau (Asplund & García Pérez 2001), revealing a potential and very significant systematic error in previous standard analyses using 1D model atmospheres. This 3D LTE correction can easily amount to $`0.6`$ dex or more at \[Fe/H\] $`=3.0`$. The final word on this topic is not written yet however, as the 3D OH analyses sofar make the crucial assumption of LTE both in the line formation and in the OH molecule formation. Existing non-LTE calculations for molecules are very rare. Asplund & García Pérez (2001) found the possibility of significant departures from LTE in the OH line formation but with little metallicity dependence using a very simplified two-level approach in 1D; the 3D case still remains to be investigated. In addition, due to the rapid cooling of the upflowing material it is not unlikely that the time available is not sufficient for equilibrium chemistry to be established, in particular for metal-poor stars with their very low temperatures in the upper layers.

\[O i\]: One could imagine that the \[O i\] line should be insensitive to the low surface temperatures in 3D models of halo stars insofar it originates from the ground level of the ionization stage where essentially all oxygen is residing. In reality though there are significant downward 3D corrections also for this line, since the line strength depends on the ratio of line and continuous opacity. The latter is dominated by H<sup>-</sup> whose number density is proportional to the free electrons which are much reduced in the low temperature gas (Nissen et al. 2002). In terms of abundance the 3D corrections amount to about $`0.2`$ dex at \[Fe/H\] $`2.0`$ for dwarfs and subgiants. Whether a similar revision is needed for giants have to await the construction of 3D models for such stars but slightly smaller corrections would not be unexpected.

O i: The most problematic of the available diagnostics in metal-poor stars is likely the O i lines. Due to their high excitation potential, the lines are formed in very deep layers where the differences between 3D and 1D models are small. Already in 1D analyses, the triplet yields higher O abundances than \[O i\] which will be aggrevated unless there are large metallicity-dependent 3D non-LTE effects on the triplet or the adopted T<sub>eff</sub>-scale needs significant revision. Preliminary 3D non-LTE calculations for the halo subgiant HD 140283 (Asplund et al., in preparation) suggest quite similar abundance corrections to the solar case, $`0.2`$ dex, which is apparently not enough (Nissen et al. 2002). The best T<sub>eff</sub>-scale for halo stars likely comes from the IR flux method, which gives too low T<sub>eff</sub> to resolve the discrepancy even with 3D models (Asplund & García Pérez 2001). The use of 3D models may have resolved the problem between the \[O i\] and the OH UV lines but clearly not all answers are in for the O i triplet.

Summary: It seems like the adoption of the new generation of 3D hydrodynamical model atmospheres holds a crucial key to the resolution of the long-standing problem of oxygen abundances in halo stars. Although the final answer has not yet emerged, the existing results together with the recent downward revision of the solar O abundance described above point to a nearly flat \[O/Fe\] trend but with a slow rise for \[Fe/H\] $`<2`$, i.e. somewhere in between the hotly contested plateau and linear trends (Nissen et al. 2002).

## 5. Concluding Remarks

The above examples illustrate that existing 1D analyses may well be hampered by significant systematic errors which can severely impact the inferred conclusions. With the advent of a new generation of telescopes and instrumentation it is clear that the dominant source of error is no longer of observational nature but rather stems from uncertainties inherent to the actual analyses. There is obviously a great need for investing also significant resources in improving the modelling aspects in order to place the findings on firmer footing. It is noteworthy that the main obstacle for achieving this is simply a lack of manpower, in spite of the existence of a very large community of astronomers relying on such model atmospheres for their applications.

#### Acknowledgments.

The work described in this contribution represent the results of successful collaborations with a large number of wonderful colleagues, in particular Carlos Allende Prieto, Mats Carlsson, Remo Collet, Ana García Pérez, Nicolas Grevesse, Dan Kiselman, David Lambert, Poul Erik Nissen, Åke Nordlund, Francesca Primas, Jacques Sauval, Bob Stein, and Regner Trampedach. Finally, I am grateful to the editors for their patience.

## References

Allende Prieto, C., Lambert, D.L., Asplund, M. 2001, ApJ, 556, L65

Allende Prieto, C., Asplund, M., García López, R.J., Lambert, D.L. 2002, ApJ, 567, 544

Asplund, M. 2000, A&A, 359, 755

Asplund, M., García Pérez, A.E. 2001, A&A, 372, 601

Asplund, M., Carlsson, M., Botnen, A.V. 2002, submitted to A&A

Asplund, M., Gustafsson, B., Kiselman, D., Eriksson, K. 1997, A&A, 318, 521

Asplund, M., Nordlund, Å., Trampedach, R., Stein, R.F. 1999, A&A, 346, L17

Asplund, M., Ludwig, H.-G., Nordlund, Å., Stein, R.F. 2000a, A&A, 359, 669

Asplund, M., Nordlund, Å., Trampedach, R., Allende Prieto, C., Stein, R.F. 2000b, A&A, 359, 729

Asplund, M., Nordlund, Å., Trampedach, R., Stein, R.F. 2000c, A&A, 359, 743

Balachandran, S.C., Bell, R.A. 1998, Nature, 392, 23

Barbuy, B. 1988, A&A, 191, 121

Boesgaard, A.M., King, J.R., Deliyannis, C.P., Vogt, S. 1999 AJ, 117, 492

Botnen, A.V., & Carlsson, M. 1999, in Numerical astrophysics, 379

Dravins, D., Nordlund, Å. 1990, A&A, 228, 203

Drawin, H.W. 1968, Z. Physik, 211, 404

Grevesse, N., Sauval, A.J., van Dishoek, E.F. 1984, A&A, 141, 10

Gustafsson, B., Bell, R.A., Eriksson, K., Nordlund, Å. 1975, A&A, 42, 407

Holweger, H., Müller, E.A. 1974, Solar Physics, 39, 19

Israelian, G., García López, R.J., Rebolo, R. 1998, ApJ, 507, 805

Kiselman, D. 1993, A&A, 275, 269

Kiselman, D. 1997, ApJ, 489, L107

Meléndez, J., Barbuy, B. 2002, ApJ, 575, 474

Mihalas, D., Däppen, W., Hummer, D.G. 1988, ApJ, 331, 815

Nissen, P.E., Primas, F., Asplund, M., Lambert, D.L. 2002, A&A, 390, 235

Nordlund, Å. 1982, A&A, 107, 1

Nordlund, Å., Dravins, D. 1990, A&A, 228, 155

Sauval, A.J., Grevesse, N., Brault, J.W., Stokes, G.M., Zander, R. 1984, ApJ, 282, 330

Stein, R.F., Nordlund, Å. 1998, ApJ, 499, 914

Uitenbroek, H. 1998, ApJ, 498, 427