text:
# Scattering phases in quantum dots: an analysis based on lattice models

## I Introduction

Phase coherence is at the heart of most phenomena studied in mesoscopic physics. However, the behavior of the electronic wave-function phase itself in an actual quantum transport device had not been studied until recent years. In the case of quantum dots, investigations have been predominantly restricted to conductance measurements  which carry no information on the transmission phase. It was not until the experiments by Yacoby et al.  that the interest in phase behavior of quantum dots really started. In these experiments a quantum dot was embedded in one of the arms of an Aharonov-Bohm (AB) ring in an attempt to analyze the transmission phase evolution as a function of the dot gate voltage. Although these experiments were the first to demonstrate the presence of a coherent component in the current through a quantum dot in the Coulomb blockade regime, they failed to give the complete evolution of the phase. This limitation was explained  as a consequence of the phase locking (or phase rigidity) that occurs in a two terminal geometry. It was shown  that the AB effect in such geometry is characterized by a parity: as a function of the AB flux the conductance exhibits either a local maximum at zero flux (positive parity) or a local minimum (negative parity). Another intriguing feature of the experimental results was the “parity conservation” over a large sequence of Coulomb blockade peaks, which reflected a similar evolution of the phase over each peak. The complete evolution of the phase could be obtained in a subsequent experiment by Schuster et al.  using a four terminal geometry. This experiment confirmed the expected evolution of the phase around the peaks and revealed that an abrupt jump of $`\pi `$ occurs in the valleys between the peaks.

Since 1995 several theoretical efforts have been devoted to explain these observations : Ref.  proposed a screening effect, Ref.  alluded to dot degeneracies, Refs.  associated the observed effect to an asymmetric deformation of the dot which leads to repeated charging of the same dot level and Ref.  pointed out some special properties of the dot states in a semi-chaotic situation (this mechanism will be further analyzed in the present work). In spite of all these efforts there is still the feeling that a more fundamental explanation is lacking. Each one of the proposed mechanisms can be criticized as relying on some particular assumptions. Only the approximate sum rule proposed recently in Ref.  is supposed to be valid in a generic (chaotic) situation. However, there is still no experimental evidence of the near resonance phase-lapse predicted by this mechanism.

On the other hand, the phase problem affects our knowledge on generic properties of scattering phases. One fundamental relation, invoked in Ref.  in connection to this problem, is the Friedel sum rule which relates the phase of the eigenvalues of the scattering matrix to the charge accumulated in the dot region. Being related to the dot charge, the Friedel phase is a continuous function of the system parameters and cannot exhibit an abrupt behavior like the one found in the experiments by Schuster et al. . However, as pointed out recently by Lee and Taniguchi and one of the authors the phase of the transmission amplitude can depart from the Friedel phase and exhibit a non-analytic behavior at the points where the modulus of the transmission vanishes. It is thus intersting to study the general conditions for the occurrence of zeroes in the transmission through a quantum dot.

The aim of this paper is to investigate the behavior of the different scattering phases in quantum dots with the help of lattice models. These type of models allow to describe a dot of arbitrary shape and to study the influence of disorder . We shall first derive the expressions for the different scattering phases in terms of Green functions. These expressions allow to relate the occurrence of zeroes in the transmission with the parity of the isolated dot wave-functions. We also study the statistical properties of the dot wave-functions in a disordered quantum dot, showing explicitly the absence of significant correlations in the chaotic case. On the other hand, as suggested in Ref. , for weak disorder one can identify certain dot levels which are much more strongly coupled to the leads than average. We shall show that in this situation the dot levels are populated in an anomalous way as a function of gate voltage. We shall finally discuss the possible role of this type of correlation effects in the phase problem.

The outline of the paper will be the following: In section II we introduce a generic lattice model for a Quantum Dot coupled to single-moded leads. In section III we derive the expressions for the different scattering phases in terms of Green functions. We discuss in particular the Friedel sum rule and the relation of the Friedel phase to the phase of the transmission amplitude. In section IV we study the conditions for the occurrence of zeroes in the transmission and show that they are independent of the strength of the coupling to the leads. The statistical properties of the isolated dot wave-functions are analyzed in section V. Finally, in section VI we discuss the role of electron correlation effects. We end the paper with some conclusions and final remarks.

## II Generic lattice model

As a model for a two dimensional quantum dot we consider a collection of $`N`$ sites on a square lattice (see Fig. 1). This model can represent a dot of arbitrary shape. The electrons in the dot are described by a tight-binding Hamiltonian with site energies $`ϵ_i`$ and a constant hopping element $`t`$ coupling nearest-neighbors only. ($`t`$ will be taken as the unit of energy). The site energies are allowed to vary following an imposed electrostatic confining potential and/or randomly in order to study the influence of disorder.

On the other hand, electron-electron interactions can be included within the constant charging energy model by adding a term $`V_{coul}=E_C(N_{dot}CV_g/e)^2`$, where $`N_{dot}`$ is the mean number of electrons in the dot, to the one-electron Hamiltonian. Its effect will be discussed in section VI.

We shall consider that the dot is coupled to electron reservoirs by two one-dimensional leads as depicted in Fig. 1. We may assume that the coupling to the left and right leads is restricted to two sites labeled by 1 and $`N`$ respectively. The first site on each lead are connected to these two sites on the dot by hopping elements $`t_L`$ and $`t_R`$ respectively. As will be discussed later, this situation can be easily generalized to the case where the first sites on the leads are connected to several sites on the dot (this multiple connection is illustrated in Fig. 1 by the dashed lines).

## III Scattering phases and Green functions

The electronic properties of lattice models are conveniently given in terms of Green functions. We need to introduce the retarded and advanced Green operators given by

$$\widehat{G}^{r,a}(\omega )=\left[\omega \widehat{H}_{dot}\widehat{\mathrm{\Sigma }}^{r,a}(\omega )\right]^1,$$
(1)
where $`\widehat{H}_{dot}`$ is the one-electron part of the isolated dot Hamiltonian and $`\widehat{\mathrm{\Sigma }}^{r,a}`$ is the (retarded, advanced) self-energy operator describing the coupling of the dot to the leads. In a Fermi-liquid like description, the self-energy should also contain terms accounting for electron-electron interactions . We shall postpone its discussion to section VI and concentrate on the one-electron properties of the quantum dot until then. When the coupling to the leads is localized at sites 1 and $`N`$ we have

$$\mathrm{\Sigma }_{l,m}^{r,a}(\omega )=\delta _{l,j_\alpha }\delta _{m,j_\alpha }t_\alpha ^2g_\alpha ^{r,a}(\omega ),$$
(2)
where $`\alpha =L,R`$, $`j_L=1,j_R=N`$ and $`g_\alpha ^{r,a}`$ are the local Green functions at the semi-infinite one-dimensional leads.

In terms of the Green operator one can express the dot total density of states $`\rho (\omega )`$ as

$$\rho (\omega )=\frac{1}{2\pi i}\text{Tr}\left[\widehat{G}^a(\omega )\widehat{G}^r(\omega )\right]=\frac{1}{\pi }\text{Im Tr}\left[\widehat{G}^a(\omega )\right].$$
(3)
We shall now discuss the connection of this quantity to the scattering matrix $`\widehat{S}(\omega )`$, defined in terms of Green functions by means of the generalized Fisher-Lee relations

$$S_{\alpha ,\beta }(\omega )=\delta _{\alpha ,\beta }2i\sqrt{\mathrm{\Gamma }_\alpha \mathrm{\Gamma }_\beta }G_{j_\alpha ,j_\beta }^r(\omega ),$$
(4)
where we have introduced the tunneling rates $`\mathrm{\Gamma }_\alpha =t_\alpha ^2\text{Im}g_\alpha ^a`$. The unitarity of $`\widehat{S}`$ can be readily shown (see appendix A). We start by defining the quantity $`\theta _F`$ as

$$\theta _F(\omega )=\text{Im}\text{ln}\text{Det}\left[\omega \widehat{H}_{dot}\widehat{\mathrm{\Sigma }}^a(\omega )\right].$$
(5)
The derivative with respect to the energy of $`\theta _F`$ is then given by

$$\frac{\theta _F}{\omega }=\text{Im}\text{Tr}\left[\widehat{G}^a\left(1\frac{\widehat{\mathrm{\Sigma }}^a}{\omega }\right)\right].$$
(6)
Thus, in the case when the energy dependence of the self-energy can be neglected one obtains the identity

$$\frac{\theta _F}{\omega }=\pi \rho (\omega ).$$
(7)
This case corresponds to a particular type of leads, having a large density of states, which can screen any deviation from charge neutrality induced by the presence of the dot. This type of leads can be called “non-polarizable leads”.

On the other hand, it can be show (see appendix B) that $`\theta _F`$ can be expressed in terms of the scattering matrix $`\widehat{S}`$ as

$$\theta _F(\omega )=\frac{1}{2i}\text{ln}\text{Det}\left[\widehat{S}(\omega )\right]$$
(8)
and thus one obtains a relation between the dot total density of states and the derivative of the scattering matrix with respect to the energy

$$\rho (\omega )=\frac{1}{2\pi i}\frac{}{\omega }\text{ln}\text{Det}\left[\widehat{S}(\omega )\right].$$
(9)
Then, by integrating this expression up to the Fermi energy, we obtain the generalized Friedel sum rule

$$N_{dot}=\frac{1}{\pi }\theta _F(E_F).$$
(10)
We see that $`\theta _F`$ is an important scattering phase. As it is related to the dot charge it should be a continuous function of the energy. It has also a simple relation to the eigenvalues of the scattering matrix. Due to the unitarity of $`\widehat{S}`$ its eingenvalues are of the form $`e^{2i\xi _{1,2}}`$ and thus $`\theta _F=\xi _1+\xi _2`$.

It should be emphasized that relation (10) holds only for the case of non-polarizable leads. For a more general case one should include also the charge induced on the leads in the Friedel sum rule. The deviation between the dot total density of states and the derivative of the Friedel phase with respect to energy has also been pointed out by Gasparian et al.  who analized the connection between densities of states and the scattering matrix for continuous models. In particular Eq. (14) in  can be written as

$$\frac{\theta _F}{\omega }=\pi \rho (\omega )\text{Im}\left(\frac{s_{L,L}+s_{R,R}}{4\omega }\right)$$
(11)
which coincides with our Eq. (6) provided that we make the approximation

$$\frac{\widehat{\mathrm{\Sigma }}}{\omega }\frac{i\text{Im}\widehat{\mathrm{\Sigma }}}{2\omega }.$$
(12)
Another scattering phase which is relevant for the interference phenomena observed in the experiments is the phase of the transmission amplitude $`\theta _t=\mathrm{arg}S_{LR}`$. In some particular cases (for instance in a one-dimensional problem) $`\theta _F`$ and $`\theta _t`$ coincide. However, as noticed recently by some authors , they are in general different. While $`\theta _F`$ is a continuous function, $`\theta _t`$ may not be defined at certain energies where the transmission vanishes. In order to be more precise, one can parametrize a general scattering matrix as

$$\widehat{S}=\left(\begin{array}{cc}ie^{i(\theta +\phi _1)}\mathrm{sin}\varphi & e^{i(\theta +\phi _2)}\mathrm{cos}\varphi \\ e^{i(\theta \phi _2)}\mathrm{cos}\varphi & ie^{i(\theta \phi _1)}\mathrm{sin}\varphi \end{array}\right).$$
(13)
with real phases $`\theta `$, $`\phi _1`$, $`\phi _2`$ and $`\varphi `$. It is then easy to show that $`\theta _F=\theta +\pi /2`$. On the other hand, when time reversal symmetry holds one has $`S_{LR}=S_{RL}`$ and thus $`\phi _2=0`$, in which case the argument of the transmission amplitude is related to $`\theta _F`$ by

$$\theta _t=\theta _F+\pi \mathrm{\Theta }(\mathrm{cos}\varphi )\frac{\pi }{2},$$
(14)
where $`\mathrm{\Theta }(x)`$ is the step function. Therefore, $`\theta _t`$ exhibits jumps of $`\pi `$ each time $`\mathrm{cos}\varphi `$ changes sign, i.e. at the points where there is a zero of the transmission. At these points the phase of the transmission amplitude deviates from the Friedel phase. Notice that the abrupt jump of $`\pi `$ of the phase of the AB oscillations between consecutive resonances is a central feature of the experimental results of Schuster et al. . We thus conclude that the study of the occurrence of zeroes of the transmission is essential to understand the experimentally observed behavior.

## IV Conditions for zeroes of the transmission

Within lattice models one can establish precise conditions for the occurrence of zeroes in the transmission amplitude through the dot. According to the Fisher-Lee relations the condition for having a zero in $`S_{LR}`$ at energy $`E_0`$ is that $`G_{1N}^r(E_0)=0`$. This is gives

$$G_{1N}^r(E_0)=\frac{C_{1N}(E_0\widehat{H}_{dot}\widehat{\mathrm{\Sigma }}^r(E_0))}{\text{Det}(E_0\widehat{H}_{dot}\widehat{\mathrm{\Sigma }}^r(E_0))}=0,$$
(15)
where $`C_{ij}(\widehat{A})`$ denotes the cofactor of the element $`i,j`$ in the matrix $`\widehat{A}`$. It is easy to see that the polynomial in the numerator is real and does not depend on the self-energy coupling of the dot to the leads at sites 1 and $`N`$. This is a direct consequence of having the coupling to the leads localized at these sites. We can thus reduce the condition for the zeroes, Eq. (15), to the simpler expression

$$C_{1N}\left(E_0\widehat{H}_{dot}\right)=0.$$
(16)
Eq. (16) clearly shows that the zeroes of the transmission are characteristic of the isolated dot structure and do not depend on the strength of the coupling to the leads. This is illustrated in Fig. 2 where the transmission for a 5x5 sites dot is shown for varying values of $`\mathrm{\Gamma }`$ within an energy range having a zero. One can observe that while the shape of the transmission varies substantially, the position of the zero is not affected.

This property allows one to relate the zeroes of the transmission to the wavefunctions of the isolated dot. In the weak coupling limit one can approximate $`G_{1N}`$ as

$$G_{1N}^a\underset{n}{}\frac{\psi _1^n\psi _N^n}{\omega \lambda _ni\mathrm{\Gamma }_L(\psi _1^n)^2i\mathrm{\Gamma }_R(\psi _N^n)^2},$$
(17)
where $`\lambda _n`$ and $`\psi _j^n`$ denote the eigenvalues and the amplitudes of the corresponding wavefunction for the isolated dot. The condition to have a zero between two consecutive eigenvalues $`\lambda _n,\lambda _{n+1}`$ is then simply given by $`\psi _1^n\psi _N^n\psi _1^{n+1}\psi _N^{n+1}>0`$.

We can now identify the sign of $`\psi _1^n\psi _N^n`$ as the parity of the corresponding dot wavefunction. By this reasoning we conclude that there should be a zero of the transmission in between dot states with the same parity regardless of the strength of the coupling to the leads.

In real systems inelastic scattering would prevent the occurrence of exact zeroes in the transmission. This situation can be described within our model by additional leads coupled to the dot as voltage probes . This effect is discussed in the next subsection.

### A simple example

A simple example which already exhibits a zero in the transmission amplitude is the case of a four sites sites dot, i.e. $`N=4`$. In this model sites 1 and 4 are the ones coupled to the leads. Sites 2 and 3 are coupled to sites 1 and 4 by hopping elements $`t`$ and we take the site energies on 1 and 4 as $`ϵ_1=ϵ_4=0`$. The transmission amplitude for such model is given by

$$S_{LR}=2i\sqrt{\mathrm{\Gamma }_L\mathrm{\Gamma }_R}\frac{2t^2(\omega \overline{ϵ})}{(\omega +i\mathrm{\Gamma }_L)(\omega +i\mathrm{\Gamma }_R)(\omega ϵ_2)(\omega ϵ_3)4t^2(\omega \overline{ϵ})(\omega +i\overline{\mathrm{\Gamma }})},$$
(18)
where $`\overline{ϵ}=(ϵ_2+ϵ_3)/2`$ and $`\overline{\mathrm{\Gamma }}=(\mathrm{\Gamma }_L+\mathrm{\Gamma }_R)/2`$.

For weak coupling and $`|ϵ_2ϵ_3|t`$ the transmission has two well resolved resonances at $`\omega ϵ_{2,3}`$. The modulus and the phase of $`S_{LR}`$ are shown in Fig. 3. As can be observed the phase of the transmission exhibits a similar evolution around each resonance. At the point $`\omega =\overline{ϵ}`$ there is a zero of the transmission and its phase exhibits an abrupt jump of $`\pi `$. From Eq. 18 it is clear that the zero is not dependent on the strength of the coupling to the leads.

Inelastic scattering can be simulated by additional voltage probes coupled to sites 2 and 3. Let us denote by $`\mathrm{\Gamma }_{inel}`$ the coupling to the additional leads. It is easy to see that the zero of the transmission now moves away from the real energy axis to $`\overline{ϵ}+i\mathrm{\Gamma }_{inel}`$. The jump in the phase of the transmission becomes smaller than $`\pi `$ and it is no longer abrupt but has a finite width given by $`\mathrm{\Gamma }_{inel}`$. The width of the phase jump can thus be taken as a measure of broadening of the dot levels due to inelastic scattering.

## V Statistical analysis of the transmission phase

According to the discussion of the previous section, the occurrence of zeroes in the transmission (and the associated jump in the transmission phase) is related to the parity of the isolated dot wavefunctions on consecutive levels. In a similar way as done by many authors for analyzing the mesoscopic conductance fluctuations of quantum dots in the Coulomb blockade regime , the parity should be studied statistically over many realizations of the dot potential.

In order to search for correlations of the parity of the dot wavefunctions we have numerically diagonalized dot lattice models up to $`30\times 30`$ sites. The confining potential is assumed to have the form of an isotropic parabola with curvature $`\alpha `$. We have investigated the influence of disorder and also the influence of different models for coupling the dot to the leads.

Figure 4a shows the spacing distribution between consecutive levels for various values of the disorder strength $`W`$ and $`\alpha =0.01`$ (in units of $`t/a^2`$, where $`a`$ is the lattice spacing). One can notice that the distribution approaches the the Wigner distribution

$`p(s)={\displaystyle \frac{\pi }{2}}s\mathrm{exp}\left[{\displaystyle \frac{\pi }{4}}s^2\right]`$

(plotted as a full line in Fig. 4a) as the disorder strength increases. This agreement suggests that our results for $`W>1`$ should be well described by random matrix theory.

For analyzing correlations in the parity of the wavefunctions we plot in Fig. 4b the probability to find consecutive levels with the same parity as a function of level spacing. As can be observed, this probability is almost constant as a function of level spacing. Although the probability is slightly larger than 0.5 for weak disorder ($`W=0.5`$) it approaches 0.5 as the disorder increases. This behavior indicates that correlations in the parity are negligible within this model.

The previous results correspond to the case where the dot is coupled to the leads at sites 1 and $`N`$, and thus the parity of the wavefunctions is given by $`\mathrm{arg}(\psi _1^n\psi _N^n)`$. This situation can be generalized to the case where many dot sites are coupled to each lead as indicated by the dotted lines in Fig. 1. We shall assume that both leads are coupled to the same number $`N_{leads}`$ of neighboring sites on each side of the dot. In this case the parity will be given by $`\mathrm{arg}(_{j_L,j_R}\psi _{j_L}^n\psi _{j_R}^n)`$, where $`j_{L,R}`$ denote the sites coupled to the left and right lead respectively. It should be emphasized that the conditions for the occurrence of zeroes discussed in section IV remains valid in this case, the first sites on the leads $`L`$ and $`R`$ playing the role of sites 1 and $`N`$.

In Fig. 5 we show the total probability to find consecutive levels with the same parity as a function of $`N_{leads}`$ and different values of the disorder strength $`W`$. For $`W=0.5`$ one can notice a slight decrease of the probability for increasing $`N_{leads}`$. The probability nevertheless remains always close to 0.5 indicating the absence of significant parity correlations even when changing the model for coupling the dot to the leads.

These results demonstrate that it is very unlikely to find a large sequence of dot levels with the same parity and consequently the difficulty to account for the experimental results of Ref.  within a one-electron model. There are, however, some peculiarities of the dot states in the limit of weak disorder which can give rise to correlation effects as discussed in the next section. These correlation effects are associated with certain dot states which are strongly coupled to the leads as discussed below.

The coupling strength $`\alpha _n=|_{j_L,j_R}\psi _{j_L}^n\psi _{j_R}^n|`$ is shown in figure 6 as a function of the level number for the case $`N_{leads}=5`$ and different values of $`W`$. In the case of extended contacts to the leads and weak disorder one can clearly distinguish states which are more strongly coupled to the leads than the average. For the case $`W=0.5`$, this states are indicated by arrows in Fig. 6. One can notice by their level number (67,79,92,106,..) that they appear in a well defined sequence, corresponding to the shell structure of states in the isotropic 2D harmonic oscillator. Contour plots of the wavefunctions with level number 67 and the next three levels are shown in Fig. 7. Although they correspond to nearly degenerate levels only the first one is strongly coupled to the leads.

## VI Effects due to Coulomb interactions

Up to now we have neglected interactions and concentrated on the one-electron properties of the dot. As commented in section II, electron-electron interactions can be included by means of the constant charging energy model. Within the traditional description of Coulomb blockade in quantum dots based on simple master equations  the main role of the charging energy is to open a gap between consecutive dot levels, dot states being populated according to the isolated dot ground state for each number of electrons. Within this picture, the evolution of the phase over a sequence of dot resonances will be fixed by the parity of the corresponding one-electron dot states. As discussed in the previous section, the parity of one-electron dot states does not exhibit significant correlations and therefore this picture is unable to account for the experimental results of Ref. .

The situation may change drastically when taking into account the finite coupling to the leads beyond lowest order perturbation theory. In this case correlation effects may lead to a population of the dot levels different to the one predicted by simple master equations. This possibility was recently pointed out by Silvestrov and Imry  and advanced as an explanation of the experimental results on the phase of the transmission through a quantum dot. The mechanism proposed by these authors would take place when there is a strongly coupled dot level followed by nearly degenerate weakly coupled levels. This situation is characteristic for the case of weak disorder discussed in the previous section.

The basic mechanism proposed by Silverstrov and Imry can be understood by considering a simple four sites problem with two (spinless) electrons. We shall assume that two of the sites correspond to the dot levels and the other two represent the left and right leads. Let us call $`ϵ_1`$ and $`ϵ_2`$ the one-electron energies of the isolated dot levels and $`\mathrm{\Delta }=ϵ_2ϵ_1>0`$ the spacing between them. The lower level is symmetrically coupled to the leads by effective hopping elements $`t_{1L}=t_{1R}=T`$ which are much larger than the ones corresponding to the upper level $`t_{2L}=t_{2R}=t`$. The energy levels on the leads sites are taken as zero. Let us analyze how the dot levels are populated as the gate voltage shift down the energy levels from $`ϵ_1,ϵ_20`$ to $`ϵ_1,ϵ_20`$.

Finding the ground state of this problem requires in principle to diagonalize the $`6\times 6`$ matrix corresponding to the system Hamiltonian in the basis $`|n_1n_2n_Ln_R>`$ where $`(n_1,n_2)`$ are the occupation numbers for the dot levels and $`(n_L,n_R)`$ correspond to the leads. The problem can, however, be reduced by a basis change in which we replace the states with one electron in the leads by its symmetric and anti-symmetric combination, i.e $`|1010>`$, $`|1001>`$, $`|0110>`$ and $`|0101>`$ are replaced by $`(|1010>\pm |1001>)/\sqrt{2}`$ and $`(|0110>\pm |0101>)/\sqrt{2}`$. In this way, the initial $`6\times 6`$ matrix is reduced into two $`3\times 3`$ blocks: one block corresponds to the empty dot states which couples to the anti-symmetric combinations and the other to the doubly occupied dot state coupled to the symmetric combinations, having the form

$`H^s=\left(\begin{array}{ccc}ϵ_1+ϵ_2+E_c& \sqrt{2}T& \sqrt{2}t\\ \sqrt{2}t& ϵ_2& 0\\ \sqrt{2}t& 0& ϵ_1\end{array}\right)H^a=\left(\begin{array}{ccc}ϵ_1& 0& \sqrt{2}T\\ 0& ϵ_2& \sqrt{2}t\\ \sqrt{2}T& \sqrt{2}t& 0\end{array}\right)`$

In the limit $`\mathrm{\Delta }0`$, i.e. $`ϵ_2ϵ_1=V_g`$ the ground state for each symmetry is given by

$`\lambda ^s={\displaystyle \frac{3V_g+E_c}{2}}\sqrt{\left({\displaystyle \frac{V_g+E_c}{2}}\right)^2+2(T^2+t^2)}`$

$`\lambda ^a={\displaystyle \frac{V_g}{2}}\sqrt{\left({\displaystyle \frac{V_g}{2}}\right)^2+2(T^2+t^2)}`$

When $`V_g0`$ the system starts in the anstisymmetric ground state and the charge of the dot levels evolve with $`V_g`$ according to

$`<n_1>={\displaystyle \frac{2T^2}{(V_g+\lambda ^a)^2+2(T^2+t^2)}}<n_2>={\displaystyle \frac{2t^2}{(V_g+\lambda ^a)^2+2(T^2+t^2)}}`$which shows that the dot levels start populating when $`V_gT`$. However, as $`tT`$, the charge in the weakly coupled level will be negligible as compared to the strongly coupled one.

At $`V_g=E_c/2`$ there is a crossing between $`\lambda ^a`$ and $`\lambda ^s`$. As a result the symmetric ground state becomes more stable and the charge on each dot level is now given by

$`<n_1>={\displaystyle \frac{2t^2+(V_g+\lambda ^s)^2}{(V_g+\lambda ^s)^2+2(T^2+t^2)}},<n_2>={\displaystyle \frac{2T^2+(V_g+\lambda ^s)^2}{(V_g+\lambda ^s)^2+2(T^2+t^2)}}`$

Thus, when crossing $`V_g=E_c/2`$ the charge on the strongly coupled state goes from $`T^2/(T^2+t^2)`$ to $`t^2/(T^2+t^2)`$ while the charge on the weakly coupled level does the opposite evolution (we neglect small corrections of order $`(T/E_c)^2`$). We see that nearly one electron is transferred from the strongly coupled level to the weakly coupled one in the middle of the charging curve, the total charge in the dot remaining constant (up to order $`(T/E_c)^2`$). The strongly coupled level is again filled when $`V_g`$ is close to $`E_c`$. The situation is similar for finite $`\mathrm{\Delta }T`$, the abrupt jump is also present but shifted to $`V_gE_c/2+\mathrm{\Delta }`$.

A more realistic description of the actual situation requires the inclusion of many sites to represent both the leads and the dot levels. One would need a large number of sites in order to simulate a continuous density of states on the leads. Using Lanczos method we are able to numerically diagonalize interacting systems with up to 16 sites and 8 electrons. The dot levels charge as a function of gate voltage for a systems with 6 sites on a line for representing each lead and 4 sites for representing the dot levels, is shown in Fig. 8. The dot levels are coupled symmetrically to the leads by hopping parameters $`t_1=0.1`$ and $`t_2=t_3=t_4=0.02`$ (all energies in units of the charging energy $`E_c`$) while the hopping parameter within the leads is taken as $`0.1`$. The weakly coupled levels lie at $`0.025`$, $`0.030`$, and $`0.035`$ respectively above the strongly coupled level. As can be observed, the strongly coupled level is repeatedly charged and discharged in a similar way as for the simple four sites model. One can notice, however, that the jumps in the charge of the strongly coupled level are not so much abrupt and become progressively less pronounced.

The next question is how this particular charging of the dot levels affects the transmission amplitude and its phase. In order to obtain the transmission through the interacting dot we map this finite size cluster into an effective non-interacting system with the same charge on each dot level for each value of the gate voltage. In this effective system we use the same hopping parameters as in the cluster calculation but replace the leads by infinite one-dimensional chains. The effective dot levels are determined self-consistently to get the correct level charges. The phase of the transmission amplitude and the transmission probability thus obtained are shown in Fig. 9. As expected, the phase follows closely the charge of the strongly coupled level. It should be remarked that this behavior is not dependent on the parity of the weakly coupled levels. On the other hand, the transmission probability exhibits broad peaks that can be associated with the repeated charging of the strongly coupled level together with narrow satellite peaks corresponding to the charging of the weakly coupled levels. These narrow peaks lie closer to the broader peaks as the cycle number increases. In the fourth cycle the narrow peak can hardly be resolved.

## VII Conclusions and final remarks

To summarize, we have analyzed the behavior of scattering phases in quantum dots using lattice models. We have first studied the definition of the different scattering phases in terms of Green functions. We have shown that abrupt jumps of $`\pi `$ in the phase of the transmission amplitude are associated with the occurrence of transmission zeroes. Then, we have shown that within lattice models and assuming single moded leads the zeroes are independent of the strength of the coupling to the leads . This property allows us to relate the occurrence of a zero with the parity of consecutive dot states in the isolated dot.

We have studied the statistical properties of the isolated dot states as a function of the disorder strength. This analysis reveals that there are no significant parity correlations between consecutive dot levels. For moderate disorder there appear states which are much more strongly coupled to the leads than the average. We have studied charging effects in this situation using exact diagonalization of small clusters. We have shown that the dot resonances as a function of gate voltage may correspond to the strongly coupled level through several cycles.

The viability of this mechanism as a possible explanation of the “phase problem” deserves further discussion. The presence of strongly coupled levels followed by many weakly coupled levels is a reminiscent of the integrable case which gradually disappears for increasing disorder (see Fig. 6). In this sense, the mechanism is not universal as it would not hold for a fully chaotic quantum dot. Thus, this mechanism could be tested experimentally by leading the dot into a fully chaotic situation (using, for instance additional gates to distort the dot shape, as in Ref. ). On the other hand, the transmission as a function of gate voltage exhibits narrow satellite peaks close to the broader Coulomb blockade peaks (see Fig. 9) which are not observed experimentally. Being so narrow these peaks could be easily washed out by thermal broadening or due to the finite bias voltage used in the measurements. Of course if certain narrow peaks are not resolved in the experiment, a number of additional scenarios might account for parity conservation, even within a non-interacting picture. Differing parity conserving scenarios can then still be distinguished according to their phase evolution pattern. It is very unlikely that a non-interacting theory can generate a a phase evolution like that of Fig. 9 which resembles closely that seen in the experiments.

One should also comment about the role of spin, not included in the present analysis. At very low temperatures spin degeneracy would lead to a Kondo effect. This effect has been recently observed experimentally  but in somewhat smaller quantum dots and would be reflected in the phase behavior as a plateau at $`\pi /2`$ in between a pair of resonances corresponding to the same dot level . At temperatures larger than the Kondo temperature spin degeneracy could at most account for the similar phase behavior over two consecutive dot resonances. Other spin effects, like Hund’s rule require almost perfect symmetry leading to orbital degeneracy and so far have been observed only in ultrasmall quantum dots .

To conclude, we point out that our analysis has been restricted to the case of single-moded leads. This is an assumption which has been implicitly been used in all the theoretical models of the problem following the characterization of the experimental set up given in Ref. . Nevertheless, the analysis of what would be measured in a multichannel situation is of interest and deserves further theoretical investigation.

###### Acknowledgements.

The authors would like to thank T. Taniguchi for fruitful discussions. One of us (ALY) also thanks discussions with F. Flores and J.A. Vergés. This work has been partially supported by the Spanish CICyT under contract PB97-0044 and the Swiss National Science Foundation. Appendix A We show in this appendix the unitarity of the scattering matrix defined from the Fisher-Lee relations. One has   $`\left(\widehat{S}^{}\widehat{S}\right)_{11}`$ $`=`$ $`|12i\mathrm{\Gamma }_LG_{11}^r|^2+4\mathrm{\Gamma }_L\mathrm{\Gamma }_R|G_{1N}^r|^2`$ (19)
$`\left(\widehat{S}^{}\widehat{S}\right)_{12}`$ $`=`$ $`2i\sqrt{\mathrm{\Gamma }_L\mathrm{\Gamma }_R}\left[(1+2i\mathrm{\Gamma }_LG_{11}^a)G_{1N}^rG_{1N}^a(12i\mathrm{\Gamma }_RG_{NN}^r)\right],`$ (20)
and similar expressions for $`\left(\widehat{S}^{}\widehat{S}\right)_{21}`$ and $`\left(\widehat{S}^{}\widehat{S}\right)_{22}`$. These elements can be rewritten as

$`\left(\widehat{S}^{}\widehat{S}\right)_{11}`$ $`=`$ $`2i\mathrm{\Gamma }_L\left[G_{11}^aG_{11}^r2i\mathrm{\Gamma }_L|G_{11}^a|^22i\mathrm{\Gamma }_R|G_{1N}^a|^2\right]+1`$ (21)
$`\left(\widehat{S}^{}\widehat{S}\right)_{12}`$ $`=`$ $`2i\sqrt{\mathrm{\Gamma }_L\mathrm{\Gamma }_R}\left[G_{1N}^rG_{1N}^a+2i\mathrm{\Gamma }_LG_{11}^aG_{1N}^r+2i\mathrm{\Gamma }_RG_{1N}^aG_{NN}^r\right].`$ (22)
It is now easy to show that the expressions between brackets vanish identically. This fact is a consequence of the identity

$$\widehat{G}^a\widehat{G}^r=\widehat{G}^a\left[\widehat{\mathrm{\Sigma }}^a\widehat{\mathrm{\Sigma }}^r\right]\widehat{G}^r$$
(23)
which follows from the definition of the Green operator. Taking the matrix elements (1,1) and (1,N) one finds

$`G_{11}^aG_{11}^r`$ $`=`$ $`2i\mathrm{\Gamma }_L|G_{11}^a|^2+2i\mathrm{\Gamma }_R|G_{1N}^a|^2`$ (24)
$`G_{1N}^aG_{1N}^r`$ $`=`$ $`2i\mathrm{\Gamma }_LG_{11}^aG_{1N}^r+2i\mathrm{\Gamma }_RG_{1N}^aG_{NN}^r`$ (25)
Appendix B

We show in this appendix that $`2i\theta _F=\text{ln}\text{Det}\left[S\right]`$. First notice that

$$\theta _F=\text{Im}\text{ln}\text{Det}\left[\omega \widehat{H}_{dot}\mathrm{\Sigma }^a\right]=\frac{1}{2i}\text{ln}\text{Det}\left[(\omega \widehat{H}_{dot}\mathrm{\Sigma }^a)\widehat{G}^r\right].$$
(27)
From the definition of the Green operator one can rewrite this expression as

$$\theta _F=\frac{1}{2i}\text{ln}\text{Det}\left[1+(\mathrm{\Sigma }^r\mathrm{\Sigma }^a)G^r\right].$$
(28)
Finally, using that $`(\mathrm{\Sigma }^r\mathrm{\Sigma }^a)_{i,j}=2i\mathrm{\Gamma }_L\delta _{i,1}\delta _{j,1}2i\mathrm{\Gamma }_R\delta _{i,N}\delta _{j,N}`$, one can easily check that

$$\theta _F=\frac{1}{2i}\text{ln}\text{Det}\left[S\right].$$
(29)


format: text

text:
# Further Adventures: Oxygen Burning in a Convective Shell

## 1. Introduction

There are many published studies of stellar convection using multidimensional hydrodynamical simulations, but few deal with convective nuclear burning occurring in stellar interior (see Deupree (1998) for a study of core hydrogen convection). Arnett (1994) (A94) and Bazàn & Arnett (1994, 1998) (BA94, BA98) have studied oxygen burning shell which is one of the last stages in a massive star presupernova evolution; we extend that work.

This is an important stage in presupernova evolution because in this convective region several phenomena take place. In or near this region: (1) most of the explosive nucleosynthesis and production of 56Ni and 57Ni occurs, (2) the “mass cut” between collapsed and ejected matter develops, and (3) mixing of different layers may happen. The standard model for treating convection in one dimensional (1D) stellar evolutionary codes, the mixing length theory (MLT), is usually used for modeling this convective region as well, even though the conditions of oxygen convective burning shell are more complicated than can be assumed for MLT to be valid (i.e., the flow is not strongly subsonic and there are both energy sources and sinks in the flow).

In Arnett (1994) this evolutionary stage was described, and the first two dimensional (2D) hydrodynamical simulation of this problem were presented. The 145 s time interval was simulated by these calculations, using the PROMETHEUS code, is much less than the duration of this stage ($`10^3`$ to $`10^5\mathrm{s}`$, see figures 10.5 and 10.6 in Arnett (1996)). During this time, convective flow was formed. Bazàn & Arnett (1994, 1998) have modeled the evolution over a longer time interval of 300 to 400 s, and noticed a mixing of composition from the neighboring stable region near the end of the simulations. In addition,

* the flow velocities were up to $`10\%20\%`$ of the sound speed,
* stellar structure was not altered significantly, and
* there were strong density and composition fluctuations in both space and time which did not reach a statistical steady state.

The main differences between this work and that of A94, BA94, and BA98 were (1) the use of a different hydrodynamic code, and (2) the simulation of a longer evolution time. The same initial model, equation of state, and nuclear reaction algorithms were used.

## 2. The numerical scheme

For this study, we used a version of the hydrodynamical code VULCAN (Livne (1993)). This code uses an algorithm that begins with a hydrodynamic time step, and is completed by a relaxation of the numerical grid, thus giving an Arbitrary Lagrangian Eulerian (ALE) scheme. In the hydrodynamic time step, the Lagrangian hydrodynamic equations in two dimensions are solved explicitly or implicitly, allowing longer time steps to be taken if the flow is strongly subsonic. The mesh relaxation phase is necessary to eliminate distortion of the cells, especially for flow with vorticity. We used a quasi-1D Lagrangian relaxation, in which each radial row of cells kept a constant mass. The code was used by Glasner & Livne (1995) for computing convective novae outbursts, and by Asida & Tuchman (1997) and Asida (2000) to simulate convection in a red giant envelope. One of the adaptation made by Asida (2000) for simularions of convection was the inclusion of Smagorinsky (1963) like sub grid scale mixing (SGSM) model, as was done in many two dimensional and three dimensional numerical studies of convection.

The equation of state and the thermonuclear reactions algorithms were essentially the same as in A94 and BA98. The equation of state was the sum of components for electrons, ions, and radiation. The nuclear reaction network used twelve species for helium, carbon, neon and oxygen burning. Neutrino cooling was included; see the above references for details. The initial model was the same as in A94 and BA98. The computational domain corresponds to a region containing the entire oxygen-burning shell in a star of $`20M_{}`$, having an initial metallicity of 0.007 (about one third solar). This shell boundaries are at radii of $`3\times 10^8\mathrm{cm}`$ and $`3\times 10^9\mathrm{cm}`$ that are equivalent to mass coordinates of $`1.4M_{}`$ and $`2.4M_{}`$.

We performed several simulations that were different in: the computational domain, spacial resolution, the initial temperature gradient and numerical parameters of the simulation. The computational domain typically includes all of the oxygen shell as well as the neighboring layers (the inner radius was located at $`2\times 10^8\mathrm{cm}`$ in some of the simulations, and the outer radius was located at $`8\times 10^{10}\mathrm{cm}`$ in other simulations). The angular extent of the wedge ran usually from 0.35 $`\pi `$ to 0.65 $`\pi `$ radians.

Typically. the oxygent shell was divided to $``$ 120 radial zones, and the spacing of zones was logarithmic in radius and linear in angular direction (i.e., $`dr=rd\theta `$) with 60 angular zones. Rotational symmetry was assumed. These characteristics are similar to the medium resolution simulations of BA98. Because of small differences in the equation of state in the initial one dimensional model and in the two dimensional code, as well as different radial zoning, we had two sets of simulations: in the first, the 1D model (with interpolation) was used as it is, and in the second it was slightly adapted, so that the temperature gradient would be superadiabatic in all zones of the convection shell.

The velocities on the inner boundary were set to be zero for the whole simulation, so that the inner core was a hard sphere. At the upper boundary there was no limitation on the velocities, but in order to eliminate mass flow out of the computational domain, an average radius was used to follow expansion or shrinking of the outer boundary. We used reflective boundary conditions on the sides (though such conditions enforce a downflow or upflow on the side boundaries, it was found in BA98 not to be important).

## 3. Results of the simulations

We present the results of one “standard” simulation with 172 radial zones and 60 angular zones. The inner boundary was at $`2\times 10^8\mathrm{cm}`$ and the outer was at $`48\times 10^8\mathrm{cm}`$. The temperature gradient was slightly super adiabatic in the oxygen shell, and no SGSM terms were used. Most of the results of the other simulations were similar and the differences are discussed mainly in the end of this section.

### 3.1. General Evolution of the Flow

In our simulations, the initial velocities were zero, and the convective flow developed as a result of the instability from round off errors. Figure 1 presents the velocity field in the beginning of the simulations for times (a) 75 s and 150 s. As we can see, the convective flow starts at the bottom of the convection region (i.e., the burning layer), and then moves up with increasing eddy size. By time 150 seconds the convective flow penetrates the upper boundary of the convection region (seen as a thick line in panel a). As a result, there is a downflow of carbon-rich material. This can be seen in Figure 2, which presents contours of carbon nucleon fraction. Panels a-e represent times of 75, 150, 300, 600 and 1200 s. We can see that the downflow (panel b) penetrates the whole convective region, and results in mixing of carbon in this region. From comparison of panels c, d and e we can see that carbon abundance becomes more uniform, as the simulation evolves.

This penetration of carbon is almost identical to that seen by BA98; see their Figure 5. The two independent hydrocodes give consistent results over the whole time spanned (400 seconds) by BA98 simulations. The small difference in the timescale for the carbon penetration is due to the small differences in the extent of the super adiabatic gradient of the initial 1D model (when we used the initial 1D model as it is, without modifying the temperature gradient, we got very similar results with a longer timescale of the penetration).

### 3.2. Carbon Enrichment

Carbon penetration, as well as other processes in the simulations, can be visualized by examining one dimensional averages of various parameters, which in principle would be equivalent to results from 1D simulations. In Figure 3 carbon abundance (nucleon fraction) is plotted as a function of stellar mass coordinate. From the initial profile (solid line) we can see the location of the oxygen shell. At 150 s (long dashed line), carbon rich material has entered the convection layer, and then mixed through the whole region. At later times, the fluctuations in the abundances of carbon decrease (compare the profiles at 300 s -short dashed line, 600 s - dotted line and 1200 s -dotted dashed line), and we have a carbon nucleon fraction of $`5\times 10^3`$.

FIG. 3.— Average carbon abundance

A change in abundance in this plot does not necessarily corresponds to mixing, because this plot represents an one dimensional average of the compositions of all the material in each radial layer. This may be revealed by a closer look at the apparent widening of the jump in carbon abundance at the upper boundary of the oxygen shell at $`m`$$`2.4M_{}`$: this jump is much wider at 150 s than it is at later times. This “anti diffusion” is possible since the widening of the interface of the shell is the outcome of two processes: a mixture of material from the two sides of the interface and large scale motions that change the shape of the interface. Thus, at later times the interface is more spherical than it was at 150 s, so that averaging the abundances over radial layers yields a sharper jump (compare Fig 2 panels b and d).

FIG. 4.— RMS fluctuations of carbon abundance

To better understand the carbon enrichment we present Figure 4 which shows the RMS fluctuations of carbon nucleon fraction ($`\sigma `$) as a function of stellar mass coordinate at several different times. The tendency toward more uniform carbon mixing is clearly seen, as was indicated by Figures 3 and 4. From this figure, we can also see, that the previously mentioned widening of the oxygen shell interface is partially due to changes in its shape since the fluctuations in carbon fraction at the boundary are relatively high, and are even higher at time 150 s.

FIG. 5.— Nuclear luminosity as a function of time.

The mixing of carbon in the burning layer causes the reaction rate to change significantly. In Figure 5, the nuclear luminosity is presented as a function of time. In the beginning there is an adjustment phase in which the energy production rate decreases to about $`5\times 10^{44}\mathrm{erg}\mathrm{s}^1`$ . This transient phase is a thermal relaxation due to the fact that the initial model is inconsistent: it does not have a two dimensional velocity field which can carry the convective energy flux that the 1D model needed. As the carbon rich material penetrates to high temperature layers, it starts to react rapidly, yielding nuclear luminosities which increase to about hundred times the initial value. In this stage, the result is mainly carbon and neon consumption. Afterward, there is a smaller decrease in nuclear luminosity (compare time 600 and 1200 s). This decrease is related to the small decrease in the average carbon abundance seen between times 600 s and 1200 s in Figure 3.

### 3.3. A Quasi-Steady State

FIG. 6.— Time average change in energy

The relaxation in the carbon abundances by time 300 s to a slowly varying state, and the small change between times 600 and 1200 s, suggest that the system may be close to a quasi - stationary state (a “steady” state on average). This one is significantly different from the steady state predicted by the one dimensional calculations used to definethe initial model. That the system is close to a steady state can be demonstrated by the following figures. Figure 6 shows the average rate of change in energy, as a function of mass coordinate (from time 400 to time 800 s). The energy produced by thermonuclear reactions (dashed line) is carried away by convection (solid line). The net change in energy is small, and balanced by changes in the gravitational energy (dotted line). This slow change is an adjustment of the structure driven by enhanced nuclear burning from the carbon ingestion. There is a net heating at the bottom, below $`1.4M_{}`$, which will be addressed in the next section.

FIG. 7.— Time average change in carbon abundance

Figure 7 shows the average rate of change in carbon nucleon fraction, as a function of mass coordinate (from time 400 s to time 800 s, and from time 800 s to time 1200 s), comparing the changes due to the divergence of the convective abundance flux (solid line) with that due to nuclear burning (long dash line). The nuclear consumption of carbon at the bottom is almost completely balanced by convective inflow; This is the burning zone. Over most of the convection region the change in carbon abundances is small. There are larger fluctuations at the top (outer) boundary of the convection zone (where the carbon abundance increases by a factor of $`20`$, see Fig. 3).

The decrease in nuclear luminosity (seen in Figure 5, after the peak at 200 s) was caused by a decrease in the net flux of carbon into the burning layer. This is easily seen in Figure 7 in the comparison of the changes in the later time interval (800 s to 1200 s - short dashed line) to the changes in the previous time interval (400 s to 800 s - solid line). The total mass of carbon in the oxygen convection shell is decreasing as a result of carbon consumption and less mixing from the upper shells. Thus our convective burning is beginning to evolve on a secular time scale by the consumption of fuel, having approached a thermal “steady” state.

in Figure 8 we present the RMS fluctuations of density (divided by the density) as a function of stellar mass coordinate for the standard simulation for several times. As the composition becomes more uniform with time, the density fluctuation decreases throughout most of the oxygen shell to about $`10^3`$ at time 1200 s. However at the edges of the convection zone, the fluctuations does not decrease and we have about $`2\times 10^2`$ at the bottom and $`7\times 10^2`$ at the top, where there is a jump in the abundances. These values are similar to those obtained by BA98.

FIG. 8.— RMS fluctuations of the density

### 3.4. Interaction with Neighboring Shells

FIG. 9.— Temperature profile at the bottom

One important interaction of the convective shell with neighboring shells is the mixing of carbon from the upper carbon rich shells, as described before. Another important interaction is heating of the nearest shells below the bottom of the convection zone: in the initial model there is a sharp decrease in temperature just below the convection zone, as convective flow starts to grow inside the lower parts of the convection zone there is some penetration of the flow to these lower temperature shells. As a result of this flow, the lower temperature shells are being heated and the higher temperature shells are being cooled. This heating is easily seen in Figure 9 where we present the temperature profile in the initial model, and at later times. Most of the heating was done by time 300 s, but even at later times, some heating exsisted, as can be seen in the left end of Figure 6. Along with this “energy mixing” there is mixing of composition in those few zones as we can see from neutron excess $`\eta `$ plot presented in Figure 10. In the initial model, there is a jump in $`\eta `$ at the boundary of the oxygen shell (corresponds to the jump in composition). Due to mixing, this jump slightly moves towards lower mass coordinate and becomes slightly wider. From comparison of Fig. 9 and 10 we see that the heating penetrates to deeper shells than composition changes.

FIG. 10.— Neutron excess $`\eta `$ at the bottom

There are two important consequences of this interaction: (1) the heating and aditional oxygen and carbon in the few zones below the 1D oxygen convection shell cause these zones to be added to the oxygen shell and, together with few zones at the bottom of it, to be part of the burning zone (we can see a broader region with temperature above $`2.2\times 10^9\mathrm{cm}`$), (2) because of the cooling of the zones at the bottom of the oxygen shell, the temperature gradient becomes less than adiabatic in those zones and the convective flow exist as a result of a super adiabatic gradient at higher zones.

FIG. 11.— Convective velocity at the bottom

Thus we see that the neighboring stable shells greatly affect the oxygen convective shell. The same convective flow that penetrates to the neighboring shells and allow the mixing of elements and energy between the shells affects further more the stable shells. In Figure 11 we present the averaged convective velocity as a function of radius at the bottom of the oxygen shell. The dashed line corresponds to the fluctuations in radial component of the velocity, while the solid line corresponds to the amplitude of the fluctuating total velocity. The plus signs present points in which the pressure varies by a factor of ten. Three regions can be identified in this plot from right to left: a region of constant velocity, a region of exponential decrease, and another region of constant velocitis but with a noticeable difference between the total and the radial velocities. (i.e. the radial velocities are much smaller on average than the tangential velocities) which is indicative of g modes. This g modes can be identified in 2D presentation of the flow (Figure 12).

FIG. 12.— Velocity field at the bottom. The thick line represents the composition jump that indicates the boundary of the oxygen shell

In another simulation we included a much larger shell above the oxygen shell. In this simulation we noticed similar flow characteristics of exponential decay (Figure 13) and tangential preference of the flow (Figure 14). From Fig. 11 and 13 we can also see that the exponential decrease at the bottom is on a scale which is about one third of the pressure scale height, while at the top it is over one pressure scale height.

FIG. 13.— Convective velocity at the top

FIG. 14.— Velocity field at the top. The thick line represents the composition jump that indicates the boundary of the oxygen shell

### 3.5. Numerical Sensitivity

We performed several simulations with different parameters as mentioned in §2. As our results are consistent with those of BA98, and since one of the main subjects in BA98 is the sensitivity of the results to numerical parameters, we would not present here all the results from our many simulations. Instead, we will focus on the conclusions from those simulations: the most sensitive feature in the results is the amount of heating below the oxygen shell. When there is more heating, the temperature gradient is less than adiabatic in a larger portion of the (bottom of) oxygen shell. As a result, the convective flow is damped, carbon can not reach the burning zone and we got less nuclear luminosity.

The opposite happened in simulations where we did not include the shell below the oxygen shell: the temperature gradient was super adiabatic so the velocities were higher (by a factor of five), more carbon entered the burning zone, and the nuclear luminosity was higher. From comparison of simulations with different resolution for this case (without the bottom shell), we found that the results are quite robust and not sensitive to the numerical resolution.

In all of our simulations the other features were essentially the same, namely: the evolution of the flow from the bottom, the penetration to neighboring shells and the enrichment of the oxygen shell by carbon.

FIG. 15.— Extra $`\eta `$ mixing at the bottom

In order to check extra mixing at scales shorter than numerical resolution, we performed a simulation with SGSM mixing terms. This simulation started from a 2D profile of the standard simulation at time 400 s and was carried out to time 800 s. The differences between this simulation and the standard simulation were: a more uniform composition at the oxygen shell - the RMS fluctuations of carbon abundance at time 800 s were similar to the results of the standard simulation at time 1200 s, and significantly more mixing at the stable bottom shell. This can be seen in Figure 15 were we plot the neutron excess $`\eta `$ for this simulation (SGSM) and the standard simulation. In this plot we can see a very small difference between the two profiles of the standard simulation (at time 400 s - solid line and at time 800 s - long dashed line). When SGSM terms were added, more mixing can be seen by time 600 s (short dashed line), and additional mixing occurs until the end of the simulation (time 800 s - dotted line).

## 4. Discussion

When we compare our results to the results of BA98 we note that the evolution of the simulations is similar for times reached previously. In both, convective flow evolves from bottom to top, with velocities exceeding $`10\%`$ of the sound speed. The flow penetrates to the carbon rich region above the one dimensional edge of the convective region, and thus causes a significant enrichment of the convective region with carbon. These higher values of carbon abundances yield a much higher nuclear luminosity than did the original oxygen burning.

However, at later times, we could identify signs of a steady state configuration in which the average abundances changed very slowly, and turbulent mixing within the convective region was able to decrease the fluctuations in the oxygen shell. However, density fluctuations near the boundaries of the convection zone (which are also composition discontinuities) did not decrease at later times, and were equal to few percents.

The most prominent feature in our simulations is the interaction of the convective shell with the neighboring shells. This interaction is due to penetration of the convective flow to the stable shells and the resulting mixing of elements and energy. These effects were present in all of our simulations, and seemed to be physically resonable and valid. However, the exact amount of mixing (especially at the bottom) depends on numerical parameters of the simulations. Moreover, since the initial 1D model did not include such overshoot, we actually simulated a transient toward a more consistent model. For example, we mentioned that the heating of the very few neighboring zones ,at the bottom of the convection shell, lowers the temperature gradient below the adiabatic value. If the convective flow is damped because of that, then the heat generation in these zones would evantually increase the temperature gradient, and convection would be restored.

The properties of the internal modes excited in the stable shells are undoubtedly affected by the computational domain and the boundary condition. Reflective waves are probably the cause for the constant level of velocities we noticed at Figure 11. However, the volume below the oxygen shell in these stars is relatively small so a finite value of the velocity is quite plausible.

An interesting point is the validity of the results as we performed 2D simulations (and not 3D). This is not an easy question to answer, however, from various comparisons of 2D and 3D simulations (see for example Kiraga et al. (1999)) it seems that the main features of the interaction of the unstable layer and its neighboring layers are similar, and the overshoot and mixingin in both cases are comparable . As we do not claim to resolve this interaction quantitatively, we think that our results are valid.

Penetration of the convective flow to neighboring stable regions is a common feature that exists in many multidimensional simulations of convection, in a variety of stellar problems. Hurlburt, Toomre, & Massaguer (1986) have studied such penetration beyond a shallow convective region in 2D, and found that it generated internal g modes in the stable region below the convection region, all the way to the bottom boundary. Two regions in this stable layer can be identified from their figure 13: close to the unstable layer there is a sharp decrease in the kinetic energy, but it does not go to zero, rather there is an almost constant level of kinetic energy throughout the rest of the stable layer. Kiraga et al. (1999) have tested penetration in both 2D and 3D, they found that the energy fluxes in the stable layer in the 3D are about two thirds of the 2D results and that the constant level is replaced with a continuas decrease when the viscosity coefficient is increased at the bottom to avoid reflection of waves.

In both these studies, ideal simplified physics was assumed. Freytag, Ludwig, & Steffen (1996) studied the structure and dynamics of shallow stellar surface convection zones, using two dimensional radiation hydrodynamic simulations with more “real” physics. They found convective motions extending “well beyond the boundary of convectively unstable region, with vertical velocities decaying exponentially with depth in the deeper parts of the lower overshoot region, as expected for linear $`g^{}`$ modes.” They suggested to approximate the average velocity at the stable region as the convective velocity at the boundary of the unstable zone multiplied by an exponential decay factor with a lapse rate that is of order of the pressure scale height.

As explained by Cox (1980, §17.8) g- modes are unstable within convection regions, and generally decrease exponentially with increasing distance from the boundaries of the oscillatory region. As Freytag, Ludwig, & Steffen (1996) point out, Landau & Lifshitz (1959) have a relevant discussion in their §34, where they discuss some properties of potential flow. If compressibility and dissipation can be neglected, a steady state potential flow which is periodic in some plane, must be damped exponentially in a direction perpendicular to that plane. Also, the shortest wavelengths will be damped fastest. Consequently, most of the overshoot is given by the largest scales which are driven at the interface, that is, the largest convective scale.

We find a striking underlying unity in these results and ours, despite the significantly different astrophysical context involved: our simulations correspond to deep shell burning convection that is driven by nuclear burning at the bottom and cooled by neutrino radiation and in which spherical properties of the domain are important; while these studies deal with envelope convection driven by photon radiation at the surface in plane parallel geometry. In both cases The flow extand beyond the formal boundaries of the convection zone, with velocities that decay exponentially in the stable shells producing internal modes.

Near the boundaries of the convection zone, where the velocities are still large, quite efficient mixing is taking place. The amount of mixing beyond that region, is not entirely clear. Freytag, Ludwig, & Steffen (1996) have used trajectories of test particles to estimate the diffusion coefficient and conclude that the diffusion coefficient corresponds directly to the exponentially decaying velocities.

Some studies of stellar evolution that used these mixing terms for hydrogen core convection (Herwig et al. (1997); Herwig et al. 1998), found that the lapse rate of decay of the diffusion should be much smaller ($``$2 percents) than the pressure scale height. They speculate that this is due to the large stability of the layers that are near the core convection zone. Yet another explanation might be true, namely the flow characters at the stable shells and the many scales involved in mixing.

Mixing is a complex phenomenon which take place in a very small length scale comparing to the length scale of the star, and since stars have a lower effective viscosity than in the simulations, and since they evolve for longer times, an accurate prediction of mixing is very difficult. Further, we note that motion does not necessarily imply mixing! In particular, oscillatory potential flow, which is irrotational,would involve “stretching” and “contracting” motions that do not give mixing, at least at lowest order.

As we can see in Figures 11-14, The flow in the stable shells is different than the flow inside the convection region, even though it is not purely potential flow (as revealed when examining the vorticity). Due to these differences in the flow, it seems plausible that the effective mixing velocity may be less than the hydrodynamic velocity, and its lapse rate different as well.

As the various shells in our model have a different composition, mixing can be directly examined from the 2D profiles and the 1D average of composition. However, as the scale of mixing is much smaller than numerical resolution, the results must be carefully examined. One way to examine the results is to compare the mixing in the standard simulation to that in the simulation with sub grid mixing model. There are small differences in the average abundances in the oxygen shell between the two simulations, and the main difference is that with SGSM terms, the profile is more uniform because of the extra mixing. The difference between the simulations is much more prominent in the lower stable shell: the width of the abundance jump is much wider when SGSM terms were used, and penetrates to deeper zones, in which the “constant level” average velocity exists.

When we used SGSM we implicitly assumed that turbulent flow exists from the scale of the mesh resolution to the scale of molecular diffusion. Since this mixing occurs in a stable region, in which internal modes have been excited, the validity of this assumption is questionable. Thus it is fair to say that further work is needed to resolve this subject of mixing beyond the convective layers.

Mazzitelli, D’Antona, & Ventura (1999) have examined the effects of overshooting predicted by full spectrum turbulence model of convection (Canuto & Mazzitelli 1991, 1992) within the context of lithium production in AGB stars. This is an interesting example of the way stellar convection may be tested. A particular aspect of their work of interest here is the suggestion that symmetric overshoot may be constrained by existing observations. We note that a general feature of numerical simmulations is up-down asymmetry. Cooling flows are narrow, faster downdrafts than the corresponding updrafts, which are broader and slower. It is suggestive that in simulations of envelope convection driven by radiative cooling at the top, there is an overshoot at the bottom (Hurlburt, Toomre, & Massaguer (1986); Freytag, Ludwig, & Steffen (1996)), while we, who have simulated shell convection driven by burning at the bottom, have found exponential overshoot at the top (and the bottom). It may be that the causes of this asymmetry can be determined, so that a complete algorithm for the overshoot maybe devised for 1D stellar evolutionary calculations.

In summary, we found that stellar evolution models should take into account penetration and an exponential decay of the velocities, and the excitation of internal modes, when simulating convective burning shells. Direct hydrodynamic simulation of stellar evolution places heavy demands on computer resources, especially for carbon and neon burning stages, which are slower. The tendency toward steady state revealed in our simulations is encouraging in that simplified algorithms may allow the construction of an improved set of one dimensional models.

We used an initial model which was derived from results of a one dimensional stellar evolution code, in which earlear stages of evolution were modeled with standard MLT. Consequently, the profile we started with is unsatisfactory, and is not consistent with our conclusion that penetration and extra mixing should be taken into account. This should be kept in mind when using directly the results of our simulations.

Enlightening discussions with Grant Bazàn, Eli Livne and Falk Herwig, and the interaction with the anonymous referee, are gratefully acknowledged. This work was supported in part by DOE grant DE-FG03-98DP00214/A001.