text:
# Nearly ideal binary communication in squeezed channels

## I Introduction

The ultimate capacity of a communication network is essentially quantum-limited, and the main concern of quantum communication is how to discriminate among quantum states that encode the relevant information . Quantum coding states are generally nonorthogonal, such that they cannot be unambiguously discriminated. As a consequence, the detection strategy should be optimized at the receiving side, in order to maximize the detection probability and/or minimize the transmission errors.

The scheme we have in mind is the following: a binary alphabet $`𝒜=\{0,1\}`$ with equal a priori probability symbols is being transmitted through a quantum communication channel. The information is encoded in two arbitrary Gaussian quantum states $`\varrho _0`$ and $`\varrho _1`$. In the following we first consider the case of pure states $`|\psi _0`$ and $`|\psi _1`$, whereas, in the second part of this letter, the analysis will be extended to the mixed-state case. Information is amplitude-keyed encoded , such that the wave functions of the two states are given by

$`\psi _0(x)=x|\psi _0`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi \sigma ^2}}}\mathrm{exp}\left[{\displaystyle \frac{x^2}{2\sigma ^2}}+if_0(x)\right]`$ (1)
$`\psi _1(x)=x|\psi _1`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi \sigma ^2}}}\mathrm{exp}\left[{\displaystyle \frac{(xa)^2}{2\sigma ^2}}+if_1(x)\right],`$ (2)
where $`f_j(x),j=1,2`$ are arbitrary phases, and $`a^+`$. Since the two states have the same a priori probability of being transmitted, the mean total energy traveling through the channel is given by $`E_T=a^2/2+(\sigma ^21/2)^2/\sigma ^2`$ (measured in unit of $`h\tau `$, $`\tau `$ being the characteristic time of the physical channel, e.g. the period for a bounded system, or the time-length of the wave-packet for a free system). The case $`\sigma ^2=1/2`$ corresponds to customary on-off coherent modulation, whereas for $`\sigma ^2<1/2`$ we are dealing with squeezed states . Although squeezing increases the total energy introduced into the channel, we will show that it can be effectively employed to improve the communication scheme, and to approach the performances of an ideal channel.

At the receiver, we consider the standard detection of the signal observable $`X`$, $`\widehat{\mu }(x)=|xx|`$, such that the output probability densities are given by $`p_0(x)=|x|\psi _0|^2=G(x;0,\sigma )`$ and $`p_1(x;a)=|x|\psi _1|^2=G(x;a,\sigma )`$, where $`G(x;a,\sigma )=(2\pi \sigma ^2)^{1/2}\mathrm{exp}[(xa)^2/2\sigma ^2]`$ is a normalized Gaussian of mean $`a`$ and variance $`\sigma ^2`$. On the basis of each measurement outcome we have to discriminate between two hypothesis: the null hypothesis $`_0`$ corresponding to the transmission of $`|\psi _0`$ (no signal), and the alternative hypothesis $`_1`$, corresponding to the transmission of $`|\psi _1`$, i. e. to the presence of the signal. The process of measurement and inference is called a decision strategy. We denote by $`Q_1`$ the power of the strategy, that is the probability of inferring the alternative hypothesis when the signal is actually present (also called the detection probability), and by $`Q_0`$ the size of the strategy, i. e. the probability of inferring the alternative hypothesis when the null hypothesis is true (also called the false-alarm probability).

In the following we employ a threshold strategy, in which the alternative hypothesis is chosen if the outcome is greater than a threshold value $`x_0`$. In order to determine the threshold value we should optimize the strategy, a goal that, in turn, requires to adopt an optimization criterion. Usually, one uses the criterion of minimizing the average cost of the decision, that is, in Bayesian terms, that of minimizing the probability of a wrong inference . Alternatively, one may accept to occasionally obtain an inconclusive inference in order to achieve error-free discrimination . Actually, these have been fruitful approaches in quantum state recognition, especially in the $`M`$-ary decision problem . However, the price of a small error probability is usually a small detection probability too, which, in turn, may imply the requirement of a high repetition rate. On the other hand, in the field of communication there exist several protocols that are robust , i.e. that may satisfactorily work also with a nonzero transmission-error rate. In this case, the main interest is that of maximizing the detection probability $`Q_1`$, while maintaining the size $`Q_0`$ to a moderated tolerable level. A decision strategy which is optimized according to such a criterion, which we will employ throughout this letter, is said to be a Neyman-Pearson (NP) optimized strategy .

## II Nearly ideal performance of a squeezed channel

The optimal NP threshold strategy for the present $`X`$-measurement is given in term of a density $`\mathrm{\Pi }(x)`$, which represents the probability of choosing the alternative hypothesis after having observed the outcome $`x`$. We have (Neyman-Pearson Lemma)

$`\mathrm{\Pi }(x)=\{\begin{array}{ccc}1& \mathrm{if}& \mathrm{\Lambda }(x)e^\kappa \\ 0& \mathrm{if}& \mathrm{\Lambda }(x)<e^\kappa \end{array}`$ (5)
where $`\mathrm{\Lambda }(x)=p_1(x;a)/p_0(x)`$ is the likelihood ratio, and $`\kappa `$ is the decision level. By varying the decision level we obtain NP strategies with different sizes. The likelihood ratio is given by $`\mathrm{\Lambda }(x)=\mathrm{exp}[(a^22ax)/2\sigma ^2]`$, and the NP strategy of Eq.(5) can be summarized as follows: the alternative hypothesis $`_1`$ is chosen if the outcome is greater than the threshold value $`x_0=(a^2+2\sigma ^2\kappa )/2a`$. The corresponding size and power are given by $`Q_0=_{}𝑑x\mathrm{\Pi }(x)p_0(x)`$ and $`Q_1=_{}𝑑x\mathrm{\Pi }(x)p_1(x;a)`$ i. e.

$`Q_0`$ $`=`$ $`{\displaystyle _{x_0}^{\mathrm{}}}𝑑xp_0(x)={\displaystyle \frac{1}{2}}\left[1\mathrm{Erf}\left({\displaystyle \frac{x_0}{\sigma \sqrt{2}}}\right)\right]`$ (6)
$`Q_1`$ $`=`$ $`{\displaystyle _{x_0}^{\mathrm{}}}𝑑xp_1(x;a)={\displaystyle \frac{1}{2}}\left[1\mathrm{Erf}\left({\displaystyle \frac{x_0a}{\sigma \sqrt{2}}}\right)\right],`$ (7)
By eliminating the decision level $`\kappa `$ between Eqs. (6) and (7) one obtains the characteristics $`Q_1(Q_0)`$

$`Q_1={\displaystyle \frac{1}{2}}\left\{1\mathrm{Erf}\left[\mathrm{InvErf}\left(12Q_0\right){\displaystyle \frac{a}{\sqrt{2}\sigma }}\right]\right\}.`$ (8)
Since the error function $`\mathrm{Erf}(x)`$ and its inverse $`\mathrm{InvErf}(x)`$ are monotone, the power at fixed size increases with the term $`a/\sqrt{2}\sigma `$. As we will see, this quantity may be enhanced by squeezing, such that for any energy $`E_T`$ squeezed channels always show larger power than coherent ones.

The value of $`E_T`$ is set by physical constraints, and a question arises about the optimal fraction of $`E_T`$ that should be employed in squeezing the channel. In fact, the energy cannot be entirely spent in squeezing, since in this case no signal amplitude is left to be discriminated. Let us define the squeezing fraction $`\gamma `$ as the fraction of the total energy $`E_T`$ that is employed to squeeze the channel. In terms of $`\gamma `$ and $`E_T`$ the amplitude and the squeezing are given by $`a=\sqrt{2E_T(1\gamma )}`$ and $`\sigma =1/2(\sqrt{\gamma E_T+2}\sqrt{\gamma E_T})`$. Using these expressions we have

$`{\displaystyle \frac{a}{\sqrt{2}\sigma }}={\displaystyle \frac{2\sqrt{E_T(1\gamma )}}{\sqrt{\gamma E_T+2}\sqrt{\gamma E_T}}}.`$ (9)
The maximum value is $`\left(a/\sqrt{2}\sigma \right)_{max}=\sqrt{E_T(E_T+2)}`$, which is reached for

$`\gamma _{\mathrm{opt}}={\displaystyle \frac{1}{2}}{\displaystyle \frac{E_T}{1+E_T}}.`$ (10)
$`\gamma _{\mathrm{opt}}`$ thus represents the optimal squeezing fraction to discriminate, according to NP criterion, amplitude-keyed signals by $`X`$ measurement. In Fig. 2 we show the characteristics $`Q_1(Q_0)`$ for optimally squeezed and coherent channels with different energies. The improvement due to squeezing is apparent.

We also notice that $`Q_1`$ is a smooth function of the squeezing fraction, which, in turn, should not be considered as a critical parameter. In facts, in order to obtain an enhancement of the strategy power, we do not need a fine tuning of $`\gamma `$. This is illustrated in Fig. 3, where a contour plot of $`Q_1`$ is shown as a function of $`Q_0`$ and $`\gamma `$. For fixed size $`Q_0`$ the power slowly varies with $`\gamma `$ in a considerably large range of values.

For a given size $`Q_0`$ the bound $`Q_1=1/2`$ defines the minimum detectable signal. As it follows from Eq. (8), this corresponds to $`a/\sqrt{2}\sigma =\mathrm{InvErf}(12Q_0)`$, and using Eq. (9) to

$`E_T^{min}={\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{InvErf}^2(12Q_0)}{1\gamma +\mathrm{InvErf}^2(12Q_0)\sqrt{2\gamma (1\gamma )}}}.`$ (11)
$`E_T^{min}`$ decreases with $`\gamma `$, i.e. squeezed channels allow one to discriminate weaker signals at given size. For small $`Q_0`$, $`E_T^{min}`$ increases quadratically for coherent channels $`E_T^{min}=y^2/2`$, and only linearly for optimally squeezed one $`E_T^{min}=1+\sqrt{1+y^2}`$, $`y`$ being the principal solution of the equation $`y\sqrt{\pi }Q_0=\mathrm{exp}(y^2)`$.

In order better to appreciate the benefit of squeezing, we compare the power $`Q_1`$ of the NP $`X`$-threshold strategy (5), with the optimal NP quantum measurement to discriminate between two pure states $`|\psi _0`$ and $`|\psi _1`$. Such an optimal measurement has been found long ago , whereas a comprehensive approach for mixed states is still lacking . For a pair of pure states the optimized measurement is given by

$`\widehat{\mu }(x|\lambda )=|\psi _1\psi _1|\lambda |\psi _0\psi _0|,`$ (12)
where $`\lambda `$ is a Lagrange multiplier which determines the decision level. The decision strategy consists in choosing the alternative hypothesis $`_1`$ for positive outcomes, and the resulting detection probability reads as follows

$`Q_1=\{\begin{array}{cc}\left[\sqrt{Q_0\omega }+\sqrt{(1Q_0)(1\omega )}\right]^2\hfill & 0Q_0\omega \hfill \\ 1\hfill & \omega Q_01\hfill \end{array},`$ (15)
where $`\omega =|\psi _0|\psi _1|^2`$ is the overlap between the two states. Notice that the ideal NP measurement has been considered to find the ultimate quantum limit to high-precision binary interferometry .

In Fig. 4 we show the power-size characteristics of the optimal strategy in comparison with that of the $`X`$-strategy for coherent and optimally squeezed channels. For squeezed channels the power increases, and approaches the optimal value already for moderate energy.

In order to summarize improvements due to squeezing we consider the mutual information between input and output $`I=_{ij}P_{ij}p_j\mathrm{log}\left[P_{ij}/(_jP_{ij}p_j)\right]`$, where $`p_0=p_1=1/2`$ are the a priori probabilities of the two symbols, and $`P_{ij}`$ is the probability of choosing hypothesis $`_i`$ when hypothesis $`_j`$ is true. In our case, $`P_{11}=Q_1`$ and $`P_{10}=Q_0`$, such that $`P_{01}=1Q_1`$ and $`P_{00}=1Q_0`$. On the left panel of Fig. 5 we show the mutual information $`I_X`$ of the $`X`$-strategy as a function of the total energy $`E_T`$ for an optimal choice of the squeezing fraction $`\gamma `$. $`I_X`$ saturates to high value already for moderate energy, showing only a weak dependence on the size of the strategy. On the right panel we show the ratio (in dB) between $`I_X`$ and the ideal value $`I_{opt}`$, corresponding to the optimal NP strategy. It results that for a squeezed channel the mutual information is approaching the ideal value for much lower energy than a coherent one. Similar plots are obtained varying the size of the strategies.

## III Fuzzy hypothesis testing and mixed Signals

So far we have considered information amplitude-keyed on pure states. However, in practice, it is more likely to deal with mixture, either because the coding stage is imperfect, or as a result of noises in the transmitter and losses in the channel. For the sake of simplicity, we consider a situation in which the null hypothesis still corresponds to coding onto the vacuum (no amplitude) state, that is $`\varrho _0=|\psi _0\psi _0|`$. On other hand, the alternative hypothesis now corresponds to coding the signal on the mixed state $`\varrho _1=𝑑bH_1(b)|\psi _b\psi _b|`$, where $`|\psi _b`$ coincides with $`|\psi _1`$ of Eq. (2) and $`H_1(b)`$ is a weight function, which will be taken of Gaussian form. The two hypothesis to be discriminated are no longer crisp, and the decision problem should be formulated in the framework of fuzzy hypothesis testing . The fuzzy null and alternative hypothesis are formulated as follows: $`_j`$ is true when a Gaussian state of amplitude $`b`$, distributed as $`H_j(b)`$, is transmitted. In our case the two membership density functions are given by $`H_0(b)=\delta (b)`$ and $`H_1(b)=(2\pi \mathrm{\Sigma }^2)^{1/2}\mathrm{exp}[(ba)^2/2\mathrm{\Sigma }^2]`$.

In order to analyze the effect of squeezing with a mixed signal we need to find the best NP $`X`$-strategy of its size. Recently, the Neyman-Pearson Lemma has been extended to to fuzzy hypothesis testing , and this allows us to solve the decision problem. The NP strategy for mixed states is a density of the form (5) with the fuzzy likelihood ratio given by

$`\mathrm{\Lambda }^\mathrm{f}(x)`$ $`=`$ $`{\displaystyle \frac{𝑑bp_1(x;b)H_1(b)}{𝑑bp_0(x)H_0(b)}}={\displaystyle \frac{𝑑bG(x;b,\sigma )H_1(b)}{G(x;0,\sigma )}}`$ (16)
$`=`$ $`\sqrt{{\displaystyle \frac{\beta ^2}{1+\beta ^2}}}\mathrm{exp}\left[{\displaystyle \frac{x^2+2ax\beta ^2a^2\beta ^2}{2\sigma ^2(1+\beta ^2)}}\right],`$ (17)
where $`\beta ^2=\sigma ^2/\mathrm{\Sigma }^2`$. The pure-state case is obtained in the limit $`\beta \mathrm{}`$.

The power-size characteristics has the same functional form (8) of the pure state case. However, for mixed signals part of the energy is degraded to noise, $`E_T=(a^2+\mathrm{\Sigma }^2)/2+(\sigma ^21/2)^2/\sigma ^2`$, such that the amplitude reads as follows $`a=\sqrt{2E_T(1\gamma )\mathrm{\Sigma }^2}`$. After inserting this expression into the term $`a/\sqrt{2}\sigma `$, and maximizing over $`\gamma `$ one obtains the optimal squeezing fraction for mixed channels

$`\gamma _{\mathrm{opt}}^\mathrm{m}={\displaystyle \frac{(2E_T\mathrm{\Sigma }^2)^2}{8E_T(1+E_T\mathrm{\Sigma }^2/2)}}.`$ (18)
As it can be easily proved from Eq. (18) $`\gamma _{\mathrm{opt}}^\mathrm{m}`$ is always smaller then $`\gamma _{\mathrm{opt}}`$ for any value of $`E_T`$, i.e. a smaller amount of squeezing can be employed in a mixed channel against a pure channel with the same energy. Correspondingly, also the power at fixed size decreases. However, squeezing a mixed channel is still extremely convenient to improve the $`X`$ strategy by comparison with a mixed coherent channel of the same energy. In order to illustrate this behavior, we define the ratio $`R=I_X^\mathrm{s}/I_X^\mathrm{c}`$ (at fixed energy and mixing parameter) between the mutual information of a squeezed and a coherent channel respectively.

In Fig. 6 we show $`R`$ as a function of $`E_T`$ for different values of the size for a weakly ($`\mathrm{\Sigma }=\sqrt{2E_T}/10`$) and a strongly ($`\mathrm{\Sigma }=2\sqrt{2E_T}/3`$) mixed channel. Notice that $`\mathrm{\Sigma }=\sqrt{2E_T}`$ is the limiting value, corresponding to a completely mixed signal with no amplitude and no squeezing. $`R`$ linearly increases for small $`E_T`$ and after a maximum of few dB (the actual height depends on the size $`Q_0`$) decreases. In the (unrealistic) limit of very high energies squeezing the channel is no longer convenient. We notice that for a strongly mixed channel such a decreasing is much slower, thus indicating that squeezing is effective in a wide range of energies. In other words, a squeezed channel is more robust against mixing of signals than a coherent one.

## IV Conclusions

In conclusion, we have shown that squeezing the channel in amplitude-keyed binary communication increases the detection probability at fixed size. We have found the optimal squeezing fraction and evaluated the mutual information for both pure and mixed signals. Optimally squeezed channels are robust against signal-mixing, and largely improve the strategy power by comparison with coherent ones, approaching the performance of the ideal receiver already for moderate signal energy.

## Acknowledgment

This work has been cosponsored by C.N.R. and N.A.T.O through grant No. 217.32 of the Advanced Fellowship Program. The author thanks people at SESAMO lab of the University of Modena and Reggio Emilia for their patience and kind hospitality during the last stage of this work.


format: text

text:
# Structure of aluminum atomic chains

## I introduction

The fabrication of the stable gold monoatomic chains suspended between two gold electrodes is one of the milestones in nanoscience. Issues brought about by this achievement are yet to be resolved: Stable chains were obtained by stretching gold nanowires; no other metal, such as Al, Cu, has been observed to form a stable monoatomic chain yet. The monoatomic chain, being an ultimate one-dimensional (1D) structure, has been a testing ground for the theories and concepts developed earlier for three-dimensional (3D) systems. For example, it is of fundamental importance to know the atomic structure in a truly 1D nanowire and how the mechanical and electronic properties change in the lower dimensionality.

The density functional theory has been successful in predicting electronic and mechanical properties of bulk metals, where each atom has 8-12 nearest neighbors depending on the crystal structure. While many neighbors in a 3D structure is a signature of the formation of metallic bonds, it is not obvious whether the “metallic” bond picture will be maintained in a monoatomic chain. In fact, for a monoatomic linear chain with one electron per atom, the dimerized state is more stable with a Peierls gap at the zone edge. The situation is expected to be more complex for the chain of aluminum atoms having $`3s^23p^1`$ valency.

The interest in metal nanowires is heightened by the observation of quantized behavior of electrical conductance at room temperature through connective necks stretching between two electrodes. Studies attempting to simulate the process of stretching by using classical molecular dynamics have shown novel atomic and mechanical properties. In particular, it was found that the 2D hexagonal or square lattice structure of atomic planes perpendicular to the axis changes to the pentagons and later to equilateral triangles when the wire is thinned down to the radius of 5-10 Å. Upon further thinning, strands (or bundles of finite atomic chains), and eventually a monoatomic chain forms at the narrowest section of the nanowire. Recently, the stability of suspended gold chains and their atomic structures have been studied extensively.

The first-principles calculations by Portal et al. showed that infinite, as well as finite gold atomic chains between two gold electrodes favor the planar zigzag geometry at a bond angle $`\alpha =`$131o. The homogenization of the charge with a depletion in the interatomic region ruled out the formation of a directional chemical bond. On the other hand, the first-principles calculations by Häkkinen et al. for a finite gold chain between two gold electrodes favored the dimerized structure. In contrast to the conclusion drawn by Portal et al., Häkkinen et al. attributed the stability of the suspended gold chain to the directional local bonding with $`spd`$ hybridization. Apparently, the stability of a finite chain depends on the strain and the atomic configuration where the chain is connected to the electrodes. In a more recent comparative study Au, Cu, Ca, K infinite chains were found to form planar zigzag structures with equilateral triangular geometry; only Au chain has a second zigzag structure with a wide bond angle $`\alpha `$=131o. Note that these atoms can be considered similar because of their $`s`$-type outermost valence orbitals. Aluminum with $`3s`$-, and $`3p`$-valence orbitals is different from Au, Cu, Ca. K. Therefore, Al is an important element for understanding the formation and stability of ultimate 1D atomic chains.

This paper presents a systematic, first-principle analysis of the binding, atomic and electronic structure of very thin Al chains. The objective is to reveal periodic linear, planar and non-planar geometries forming stable structures. An emphasis is placed on the planar structures forming zigzag chains. Some of the Al chain structures are compared with the corresponding structures of Au chains. It is found that by going from bulk Al to a chain structure the character of bonding changes and acquires directionality. The higher the coordination of individual atoms, the stronger is the binding energy. It is hoped that the present analysis will contribute to the understanding of atomic structure and related physical properties (such as electrical and thermal conductance, elascticity etc.) of infinite and finite atomic chains.

## II method

First-principle calculations were carried out within the density functional theory. Al and Au chains are treated within the supercell geometry. To minimize the interchain interaction the distance between the

chains is taken to be 20 Å. The wave functions are expressed by plane waves with the cutoff energy, $`|𝐤+𝐆|^2`$275 eV. The Brillouin zone (BZ) integration is performed within Monkhorst-Pack scheme using (1$`\times `$1$`\times `$40) $`𝐤`$-points. The convergence with respect to the energy cutoff and number of k-points were tested. Ionic potentials are represented by ultra soft Vanderbilt type pseudopotentials and results are obtained by Generalized Gradient Approximation for fully relaxed atomic structures. Preconditioned conjugate gradient method is used for wave function optimization. Since ionic relaxations are carried out by the conjugate gradient method, the optimized (fully relaxed) structures obtained in this study are stable structures. In certain cases the stability of a structure is tested by calculating the total energy while the atoms are displaced in special directions. Numerical calculations are performed by using VASP code. The $`z`$-axis is taken along the chain axis, and $`y`$-axis ($`x`$-axis) is perpendicular to (in) the plane of zigzag structure.

## III Results

### A Optimized structures and cohesive energies

The variation the total energy $`E_T`$, of the atomic Al chain calculated for the fully relaxed linear, planar (zigzag and ladder), and non-planar (cross) structures is shown in Fig.1. The geometries of these structures and their relevant structural parameters are shown by insets. Since the total energies are given with respect to the energy of the free Al atom, the cohesive energy $`E_C=E_T`$. The zigzag geometry displays two minima; one occurs at $`s=`$1.26 Å and has cohesive energy $`E_C=`$2.65 eV/atom; other has shallow minimum and occurs at $`s=`$2.37 Å with cohesive energy $`E_C=`$1.92 eV/atom. The high cohesive energy zigzag structure (specified as $`T`$) having the bond length $`d=`$2.51 Å, and the bond angle $`\alpha `$60o forms equilateral triangles. This geometry allows for four nearest neighbors, which is less than the six nearest neighbors occurring in the Al(111) atomic plane and twelve nearest neighbors in the close packed bulk metal. The equilateral triangular geometry can also be viewed as if two parallel linear chains with an interchain distance of 2.17 Å are displaced by $`d/2`$ along the chain axis ($`z`$-direction). This is reminiscent of the hollow site registry of 2D atomic planes which usually increases the cohesive energy.

The low cohesive energy zigzag structure (specified as $`W`$) has $`d=`$2.53 Å and wide bond angle $`\alpha `$139o, and allows for only two nearest neighbors with bonds slightly larger than those of the $`T`$-structure. We also found that the cohesive energy decreases if an Al atom is displaced perpendicular to the zigzag plane. Therefore, both zigzag structures are planar. The minimum energy of the linear structure ($`\alpha =`$ 180o and denoted as $`L`$) has relatively short bond length, $`d=s=`$2.41 Å. It is $``$0.5 eV above the minimum energy of the $`W`$-structure and has cohesive energy $`E_C=`$1.87 eV/atom.

Two linear chains can form a ladder structure which allows for three nearest neighbors with $`\alpha =90^o`$ and $`E_C`$2.4 eV/atom intermediate to the $`T`$-, and $`W`$-structures. The cohesive energy is further increased to 2.5 eV/atom when the separation between chains is sligtly increasesd. This way two strands (specified as $`S`$-structure) form, which are held in place by the uniaxial stress between two electrodes.

Non-planar cross structure (specified as $`C`$-structure) has four atoms which form two perpendicular dumbbells (A and B) in the unit cell. The lengths of these dumbells are different ( A: 2.8 Å and B: 4.15 Å) and the chain is made by the ABABA… sequence of these dumbbells. Al atoms in A has five non-planar bonds, and those in B have four bonds of $``$2.8 Å. The cohesive energy of this structure is calculated to be 3.04 eV/atom. Since the atoms of the A-dumbbells are bound to the nearest five atoms forming equilateral triangles in different planes, and those of the B-dumbbels have four non-planar bonds, this cohesive energy is highest among the 1D structures described in Fig. 1. The $`C`$-structure was revealed first in the extensive analysis of Gülseren et al. by using empirical glue potential. The overall features of the stable $`C`$-structure determined by these calculations are confirmed here, but the structural parameters are more accurately determined by the present first principle calculations.

It is worth noting that the Au atomic chain also forms two different zigzag structures similar to those of Al; but Cu, Ca and K do not. Our calculated values for $`s`$, $`d`$, $`E_C`$ are respectively, 1.36 Å, 2.71 Å, 2.23 eV/atom for the $`T`$-structure; 2.33 Å, 2.56 Å, 1.90 eV/atom for the $`W`$-structure; 2.59 Å, 2.59 Å, 1.68 eV/atom for the $`L`$-structure of Au chain. The nearest neighbor distance of the $`T`$-structure of the Au chain is reduced only $``$ 6 $`\%`$ from that of the bulk. Is this puzzling similarity of 1D atomic structures of Al and Au (despite their dissimilar valencies) only a coincidence? We now address this issue.

We calculated the cohesive energy of bulk Al (Au), $`E_C=`$3.67 (3.20) eV/atom at the nearest neighbor distance $`d=`$2.86 (2.96) Å (or lattice parameter $`a=`$4.04 (4.18) Å). The energetics of 1D and bulk structures are compared in Table I. Simple arguments based on the counting of nearest neighbor couplings would suggest a relatively small cohesive energy, e.g., $``$ 1.3 eV for the $`T`$-structure. On the contrary, 1D structures studied here have cohesive energies higher than one can estimate by comparing their coordination numbers with that of bulk. Apparently, the bonds in 1D structures become stronger. In fact, it was found previously that the linear Al chain has a Young’s modulus stronger than bulk. Recent scanning tunneling microscope studies revealed that the bond strength of the Au nanowire is about twice that of a bulk metallic bond.

### B Charge density analysis

Figure 2 shows the charge density contour plots of bulk, $`L`$-, $`W`$-, and $`T`$-structures. In contrast to uniform metallic charge density of bulk, the bonding acquires directionality in 1D structures of Al. For the $`L`$-structure the charge is accumulated between atoms forming a directional bond, and is mainly due to the $`\sigma `$ states (formed by $`3s+3p_z`$ orbitals) and partly due to $`\pi `$ states (formed by $`3p_x`$ and $`3p_y`$ orbitals perpendicular to the chain axis). The calculated charge distribution suggests that directional “covalent” bonds are responsible for the bonding. This situation is maintained in the zigzag $`W`$-structure except for a slight distortion of the bond charge. Actually, the $`W`$-structure with wide bond angle is not dramatically different from the $`L`$-structure. In the $`T`$-structure that forms equilateral triangles, the charge density is apparently different from that of the $`W`$-structure. We see a continuous (connected) region of high charge density between double atomic chains. However, this is nothing but the overlap of charges of four bonds emerging from each chain atom, and is confirmed by the contour plot of an individual bond charge in a plane perpendicular to the zigzag ($`xz`$) plane and passing through an Al–Al bond. We also notice that the charge becomes slightly delocalized by going from $`L`$ to $`T`$-structure. These charge distributions of the Al–Al bond described above is different from the corresponding charge distribution of Au zigzag structures shown in Fig. 3. Clearly, there are no directional bonds in the Au chain; valence charge is delocalized. This finding is in confirmity with the results of Portal et al. on the infinite Au chain. On the other hand, Häkkinen et al. deduced directional bonding with $`spd`$ hybridization in finite Au chains between two Au electrodes by performing similar type of pseudopotential plane wave calculations.

### C Electronic structure

A comparative analysis of the electronic band structure of Al monoatomic chains illustrated in Fig. 4 provides further insight into the stability and character of bonding. The band structure of the $`L`$-structure is folded for the sake of comparison with the zigzag structures. Two filled $`\sigma `$ bands arise from the $`3s+3p_z`$ valence orbitals and make the bond charge shown in Fig. 2a. Because of the linear geometry $`3p_x`$ and $`3p_y`$ are equivalent, and give rise to doubly degenerate $`\pi `$ band crossing the Fermi level. As pointed out by Peierls, a one-dimensional metal with a partly filled band will distort away from a regular chain structure to lower its energy. According to the above analysis, a linear chain of uniformly spaced Al atoms with spacing $`s=d`$ has a quarter-filled band which crosses the Fermi level at $`k_z=\pm \pi /4d`$. A distorted unit cell $`4d`$ in length will cause this point to coincide with the edge of the Fermi distribution. The change in the crystal potential due to the $`4d`$ distortion will open up a Peierls gap at the reduced zone-edge lowering the total energy. In practice, the gain in energy due to such Peierls distortions is rather small even for a $`2d`$ dostortion (dimerization) and is likely to be below computational error for the $`4d`$ distortion here. Thus, although the linear chain of Al atoms is unstable, in principle, the effect of such a distortion on cohesive energy is clearly negligible.

The symmetry between $`3p_x`$ and $`3p_y`$ orbitals is broken in the zigzag structure, and hence the $`\pi `$-band is split. Apart from this band splitting and slight rise of bands, the overall form of the energy band structure is maintained in the $`W`$-structure. The $`W`$-structure is, however, more stable than the $`L`$-structure because of its relatively stronger electronic screening. In the $`T`$-structure the split $`\pi `$-bands are lowered, and the form of the $`\sigma `$-bands undergo a significant change due to the equilateral triangular geometry. Despite slight delocalization of charge, the total energy of the $`T`$-structure is lower than the $`W`$-structure. The relative stability originates from the increased number of nearest neighbors. It is noted that the bands in all chain structures can be considered similar as far as Fermi level crossing of the bands is concerned.

In contrast to the Al chains described above, the energy band structure of Au undergoes significant changes in different structures near the Fermi energy. For example, for the linear structure one band crosses the Fermi level near the $`X`$-point of the BZ, another band at the $`\mathrm{\Gamma }`$-point is very close to the Fermi level. For the $`W`$-structure, two bands cross at the Fermi level and at the zone boundary with negligible Peierls distortion gap and the rest of the bands are lowered. The lowering of the state density at the Fermi level stabilizes the zigzag structure relative to the linear structure. In the $`T`$-structure two bands cross the Fermi level. Inspite of these changes the character of the bonding remains essentially metallic in the chain structures.

## IV Discussion and conclusions

The analysis of the above results and comparison with Au chain reveal that the metallic bond of bulk Al changes to a directional covalent bond in the 1D monoatomic chain. The metallicity is ensured by the $`\pi `$-states. For that reason, our efforts of calculating the bands of the Al chains with tight binding method by using the bulk parameters have not been too successful. This suggests that the transferability of energy parameters fitted to bulk is not satisfactory for the 1D ($`T`$-, $`W`$\- and $`L`$-) structures.

The $`T`$ structure with two parallel linear chains can be viewd as the 1D analog of the 3D close packing. In this respect, the $`T`$-structure may be considered in a different class and as a precursor of the 2D hexagonal lattice. Adding one more parallel chain in registry with the quasi 1D $`T`$-structure, one starts to build the hexagons, where 2/3 of the atoms have four and 1/3 of the atoms have six nearest neighbors. As a natural extention of these arguments, another intermediary, quasi 1D structure, for example, is a ladder structure which consists of two parallel linear chains forming a row of squares with a lattice constant of $`d`$ and allowing for three nearest neighbors. This metastable structure is a 1D analog of the top site registry of 2D atomic planes. Our calculations show that the cohesive energy of the ladder structure is increased when the distance between two chains increases, so that the chain turns to two strands. The cohesive energy of the strands is found between the $`T`$\- and $`W`$-structure. By going from planar to non-planar geometry the cohesive energy further increases The present work suggests that the Born-Oppenheimer surface for these quasi 1D structures is rather complex, and generally $`E_C`$ increases with increasing coordination number and decreasing bond angle. As clarified in Sec. II, the stable structures correspond to the local minima on the Born-Oppenheimer surface and are expected to be vibrationally stable at least at low temperatures. The 1D $`T`$-structure found for Al, Au, Cu, Ca and K appears to be common to metals, in a way an intermediate structure between a truly 1D and 2D structures.

We also note that the linear and zigzag structures of Al have two bands crossing the Fermi energy. Calculations using the Green’s function method yield one conduction channel for each band of uniform chain crossing the Fermi level, and hence the ballistic conductance of the $`T`$-, $`W`$-, and $`L`$-structures $`G=2(2e^2/h)`$. This value for the conductance arises from the fact that the channel capacity or the maximum conductance per channel is $`2e^2/h`$. It is straightforward to motivate the maximum value by appealing to the Heisenberg’s uncertainity principle. Recalling that conductance $`G=\mathrm{\Delta }I/\mathrm{\Delta }V`$, and $`\mathrm{\Delta }I=\mathrm{\Delta }Q/\mathrm{\Delta }t`$, then for a single channel in extreme quantum limit $`\mathrm{\Delta }Q=e`$. One can readily write $`G=e^2/\mathrm{\Delta }E\mathrm{\Delta }t`$. Now invoking the uncertainity principle, $`\mathrm{\Delta }E\mathrm{\Delta }th`$, one finally obtains $`G2e^2/h`$. Here the factor of two is due to spin. The maximum conductance per channel can never be greater than $`2e^2/h`$. We note that the value of the ballistic conductance for the infinite chain (i.e. $`4e^2/h`$) is at variance with the experimental results yielding only $`G2e^2/h`$ for the finite Al chain. The discrepancy was explained by the fact that the electronic states are modified due to the finite size of the chain and the atomic configuration where the chain is coupled to the electrodes.

It is important to remark that creation of metallic overlayers on semiconductors is required in chip technology. Thus if one could place these quasi 1D structures on semiconductor surfaces without losing the metallic behavior of the chains, the technology would be considerably enhanced. Unfortunately, when one examines the deposition of monolayers of metals like Al, Au and Ga on Si, the lowest energy configuration turns out to be semiconducting in nature. A metastable state, in which Al forms a metallic zigzag structure on the Si(100) surface, has been reported. It remains to be seen how feasible it is to fabricate such a structure.

In summary, we have found that a zigzag chain of aluminum in triangular configuration is most stable among the planar structures we studied. The structural results for planar geometries are similar to gold but bonding is different. A new metastable ladder structure intermediate to distorted linear and triangular structure is also reported. The metallicity has its origin in the $`\pi `$-bands. The stabilization of these metallic monoatomic chains on semiconductors remain an experimental challenge.

## V acknowledgement

We thank Dr Oguz Gülseren for his stimulating discussions.


text:
# 1 Signal cross section versus bottom squark mass for 𝜆''₃₁₃=1.

## Acknowledgments

This work was supported in part by DOE grant No. DE-FG02-94ER40817.


text:
# A New Code for Nonradial Stellar Pulsations and its Application to Low - Mass, Helium White Dwarfs

## 1 INTRODUCTION

At present, nonradial stellar pulsations constitute a valuable tool in extracting information about the properties and internal structure of objects that undergo such kind of oscillations. As observational techniques are refined, a growing variety of stellar objects is discovered that undergo nonradial oscillations. The most extensively studied case is our Sun, which has been shown to vibrate in a large amount of modes, investigation of which has been so powerful in revealing clues of its internal structure that it is known as the discipline of Helioseismology. Also, it has been established that different kinds of stellar objects which are located in a variety of places in the HR diagram, covering several evolutionary stages, undergo nonradial pulsations, e.g.: Ap, WR, $`\delta `$ Scuti, $`\beta `$ Cephei, and variable white dwarfs (hereafter WDs): DAV, DBV, DOV and PNNV. The study of oscillations in such variable stars, and the subsequent comparison and fitting with theoretical pulsation patterns is now known as Asteroseismology. For a general introduction to this topic see, e.g., Cox (1980); Unno, et al. (1989); Brown & Gilliland (1994); and Gaustchy & Saio (1995; 1996).

From the observational point of view (with the obvious exception of our Sun), WDs represent the best established kind of nonradial pulsators. WDs are the final fate for low and intermediate mass stars. During its cooling, they go across very narrow instability strips inside which suffer nonradial pulsations in the branch of g- (gravity) modes with periods between 100 and 2000 sec and amplitudes up to 0.3 magnitudes.

As WDs cools down, the outermost layers main constituent (hydrogen in DAVs, helium in DBVs, and probably carbon or oxygen in DOVs and PNNVs; see e.g. Brown & Gilliland 1994 and Gaustchy & Saio 1995; 1996) gets partial ionization conditions which forces the development of an outer convection zone (OCZ). In the bottom of such zone the instabilities are originated ($`\kappa `$ mechanism of overstability). This gives rise to the blue edge of the instability strip (see, e.g., Dziembowski & Koester 1981 and Winget, et al. 1982; see also Tassoul, Fontaine & Winget 1990, hereafter TFW)111Recently Goldreich & Wu (1999) have proposed the “convective driving mechanism” as the responsible for overstability of g-modes in DA WDs.. The physical reason for the appearance of a red edge is still under debate.

At present, there is a general consensus that variable WDs are very interesting objects for asteroseismological studies. Their very simple internal structure allows us to predict theoretically the modes with a very high degree of detail and sophistication. Also, they have a very rich spectrum of periods which may give us information about stellar masses (by means of measuring the mean period spacing), core composition (by measuring the rate of change of the period $`\dot{P}`$), mass of the surface helium and hydrogen layers (if present) (by means of the departures from uniform period spacing and mode trapping), velocity of rotation and strength of the magnetic field (studying the fine structure of a multiplet), etc. (see, e.g., Bradley 1996). So, it is not surprising that in the last years DAVs and DBVs as well as DOVs have been the preferred target for the network called “Whole Earth Telescope” (WET). WET observations have been of an unprecedent quality which in some cases (see, e.g., Winget et al. 1991; Kawaler 1993; Bradley & Winget 1994) allowed to apply the powerful tools above mentioned.

In view of the growing importance of this line of research, in our Observatory we have begun the study of nonradial oscillations. To this end, we have developed a code to compute eigenfunctions and eigenfrequencies for spherical objects in the linear, adiabatic approximation. The main purpose of the present paper is to describe the numerical techniques employed in our program to search for and to compute the modes and its application to the case of a 0.3 $`M_{}`$ pure helium WD. This code is based on a modification of the general Newton - Raphson technique presented in Kippenhahn, Weigert & Hofmeister (1967) to solve the set of difference equations which represent the differential equations of linear, adiabatic, nonradial pulsations of spherical stars.

In order to test the code we have applied it to the case of polytropic spheres whose eigenfrequencies have been computed with very high accuracy (Christensen - Dalsgaard & Mullan 1994). Also, we have coupled our pulsation code to an stellar evolution code. We have done it in order to be able to follow the changes in the pulsational spectrum of a given model during its evolution with a number of steps as large as it would be needed by specific studies. In making such computation we selected an evolutionary time step short enough such that eigenvalues and eigenfunctions corresponding to a given model are used as a initial approximation to those of the next one. This is very useful for a relaxation code like this and is intended to work as it is the case of the Heneyey technique in stellar evolution.

The rest of the paper is organized as follows: Section 2 is devoted to present the differential equations of linear, adiabatic, nonradial pulsations we have to solve. In Section 3 a detailed description of the general Newton - Raphson technique we employ for computing the eigenmodes is presented. Section 4, presents the results for the case of polytropic spheres, and in Section 5 we describe the work we have done with a realistic helium WD model. Finally, in Section 6 summarize our results.

## 2 The Differential Equations of Nonradial Pulsations

The differential equations that govern linear, nonradial pulsations of spherical stars in the adiabatic approximation are (Unno, et al. 1989):

$$x\frac{dy_1}{dx}=\left(V_g3\right)y_1+\left[\frac{\mathrm{}(\mathrm{}+1)}{C_1\omega ^2}V_g\right]y_2+V_gy_3,$$
(1)
$$x\frac{dy_2}{dx}=\left(C_1\omega ^2A^{}\right)y_1+\left(A^{}U+1\right)y_2A^{}y_3,$$
(2)
$$x\frac{dy_3}{dx}=\left(1U\right)y_3+y_4,$$
(3)
$$x\frac{dy_4}{dx}=UA^{}y_1+UV_gy_2+\left[\mathrm{}(\mathrm{}+1)UV_g\right]y_3Uy_4.$$
(4)
$`x`$ is the independent variable defined as $`x=r/R`$ where $`r`$ is the radial coordinate and $`R`$ is the stellar radius. The dependent variables are defined as

$$y_1=\frac{\xi _r}{r},$$
(5)
$$y_2=\frac{1}{gr}\left(\frac{p^{^{}}}{\rho }+\mathrm{\Phi }^{^{}}\right),$$
(6)
$$y_3=\frac{1}{gr}\mathrm{\Phi }^{^{}},$$
(7)
and

$$y_4=\frac{1}{g}\frac{d\mathrm{\Phi }^{^{}}}{dr},$$
(8)
$$\omega ^2=\frac{\sigma ^2R^3}{GM_{}}.$$
(9)
$`\xi _r`$ represents the radial displacement of the fluid and $`p^{^{}}`$ and $`\mathrm{\Phi }^{^{}}`$ are the Eulerian variations of the equilibrium values of the pressure $`p`$ and the gravitational potential $`\mathrm{\Phi }`$ respectively. $`g`$ is the local acceleration of gravity and $`M_{}`$ is the stellar mass. $`\omega ^2`$ is the dimensionless square of the angular frequency of oscillation, $`\sigma ^2`$.

Among the coefficients of Eqs. (1-4) there appears the harmonic degree $`\mathrm{}`$, corresponding to the spherical harmonic $`Y_{\mathrm{},m}(\theta ,\varphi )`$ which describes the angular dependence of the oscillation pattern. The adimensional quantities $`V_g`$, $`U`$, $`C_1`$, and $`A^{}`$, inherent to the non - perturbed model, are defined as:

$$V_g=\frac{V}{\mathrm{\Gamma }_1}=\frac{1}{\mathrm{\Gamma }_1}\frac{d\mathrm{ln}p}{d\mathrm{ln}r}=\frac{gr}{c^2}=\frac{gr\rho }{\mathrm{\Gamma }_1p},$$
(10)
$$U=\frac{d\mathrm{ln}M_r}{d\mathrm{ln}r}=\frac{4\pi \rho r^3}{M_r},$$
(11)
$$C_1=\left(\frac{r}{R}\right)^3\frac{M_{}}{M_r},$$
(12)
$$A^{}=\frac{r}{g}N^2=r\left(\frac{1}{\mathrm{\Gamma }_1}\frac{d\mathrm{ln}p}{dr}\frac{d\mathrm{ln}\rho }{dr}\right).$$
(13)
Here $`\mathrm{\Gamma }_1`$ is the first adiabatic index, $`c^2`$ is the square of the local velocity of sound, $`M_r`$ is the mass contained in a (non - perturbed) sphere of radius $`r`$, and $`N^2`$ is the Brunt \- Väisälä frequency.

Eqs. (1-4) together with the adequate boundary conditions conform the eigenvalue problem to be solved.

## 3 The Numerical Code

In order to write the pulsation equations in a form adequate for numerical calculations we divide our non - perturbed model in a number of spherical, concentric shells. This is equivalent to divide the dominion of the independent variable in $`N`$ mesh points ($`N1`$ shells) non - necessarily equidistantly spaced $`x_j,j=1,\mathrm{},N`$. In our treatment we define $`x_1=1`$ as the surface and $`x_N=0`$ as the central point of the model.

Now, we replace Eqs. (1-4) by difference equations. The system of equations may be rewritten as

$$\frac{dy_i}{dx}=f_i(y_1,y_2,y_3,y_4,\lambda ),i=1,2,3,4,$$
(14)
where $`\lambda =\omega ^2`$. In finite differences, Eq. (14) may be written as222We have also employed finite differences of the form
$$\frac{[y_i]_{j+1}[y_i]_j}{x_{j+1}x_j}=\frac{1}{2}\left(f_i([y_1,y_2,y_3,y_4]_j;\lambda )+f_i([y_1,y_2,y_3,y_4]_{j+1};\lambda )\right);i=1,2,3,4;j=1,\mathrm{},N1.$$
Applying such form for the difference equations we have found results negligibly different from those presented in this paper.

$$\frac{[y_i]_{j+1}[y_i]_j}{x_{j+1}x_j}=f_i([y_1,y_2,y_3,y_4]_{j+\frac{1}{2}};\lambda );i=1,2,3,4;j=1,\mathrm{},N1$$
(15)
where, for any quantity $`\psi `$,

$$[\psi ]_{j+\frac{1}{2}}=\frac{[\psi ]_j+[\psi ]_{j+1}}{2},$$
(16)
being $`[\psi ]_j`$ the value of the quantity $`\psi `$ at the meshpoint $`j`$. The outer boundary conditions are given by333The first two of these equations are the so called “zero boundary condition”, which we shall employ in the treatment of polytropes. In the case of realistic WD models we have replaced the first one by Eq. (29) (see below). (Unno, et al. 1989)

| $`[y_1]_1[y_2]_1+[y_3]_1=0,`$ |
| --- |
| $`(\mathrm{}+1)[y_3]_1+[y_4]_1=0,`$ |
| $`[y_1]_1=1.`$ |

(17)
The last equation is the normalization condition, usually employed in previous works. The central boundary conditions are given by (Unno, et al. 1989)

| $`[y_1]_N\frac{[C_1]_N\lambda }{\mathrm{}}[y_2]_N=0,`$ |
| --- |
| $`\mathrm{}[y_3]_N[y_4]_N=0.`$ |

(18)
In order to solve this system of difference equations we shall use the general Newton - Raphson technique closely following the formulation presented in Kippenhahn, et al. (1967) for the case of stellar evolution. We shall assume we are able to get an approximate solution for the system and we want to improve it iteratively. In the case that this initial solution is not far from the actual solution of the system, we may expand the equations up to first order in the corrections for the values of the eigenfunctions at each meshpoint and also for the eigenvalue (square of the eigenfrequency). In this standard way we get a linear system of equations in the corrections that must be solved.

In condensed way, the algebraic system of equations for the first order corrections may be expressed as

$$\frac{B_k}{[y_1]_1}\delta [y_1]_1+\mathrm{}+\frac{B_k}{[y_4]_1}\delta [y_4]_1+\frac{B_k}{\lambda }\delta \lambda =B_k;k=1,2,3,$$
(19)
| $`\frac{G_i^j}{[y_1]_j}\delta [y_1]_j+\mathrm{}+\frac{G_i^j}{[y_4]_j}\delta [y_4]_j+\frac{G_i^j}{[y_1]_{j+1}}\delta [y_1]_{j+1}+\mathrm{}+\frac{G_i^j}{[y_4]_{j+1}}\delta [y_4]_{j+1}+`$ |
| --- |
| $`\frac{G_i^j}{\lambda }\delta \lambda =G_i^j;i=1,2,3,4;j=1,2,\mathrm{},N1,`$ |

(20)
$$\frac{C_m}{[y_1]_N}\delta [y_1]_N+\mathrm{}+\frac{C_m}{[y_4]_N}\delta [y_4]_N+\frac{C_m}{\lambda }\delta \lambda =C_m;m=1,2,$$
(21)
where $`\delta [y_i]_j`$ are small corrections to the eigenfunction $`y_i`$ at the $`j`$ mesh point, and $`\delta \lambda `$ stand for the correction to the eigenvalue $`\lambda `$. $`B_k`$, $`G_i^j`$, and $`C_m`$ are the values of the difference equations when applied the solution to be iteratively improved. Obviously, all of them must be zero when evaluated with their exact solution.

Now, we have to invert a big matrix, that has non - zero elements only in blocks located on the diagonal and in the last column (because of the derivative with respect to the eigenvalue). Notice that this is an important difference compared with the case of stellar evolution. Thus, in handling the big matrix we need a specific algorithm to solve for the corrections.

Evaluating Eq. (19) for $`k=1,2,3`$ and Eq. (20) for $`i=1,2,3,4`$ and $`j=1`$, we may write the first block of such matrix as

$$\left[\begin{array}{ccccccc}\frac{B_1}{[y_1]_1}& \mathrm{}& \mathrm{}& \frac{B_1}{[y_4]_1}& 0& 0& 0\\ \frac{B_2}{[y_1]_1}& \mathrm{}& \mathrm{}& \frac{B_2}{[y_4]_1}& 0& 0& 0\\ \frac{B_3}{[y_1]_1}& \mathrm{}& \mathrm{}& \frac{B_3}{[y_4]_1}& 0& 0& 0\\ \frac{G_1^1}{[y_1]_1}& \mathrm{}& \mathrm{}& \frac{G_1^1}{[y_4]_1}& \frac{G_1^1}{[y_1]_2}& \frac{G_1^1}{[y_2]_2}& \frac{G_1^1}{[y_3]_2}\\ \mathrm{}& \mathrm{}& & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \frac{G_4^1}{[y_1]_1}& \mathrm{}& \mathrm{}& \frac{G_4^1}{[y_4]_1}& \frac{G_4^1}{[y_1]_2}& \frac{G_4^1}{[y_2]_2}& \frac{G_4^1}{[y_3]_2}\end{array}\right]\left[\begin{array}{c}\delta [y_1]_1\\ \delta [y_2]_1\\ \delta [y_3]_1\\ \delta [y_4]_1\\ \delta [y_1]_2\\ \delta [y_2]_2\\ \delta [y_3]_2\end{array}\right]=\left[\begin{array}{ccc}0& \frac{B_1}{\lambda }& B_1\\ 0& \frac{B_2}{\lambda }& B_2\\ 0& \frac{B_3}{\lambda }& B_3\\ \frac{G_1^1}{[y_4]_2}& \frac{G_1^1}{\lambda }& G_1^1\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \frac{G_4^1}{[y_4]_2}& \frac{G_4^1}{\lambda }& G_4^1\end{array}\right]\left[\begin{array}{c}\delta [y_4]_2\\ \delta \lambda \\ 1\end{array}\right].$$
(22)
Now, we define vectors $`U_s,V_s,W_s,(s=1,\mathrm{},4N5)`$ such that we may write

$$\left[\begin{array}{c}\delta [y_1]_1\\ \delta [y_2]_1\\ \delta [y_3]_1\\ \delta [y_4]_1\\ \delta [y_1]_2\\ \delta [y_2]_2\\ \delta [y_3]_2\end{array}\right]=\left[\begin{array}{ccc}U_1& V_1& W_1\\ U_2& V_2& W_2\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}\\ U_7& V_7& W_7\end{array}\right]\left[\begin{array}{c}\delta [y_4]_2\\ \delta \lambda \\ 1\end{array}\right]$$
(23)
for this block. The coefficients of this vectors may be calculated easily solving

$$\left[\begin{array}{ccccccc}\frac{B_1}{[y_1]_1}& \mathrm{}& \mathrm{}& \frac{B_1}{[y_4]_1}& 0& 0& 0\\ \frac{B_2}{[y_1]_1}& \mathrm{}& \mathrm{}& \frac{B_2}{[y_4]_1}& 0& 0& 0\\ \frac{B_3}{[y_1]_1}& \mathrm{}& \mathrm{}& \frac{B_3}{[y_4]_1}& 0& 0& 0\\ \frac{G_1^1}{[y_1]_1}& \mathrm{}& \mathrm{}& \frac{G_1^1}{[y_4]_1}& \frac{G_1^1}{[y_1]_2}& \frac{G_1^1}{[y_2]_2}& \frac{G_1^1}{[y_3]_2}\\ \mathrm{}& \mathrm{}& & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \frac{G_4^1}{[y_1]_1}& \mathrm{}& \mathrm{}& \frac{G_4^1}{[y_4]_1}& \frac{G_4^1}{[y_1]_2}& \frac{G_4^1}{[y_2]_2}& \frac{G_4^1}{[y_3]_2}\end{array}\right]\left[\begin{array}{ccc}U_1& V_1& W_1\\ U_2& V_2& W_2\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}\\ U_7& V_7& W_7\end{array}\right]=\left[\begin{array}{ccc}0& \frac{B_1}{\lambda }& B_1\\ 0& \frac{B_2}{\lambda }& B_2\\ 0& \frac{B_3}{\lambda }& B_3\\ \frac{G_1^1}{[y_4]_2}& \frac{G_1^1}{\lambda }& G_1^1\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \frac{G_4^1}{[y_4]_2}& \frac{G_4^1}{\lambda }& G_4^1\end{array}\right].$$
(24)
From Eq. (20), and after short manipulation, we may write an arbitrary block of the big matrix (except the central) in the form

$$\left[\begin{array}{c}\delta [y_4]_j\\ \delta [y_1]_{j+1}\\ \delta [y_2]_{j+1}\\ \delta [y_3]_{j+1}\end{array}\right]=\left[\begin{array}{ccc}U_{4j}& V_{4j}& W_{4j}\\ U_{4j+1}& V_{4j+1}& W_{4j+1}\\ U_{4j+2}& V_{4j+2}& W_{4j+2}\\ U_{4j+3}& V_{4j+3}& W_{4j+3}\end{array}\right]\left[\begin{array}{c}\delta [y_4]_{j+1}\\ \delta \lambda \\ 1\end{array}\right],$$
(25)
where $`j=2,\mathrm{},N2`$. Now, the coefficients of vectors $`U_s,V_s,W_s`$ may be evaluated solving the system

$$\left[\begin{array}{cccc}\alpha _1^j& \frac{G_1^j}{[y_1]_{j+1}}& \frac{G_1^j}{[y_2]_{j+1}}& \frac{G_1^j}{[y_3]_{j+1}}\\ \alpha _2^j& \frac{G_2^j}{[y_1]_{j+1}}& \frac{G_2^j}{[y_2]_{j+1}}& \frac{G_2^j}{[y_3]_{j+1}}\\ \alpha _3^j& \frac{G_3^j}{[y_1]_{j+1}}& \frac{G_3^j}{[y_2]_{j+1}}& \frac{G_3^j}{[y_3]_{j+1}}\\ \alpha _4^j& \frac{G_4^j}{[y_1]_{j+1}}& \frac{G_4^j}{[y_2]_{j+1}}& \frac{G_4^j}{[y_3]_{j+1}}\end{array}\right]\left[\begin{array}{ccc}U_{4j}& V_{4j}& W_{4j}\\ U_{4j+1}& V_{4j+1}& W_{4j+1}\\ U_{4j+2}& V_{4j+2}& W_{4j+2}\\ U_{4j+3}& V_{4j+3}& W_{4j+3}\end{array}\right]=\left[\begin{array}{ccc}\frac{G_1^j}{[y_4]_{j+1}}& \beta _1^j& \gamma _1^j\\ \frac{G_2^j}{[y_4]_{j+1}}& \beta _2^j& \gamma _2^j\\ \frac{G_3^j}{[y_4]_{j+1}}& \beta _3^j& \gamma _3^j\\ \frac{G_4^j}{[y_4]_{j+1}}& \beta _4^j& \gamma _4^j\end{array}\right],$$
(26)
where quantities $`\alpha _i^j,\beta _i^j,\gamma _i^j`$ are defined as

| $`\alpha _i^j`$ | = | $`\frac{G_i^j}{[y_4]_j}`$ | + | $`U_{4j3}\frac{G_i^j}{[y_1]_j}`$ | + | $`U_{4j2}\frac{G_i^j}{[y_2]_j}`$ | + | $`U_{4j1}\frac{G_i^j}{[y_3]_j},`$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  |  |  |  |  |  |  |  |  |
| $`\beta _i^j`$ | = | $`\frac{G_i^j}{\lambda }`$ | + | $`V_{4j3}\frac{G_i^j}{[y_1]_j}`$ | + | $`V_{4j2}\frac{G_i^j}{[y_2]_j}`$ | + | $`V_{4j1}\frac{G_i^j}{[y_3]_j},`$ |
|  |  |  |  |  |  |  |  |  |
| $`\gamma _i^j`$ | = | $`G_i^j`$ | + | $`W_{4j3}\frac{G_i^j}{[y_1]_j}`$ | + | $`W_{4j2}\frac{G_i^j}{[y_2]_j}`$ | + | $`W_{4j1}\frac{G_i^j}{[y_3]_j}.`$ |
|  |  |  |  |  |  |  |  |  |

(27)
with $`i=1,2,3,4`$.

Finally, in the central point of the model, we have $`j=N1`$ in Eq. (20) and $`m=1,2`$ in (21). It is easy to show that the last block of the big matrix may be written as

$$\left[\begin{array}{cccccc}\alpha _1^{N1}& \frac{G_1^{N1}}{[y_1]_N}& \frac{G_1^{N1}}{[y_2]_N}& \frac{G_1^{N1}}{[y_3]_N}& \frac{G_1^{N1}}{[y_4]_N}& \beta _1^{N1}\\ \alpha _2^{N1}& \frac{G_2^{N1}}{[y_1]_N}& \frac{G_2^{N1}}{[y_2]_N}& \frac{G_2^{N1}}{[y_3]_N}& \frac{G_2^{N1}}{[y_4]_N}& \beta _2^{N1}\\ \alpha _3^{N1}& \frac{G_3^{N1}}{[y_1]_N}& \frac{G_3^{N1}}{[y_2]_N}& \frac{G_3^{N1}}{[y_3]_N}& \frac{G_3^{N1}}{[y_4]_N}& \beta _3^{N1}\\ \alpha _4^{N1}& \frac{G_4^{N1}}{[y_1]_N}& \frac{G_4^{N1}}{[y_2]_N}& \frac{G_4^{N1}}{[y_3]_N}& \frac{G_4^{N1}}{[y_4]_N}& \beta _4^{N1}\\ 0& \frac{C_1}{[y_1]_N}& \frac{C_1}{[y_2]_N}& \frac{C_1}{[y_3]_N}& \frac{C_1}{[y_4]_N}& \frac{C_1}{\lambda }\\ 0& \frac{C_2}{[y_1]_N}& \frac{C_2}{[y_2]_N}& \frac{C_2}{[y_3]_N}& \frac{C_2}{[y_4]_N}& \frac{C_2}{\lambda }\end{array}\right]\left[\begin{array}{c}\delta [y_4]_{N1}\\ \delta [y_1]_N\\ \delta [y_2]_N\\ \delta [y_3]_N\\ \delta [y_4]_N\\ \delta \lambda \end{array}\right]=\left[\begin{array}{c}\gamma _1^{N1}\\ \gamma _2^{N1}\\ \gamma _3^{N1}\\ \gamma _4^{N1}\\ C_1\\ C_2\end{array}\right].$$
(28)
In this case the quantities $`\alpha _i^{N1}`$, $`\beta _i^{N1}`$, $`\gamma _i^{N1}`$ are evaluated from Eqs. (27) at $`j=N1`$.

We note that Eq. (28) can be solved in order to obtain the corrections for eigenfunctions $`y_1`$, $`y_2`$, $`y_3`$ and $`y_4`$ in the central point of the object, and for the eigenvalue $`\lambda `$. Also, it must be noted that the correction in eigenfunction $`y_4`$ pertaining to the immediate outer mesh point of grid is obtained. This correction, namely $`\delta [y_4]_{N1}`$, serve as “coupling” between the points $`N`$ and $`N1`$. In fact, we may use Eq. (25) with $`j=N2`$ in order to obtain the rest of corrections belonging to the eigenfunctions at the mesh point $`N1`$. Next, the employment of this procedure for successive downwards values of $`j`$ (in Eq. 25) using $`\delta [y_4]_{j+1}`$ as coupling between the quantities belonging to consecutive meshpoints, as well as Eq. (23), leads to find the corrections for the eigenvalue and eigenfunctions for the complete model. These corrections are applied to the initial solution and the algorithm is employed iteratively up to the stage at which all the (absolute value) of the corrections are smaller than some previously adopted value. At this point, the complete set of difference equations that represent Eqs. (1-4) has been solved.

In order to search for the first approximation to the eigenfunctions and the eigenvalue of a mode, we have applied the method of the discriminant presented in Unno et al. (1989). Specifically, the expression adopted reads $`D(\omega )=(\mathrm{}+1)[y_3]_1+[y_4]_1`$. Notice that $`D(\omega )=0`$ corresponds exactly to the second outer boundary condition given in Eqs. (17). We refer the reader to that book for more details.

Obviously, for applying this technique we need to define a grid of meshpoints. In getting an appropriate distribution of meshpoints we have employed a simple recipe. As a first approximation to the solution of the oscillation equations we have computed the eigenfunctions and eigenvalue taking the points at which the equilibrium model is defined. After convergence, we have redistributed the grid asking for the same criterium usually employed in stellar evolution (Kippenhahn, et al. 1967): that the relative variation of each eigenfunction inside a zone must be below some prescribed limit (if the eigenfunctions are below some other very small, prescribed value we have asked for absolute differences in place of relative ones). If necessary, our program add or remove meshpoints, using spline interpolation over quantities of non - perturbed model, and lineal interpolation in eigenfunctions, since these will be improved after in the general Newton - Raphson stage.

## 4 Application to Polytropic Spheres

As a first test for our code, we have applied it to the well - known case of nonradial pulsations of polytropic spheres. Polytropes are particularly adequate for the purpose of testing, because in this case we separate the inaccuracies of the non - perturbed model from those inherent to the nonradial pulsation code. In particular we have computed the p- (pressure), f- (fundamental) and g- (gravity) modes of polytropes with indices $`n=1.5,2,2.5,3,3.5`$ and $`4`$ for the case of $`\mathrm{}=2,3`$ and $`4`$.

Because in this case $`V_g`$ and $`A^{}`$ are divergent at surface (see Eqs. 10 and 13), we have integrated the structure by means of an accurate Runge - Kutta technique (Press, et al. 1992) and most of the meshpoints were located near surface. Here we assumed $`\mathrm{\Gamma }_1=\frac{5}{3}`$.

To compare our eigenvalues with those available in the literature we shall employ the work of Christensen - Dalsgaard & Mullan (1994) which presents accurate eigenfrequency tabulations for a variety of polytropic spheres. For such purpose we have plotted in Fig. 1 the relative difference $`|\omega ^2\omega _{ChD}^2|/\omega _{ChD}^2`$ for the cases of p-modes for $`n=1.5,3`$ and $`\mathrm{}=2,3`$, g-modes for $`n=3`$ and $`\mathrm{}=2,3,4`$, and also for p-modes for $`n=4`$ and $`\mathrm{}=2,3`$. The comparison indicates that the higher radial order the larger the relative differences are, without any dependency on $`\mathrm{}`$ value. The absolute value of relative differences found in the eigenvalues are $`4\times 10^4`$, except in the $`n=3`$ g-modes, for which the difference is $`2\times 10^3`$. As Christensen - Dalsgaard & Mullan (1994) state that their computations are accurate to $`10^8`$ we conclude that our code is able to compute the eigenvalues of polytropes with a precision of $`10^3`$ which is enough for our purposes.

## 5 Application to Realistic Models: Nonradial g-Modes in Helium White Dwarfs

With the aim of investigating the behaviour of the code when it is applied to realistic models, in this section we shall consider nonradial g-modes in evolutionary models of low mass, pure helium WDs. This election is due mainly because of two reasons. First, at present we have available a detailed and updated code we have employed in the calculation of the evolution of such stars (see, e.g., Althaus & Benvenuto 1997; Benvenuto & Althaus 1998). Such evolutionary code is fully described in the above-cited papers. Second, the computation of nonradial modes in pure helium WDs provide us an important way for testing the pulsation code when applied to realistic models. At present, to our notice there is no study available in the literature on the pulsation properties of these objects. In evaluating our results we shall pay special attention to the asymptotic behaviour of periods of oscillation, comparing in particular the spacing of adjacent periods (for the same $`\mathrm{}`$) with that predicted by asymptotic theory of Tassoul (1980) in stars with homogeneous composition. This theory has been employed as check for the nonradial spectrum of polytropes (Mullan & Ulrich 1988; Mullan 1989) in the Cowling approximation.

### 5.1 Outline of Pulsation Calculations

One possible way for computing nonradial oscillations modes in evolutionary models is to first compute such structures with a evolutionary code, and then to choose a subset of these models (usually belonging to a predetermined interval in $`T_{eff}`$, called hereafter $`T_{eff}`$-strip) for pulsation analysis. This procedure has been employed in most of the studies of adiabatic pulsations in WDs to date (see, e.g., TFW, Bradley, Winget & Wood 1993, Bradley 1996, Brassard et al. 1992a, 1992b).

Here, we present an alternative way for computing nonradial modes in evolutionary sequences. The basic idea is very direct: if the time step in the evolutionary sequence is short enough, eigenvalues and eigenfunctions corresponding to a given model should strongly resemble those corresponding to the previous one. Then, if we couple the pulsation code to the evolutionary code, it is feasible to avoid the scan in $`\omega ^2`$ for each model: the search of modes is carried out only for the first model inside the relevant $`T_{eff}`$-strip. The modes are computed and stored and serve as a trial solution for the next model and so on. In this way it is possible economize CPU time, and (more importantly) it is possible to follow the changes in the pulsational spectrum due to the evolution of the stellar structure in a continuous fashion.

### 5.2 Details of Calculations

Let us briefly describe how our pulsation and evolutionary codes work together. Firstly, $`T_{eff}`$-strip is chosen, as well as the frequency window to be scanned. The evolutionary code computes the model up to the moment when the model reaches the hot edge of the $`T_{eff}`$-strip. Then, the program calls the pulsation routine beginning the scan for eigenfrequencies, as described in §3. When a mode is found, the code generates an approximate solution for $`y_i;i=1,\mathrm{},4`$ and $`\omega ^2`$, which is improved iteratively. Then, such solution is tested and, if necessary, meshpoints are redistributed as outlined in §3, and iterated to convergence. Before improving the solution, the eigenmodes have been stored on the original grid in a common block, with the aim of being employed later as an approximate solution for the next stellar model of the sequence. When the computation of the mode has been finished, the code begins to look for eigenmodes again, and the process is repeated until the given frequency window is fully covered. Thus, we have finished the computation of all eigenmodes of the first model belonging to the $`T_{eff}`$-strip and the structure of each one is now in the computer memory. Then, the evolutionary code generates the next stellar model and the code calls pulsation routine again, but in this occasion the search for modes is skipped. Instead of this, stored eigenmodes are taken as input to the general Newton - Raphson scheme to approximate the modes of this subsequent stellar model. Then, the solution is iterated, and so on. The whole procedure is automatically repeated for all evolutionary models inside the $`T_{eff}`$-strip. When the star gets outside the $`T_{eff}`$-strip, computation is finished.

In terms of CPU, a large fraction of the running time is spent in the search for modes, but this is executed only once. The process of relaxation and improvement of the solution is faster, and the building of each stellar model is almost instantaneous.

### 5.3 Helium White Dwarfs Models

As mentioned before, we have choose models of helium WDs for checking the efficiency our pulsation code in realistic stellar configurations. To be specific, we computed nonradial g-modes for a helium WD model with mass $`M_{}=0.3M_{}`$, in the range of $`8000KT_{eff}25000K`$. The complete sequence comprise 216 models. Convection, present in the outermost layers, is treated with the ML3 version of Mixing Length Theory.

### 5.4 Results

Since one of our interests here is to evaluate the asymptotic behaviour of eigenperiods, we have computed nonradial g-modes for $`\mathrm{}=1,2`$ and $`3`$ with radial order from $`k=1`$ to $`56`$. The resulting set of modes is adequate for our purpose, since it covers a broad range in periods, enough for enabling the comparison with the predictions of the asymptotic theory of Tassoul (1980). Here we employ the mechanical external boundary condition given in Unno, et al. (1989)

$$y_1\left[1+\left(\frac{\mathrm{}(\mathrm{}+1)}{\omega ^2}4\omega ^2\right)\frac{1}{V}\right]y_2+y_3\left[1+\left(\frac{\mathrm{}(\mathrm{}+1)}{\omega ^2}\mathrm{}1\right)\frac{1}{V}\right]=0,$$
(29)
which replaces the first of Eqs. (17).

For each computed mode, the quantities of interest are the period $`P_k`$, eigenfunctions $`y_i;i=1,\mathrm{},4`$; kinetic energy K.E., and first order rotation splitting coefficient, $`C_{\mathrm{},k}`$. In addition, we compute the variational period, $`P_k^v`$, and the weight function, WF, in the form given by Kawaler, Hansen & Winget (1985). Also, for each non-perturbed model we obtain the asymptotic spacing of periods (for the same $`\mathrm{}`$), $`\mathrm{\Delta }P_A`$, given by (TFW, Tassoul 1980)

$$\mathrm{\Delta }P_A=\frac{P_0}{\sqrt{\mathrm{}(\mathrm{}+1)}},$$
(30)
where $`P_0`$ is defined as

$$P_0=2\pi ^2\left[_0^1\frac{N}{x}𝑑x\right]^1.$$
(31)
As in TFW, for computing $`P_0`$ we have ignored the presence of an OCZ, integrating from the center to the surface of model but setting $`N=0`$ in the OCZ where $`N^2<0`$. In this approximation, we have overestimated the value of $`\mathrm{\Delta }P_A`$, but since the thickness of the OCZ (in the radial coordinate) is very small, this is a good approximation.

The results obtained for $`\mathrm{}=1,2,3`$ are qualitatively very similar. Although pulsations have not yet been observed in these stars, in other pulsating WDs the modes with $`\mathrm{}=1`$ dominate the observations, thus we shall concentrate in $`\mathrm{}=1`$.

In Fig. 2 we have plotted the eigenfunction $`y_1`$ for the modes g1, $`\mathrm{}`$, g5 with $`\mathrm{}=1,2`$, for a model with $`T_{eff}=11900`$ K. From Fig. 2 it is evident that $`y_1`$ has large amplitudes in the whole star, especially in the core. This feature is strongly emphasized when we inspect modes with increasing $`k`$. This is a remarkable difference in the behaviour of $`y_1`$ comparing with the case of oscillating DA and DB WDs. In such WDs $`y_1`$ has a lower central amplitude (for modes with $`\mathrm{}=2`$ in a $`0.6M_{}`$ DA WD model see, e.g., Fig. 14a and 16a of Brassard et al. 1992b).

Fig. 3 displays the normalized weight function (WF) for g1,$`\mathrm{}`$,g5 ($`\mathrm{}=1`$) of the same model. WF provides a measure of the relative contribution from different zones in the model to the period formation. The plot indicates that period formation happens mainly in external region of star, where WF is large. However there are important contributions from central locations ($`x0.2`$), which is somewhat different to the WF corresponding to the case of DA and DB WDs, for which WF displays appreciate values only in the envelope (in the “normal” modes; see Figs. 15 and 17 of Brassard et al. 1992b).

Let us discuss our results concerning periods ($`P_k`$) and kinetic energies (K.E.). The periods here computed show the expected trend in g-modes, with values increasing with the order of modes $`k`$. Since our calculations covers a wide $`T_{eff}`$-strip, it is possible infer the changes of periods during WD cooling.

It is worth mentioning that we have computed the Brunt - Väisälä frequency as in TFW (the Ledoux term $`B`$ is zero for the whole of model, since it is chemically homogeneous). We have found that $`N^2`$ decreases as the model cools down, a feature that is very pronounced in the core. It is due mainly to the increasing degeneracy in that region, decoupling more and more pressure from temperature ($`\chi _T`$ drops). As consequence eigenperiods increase monotonically when $`T_{eff}`$ drops below $`20000`$ K, as clearly shown in Fig. 4A.

Next, we consider the behaviour of the kinetic energy (K.E.) of the modes. We have found that, as expected in a chemically homogeneous star, K.E. is a smooth function of $`k`$ and thus of the period. In Fig. 5 we show $`\mathrm{log}`$(K.E.) vs. $`T_{eff}`$ for the models included in Fig. 4. For hotter models it is clear that the more energetic modes are those with lower order, because they penetrate deep in the star, where density is high. However, when $`T_{eff}`$ drops below $`20000`$ K the energy of high order modes (which are mainly concentrated in the outer envelope of star) strongly increase. The explanation of this effect (see Brassard et al. 1992b in context of a DA WD star) resides in the fact that at such $`T_{eff}`$ value the OCZ of the object suddenly gets thicker (see the dotted line in Fig. 17 of Althaus & Benvenuto 1997). The high order modes that oscillates in the envelope feel gradually the presence of convection as the star cools down. Since that in OCZs g-modes become evanescent, such modes are forced to get larger amplitudes somewhat below the bottom of the OCZ where the density is larger. Since K.E. is proportional to the integral of square of displacement, weighted by $`\rho `$, these modes oscillate with larger energies. In sharp contrast to the behaviour of high $`k`$ modes, low order modes are rather insensitive to the thickening of the OCZ.

Now let us discuss the period spacing of the modes. To our knowledge there is no study available in the literature devoted to nonradial pulsations of helium WDs. This situation does not allow us to perform a direct comparison of our results. Nevertheless we can have a good idea related to the reliability of our pulsation code examining the period spacing of consecutive modes of the same $`\mathrm{}`$. Tassoul (1980) predicted that the asymptotic period spacing for g-modes for a chemically homogeneous object in the adiabatic approximations is a constant value $`\mathrm{\Delta }P_A`$ (see Eq. 30). In Fig. 6 we show the forward period spacing $`\mathrm{\Delta }P_k=P_{k+1}P_k`$ vs. radial order $`k`$ ($`\mathrm{}=1`$) for models with $`T_{eff}=13200,11800`$ and $`9400`$ K. As reference we also include in this figure the corresponding asymptotic spacing predicted by Tassoul (1980) treatment. From Fig. 6 it is clear that the trend of the numerical solution is the correct one.

In order to get an estimation of the accuracy of our treatment, we have calculated the variational eigenperiods ($`P_k^v`$) for each mode in the formalism presented by Kawaler, et al. (1985). The comparison of variational and numerical eigenperiods for the values of $`\mathrm{}`$ here considered gives a difference between these treatments which is lower than $`1\%`$ (see Fig. 7 for the case of $`\mathrm{}=1`$).

By the way, let us comment on that we have also tested our code running it on the same models for which oscillation modes were previously computed. Specifically we have employed for such test two carbon - oxygen DA WD models of 0.5 $`M_{}`$ and 0.85 $`M_{}`$ whose structure was computed with the WDEC evolutionary code and its vibrational properties were previously analyzed (Bradley 1996). In this case $`M_H/M_{}=10^4`$; $`M_{He}/M_{}=10^2`$ and $`T_{eff}12500`$K. This is a valuable test for our pulsation code because we could confront our results with those of other researchers working on detailed models. Also the fact that these models were computed with the same evolutionary code is valuable because differences in the resulting periods are due only to the pulsation code and not to differences in e.g. between the treatment of the equation of state in our evolutionary code and in WDEC. Considering a large amount of modes, the differences between the previously computed periods and those we have found has been below $`0.1\%`$. Thus, we judge that our code produce results accurate enough for our purposes. This comparison has been made possible by Paul Bradley who provided us with the DA WD models and the periods of many its modes.

## 6 Summary and Conclusions

In this paper we have presented a general code intended for calculating linear, adiabatic, nonradial eigenfunctions and eigenfrequencies corresponding to spherically symmetric stars by means of a finite difference scheme. This is a very simple strategy based on a slight modification of the general Newton - Raphson method for computing stellar evolution presented by Kippenhahn, Weigert & Hofmeister (1967) because of the different structure of the matrix to be inverted in finding the iterative corrections.

As a first test to this code we have applied it to the simplest case: polytropic spheres. As there exist very precise tabulations of the eigenfrequencies for such configurations, we have been able to perform a quantitative comparison of our results with the best available in the literature (Christensen - Dalsgaard & Mullan 1994). We have found that our code is able to compute polytropic eigenfrequencies with an accuracy of $`10^3`$ or much better for low $`k`$ values.

Then, we have coupled this nonradial oscillation code to our evolutionary code and applied it to the computation of the pulsation spectrum of pure, low mass helium white dwarfs (WDs). Specifically we considered the case of an object with a mass of 0.3 $`M_{}`$ for a wide range of effective temperatures ($`T_{eff}`$-strip). In making such computation we have employed a new method. Taking a short - enough time step in the evolutionary sequence, eigenvalues and eigenfunctions corresponding to a given model strongly resemble those corresponding to the previous one. Then, we scan for modes in the first model inside the $`T_{eff}`$-strip. Eigenmodes are computed and stored and serve as a trial solution for the next model and so on. In this way it is possible to follow the changes in the the pulsational spectrum due to the evolution of the stellar structure in a continuous and fairly detailed fashion (see Fig. 4).

Regarding the computation of g-modes for the 0.3 $`M_{}`$ pure helium WD model we have presented eigenfunctions and weight function (WF) for low order $`k`$ modes, and periods, period spacing, kinetic energy (K.E.) and variational periods for a wide range of modes. These quantities were computed for the case of $`\mathrm{}=1,2`$ and $`3`$ but in Figs. 3 \- 7 we have presented only the case $`\mathrm{}=1`$ for brevity. To our notice this is the first effort in studying the pulsation properties of helium WDs.

Our experience with realistic helium WD models indicates that our numerical code for pulsations work nicely and in particular the solutions we have found have an asymptotic period spacing in good agreement with the predictions of the analytical theory of Tassoul (1980). This, and the fact that comparison with more complex WD models with layers of different internal composition show us that our numerical tool for computing eigenmodes work adequately, especially regarding the coupling we have performed between it and our evolutionary code. We plan in the near future to apply this new code to the study of pulsation of WDs in general.

The authors would like to warmly acknowledge Paul Bradley for his kindness in providing us with numerical models of carbon - oxygen DA WDs together with his pulsation calculations and also for the time he spent in doing so. This allowed us to perform key a test to our code, such that it made possible to us to be confident with the results the code produces. We also want to thank our referee D. Winget for suggestions that allowed us to improve the original version of this paper.


text:
# References

MORPHOLOGY OF CIRCUMSTELLAR ENVIRONMENT AND
SOME CHARACTERISTICS OF CIRCUMSTELLAR SHELLS OF STARS
WITH THE R CORONAE BOREALIS VARIABILITY

A.E.Rosenbush

© 2000

Main Astronomical Observatory of National Academy of Sciences of Ukraine
03680, Golosiiv, Kyiv-127, Ukraine
e-mail: mijush@mao.kiev.ua

The well-known light minima of stars with the R Coronae Borealis variability are caused by the formation of an additional circumstellar dust shell, the screening shell, inside the permanent shell.

Under the assumption of uniform distribution of matter in the circumstellar environment we estimated the optical thickness of the permanent gas–and–dust shell at 0.2–0.7, and its geometrical thickness is no less than 0.4 of its own radius. The wavelength dependence of extinction is close to neutral.

From spectral observations of R CrB itself in the 1985 minimum we traced the transformation of the stellar linear and molecular absorption spectrum to the emission spectrum and established that the fast variation of the U–B colour index by –0.6 in the light decline was caused purely by a change of the spectrum type. The spectrum transformation causes an increase of star brightness in the U, B, and V bands by about 1.4, 0.75, and 0.75 mags, correspondingly.

It is suggested that a high–velocity ($`>`$200 km/s) matter stream through the circumstellar environment is the cause of the excitation of the emissions observed during light minima when the photospheric flux is weakening.

INTRODUCTION

The set of unique characteristics of the stars with the R Coronae Borealis (RCB) type variability keeps constantly the high interest to them from observers and theorists. For all this, side by side with the progress in observations, negligibility in the theory is noted.

The principal variability – deep light declines up to 8 mag which last hundreds of days and are caused by dust condensation on the line of sight – remains the phenomenon with many questions, including even the lack of a simple geometric model of the phenomenon .

The basis of the modern model is the cloud structure of the circumstellar permanent dust shell  when one of the dust clouds in this shell forms on the line of sight, and, as a result, a visual light minimum is observed. This explains the lack of anticorrelation between variations of visual and infrared brightness during a minimum . The lack of clear influence of the permanent shell, in which up to 40 percent of the bolometric luminosity of star is transformed, on apparent characteristics of the star at a light maximum supports the idea of its cloud structure. The difficulties of this model in simulating light minima were well known .

As the result of our studies, we arrived at the understanding that it is necessary to investigate not variables of the R Coronae Borealis type but the R Coronae Borealis phenomenon observed in various objects, for instance, in novae. It is need the high mass loss rate and the enhanced abundance of carbon. The hydrogen is necessary for the phenomenon. Hydrogen deficiency is a consequence of the phenomenon, and its full exhaustion means the cessation of the variability. This view helps us to understand many sides of the phenomenon. At a time we concluded that the 1972 rejection of the model of a spherical dust shell for the interpretation of the R Coronae Borealis phenomenon  was not sufficiently justified: not all sides of the phenomenon were known and were taken into consideration. Below we consider some principal consequences from the assumption of the uniform distribution of matter in the circumstellar environs.

STRUCTURE OF CIRCUMSTELLAR ENVIRONS

Three spherical shells can exist in the environs of a star: two constant shells, the permanent and the fossil (the latter is not connected directly with current light minima), and a temporal, screening one. Both shells are formed in the outflow of mass with a high loss rate estimated at about $`10^6`$ mass of the Sun per year.

The radius of permanent shell was estimated at 26 radii of the star or 2340 radii of the Sun . The formation of this shell can occur in the framework of the homogeneous dust nucleation theory by Fadeyev , which gives the dust condensation zone radius close to this value. The mass of the shell is about $`10^6`$ mass of the Sun.

Sometimes conditions can arise for dust condensation very close to the star, i.e., a screening shell with a radius of 4-10 radii of star is formed . We estimated the full mass of this shell at about $`10^7`$ mass of the Sun. The screening shell has its own infrared excess which compensates the weakening of the permanent shell luminosity. The effective temperature of screening shell during its maximum optical thickness drops down to 700-900 K, which is equal to the permanent shell temperature.

OPTICAL THICKNESS OF PERMANENT SHELL

It is known that the permanent shell re–radiates up to 42 percent of the bolometric luminosity of the star itself . If we assume that the shell only re–radiates, we have an indirect estimate of optical thickness of this shell at about 0.7.

FG Sge allows a more direct estimation. Since 1992 it shows the RCB type variability that allows us to compare stellar parameters before and after 1992. Since 1992 FG Sge was at the state of maximal brightness, i.e., it had only the permanent shell, in the end of 1997 - the beginning of 1998, in 1999 and 2000. The difference of average brightness in these periods from the level before 1992 may be due to the extinction in the permanent shell. Our estimate for this is about 0.26m. It means that the permanent shell absorbs nearly 20 percent of the bolometric luminosity of the star itself. From published data we can estimate that the bolometric luminosity of shell is nearly 25 percent of the bolometric luminosity of the star itself.

SIZE OF DUST PARTICLES

RY Sgr has a number of observations in the UV, where dust particles primarily show an evidence of their existence after condensation.

From the graphic and tabulated data of , in addition to the author’s results, we can obtain the following data.

Firstly, the light minima were not seen in the UV. The decline of the 240 nm flux during a minimum occurs probably only up to the weakening of visual brightness by nearly 2m. The following rise of extinction in the visual is not accompainied by a rise of extinction in the UV. According to , this may occur by the change of a dust particle size from 0.01 micron to 0.3 micron. The particle density does not increase, otherwise the extinction should continue to increase. From this we can draw the unexpected and interesting conclusion that the light minimum may be considered only as a consequence of a change of the dust particle sizes but not their density, and the dust nuclei are constantly present in the stellar atmosphere.

Secondly, the amplitude of the UV brightness pulsations outside minima is a factor of 1.8 higher, on the average, than in the visual pulsations. The higher amplitude of the UV pulsations can be interpreted as the consequence of condensation of dust particles, their growth up to radii of 0.003 micron and disruption in every pulsation.

The dust particles move away from the star, dragging the gas of the stellar atmosphere.

OUTFLOW VELOCITY AND ORIGIN OF EMISSION LINES

The kinetic energy of helium atoms which are the main gas component of the stellar atmosphere is as high as 800 eV and more at a 200 km/s velocity. This energy is sufficient to excite emission by the interaction with low–velocity atoms. It is known that emissions of C IV\] 155.0 nm, C II 133.5 nm, etc. are observed beyond minima, i.e., we have a clear indication on their origin in the permanent shell. If there is no dust, there are no UV emissions mentioned above. This is observed in XX Cam and HD 182040.

This transformation of energy occurs on the inner boundary of permanent shell on a level of 10–26 star radii at a density of about $`10^{14}`$$`10^{15}`$ g/cm or a helium density of about $`10^9`$ cm-3.

Thus we can identify the region of the origin of sharp emission lines or the ”chromosphere” of the star with the zone of dust condensation or with the permanent and screening dust shells.

During light minima the absorption spectrum disappears and the emission spectrum appears. Such transformation determines almost completely the behaviour of U–B and B–V indices observed during a visual minimum.

VARIATIONS OF COLOUR INDICES AND SPECTRUM IN A LIGHT MINIMUM

Based on the R CrB minimum in 1985, for which our extensive spectral observations are available and a detailed photometry is published, we studied the effect of the disappearance of absorption lines and appearance of emission lines on the colour indices of the broad–band photometric U, B, and V system.

We obtained that the key points of spectral changes are three levels of brightness weakening.

I – 1.4–1.6m – no considerable changes;

II – 2.9–3.4m – the disappearance of absorptions and the appearance of emissions;

III – 4.3m – the maximum of emission intensity and its following sharp drop.

It is necessary to compare these levels with the light and colour curves. After the light fading to the first level a sharp decrease of U–B began. At the second level, a sharp increase of B–V occurred. At the third level, the brightness is minimal, it is the reverse point of the light curve and a sharp extremum of U–B; the B–V behaviour did not change.

Thus we make the opposite conclusion that the colour change is due mainly to the change from the absorption spectrum to the emission, and all variations in U–B from +0.3 to –0.3m in 1985 were related to this cause only. The B–V variation did not exceed –0.1m. The estimation of spectral transformation on our spectrograms for the 1985 minimum in the B band gives a value of line contribution of about 0.75m. Thus one can say that the lines contribute about 0.75m in the B band, and their contribution in the V band was slightly smaller. The contribution in the U band was well above, and it reached 1.4m. The light minima are distinguished by the intensity of emission lines. Therefore, the colour variations were dissimilar but the cause of variations is the same.

This moves away one of the principal objections against the uniform distribution of dust in the circumstellar shells.

The well–known loops of a star on the ”colour–brightness” diagram during its light minimum should be interpreted with this circumstance taken into consideration.

CONCLUSION

As an addition to the above, we present a possible version of the development of the RCB type variability in FG Sge as a consequence of formation of an uniform circumstellar dust shell.

Years before 1992 FG Sge had a small infrared excess appropriate to of a low dust mass of 4\*$`10^{10}`$ mass of the Sun. The IR excess was constantly rising, and at some stage of its increase the first RCB minimum occurred. This behaviour is very similar to the RCB phenomenon in the DQ Her novae , when the IR excess appeared well before the beginning of the RCB–like visual decline at the transition stage of a nova outburst. I.e., in FG Sge the first RCB minimum was also connected with the formation of a screening shell. At the same time there was no permanent gas–and–dust shell, as no characteristic emissions were observed during the first minimum, but only the blue-shifted circumstellar absorption D Na I lines. It is unlikely that the first dust cloud as a cause of a light fading is formed at once on the line of sight. If this were so, the ”chromospheric” emissions would necessarily be observed in the 1992 minimum as they were observed in following ones. This supports our assumption about the origin of emissions in the permanent shell.

Since 1972 we have again all grounds to use the dust shell model to interpret the R Coronae Borealis phenomenon. The lack of obvious progress in the theory is connected with the one–sided approach to this ”unique” variability. It is a widespread phenomenon, and it is more correct to study the RCB phenomenon, observing various objects, including novae at the transition stage of outburst.

The full text will be published as two papers in the ”Astrofizika/Astrophysics” journal in 2000 and 2001.

REFERENCES