text:
# Evolution of the structural and bonding properties of Aluminum-Lithium clusters

## I Introduction

The discovery of C60 buckminsterfullerene , observation of magicity in metallic clusters , unusual thermodynamical properties of clusters of atoms like Na , Sn  etc., has given rise to an explosive growth of research in the field of cluster physics. Extensive theoretical and experimental  work has been carried out to understand the physical and chemical properties of clusters like size-evolutionary pattern in cluster geometry, thermodynamics, vibrational and rotational properties, optical properties, electronic structure, and bonding, chemical reactivity, as well as the emergence of the bulk properties with the increase in cluster size.

Recently observed aromaticity  and antiaromaticity  in inorganic molecules, and in particular Al-Li clusters, has stimulated further research in these systems. These properties are well-known and important in organic chemistry . However, recently the first experimental and theoretical evidence of aromaticity in all-metal system viz. Al$`{}_{}{}^{2}{}_{4}{}^{}`$ and $`M`$Al$`{}_{}{}^{}{}_{4}{}^{}`$ ($`M`$ = Li, Na, Cu) was reported by Wang and coworkers . They found the Al$`{}_{}{}^{2}{}_{4}{}^{}`$ dianion to be square-planar possessing two $`\pi `$-electrons, thus conforming the structural criterion for aromaticity. The concept of antiaromaticity, introduced by Breslow  et al., corresponds to the destabilization seen in the cyclic systems with $`4n\pi `$-electrons. Although such molecules have not yet been observed except in organic chemistry, calculations reported by Shetty  et al. and others  indicate the Al4Li4 cluster to be antiaromatic.

Aluminum-lithium clusters, which we study in this paper, are also interesting due to their unusual structural and bonding properties that are quite different from those of the pure clusters of the constituent elements. Both pure Al and Li clusters have been extensively studied using the spherical jellium model (SJM) . However, it is observed that although the SJM has been quite successful in describing the gross electronic structure and the stability of alkali metal atom clusters , the Al clusters present an interesting contrast. Rao and Jena  have reported a comprehensive study of the aluminum cluster using ab initio density functional theory. They found that in small Al clusters, the effective valency of Al is one  due to a large energy gap of about 5eV separating the $`3s^2`$ and $`3p^1`$ orbitals. However, Rao  et al. as well as Dhavale  et al. did not find any signature of such monovalent behavior, contrary to the expectations. Investigation by Cheng  et al. shows that a single Al atom in Li clusters introduces a localized bond between the Al impurity and the Li host. They also found a magic cluster AlLi5, and suggested that AlnLi5n can aggregate at least for some values of $`n2`$ and may exhibit properties characteristic to an assembly of such clusters from the AlLi5 subunits. Akola  et al., however, found that this idea does not apply after $`n=2`$. Further, their investigation on bonding in Li-rich Al-Li clusters show a substantial charge transfer from Li atom and nearby Al atom, strengthening the ionic Al-Li bond, while the Al-Al bond gained a more covalent nature. Such a charge transfer has been observed by Kumar  in several mixed Al-Li clusters. He found a layered Al10Li8 compound to be magic with electronic and geometric shell closing. This was the first instance where a shell closure in $`sp`$ bonded metal cluster was found to occur at 38 valence electrons. Note that the magic numbers for alkali metal atom clusters  are 2, 8, 18, 20, 40,…, whereas, some of the magic Al clusters are Al$`{}_{}{}^{+}{}_{7}{}^{}`$, Al$`{}_{}{}^{}{}_{7}{}^{}`$, Al$`{}_{}{}^{}{}_{11}{}^{}`$, Al$`{}_{}{}^{}{}_{13}{}^{}`$, etc .

The Al-Li bulk is especially interesting due to the fact that it forms a stable alloy at over the wide range of concentration. However, the most stable intermetallic B32 phase corresponds to the 50:50 concentration . In this B32 phase, a mixture of covalent and ionic bonding is seen due to a remarkable charge transfer from the Li atoms to the Al atoms. Moreover, the close resemblance of the density of states (DOS) (without any band gap)  to that of the covalently bonded diamond structure indicates a tetravalent behavior of Al atom. Thus, based on these observations, one can expect that the behavior of Al in AlnLin clusters would be similar to that of the tetravalent atom Si , Ge , Sn, and Pb . Indeed, our present work shows that in certain clusters, structure and bonding of Aln is similar to that of the group IV A clusters. Moreover, earlier work on the heterogeneous Al-Li clusters has focussed on specific aspects like systematics in the geometry, stability, shell closure, magic behavior etc. Bonding in such clusters was discussed on the basis of total charge density as well as the difference of the self-consistent charge density, $`\rho _{scf}`$, and the superimposed atomic charge densities $`\rho _{superimposed}`$ of the constituent atoms. On the contrary, analysis of the molecular orbitals (MOs) has revealed some unusual properties like aromaticity and antiaromaticity in all-metal Al-Li clusters. In our earlier work, we have examined some of the issues concerning the geometry and the stability of various Al-Li clusters viz. AlnLin ($`n`$ =, 1-10, 13) , AlnLi7 ($`n`$=1-7) , AlLin ($`n`$=1,8)  and Al13Lin ($`n`$=1-4,10,19,20,21) . We have found that the geometries of these clusters were dictated by the geometry of the core Al cluster, enclosed in the Li cage. However, since these calculations were performed by the density based molecular dynamics (DBMD), an analyses of the bonding and electronic structure could not be done. In the present work, we study the systematics of the geometries, energetics, electronic structure and the bonding properties in various Al-Li clusters viz. AlnLin ($`n`$=1-11), Al$`{}_{}{}^{}{}_{2}{}^{}`$, Al$`{}_{}{}^{2}{}_{2}{}^{}`$, Al2Li, Al2Li-, and Al6Li8 using Kohn-Sham formulation of the density functional theory (DFT) within the pseudopotential and the generalized gradient approximation (GGA). The bonding in these clusters is analyzed via the electron localization function (ELF)  and the molecular orbitals. In section II we shall describe the computational details, followed by a discussion of the results in section III.

## II Computational Details

The ground state geometries as well as other low-lying structures were obtained in two stages. In the first stage, Born-Oppenheimer molecular dynamics (BOMD)  based on Kohn-Sham (KS)  formulation of density functional theory (DFT) was employed. The computer code used for this purpose was developed in our own group. The total energy during each of the molecular dynamics step was minimized using damped equation of motion . The calculations were performed using norm-conserving pseudopotentials of Bachelet  et al., in the Kleinman-Bylander  form with $`s`$ part treated as nonlocal. The exchange-correlation potential was calculated using the local density approximation (LDA) given by Ceperley-Alder . A cubic supercell of length $`40a.u.`$ with an energy cutoff of $`17Rydberg`$ was found to provide sufficient convergence of the total energy. Increasing the energy cutoff did not improve the total electronic energy significantly.

The damped equation of motion scheme  permits use of a fairly moderate time step $``$100a.u. Starting from a random configuration, the AlnLin clusters were heated to 1400-1700$`K`$, and allowed to span the phase space for a few thousand iterations. In our previous investigations  on AlnLin clusters, it was found that the Li atoms segregate at the surface with the Al atoms forming an inner core. In order to avoid any such segregation during the molecular dynamics run, we have interchanged the Al and the Li atoms. This has ensured that the cluster visits all its local minima as well as the global minimum. At the end of each ionic displacement, the norm of the eigenstates defined as $`h\psi _iϵ_i\psi _i^2`$ (where $`ϵ_i`$ is an eigenvalue corresponding to the eigenstate $`\psi _i`$ of the Kohn-Sham Hamiltonian $`h`$) was within the range of $`10^410^7`$a.u.

In the second stage, various low-lying structures were obtained by conjugate gradient and/or steepest descent  method starting from various suitable configurations during the molecular dynamics run. These configurations were then optimized using ultrasoft pseudopotentials  within the generalized gradient approximation (GGA) implemented in the VASP  package. The Perdew-Wang  exchange-correlation potential for GGA has been used. The size of the simulation cell was varied according to the cluster studied (see Table-I).

Table - I

Size of the supercell (in Å) for various clusters, where $`n`$ is the total number of atoms in the clusters.

|  | $`n`$ | Simulation Cell (Å) |  |
| --- | --- | --- | --- |
|  | 2,4,6 | 16x16x16 |  |
|  | 8,10,12 | 18x18x18 |  |
|  | 14,16.18,20 | 20x20x20 |  |

The geometries were optimized with a kinetic energy cutoff of the order of $``$ 12 $`Rydbergs`$. This energy cutoff was sufficient for the convergence of the total energy. Increasing the energy cutoff did not improve the total electronic energy significantly. The structures were considered to have converged when the forces on each ion was less than 0.01eV/Å with a convergence in the total energy within the range of $`10^410^6`$eV. In general, we find that there are many isomeric structures nearly degenerate to the lowest energy state. Few of the structures can be obtained by interchanging the Al and the Li atoms or by rearranging the position of Li atoms. In the present work, we discuss only a few geometries corresponding to the lowest energies.

The nature of the bonding has been investigated using the electron localization function (ELF)  along with the molecular orbitals (MO) as well as the difference charge density. The ELFs have been found to be useful for elucidating the bonding characteristics of a variety of systems, especially in conjunction with the charge density. The value of the ELF lies between 0 and 1, where 1 represents a perfect localization of the valence charge density. The difference charge density is the difference between the self-consistent charge density, $`\rho _{scf}`$, and the superimposed atomic charge densities, $`\rho _{superimposed}`$, of the constituent atoms.

## III Results and Discussion

### A Structure, energetics and stability

The ground state geometries and some of the low energy structures of AlnLin, and Al6Li8 clusters are shown in figures 1 ($`n`$=2-6, and Al6Li8), and 2 ($`n`$=7-11). The ground state geometries are divided into two broad categories, for (i) $`n`$=2-4 (figures 1a(i)-1c(i)), and (ii) $`n`$=7-9, and 11 (figures 2a(i)-2c(i), 2e(i)), depending on the growth pattern. In the first category a quinted roof-like structure of Al2Li2 is seen (a quinted roof structure is a rhombus with a bend at its diagonal). The second category shows a pentagonal ring which is a signature of an icosahedral growth. The ground state geometries of $`n`$=5, 6, and 10 clusters (figures 1d, 1e, 2d) show neither a quinted roof nor a pentagonal ring. In what follows, we discuss the structure and stability of these clusters. The stability is discussed via the binding energy (Eb), the dissociation energy ($`\mathrm{\Delta }`$E) and the second difference in the total energy with respect to a single Al-Li pair ($`\mathrm{\Delta }^2`$E), and the energy gap between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO). These quantities are plotted in figure 3a (Eb), 3b ($`\mathrm{\Delta }`$E, $`\mathrm{\Delta }^2`$E), and 3c (HOMO-LUMO gap) as a function of the number of Al atoms in the cluster and are defined as

$`\mathrm{E}_b[\mathrm{Al}_n\mathrm{Li}_n]=\left[\mathrm{E}[\mathrm{Al}_n\mathrm{Li}_n]+n\left(\mathrm{E}[\mathrm{Al}]+\mathrm{E}[\mathrm{Li}]\right)\right]/2n`$ (1)
$`\mathrm{\Delta }\mathrm{E}[\mathrm{Al}_n\mathrm{Li}_n]=\left[\mathrm{E}[\mathrm{Al}_n\mathrm{Li}_n]\mathrm{E}[\mathrm{Al}_{n1}\mathrm{Li}_{n1}]\right]\mathrm{E}[\mathrm{Al}\mathrm{Li}]`$ (2)
$`\mathrm{\Delta }^2\mathrm{E}[\mathrm{Al}_n\mathrm{Li}_n]=\left[\mathrm{E}[\mathrm{Al}_{n1}\mathrm{Li}_{n1}]2\mathrm{E}[\mathrm{Al}_n\mathrm{Li}_n]\right]\mathrm{E}[\mathrm{Al}_{n+1}\mathrm{Li}_{n+1}]`$ (3)
Let us recall that a pronounced maximum in the dissociation energy $`\mathrm{\Delta }`$E, along with a corresponding minimum in $`\mathrm{\Delta }^2`$E, signifies stability. In addition, in figure 3d, we have plotted the smallest bond distances of Al-Al, Al-Li and Li-Li as a function of the number of Al atoms in the cluster.

The ground state geometry of Al2Li2 is a bent rhombus (figure 1a(i)). The planar geometry (figure 1a(ii)) is the low-lying structure of the Al2Li2. This cluster is one of the most stable clusters of AlnLin. A substantial rise in the binding energy (figure 3a) is seen when two Al-Li dimers combine to form Al2Li2 cluster. Moreover, a peak in the dissociation energy ($`\mathrm{\Delta }`$E) (figure 3b), a minima in the second difference in energy ($`\mathrm{\Delta }^2`$E) (figure 3b) and a rather large HOMO-LUMO gap (figure 3c) also signifies a relatively high stability of this cluster. Interestingly, the planar Al2Li2 structure, with some bending, has been observed by Kumar  in the magic cluster Al10Li8. However, because of his restricted nature of search, he did not find this 8 valence electron system to have particular stability.

A similar quinted roof structure is also seen in the ground state structures of the clusters Al3Li3 and Al4Li4(figures 1b(i) and 1c(i)). However, the Al3Li3 cluster is not energetically as stable as Al2Li2, as seen from figures 3a, 3b, and 3c. The lowest energy structure of Al4Li4 shows two units of Al2Li2 (figure 1c(i)). These two units are arranged so as to form a rectangle out of the four Al atoms. Indeed, it is this rectangular shape that gives rise to an antiaromatic nature in this cluster . This structure has been extensively studied by us  and others . The two low-lying structures of Al4Li4 (figure 1c(ii) and 1c(iii)) show a square-plane and a quinted roof of four Al atoms, respectively.

In Al5Li5 cluster, a square-plane of Al4 is seen, with the fifth Al atom capping this plane forming a pyramid (figure 1d). Such a square-planar geometry of Al4 has been observed in the all-metal aromatic compounds Al$`{}_{}{}^{2}{}_{4}{}^{}`$ and $`M`$Al$`{}_{}{}^{}{}_{4}{}^{}`$ ($`M`$ = Li, Na, and Cu) . Al5Li5 shows the first three dimensional structure of Aln in AlnLin clusters. Recall that the appearance of such three dimensional geometry in pure Al clusters  is seen for Al6. The Al5Li5 cluster is the most stable among the clusters studied. This is clear from the binding energy, $`\mathrm{\Delta }`$E, $`\mathrm{\Delta }^2`$E, and the HOMO-LUMO gap plots shown in figures 3a-3c, respectively. An interesting aspect of the Al5Li5 cluster is the formation of a plane composed of four Li atoms with one Al atom at the center, akin to the face of a fcc cell. Similar structure for Al5Na5 was predicted by Dhavale  et al. A complete fcc structure can be formed by fourteen atoms: six Al atoms at the face centered sites and eight Li atoms at the vertices of the cube. Such a structure was first studied by Shah  et al. using DBMD. They found this structure to be the lowest energy structure of Al6Li8. However, our investigation shows that this fcc structure, with slight distortion, is 1.11eV higher than that the most stable structure. The lowest energy structure and the fcc structure of this cluster is shown in figure 1f(i) and 1f(ii), respectively. It can be noted that these two structures are nearly similar, the main difference being the elongation of one of the Al4 square-planes in the fcc structure into a rectangle. It is interesting to note that the series Al$`{}_{}{}^{2}{}_{4}{}^{}`$, Li-Al$`{}_{}{}^{}{}_{4}{}^{}`$, Al5Li5, and Al6Li8 shows a three dimensional octahedral growth of Al cluster surrounded by Li atoms. This results in an effective delocalization of the electron density (see section III B).

The Al6Li6 cluster also shows a geometry similar to that of Al5Li5, with distortions in the Al4 square-plane (figure 1e(i)). In the first low energy structure (figure 1e(ii)), a tetrahedra of Al4 is seen, while the other low lying structure (figure 1e(ii)) shows an octahedra of Al6. The octahedral structure is the ground state geometry of pure Al6 cluster .

The appearance of a pentagonal ring, which is a precursor to the icosahedral growth, is seen for $`n=7`$ (Al13 is a distorted icosahedron in its ground state ). The ground state geometry of Al7Li7 shows a pentagonal bipyramid which is capped by the Li atoms (figure 2a(i)). Appearance of a similar pentagonal ring in pure Al clusters (in neutral as well as singly charged form) is seen for Al9  as well as in tetravalent atom clusters viz. Sn, and Pb  for seven atom clusters. Incidentally, these clusters are isoelectronic to Al7Li7 with 28 valence electrons. The low energy structure of Al7Li7 (figure 2a(ii)) as well as the ground state geometries of $`n=8`$, 9, and 11 clusters (figure 2b(i), 2c(i) and 2e(i), respectively) also shows such pentagonal rings. However, the lowest energy structure of Al10Li10 cluster (figure 2d(i)) is completely different from that of the Al7Li7 to Al9Li9 and Al11Li11, in that, no signature of the pentagonal ring is seen. In fact, the geometry is similar to that of Sn10, i.e. a tetracapped trigonal prism (TTP)  with some distortions. This TTP structure is enclosed in a Li cage resulting in its distortion. It also shows magical behavior, since it exhibits a peak in the dissociation energy ($`\mathrm{\Delta }`$E), a minima in the second difference in energy ($`\mathrm{\Delta }^2`$E), and a large HOMO-LUMO gap (figure 3a-3c, respectively).

We have performed an analysis of the interatomic bond distances in order to understand the mixing and segregation behavior of the Al and the Li atoms. In figure 3d, we plot the smallest bond distances of Al-Al, Al-Li and Li-Li as a function of number of Al atoms in the AlnLin cluster. We find that in smaller clusters ($`n4`$) the Al-Al bond is dominant. However, from $`n=5`$ onwards, the Al-Li bond distance is comparable to that of Al-Al. Recall that at $`n=5`$, the first three dimensional Aln structure emerges, thereby increasing the surface area. This helps the Li atoms to maximize the number of Al-Li bonds and thereby move closer to the Al atoms.

It is well known that the density based method is not as accurate as the Kohn-Sham method, nonetheless, it has been extensively used to investigate the geometry, and stability of various Al-Li clusters . Here, we make some pertinent comments on the difference in the geometry of the AlnLin clusters obtained by these two methods. Our results show that the geometries of some of the clusters were in good agreement with those predicted by DBMD. For instance, the structure of Al2Li2 was found to be a quinted roof by both the methods (see ref. VI for DBMD results). The geometries of other clusters like Al7Li7 to Al9Li9 show some substitutional disorder. However, there are certain clusters where the geometries were quite different by the KS method. For instance, the structures of, Al3Li3, predicted by DBMD, is a capped trigonal bipyramid . When optimized by KS, this structure gets Jahn-Teller distorted into a structure with a quinted roof of Al2Li2. The evolution in the geometry by DBMD method shows an earlier appearance of three dimensional structure of Aln at $`n=4`$ (Al4 is tetrahedra in Al4Li4), as compared to the pyramidal structure of Al5 in Al5Li5 by KS method.

### B Bonding

#### 1 Al2Li2

As noted above, Al2Li2 is one of the most stable cluster. This cluster has 8 valence electrons, which corresponds to a closed shell in the jellium sense. However, the behavior of the MOs as well as the eigenvalue spectrum does not resemble to that of the SJM. Kumar  found that the bonding between Al-Li in the planar Al2Li2 structure is ionic. On the contrary, our investigation on the bent rhombus structure does not show signature of any significant charge transfer from the Li atoms to the Al atoms. However, an analysis of the difference charge density shows that, a charge transfer from the Li atoms to the bonding region between the Al and Li atoms takes place (figure not shown). In order to get better insight of the bonding in Al2Li2, we have studied various clusters viz. Al2, Al$`{}_{}{}^{}{}_{2}{}^{}`$, Al$`{}_{}{}^{2}{}_{2}{}^{}`$, Al2Li, Al2Li-, and Al2Li2. This series represent a growth of Al2Li2 from the Al2 dimer. In figure 4, we show the eigenvalue spectrum, and the Al-Al bond distance for the optimized geometries of these clusters. The eigenvalue spectrum of Al2 (figure 4a), in neutral, singly and doubly charged form, shows a triply degenerate HOMO state. These states are partially filled two $`\pi `$ and one $`\sigma `$ bond (figure not shown). Successive substitution of each electron in the charged dimer, by Li atoms splits the HOMO state reducing the degeneracy. The first Li atom forms a $`\sigma `$ bond with the $`p_yp_y\pi `$ orbital of the Al2 dimer (isodensity surfaces not shown). The second Li atom further splits the HOMO state introducing a substantial gap of 0.68eV. However, in this case, a hybridization of the $`s`$ orbital of Li atom and the $`\pi `$ orbital of Al2 takes place. In figure 5, we show the second HOMO and the HOMO state of the Al2Li2 cluster. Clearly, these states shows the $`\pi `$ bonded orbitals of Al2 hybridized with the $`s`$ orbitals of the Li atoms. In figure 4b, we plot the Al-Al bond distance in these clusters. A contraction of the Al-Al bond upon addition of an electron or a Li atom is seen. The contraction of the bond on addition of a Li atom being more than that due to the addition of an electron. Thus, the hybridization of the $`s`$ orbitals of Li with the $`\pi `$ orbitals of Al2 enhances its stability.

#### 2 Al4Li4

The bonding in Al4Li4 cluster has been extensively studied by us  and others . This cluster has been found to be antiaromatic with four $`\pi `$-electrons. For the sake of completeness, we discuss the bonding in this cluster. This cluster, as discussed earlier, is composed of two Al2Li2 units arranged edge-to-edge to form a Al4 rectangle. This rectangular structure leads to antiaromaticity. In this cluster $`sp^2`$ hybridization of Al takes place leaving one empty unhybridized $`p`$ orbital. The valence electron of each of the four Li atoms is then transfered to this empty $`p_z`$ orbital, thus providing four $`\pi `$-electrons for antiaromaticity. Interestingly, its first low-lying structure (figure 1c(ii)) can also be considered as a candidate for antiaromaticity. An analysis of the MOs shows that this cluster also has four $`\pi `$-electrons similar to that of the lowest energy structure, thereby conforming the structural and electron count criteria for antiaromaticity. However, a detailed analysis of the magnetic field induced ring currents has to be done in order to understand the antiaromaticity in this structure.

#### 3 Al5Li5 and Al6Li8

As discussed in section III A, the structures of Al5Li5 and Al6Li8 (in fcc geometry) clusters show a square-plane of Al4. Such Al4 square-plane in Al$`{}_{}{}^{2}{}_{4}{}^{}`$ and Li-Al$`{}_{}{}^{}{}_{4}{}^{}`$ has lead to aromaticity in these clusters . Moreover, these clusters (Al$`{}_{}{}^{2}{}_{4}{}^{}`$, Li-Al$`{}_{}{}^{}{}_{4}{}^{}`$, Al5Li5, and Al6Li8) show a three dimensional octahedral growth of the Al cluster. It is, therefore, interesting to discuss the change in bonding, and hence the effect on aromaticity in these clusters. Recall that Al$`{}_{}{}^{2}{}_{4}{}^{}`$ and Li-Al$`{}_{}{}^{}{}_{4}{}^{}`$ are aromatic clusters with two $`\pi `$-electrons . The HOMO state in these clusters are completely delocalized $`\pi `$ orbitals. There are two more delocalized $`\sigma `$ bonds: one composed of radial $`p`$ orbitals, and the other composed of perpendicular $`p`$ orbitals. Due to the presence of a Al4 square-plane, the bonding in Al5Li5 and Al6Li8 clusters is expected to be similar to that in Al$`{}_{}{}^{2}{}_{4}{}^{}`$. In that case, these clusters would have two $`\pi `$-electrons similar to that of Al$`{}_{}{}^{2}{}_{4}{}^{}`$ and Li-Al$`{}_{}{}^{}{}_{4}{}^{}`$, and show some signature of aromaticity. In figure 7a and 7b, we show the HOMO states of Al5Li5 and Al6Li8. Interestingly, the HOMO state of Al5Li5 does show such a delocalized $`\pi `$ bond between the Al4 square (figure 7a). However, a lone electron is seen at the fifth Al atom that caps this square. This reduces the number of $`\pi `$ electrons count to one, thereby losing the aromaticity. The HOMO state of Al6Li8, on the other hand, is quite different than that seen in the previous clusters (see figure 7b). It shows a localized bond composed of $`p`$ orbitals of the Al atoms. In figure 6, we show the eigenvalue spectrum of Al6Li8, along with that of the other AlnLin ($`n=`$1-11) and Sn10 clusters. It is clear from this that the eigenvalue spectrum of Al6Li8 shows nearly a jellium-like structure. However, an analysis of the behavior of the molecular orbitals (figure not shown) of Al6Li8 shows that the lower six states are jellium-like, whereas, the higher states show localized bonds arising out of $`p`$ orbitals of the Al atoms. The isodensity surfaces for the 6th state of Al5Li5 and the 10th state of Al6Li8 are shown in figure 7c, and 7d, respectively. These are $`\sigma `$-bonded states, composed of radial $`p`$ orbitals and are analogous to the radial $`p`$ orbitals of Al$`{}_{}{}^{2}{}_{4}{}^{}`$, as noted by Wang et al . Since, these bonds involve all the Al atoms, they exhibit a substantial delocalization.

#### 4 AlnLin ($`n`$=7-11)

The bonding in these clusters (AlnLin, $`n=711`$), is quite different from that of the smaller clusters. The analysis of the behavior of the MOs shows that the eigenvalue spectrum can be classified into two groups: the lower $`n`$ states are jellium-like, whereas, the higher states forms a complex band of localized bonds arising out of the interaction of the $`p`$ electrons of each of the Al atoms. The lower states can be identified as 1s, 1p, 1d, 2s,…, in conformation to the SJM. Typical representative isodensity plots of these states are shown in figure 8. They represent the 1s, 1p, 1d and 2s states of the jellium composed of the Al-$`3s`$ orbitals. In figure 6, we show the eigenvalue spectrum of AlnLin ($`n=`$1-11) and Sn10 clusters. This figure shows a clear energy separation of lower jellium-like $`n`$ states of the clusters with $`n>6`$, from the higher orbitals. In order to study this more extensively, we have performed a spherical harmonics analysis of the KS orbitals . We find that the lower $`n`$ orbitals have a predominant $`s`$ character (75-85%) indicating that these state are arising out of the $`3s`$ orbitals of the Al atoms, which hybridizes to form a jellium. However, this description cannot be extended to the higher states. The higher states are composed of Al-$`3p`$ orbitals.

The localization characteristics in the bonding in the Al-Li clusters is also analyzed via the electron localization function plots. Previous investigations show that even in large Al-Li clusters, the Al-Al bond is covalent. However, our earlier work on the structural properties shows a clustering of Al atoms enclosed in Li cage. This clustering could lead to a delocalization of the electron density. In order to study this, we have analyzed the electron localization function in the AlnLin clusters. In figure 9, we show the isovalued surface of the ELF for the clusters Al2Li2, Al5Li5, and Al10Li10, at the values 0.85, 0.8, and 0.72, respectively. These plots shows a localization of the electron density along the Al-Al bond (increasing the value of ELFs does not show any bond between Al-Al). It has been noted by Silvi et al , that an ELF value of about 0.7 or greater is a indication of a localized bond in that region. Thus, the Al-Al bonds in these clusters are predominantly covalent, in accordance with the previous observation. However, as the cluster grows in size, the degree of localization of the electron density in these clusters is reduced. Further, the localization in Al5Li5 is only in the plane forming the square (for out of plane, the Al-Al bond is seen for the ELF value of 0.7). A similar behavior of the ELF for Al7Li7 (not shown in figure) is seen, in that, a localized bond along the pentagon (ELF value = 0.78) and delocalized bond in the perpendicular direction (ELF value = 0.6) is seen.

## IV Tetravalent behavior of Al in AlnLin clusters

Previous investigations on mixed Al-Li clusters , have shown a substantial charge transfer from the Li atom to the Al atom. Due to such charge transfer, Al behaves as a tetravalent atom in mixed Al-Li clusters with nearly 50:50 percent concentration. In this section, we present the evidence of such behavior. The analysis of the geometries of Aln in AlnLin shows remarkable similarities to those of the clusters of the tetravalent atoms viz., Si, Ge, Sn and Pb. The transition from planar to three dimensional geometry at $`n=5`$, the appearance of a pentagonal ring at $`n=7`$, and the formation of the distorted TTP structure of 10 atom cluster are some of the features of the clusters of tetravalent atoms. Such behavior of Aln in AlnLin clusters is also observed. Moreover, the eigenvalue spectrum, in figure 6, of the Al10Li10 show remarkable similarity to that of the Sn10 cluster. The bonding between Al-Al in these clusters is covalent, which is similar to that in the clusters of tetravalent atoms. Finally, the resemblance of the DOS of the Al-Li in its most stable bulk phase i.e. B32, to that of diamond  also indicates a tetravalent behavior of the Al atom in Al-Li systems with 50:50 concentration.

## V Conclusion

In the present work, we have reported the systematic investigation of the geometries, energetics, electronic structure and bonding in various Al-Li clusters viz. AlnLin ($`n`$=1-11), Al$`{}_{}{}^{}{}_{2}{}^{}`$, Al$`{}_{}{}^{2}{}_{2}{}^{}`$, Al2Li, Al2Li-, and Al6Li8 using BOMD within the framework of density functional theory. The bonding in these clusters was discussed via the electron localization function as well as the behavior of the molecular orbitals. The growth pattern is divided into two broad categories: first consisting of a bent rhombus of Al2Li2 ($`n`$=2-4) and second consisting of a pentagonal ring ($`n`$=7-9,11). We find that the 8, 20, and 40 valence electron systems to be magic, exhibiting a peak in the dissociation energy, a minima in the second difference in energy, and a large HOMO-LUMO gap. The structural transition of Aln in AlnLin, from two dimension to three dimension increases the surface area, thereby helping the Li atoms to maximize the number of Al-Li bonds, as a result of which the Li atoms move closer to the Al atoms. In clusters with more than six Al atoms, a charge transfer from Li atoms to Al atoms makes Al behave as a tetravalent atom like Si, Ge, Sn, and Pb in AlnLin clusters. These negatively charged Aln structure is then stabilized by the positive Li environment. The behavior of the MOs in these clusters ($`n>6`$) can be divided two groups: a lower group of jellium-like states arising out of the $`s`$ electrons of the Al atoms and a higher group of localized bonds formed through the interaction of the $`p`$ electrons on each of the Al atoms. The formation of the three dimensional Aln structure destroys the aromatic and antiaromatic nature in these clusters.

## VI Acknowledgement

SC gratefully acknowledges the financial support of CSIR (New Delhi). VVP wishes to thank NSERC of Canada and the University of Pune for partially supporting this research.

List of Figures

1. The ground state and the low-lying geometries of the AlnLin ($`n`$=2-6) and Al6Li8 clusters. The black circles represent Al atoms and the white circles represent the Li atoms. The lowest energy structure is represented by (i).
2. The ground state and the low-lying geometries of the AlnLin ($`n`$=7-11) clusters. The black circle represents Al atoms and the white circles represent the Li atoms. The lowest energy structure is represented by (i).
3. (a) The binding energy per atom (in eV/atom) of the AlnLin cluster ($`n`$=1-11). (b) The dissociation energy ($`\mathrm{\Delta }`$E) and the second difference in the total energy ($`\mathrm{\Delta }^2`$E) of the AlnLin clusters ($`n`$=2-10) with respect to a Al-Li dimer. (c) The HOMO-LUMO gap (in eV) of the AlnLin cluster ($`n`$=1-11). (d) The minimum interatomic distance Al-Al, Al-Li and Li-Li (in Å) of the AlnLin cluster ($`n`$=1-11).
4. (a) The eigenvalue spectrum (in eV). The continuous lines correspond to occupied states and dashed lines correspond to the empty states. The numbers on the right indicate the degeneracy of the states. The numbers on the left indicate the occupancies of the states. All other occupied states have two electrons. (b) The Al-Al bond distance (in Å).
5. The isodensity surface of the second HOMO and the HOMO state of Al2Li2 at 1/5th of its maximum value. The black circles represent Al atoms and the gray circles represent the Li atoms.
6. The eigenvalue spectrum (in eV) of AlnLin ($`n`$=1-11), Al6Li8, and Sn10 clusters. The continuous lines correspond to the occupied states and the dotted lines correspond to the empty states. The numbers on right indicate the degeneracy of the states. All states are doubly occupied. The thick dashed line shows the $`n^{\mathrm{th}}`$ state of these clusters, where $`n`$ is the number of Al atoms. For Sn10 the thick dashed line corresponds to the 10th state.
7. (a) The isodensity surface of the HOMO state of Al5Li5, at 1/5th of its maximum. (b) The isodensity surface of the HOMO state of Al6Li8, at 1/5th of its maximum. (c) The isodensity surface of the $`6^{\mathrm{th}}`$ state of Al5Li5, at 1/5th of its maximum. (d) The isodensity surface of the $`10^{\mathrm{th}}`$ state of Al6Li8, at 1/8th of its maximum. In all the figures, the black circles represent Al atoms and the gray circles represent the Li atoms.
8. The isodensity surface of various orbitals. These are: (a) 1st orbital of Al7Li7 at 1/5th of its maximum value, (b) 3rd orbital of Al8Li8 at 1/5th of its maximum value, (c) 8th orbital of Al9Li9 at 1/10th of its maximum value, and (d) 10th orbital of Al10Li10 at 1/5th of its maximum value, representing the 1s, 1p, 1d, and 2s state of the jellium composed of the Al-3$`s`$ orbital. In all the figures, the black circles represent Al atoms and the gray circles represent the Li atoms.
9. The isodensity surface of the electron localization function (ELF) of: (a) Al2Li2 at the value 0.85, (b) Al5Li5 at the value 0.80, and (c) Al10Li10 at the value 0.72. In all the figures, the black circles represent Al atoms and the gray circles represent the Li atoms.


format: text

text:
# Timing and spectral studies of LMC X-4 in high and low states with Beppo-SAX : Detection of pulsations in the soft spectral component

## 1. Introduction

LMC X-4 is an eclipsing, accretion-powered, binary X-ray pulsar in the Large Magellanic Cloud orbiting around a 20 M O7 III-V companion with an orbital period of $``$1.4 day (Kelley et al. 1983, Ilovaisky et al. 1984) and pulse period of $``$13.5 s (White (1978); Li, Rappaport & Epstein (1978); Lang et al. (1981)). X-ray eclipses with a 1d.4 recurring period were discovered by Li et al. (1978) and White (1978). X-ray pulsations with a period of 13.5 s were first detected in LMC X-4 from SAS 3 observations in 1976 and 1977 only during the infrequent flaring events by Kelley et al. (1983), who also derived the orbital solution from pulse arrival delay measurements. During a high state in 1983, EXOSAT observations detected the pulsations both during flaring and non-flaring episodes (Pietsch et al. (1985)). Using subsequent pulse timing measurements with GINGA, ROSAT (Levine et al. (1991); Woo et al. (1996)) and RXTE (Levine et al. (2000)), an orbital period decay with a time scale of $`10^6`$ yr-1 was measured, which is similar to the period decay rate in other high mass accreting binary pulsars (Cen X-3, SMC X-1 etc.).

Analogous to the well known X-ray pulsar Her X-1, LMC X-4 exhibits a periodic long term intensity variation at 30.5 day (Lang et al. (1981)). The super-orbital period has recently been found to be decreasing (Paul & Kitamoto 2002). In the case of LMC X-4, the X-ray flux varies by two orders of magnitude between the high and low intensity states of the super-orbital period. Varying obscuration by the accretion disk due to precession provides a good explanation for the long term periodic intensity variations. The 13.5 s X-ray pulsation is not yet detected during the low state of the super-orbital period indicating that most of the radiation observed in low state is probably reprocessed emission from the stellar wind of the companion star. The X-ray spectrum during eclipse of the neutron star, which consists only of scattered/reflected radiation was found to be nearly identical in shape and strength throughout the super-orbital period, supporting the disk obscuration hypothesis (Woo et al. (1995)). However, observations with the RXTE (Heindl et al. (1999)) did not detect any significant change in the absorption column density in different phases of the super-orbital period. One possibility that was pointed out is that in low state, part of the X-ray emission from the central source is completely blocked by the disk and the rest is observed with the Galactic absorption. It should be noted that, for sensitive measurements of absorption column density, observations with good soft X-ray spectral capability are required.

The X-ray spectrum of LMC X-4 consists of, in addition to a hard power-law component, a soft excess below 1 keV and an iron emission line at 6.4 keV (Woo et al. 1996). Iron K shell emission lines in X-ray pulsars are usually narrow, but the equivalent width can sometimes be as large as $``$1 keV or more. These lines are produced by illumination of neutral or partially ionized material in accretion disk, stellar wind of high mass companion star, material in the form of circumstellar shell, material in the line of sight, or accretion column. Expectedly, the iron emission line parameters are often found to be correlated to the continuum hard X-ray flux, absorption column density etc. Using a large number of RXTE/PCA observations of LMC X-4 made in different phases of the super-orbital period, a correlation was found between the line flux and the continuum flux, along with highly variable nature of the equivalent width of the line during low state (Naik & Paul (2003)).

It is known that X-ray pulsars those do not suffer from strong absorption by material in the line of sight show presence of soft X-ray excess above the extrapolated hard power-law component (Paul et al. (2002), and references therein). In some pulsars, the soft excess is also found to pulsate, often with a phase difference with respect to the hard component. Origin of the pulsating soft component is not yet understood clearly, at least for the bright pulsars in LMC and SMC. In LMC X-4, a small Galactic absorption column density and lack of local absorbing material allows us to probe the soft X-ray emission with some detail. Observations of LMC X-4 with ASCA detected the soft excess but were found to be inconclusive regarding the nature of the soft component (Paul et al. (2002)). The latter is due to a small pulse fraction of this pulsar and also to the fact that the ASCA spectrometers are not sensitive at energies where the soft excess dominates the power-law component in LMC X-4.

We have carried out detailed timing and spectral analysis of two observations of LMC X-4 with the narrow field instruments of Beppo-SAX during one high and one low intensity states of LMC X-4. We aimed to detect pulsations in the Beppo-SAX observation carried out in the high state and to determine the mid-eclipse time during this observation. From pulse phase averaged spectral studies in the two states, we investigate relations of the iron line parameters with the continuum flux. To investigate the nature of the soft excess, pulse-phase-resolved spectral analysis has been carried out for the high state observation. In the subsequent sections we give details of the observations, the results obtained from the timing and spectral analysis, followed by a discussion on the results obtained from these two Beppo-SAX observations.

## 2. Observations

The 2–12 keV X-ray light curve of LMC X-4 taken with the RXTE-ASM for last six years, is shown in Figure 1, folded at the present super-orbital period of 30.276 days (Paul & Kitamoto 2002). Two observations of LMC X-4 were made with the Beppo-SAX narrow field instruments during 1997 March 13 10:48 to March 15 08:17 (UT) and 1998 October 20 22:39 to October 22 08:05 (UT) in the low and high states of the source respectively. The durations of these two observations with respect to the phase of the long period are also shown in Figure 1. We have used archival data from these two observations in the present work.

For the present study, we have used data from the Low and Medium Energy Concentrator Spectrometers LECS and MECS and the Phoswich Detector System PDS on-board Beppo-SAX satellite. In the 1997 observation, the hard X-ray instrument PDS with a larger field of view contained mostly photons from a nearby transient source EXO 053109–6609.2, which is only 0.5 degree away from LMC X-4 (Burderi et al. 1998). The broad band spectrum (0.1–100 keV) from the 1998 Beppo-SAX observation has been reported by La Barbera et al. (2001). To obtain energy resolved pulse profiles in hard X-rays, we have used PDS data from the 1998 observation.

The MECS consists of two grazing incidence telescopes with imaging gas scintillation proportional counters in their focal planes. The LECS uses an identical concentrator system as the MECS, but utilizes an ultra-thin entrance window and a drift less configuration to extend the low-energy response to 0.1 keV. Time resolution of the instruments during these observations was 15.25 $`\mu `$s and energy resolutions of LECS and MECS are 25% at 0.6 keV and 8% at 6 keV respectively. The PDS consists of a square array of four independent NaI(Tl)/CsI(Na) phoswich scintillation detectors. The energy resolution of PDS instrument is 15% at 60 keV. The time resolution of the PDS instruments during 1998 observation was 0.25 ms. For a detailed description of the Beppo-SAX mission, we refer to Boella et al. (1997).

## 3. Timing Analysis

The light curves of LMC X-4 during these two Beppo-SAX observations obtained with the LECS and MECS detectors are shown in Figure 2 with respect to the orbital phase (mid-eclipse time = orbital phase 0.5). During the 1997 observation, which was in a low state, a smooth intensity variation is observed throughout the orbital period and the eclipse is not clearly seen. Near the end of the low state observation, an increase in count rate is noticed which is similar to the flares observed earlier in the low state with ROSAT (Woo et al. 1995). Similar to the other high state observations of LMC X-4, the 1998 observation clearly shows an eclipse during which the X-ray intensity decreases sharply to a very low level.

For timing analysis, the arrival times of the photons were first converted to the same at the solar system barycenter. Light curves with time resolution of 0.25 s were extracted from circular regions of radius $`4\mathrm{}`$ around the source. The pulses from LMC X-4 loose coherence within a relatively short time scale of few thousand seconds due to a short orbital period of 1.4 days while the semi-amplitude of the arrival time delay due to orbital motion is 26.3 s. Additionally, the relatively small effective area of the Beppo-SAX telescopes and small pulse fraction of LMC X-4 makes it difficult to detect pulse arrival times from small segments of the light curve. Therefore, to detect the pulsations and also the mid-eclipse times during these observations, we have first corrected the light curve for the binary motion. The orbital parameters, namely the semi-amplitude was taken to be 26.31 s and the mid-eclipse time was derived from the quadratic solution given by Levine et al. (2000). The binary arrival time delay correction was applied with a range of trial mid-eclipse times around the value extrapolated from Levine et al. (2000). Pulse folding and $`\chi ^2`$ maximization method was applied to all the corrected light curves. The distribution of maximum $`\chi ^2`$ against the trial mid-eclipse times obtained from each pulse folding analysis is shown in Figure 3 for the high state observation in 1998. The maximum $`\chi ^2`$ distribution has a Gaussian profile around the expected value of the mid-eclipse time, the center of which gives the correct value. Using the same method we did not detect any pulsations in the low state. The result of the pulse folding analysis corresponding to the peak of the curve in Figure 3 is shown in Figure 4, which clearly shows the detection of pulsations in this observation contrary to La Barbera et al. (2001). In the high state, the pulsations were in fact detected independently in light curves of the LECS, MECS and PDS detectors. The pulse profiles obtained from the high state light curves of the LECS, MECS, and PDS in different energy bands are shown in Figure 5. From this analysis, we have derived the pulse period of LMC X-4 to be 13.50260(12) s and the mid-eclipse time MJD 51106.6399(25) corresponding to orbit number 5877.

Combining the new determination of mid-eclipse time with the previous measurements (Levine et al. 2000 and references therein, Paul et al. 2002) we fitted a second-order polynomial to derive the orbital decay rate. Residuals of the orbital mid-eclipse time history of LMC X-4 relative to a constant period are shown in Figure 6. The solid line in the figure shows the quadratic nature of the mid-eclipse time history. We derive a period decay rate of $`\frac{\dot{\mathrm{P}}_{\mathrm{orb}}}{\mathrm{P}_{\mathrm{orb}}}=(9.89\pm 0.05)\times 10^7`$ yr-1.

The pulse profiles in different energy bands are shown in Figure 5. In the low energy band (0.1–1.0 keV of LECS, top panel), it is nearly sinusoidal and has a complex structure in the energy band of 4.0-10.0 keV (MECS, bottom panel) with multiple dips superposed on a smooth sinusoidal profile. A phase difference between the sinusoidal pulse components of the two energy bands is also visible in the figure. At intermediate energies (1.0–4.0 keV, both LECS and MECS), the pulse profile is a mixture of the above two. The pulse fraction is less than 10% and these features are similar to the known pulsation properties of LMC X-4. Though the pulsations are seen in the PDS light curves in 15.0–60 keV energy band (top-right three panels of Figure 5), the pulse profiles do not show the complex dipping feature, as seen in the MECS profiles, and are similar to those obtained from the LECS light curves with certain phase difference. The light curve above 60 keV is mainly background dominated and we did not detect pulsations in the 60-200 keV range. The transient X-ray pulsar EXO 0531-66 which has a pulse period of 13.7 s, close to that of LMC X-4 was in the field of view of the PDS instrument. This source was detected in the LECS and MECS instruments at an intensity level of about 0.2% of LMC X-4. From a period search of the PDS light curve in a range covering the pulse period of both LMC X-4 and EXO 0531-66, we have verified that contamination of the PDS folded light curve of LCM X-4 was minimal.

## 4. Pulse phase averaged spectroscopy

For spectral analysis, we have extracted LECS spectra from regions of radius $`6\mathrm{}`$ centered on the object (the object was at the center of the field of view of both the instruments). The combined MECS source counts (MECS 1+2+3 in the 1997 observation and MECS 2+3 in the 1998 observation) were extracted from circular regions with a $`4\mathrm{}`$ radius. For spectral fitting, the September 1997 LECS and MECS1 response matrices were used. Background spectra for both LECS and MECS instruments were extracted from appropriate source-free regions of the field of view with extraction region on the detector similar to the source extraction regions. Some rebinning was done to allow the use of $`\chi ^2`$-statistic. Events were selected in the energy ranges 0.1–4.0 keV for LECS and 1.65–10.0 keV for MECS where the instrument responses are well determined. Combined spectra from the LECS and MECS detectors, after appropriate background subtraction, were fitted simultaneously. All the spectral parameters, other than the relative normalization, were tied to be the same for both detectors and the minimum value of hydrogen column density $`N_H`$ was set at the value of Galactic column density in the source direction.

### 4.1. Low state

In the low state observation of 1997, the source count rate is about two orders of magnitude lower than the high state. Therefore, spectral fitting required significant rebinning. The spectrum, when fitted to a single power-law model with line of sight absorption, showed significant soft excess below 1 keV. Addition of a soft component in the model improves the spectral fit; we have tried two different components such as black-body emission and bremsstrahlung component to fit the soft excess. Results of the spectral fits are given in Table 1. From the spectral fitting, we are unable to distinguish between these two models of the soft component. The low state photon spectrum along with a spectral model comprising of three components, a hard power-law, a Gaussian emission line and a bremsstrahlung emission are shown in Figure 7. Two notable features of the low state spectrum are that the power-law component is very hard with a photon index of $``$ 0.1, and the iron line equivalent width is large, $``$ 1.3 keV. Both of these have been observed earlier with the RXTE (Naik & Paul (2003)). The line of sight absorption during low intensity state is found to be $``$ 5 $`\times `$ 1020 atoms cm-2, similar to the Galactic column density, and the characteristic temperature of the soft component is lower than that in the high state.

### 4.2. High state

A power-law fit to the high state spectrum shows a large soft X-ray excess. Similar to the low state, we tried to fit the soft excess with different single components : black-body, thermal bremsstrahlung and soft power-law, none gives a satisfactory fit. Finally we found that a combination of any two of these components for the soft excess improves the spectral fit to some extent. The best-fit model is the one comprising a hard power-law of photon index 0.65, an iron Kα emission line of equivalent width 240 eV, a black-body component of temperature 0.15 keV and a soft power-law with photon index $``$ 3. The LECS and MECS count rate spectra are shown in Figure 8 along with contributions of individual components. Irrespective of the model chosen to fit the soft excess, the soft component starts dominating the spectrum at energies below 0.8 keV.

## 5. Pulse phase resolved spectroscopy

Since we have not detected any pulsations in the low state observation of 1997, the pulse phase resolved spectroscopy is done only on the 1998 data set to understand the nature of the soft component. The photon arrival times in the LECS and MECS event files were corrected for the solar system barycenter and for the arrival time delays due to orbital motion. Following this, spectra were accumulated into 16 pulse phases by applying phase filtering in the FTOOLS task XSELECT. As in the case of phase-averaged spectroscopy, the background spectra were extracted from source free regions in the event files and appropriate response files were used for the spectral fitting.

Each pulse phase resolved spectra has much inferior signal to noise ratio compared to the phase averaged spectrum. Therefore, sometimes it is not possible to constrain all the model parameters, specially if a complicated spectral model is used. Since our aim is to investigate nature of the soft excess (whether pulsating or not) we used only two models, in which a single component (either a soft power-law or a bremsstrahlung) is used for the soft excess and a relatively low reduced $`\chi ^2`$ is obtained. For the phase resolved spectra, the iron-line energy, line-width and $`N_H`$ were fixed to their phase-averaged values and all the other spectral parameters were allowed to vary. The continuum flux and the fluxes of the soft and hard components in 0.1 – 10.0 keV energy range were estimated for all the 16 phase-resolved spectra. The modulation in the X-ray flux for the hard and soft spectral components and the total flux are shown in Figure 9 along with the 1$`\sigma `$ error estimates. Pulse phase resolved spectral analysis shows that modulation of the hard power-law flux is very similar to the pulse profile at higher energies. A pulsating nature of the soft-spectral component is clearly detected irrespective of the spectral model used. The pulsating soft component has a nearly sinusoidal profile, dissimilar to the complex profile seen at higher energies, with a certain phase difference with the hard component. These properties are very similar to what is seen in Her X-1 (Endo et al. (2000)) and SMC X-1 (Paul et al. (2002)).

## 6. Discussion

### 6.1. Orbital evolution

There are several ways to measure orbital ephemeris of X-ray binaries from X-ray observations. Pulse frequency modulation or arrival time delay measurement of X-ray pulsars gives very accurate measurement of mid-eclipse time (T$`_{\frac{\pi }{2}}`$) of binary X-ray pulsars and has so far been used to measure the orbital parameters of a large number of binary X-ray pulsars. A quadratic nature of the mid-eclipse time history during the last two/three decades unambiguously established an orbital evolution with a time scale of about $`10^510^6`$ yr in several binary X-ray pulsars (Cen X-3, Nagase et al. 1992; LMC X-4, Levine et al. 2000; SMC X-1, Wojdowski et al. 1998 etc.). In case of low mass X-ray binary pulsars, the orbital evolution is mainly due to conservative mass transfer from the companion to the neutron star whereas in case of high mass X-ray binaries, it is either due to mass loss from the system and/or strong tidal interaction between the two stars. The measurement of one new mid-eclipse times of LMC X-4 with Beppo-SAX and two with ASCA (Paul et al. (2002)) allows us to determine the orbital period derivative of the source with better accuracy. It has been pointed out before that tidal interaction is probably the dominating effect in orbital evolution of LMC X-4 ( Levine et al. 2000).

### 6.2. The iron emission line

Spectral analysis of the high and low states of LMC X-4 showed significant differences. Presence of a prominent iron emission line can be seen in the low intensity state of the source (Figure 7). While the continuum flux in 0.1–10 keV band in the two states differ by a factor of $``$ 60, the iron line flux differs only by a factor of $``$ 12. This results in a much larger equivalent width of 1.26 keV in the low state compared to 200 eV in the high state. In view of an obscuring precessing accretion disk model for the super-orbital period, it is interesting to note that the column density is very similar in the two states and are close to the Galactic value. It was suggested (Heindl et al. (1999)) that if the low state X-rays are due to Compton scattering by circumstellar material, significant obscuration by the accretion disk in the low state may still be true. Partial coverage of the hard X-ray emission could be a possible explanation for the low state. However, we found that addition of a partial covering model component to the hard power-law does not fit the low state spectrum.

The spectral results obtained from present work are consistent with those obtained from the RXTE-PCA observations (Naik & Paul (2003)). However, a much better sensitivity for weak sources and better energy resolution of Beppo-SAX LECS and MECS compared to RXTE-PCA instruments gives improved confidence in the equivalent width measurements of the iron emission line. A smaller iron line flux in the low state indicates that most of the line emission is produced in a region comparable to or smaller than the size of the obscuring material, probably the accretion disk. But an increased equivalent width of the iron line in low state also indicates that a part of the emission line must originate in a region further away from the neutron star. It is known that the HMXB pulsars often show several iron line components from different ionization species, probably produced in different regions of the system. The best example is the case of Cen X-3 where at least three components with different variability characteristics were resolved (Ebisawa et al. 1996). The flux evolution of different line components during eclipse egress suggests that the 6.4 keV line is emitted close to the neutron star, while the other components are probably emitted in a more extended highly ionized plasma (Nagase et al. 1992, Ebisawa et al. 1996). In the case of LMC X-4 also, La Barbera et al. (2001) claimed detection of more than one component with one of them being at 6.1 keV. High spectral resolution observations with Chandra or future ASTRO-E mission will allow one to find out whether these two species of iron have different ionization state. High spectral resolution observations in different phases of the super-orbital and orbital periods of LMC X-4, including during the eclipse, will finally help to settle the issue regarding the origin of the iron emission line in LMC X-4.

### 6.3. Pulsations of the soft excess and its origin

A soft excess above the hard power-law component is now known to be present in several accreting pulsars and the soft component is also known to be pulsating in some of these sources. The soft excess has been modeled with several different types of emission, but no single model is applicable to all sources with soft excess. A black-body type pulsating soft excess can describe the emission from Her X-1 well (Endo et al. 2000), however, it runs into problem in case of the more luminous sources like SMC X-1 and LMC X-4 (Paul et. al. 2002). In some sources, a bremsstrahlung type of soft emission describes the spectrum well. However, since the emission region has to be very large and the cooling time scale is also very large, such an emission is not expected to pulsate. The soft excess in LMC X-4 is more complex and more than one component is needed to model the soft excess. However, when a bremsstrahlung component is included, it produces most of the soft excess and the pulsating soft excess requires this component to pulsate. LMC X-4 has a very small pulse fraction, and the soft excess peaks at energies lower than the ASCA threshold. As a result, ASCA observations of LMC X-4 were inconclusive regarding the pulsations of the soft excess in LMC X-4. In the present work with Beppo-SAX, we have obtained definite detection of pulsations in the soft excess irrespective of the spectral model used.

In the low state, the LECS and MECS light curves (Figure 2) show large amplitude but gradual intensity variation over the entire binary orbit. The sharp eclipse, that is easily detectable in the high state, is not clear in the low state at all. This indicates that in low state, the observed X-rays are due to reprocessing from a scattering region that has size comparable to that of the companion star. It is possible that in different orbital phases part of the scattering region is blocked from view by the companion star. The amplitude of variation, i.e, ratio between the maximum and minimum in the low state binary light curve is larger in MECS compared to LECS. This is possible if the soft excess is produced in a larger region and the hard power-law component originated near the neutron star is heavily absorbed with the absorption highly varying with orbital phase. Using a radius of 8.1 R of the companion star (Woo 1993) and a binary separation of 26.31 lt-s (Levine et al. 1991), the solid angle subtended by the companion star at the neutron star is calculated to be $``$4$`\pi `$/7. In the present work, the soft excess flux in the low and high states of LMC X-4 are measured to be about about 18% of the hard power-law flux. Therefore, one probable site for production of the soft excess is the surface of the companion star, which can work as a reprocessing agent. Different visibility of the stellar surface that is irradiated by the neutron star can cause the observed orbital soft X-ray intensity modulation in the low state.

## Acknowledgments

We thank the referee for her/his valuable suggestions that helped us to improve the content of this paper and also regarding analysis of the PDS data. We thank A. R. Rao for some valuable discussions. The Beppo-SAX satellite is a joint Italian and Dutch program. We thank the staff members of Beppo-SAX Science Data Center and RXTE/ASM group for making the data public.


text:
# Universality in the merging dynamics of parametric active contours: a study in MRI-based lung segmentation

## Abstract

Measurement of lung ventilation is one of the most reliable techniques of diagnosing pulmonary diseases. The time consuming and bias prone traditional methods using hyperpolarized H3He and 1H magnetic resonance imageries have recently been improved by an automated technique based on multiple active contour evolution. Mapping results from an equivalent thermodynamic model, here we analyse the fundamental dynamics orchestrating the active contour (AC) method. We show that the numerical method is inherently connected to the universal scaling behavior of a classical nucleation-like dynamics. The favorable comparison of the exponent values with the theoretical model render further credentials to our claim.

PACS: 87.10. +e, 87.68. +z, 82.60. Nh

Ventilation analysis is an authentic way of diagnosing lung airway diseases. The ratio of the volume of ventilated (functional) portions of lungs to the total lung volume is known as lung ventilation, which is used in validating pulmonary drugs Ray . The process involves two complementary magnetic resonance (MR) imaging modalities, the hyperpolarized helium-3 ($`\mathrm{H}^3`$He) imagery and the proton (1H) imagery. Lung functionality, including the volume of ventilated lungs, can be obtained from the former modality while lung anatomic details, including total lung volume, are accessed through the latter Ray . Since manual investigation of the MR imagery to compute lung ventilation is extremely time consuming, an active contour (AC) or snake based automated method has been proposed Ray ; Xu  to compute the total lung volume from proton MR imagery on a 2-D slice-by-slice basis Ray . A snake is defined as a massless 2-D thin string (closed or open) that can move on the image domain driven by two types of forces, internal elastic forces and external image forces Kass . Under the influence of these two forces an initial contour (snake) clings to image edges and delineates an object. A snake always gives continuous edges unlike any traditional edge detector (such as the Canny method Canny ), thereby eliminating any post-processing steps to connect the detected broken edges. These two properties are particularly useful when the object outline is broken and noisy as in most of the 1H MR imagery. The snake method of finding the object boundary relies on the initial snake placement inside the image. If a small initial snake is placed inside a lung cavity on the MR image, while growing, the snake may be stopped by the associated numerical artifacts and may not capture the actual lung outline Ray . Starting with a larger snake may result in missing the lung cavity completely. A possible solution is to start with multiple non-overlapping small snakes inside the lung cavity and evolve (grow) them until they merge with each other and capture the cavity outline Ray . During such a process, the growing snakes merge with each other into a single contour. This automatic merging of non-overlapping snakes is characterized by certain attributes: a) during the evolution process no two snakes overlap with each other, b) every snake stops evolving at the object edge as a single snake does during its course of evolution, and c) growing convex shaped snakes (e.g., circular or rectangular snakes) inside a convex object recovers the object boundary. Although merging of multiple snakes is experimentally verified in a multitude of cases, a concrete theoretical understanding of this merging snake approach remains a challenge.

The aim of this letter is to bridge this gap by showing that the underlying principle of multiple snake-merging is governed by a universal power law behavior which originates from an inherent nucleating structure. From a biomedical engineering perspective, the study of the efficacy of the lung cavity delineation method is crucial for robust clinical application. The lung cavity segmentation by merging snake method involves three steps: a) initially small non-overlapping contours are placed inside the lung cavity, b) generalized gradient vector flow (GGVF) fields Xu  are computed with a Dirichlet boundary condition on the initial circular snakes Ray , and c) all the snakes are evolved simultaneously and independently of each other with the GGVF force field as the external force for the snakes. This automated lung cavity segmentation is attractive for a number of reasons. While other merging snake algorithms, such as the one proposed by McInerney and Terzopoulos McInerney  is computationally non-trivial compared to the original snake evolution algorithm of Kass et al. Kass , the merging snake algorithm by Ray et al. maintains the same computational simplicity of Kass et al. algorithm. Also, based on the position of an initial snake, the rigidity parameters of the snakes can be varied, so that on one hand delicate high curvature features, such as costophrenic angles, can be accurately captured, while on the other hand, snakes can be made sufficiently stiff in order to avoid capturing artifacts. Fig. 1 shows multiple snake initialization, evolution, merging and delineation of lung boundary by the Ray et al l method.

In order to map the merging snake scenario to statistical thermodynamical systems, we first provide the mathematical background for a parametric active contour or snake. A snake is a curve $`𝐂(s)=(p(s),q(s))`$ defined by the parameter $`s[0,1]`$. The snake is evolved in such a way that it minimizes the energy functional Xu

$$E_s=_0^1(\frac{1}{2}\{\alpha |𝐂^{}(s)|^2+\beta |𝐂^{\prime \prime }(s)|^2\}+E_{\mathrm{ext}}[𝐂(s)])𝑑s$$
(1)
where the first two terms give the internal energy of the snake ($`\alpha ,\beta 0`$) and $`E_{\mathrm{ext}}`$ represents the external energy added to the system. $`𝐂^{}`$ and $`𝐂^{\prime \prime }`$ are the first and second derivatives of the snake with respect to $`s`$. An example of external force for the snake is the gradient force: $`E_{\mathrm{ext}}=|I(x,y)|^2`$, where $`I(x,y)`$ is the image pixel intensity. In general, in the presence of an external force $`𝐕(𝐂(s,t))`$, the time dynamics of such a snake is given by Ray

$$\frac{𝐂(s,t)}{t}=\alpha 𝐂^{\prime \prime }(s,t)\beta 𝐂^{\prime \prime \prime \prime }(s,t)+𝐕(𝐂(s,t)).$$
(2)
Note that the snake in eqn. (2) is now a function of time $`t`$ as well. The stationary solution of eqn. (2) corresponds to a snake that minimizes the energy functional. Ray et al proposed an external force $`V(x,y)`$ obtained by solving the following partial differential equation (PDE) applying Dirichlet boundary condition:

$$g(|f|)^2𝐕(1g(|f|))(𝐕f)=\mathrm{𝟎},$$
(3)
where $`g(\alpha )=\mathrm{exp}(K\alpha )`$ and $`f(x,y)=E_{\mathrm{ext}}(x,y)`$, $`K`$ being a tunable parameter controlling the smoothness of the external snake force field. Here we study the resultant dynamics due to the evolution of a number of such snakes, defined by the above system of forces.

When multiple snakes are evolved inside the desired closed boundary, the growth algorithm confirms that they all finally merge into a single snake after a finite time. This merging of snakes is basically a topological effect and naively the phenomenology reminds one of ”nucleation” as seen in classical first-order thermodynamic systems nucleation . In our case, we allow a finite number of GGVF snakes, each with a finite starting radius, to evolve in a 2-D plane and then numerically evaluate a few measurables – the nucleation time (NT), the bounded area (BA) after nucleation has occurred, the critical radius (CR) at the time of nucleation and also the nucleation rate (NR), all as functions of snake evolution time. The respective quantities are defined as follows – nucleation time is the time required for all the snakes to merge together as a single unit, bounded area is the sum of the areas of the initial snakes before complete merging and the area under the single snake after merging, critical radius is the radius of curvature of the merged structure and nucleation rate is the ratio of the number of snakes to the bounding area, before the merging has actually taken place. One should note that by complete merging we refer to the critical phase when all the initial snakes merge together for the first time.

We perform two numerical experiments with merging. In the first experiment, we start with a circular binary image of radius $`R`$ containing $`N`$ number of circular snakes, each of radius $`r<R`$, randomly distributed inside. As described before, the initial snakes are driven by GGVF forces and they maintain a non-overlapping dynamics. In another numerical experiment, we vary the radii of the smaller circles and later also vary their total number. Both numerics are repeated with varying sizes of the initial domain $`R`$. The enclosed figures 3 and 4 show the variation of the bounding area and the critical radius against time respectively, in non-dimensionalized units, in loglog plots.

We find that all the graphs show excellent scaling in a broad domain, with the nonlinear zones referring to the saturation limits in each case. In the figures shown, the bounding image radius is 80 units while the corresponding starting snake radius is 5 units. We have simulated using different number of initial snakes, indicated in the legend of the figures. The plots show power law variations of each of the measured quantities with time, i. e. $`yx^\alpha `$ and the two relevant exponents, that of critical radius and bounding area, have the values $`\alpha _{\mathrm{BA}}0.91.0,\alpha _{\mathrm{CR}}=0.25`$. We have used multiple combinations of the parameter values but the exponents remain universal, that is they invariably follow a power law statistics.

To analyse the results further, we begin with the hypothesis that the phenomenology of the whole process, that of various snakes merging together under certain boundary conditions in the presence of suitably defined external forces, is fundamentally similar to that of a domain growth process as encountered in the realms of classical nucleation. In what follows, we try to evaluate the time rate of evolution of a system of snakes before and after all the snakes have completely merged with each other. We start from a model of solutes trying to nucleate in a solution and then compare the results with the merging mechanism of snakes we find in the AC method.

Our starting description is that of the diffusive growth of two dimensional spherical droplets (circles, in 2-D) in a Lifshitz-Slyozov nucleation ; LStheory ; binderstauffer  type of continuum theory. We start with the thermodynamic definition of the work $`\mathrm{\Delta }W`$ required to form a nucleus from an aggregate of solute particles

$$\mathrm{\Delta }W=\delta (ET_0S+P_0V)$$
(4)
where $`\delta (ET_0S)`$ is the increment in the free energy of the system and $`\delta (P_0V)`$ is the associated external work done. Optimizing this total work factor, one can show nucleation  that the critical radius of a nucleating system, this being the rate of merger of snakes in AC method, is given by

$$R_{\mathrm{CR}}=\frac{\beta s^{}}{\mu ^{}(P)\mu _0^{}(P)}$$
(5)
Here $`\mu _0^{}`$ and $`s^{}`$ are the chemical potential and molecular surface area of the nucleus, $`\beta `$ is the surface tension and $`\mu ^{}`$ is the chemical potential of the solute in the solution. As more and more solutes start nucleating i. e. snakes merge, the solution approaches the critical limit of supersaturation and in the steady state, the critical radius grows as

$$\frac{dR(t)}{dt}=D(\frac{c}{r})|_{r=R}$$
(6)
where $`c(R)`$ represents the spherically symmetric concentration distribution of snakes around a nucleus of radius $`a`$, $`D`$ being the diffusion coefficient. Following the treatment of Lifshitz-Slyozov nucleation ; LStheory , one can show that if $`R_{\mathrm{CR}}(0)`$ is the critical radius at the beginning of the merging process, then for a diffusive nature of relaxation, the radius follows a dynamics

$$\frac{dR(t)}{dt}=\frac{R_{\mathrm{CR}}^2(0)}{R(t)}(\frac{1}{R_{\mathrm{CR}}(t)}\frac{1}{R(t)})$$
(7)
We now define the dimensionless quantities $`x(t)=\frac{R_{\mathrm{CR}}(t)}{R_{\mathrm{CR}}(0)}`$, $`u(t)=\frac{R(t)}{R_{\mathrm{CR}}(t)}`$ and $`\tau =2\mathrm{log}x(t)`$. The last one increases monotonically from 0 to $`\mathrm{}`$ as the time $`t`$ increases likewise. Combining (6) and (7) in 2D, the dynamics is given by

$$\frac{du^2(\tau )}{d\tau }=\gamma (1\frac{1}{u})u^2$$
(8)
where $`\gamma =\frac{dt}{R_{\mathrm{CR}}(0)dx}>0`$. A linear stability analysis of the above nonlinear equation gives a fixed point at $`\gamma _0=\frac{27}{4}`$. We are interested in the dynamics around an $`ϵ`$-neighborhood of this point $`\gamma _0`$. If $`\gamma (\tau )=\frac{27}{4}[1ϵ^2(\tau )]`$ ($`ϵ0`$ as $`\gamma \mathrm{}`$), then near the critical point ($`u_0=(\gamma /2)^{1/3}=3/2`$) the merging snakes follow a dynamics defined by

$$\frac{du}{d\tau }=2(u\frac{3}{2})^2\frac{3}{4}ϵ^2$$
(9)
The time variation of the merging nuclei (snakes) is

$$x(t)=\frac{4}{27}\frac{\sqrt{t}}{R_{\mathrm{CR}}(0)}$$
(10)
Thus the bounding area ($`x^2(t)`$) of the merging snakes grow at the rate of $`t`$ which is roughly speaking our numerical estimate (0.9-1.0) also (Fig. 3). However, the above analysis is only valid before complete merging of the snakes has occurred. In the merged phase when the resultant asymmetrical structure continues growing finally to coalesce with the bounding image, the system dynamics is modified. To analyse the situation, we start from eqn.(7). The equation may alternatively be represented as

$$\frac{dR(t)}{dt}=\frac{D\delta c(R(t))}{R(t)}\frac{1}{T}\frac{1}{R^2(t)}.$$
(11)
$`T`$ in the above equation is the temperature of the solution in the equivalent thermodynamic problem which in our case, is a measure of the average surface energy $`E`$, $`ET`$, of the system. Evidently, in 2D, $`E`$ goes as $`2\pi R\beta `$. Plugging these values in the above equation, we see that after complete merging of the snakes has taken place, the effective radius of curvature of the resultant structure grows as $`R_{\mathrm{CR}}(t)t^{1/4}`$ (Fig. 4). Once again our numerical result is in exact harmony with the theory.

Our above analysis, both numerical and analytical clearly suggests that below the apparently simplistic level of the GGVF application, the system dynamics has a more fundamental symmetry. This symmetry comes from the fact that the GGVF method lies in the same universality class as that of a classical nucleation model. This first-order critical response of the system is what provides the subtlety of the underlying physics. There is, however, a notable shortcoming with the GGVF technique, that it does not allow us the liberty of starting with an initial condition at an arbitrary location. If the snake starts at a position in which a major portion of the initial snake is outside the desired boundary or vice versa, then the snake driven by GGVF will not converge to the actual boundary. Although it can be shown by the Reed-Simon’s theorem reedsimon  that a convex set (a circle, say) growing within a larger convex set (a larger circle or rectangle, say) will always merge with the outer boundary under the action of isotropic driving forces, to the best of our knowledge, no such mathematical lemma exists for a convex set growing in a concave set or vice versa.

In summary, our achievement in the present paper has been to analyze the physical foundation of the GGVF merging technique and to provide answer to the rather puzzling question as to why it works so accurately. In the process, we have shown that the answer lies in the general scaling behavior of the underlying nucleation dynamics, defined by proper scaling laws. This clearly indicates that an active contour system is in the same universality class as the nucleation model we considered. Our results, in unison with probable biomedical applications, are expected to inspire further studies in the understanding of lung-based diseases.


text:
# What’s in a Name?

## I Introduction

Despite its long tradition in mathematics and computer science, graph theory Bollobás (2002) has reached great popularity only recently through innovative research in the novel area which became known as complex networks Albert and Barabási (2002). By integrating theoretical principles, especially from statistical mechanics, with experimental and simulated data, recent investigations have shown that several important natural phenomena such as infectious diseases, ecological systems, protein folding, society and the internet, are characterized by scale-free and/or small-world behavior as far as their connectivity is concerned Albert and Barabási (2002). In particular, studies modeling linguistics in terms of complex networks have indicated that several aspects of human language, such as word proximity Cancho and Solé (2001) and synonyms Albert and Barabási (2002); Motter et al. (2002), are at least partially characterized by power law behavior. As language corresponds to one of the most essential manifestations of the human brain, such findings can be taken as an indication that that complex structure, or at least its portions more closely related to language and associations, may also be intrinsically organized according to power laws and scale free behavior Motter et al. (2002). As the conscious and predominantly sequential flow of ideas in humans, James’s *fringe of consciousness* James (1955), is closely related to the externalization of ideas through language, it is also possible that the scale free properties found in linguistic structures can also be an intrinsic property of consciousness.

The current work aims at investigating such possibilities through a psychophysical experiment involving human subjects to associate words from the English language. By understanding the presented words and associations as graph nodes and edges, respectively, it is possible to perform a quantitative analysis of the digraph connections by considering statistics (average and standard deviation) of network measurements such as the node degree, the average length, and the clustering coefficient.

The current article starts by describing the experimental approach and proceeds by analysing and discussing the respectively obtained data.

## II The Psychophysical Experiment

Along the last decades, psychophysics has establishing itself as an important area in psychology and neuroscience, providing invaluable means for quantifying perception. Provided the experiments are carefully devised and conducted, objective and relatively precise information can be obtained about the dynamics of perception. As in physics, the experiments have to be planned and performed while most factors likely to influence the investigated phenomena are kept constant. The popularization of personal computers has motivated the use of such machines for automating of psychophysical experiments, accounting for enhance repeatability and standardization.

In this work, a program was developed in SCILAB with the specific finality of investigating associations between words from the English language. Starting with the word ’sun’, the subject is prompted to associate a subsequent word. No specific instructions are given regarding the type of association, except that special characters, plurals and verb conjugations are to be completely avoided. There is no time limit for providing the new word, and the experiment can be broken into several sections, while collecting all the obtained data stored into files. An illustration of the first steps of the experiment is provided below, where the words supplied by the subject are represented in italic:

* sun $``$ *desert*
* desert $``$ *pyramid*
* sun $``$ *gold*
* pyramid $``$ *triangle*
* pyramid $``$ *desert*
* triangle$``$ *square*
* $`\mathrm{}`$

Observe that the only predefined word is that presented first, all the others being subsequently defined by the subject. The presented and suggested words are henceforth referred to as *presented word* and *input word*, respectively. After each new word is input, its presence in the current list of words is verified, the word being included otherwise. Each word is treated as a graph node, and each pair of words is understood as a graph edge $`(`$presented word, input word$`)`$. The therefore obtained direct graph (i.e. a digraph), with the frequency of each association treated as the weight of the respective edge, provides an interesting formal representation of the word associations. The whole sequence of presented and input words is recorded for further analysis. In order to guarantee the words to be presented in a uniform fashion, in the sense that each word is presented about as many times as the others, a density probability function $`p(w)`$ describing the number of times each word is presented is kept all times. The presented words are drawn from the complemented density function, i.e. $`max\{p(w)\}p(w)`$, so that the less frequently presented words have higher likelihood to be chosen, leading to a levelling effect. The experiment terminates after a pre-defined number of words are presented, and the more recent input words, which have consequently been presented only a few times, are excluded from the data and respective network.

## III Results

The above experiment was performed with a single subject along a whole week, totaling 305 different words from which 250 words were chosen (the remainder, more recently input words, were discarded for the sake of enhanced uniformity). A total of 1930 associations were recorded. The types of the input words is given in Table 1, and Table 2 shows the frequency of types of associations. Figure 1 presents the population of each input word, and Figure 2 depicts the occurrence of new words along the presentation stages identified by $`i`$. Figure 3 gives the histogram of repeated associations. The average and standard deviation of the node degree $`k`$, clustering coefficient $`C`$ and average length $`\mathrm{}`$ are presented in Table 3. Figure 4 shows the histogram of equal words apart by specific lags along the presentation sequence. For instance, the ordinate value at lag 100 indicates the number of equal words distant of lag along the sequence. For the sake of enhanced uniformity, the sequence is considered up to its total length minus the maximum lag value. The loglog curves of the cumulative output and input node degree (recall that we are dealing with a digraph) are presented in Figures 5 and 6.

## IV Discussion

The several interesting trends and phenomena identified by analysis of the experimental data are characterized and briefly discussed in the following:

*A. Attractor formation:* As shown in Figure 2, the number of new words input by the user tended to diminish, reaching a near equilibrium state where very few new words are likely to be added. This suggests an attraction basin defined by the initial word.

*B. Word density asymetry:* As illustrated in Figure 1, the subject tended to enter some words more often than others. Particularly, there was a generalized preference for adjectives such as *good* and *long*, among others. This is hardly surprising, as adjectives are more immediately applicable to several words.

*C. Edges asymmetry:* As clearly seem from Table 2, not every association is reciprocal, i.e. the existence of an edge $`(i,j)`$ does not necessarily implies the presence of $`(j,i)`$. Examples of such asymmetric cases obtained in the considered experiment include $`(sky,blue)`$ and $`(blue,sky)`$. To some extent, such asymmetries are observed in cases involving a more common word followed by a less common one, such as a general adjective and a specific noun. Another characteristic that has been verified from the experiment is the tendency of the associations to correspond to synonyms and antonyms, especially regarding pairs of adjectives.

*D. Wide variation of node degree:* The high standard deviation obtained for the node degree indicates that the number of associations induced by each presented word varies considerably. It is possible that more common words which usually appear connected to several other words, such as adjectives, tend to favor higher number of associations.

*E. Associations asymmetry:* One of the clearest results deriving from the reported experiment was the fact that some associations tended to be much more stable than others, in the sense that they were more systematically repeated and yielded a smaller number of variations. Examples of such cases include $`(bread,butter)`$ and $`(pecker,wood)`$. This property seems to be connected to the node degree wide dispersion, in the sense that association pairs involving at least one word characterized by higher node degree tended to favor a higher number of different associations.

*F. Context biasing:* As is clear from Figure 4, the choice of a word by the subject tended to be influenced by those more recently input. The memory effect seems to disappear for lags higher than 250 presentations.

*G. Small-world features* The relatively low average length shown in Table 3 suggests that the obtained association follows the small-world paradigm. This is an immediate consequence of the fact that the experiment inherently targets word associations.

*H. Power-law features:* As could be expected, the output and input degree distributions shown in Figures 5 and 6 resulted markedly different, with the latter being more compatible with power-law scaling, especially at the initial portion of the curve.

As the limited power-law trend corresponds to the possibly more complex and interesting identified features, it is further discussed in the following, including a possible explanation. Although difficult to be defined, the conscious portion of thinking is a predominantly sequential process. While solving a problem, or just relaxing, the flow of ideas and concepts takes place as a sequence of ideas associated in some way which is highly dependent on the context defined by the more recent thoughts. To a large extent, the successive ideas are characterized by some strong or weak association. For instance, after thinking about the sky, next possible ideas are likely to be blue, air, sun or clouds. Therefore, at least part of the flow of thoughts can be thought in terms of a Markovian system. At the same time, memories are often related to associations Humphreys et al. (1989). From the computational point of view, it is possible to enhance the storage potential by organizing the stored concepts in a hierarchical fashion, so that the description of new concepts at lower hierarchical levels can include only the features not covered by the upper levels. For instance, the description of a cat can be derived from that of mammals, including only those characteristics that are intrinsic to cats (see Figure 7). This concept of inheritance leads naturally to associations between concepts and ideas, even between two non-adjacent hierarchical levels, a phenomenon that can be related to Hofstadter ’strange loops’ Hofstadter (1999). Though additional features are certainly incorporated into the brain dynamics, such hierarchical and associative schemes lead to the interesting situation where several concepts end up associated, even if indirectly, to those in the upper hierarchies. Consequently, the concepts tend to become more and more associated as one moves from the lower to the upper hierarchical levels, possibly leading to a rich gets richer scheme, and hence to scale free organization.

Given that the adjacency matrix of the obtained graph can be immediately understood as the transition matrix of a Markovian systems, it is possible to use Monte Carlo simulation in order to produce sequences of associated words, such as that illustrated below. As the context is limited to one association level, such sequences are characterized by subsequent repetitions of words.

> horse, brown, bear, brown, sugar, sweet,
> good, earth, land, good, well, good, time,
> out, sun, hot, water, cold, water, cold,
> wool, sheep, four, clock, six, tea, leaf,
> thin, sheet, wide, field

## V Concluding Remarks

The present work has illustrated how complex network and graph theory concepts and formalisms can be applied to characterize human cognitive activities, namely the association of words. While previous related works such as Motter et al.Motter et al. (2002) investigated word associations through the use of static databases, the current approach considered psychophysical experiments. The main differences implied by such an approach are the fact that the importance of associations can be inferred from the respective frequencies. In addition, a random element is implied by the fact that the user is likely to vary the chosen associations while affected by the context established by the presentation sequence.

Although limited to a single subject, the obtained experimental results led to a series of interesting findings, including the identification of attraction basins, context biasing, association asymmetries, small-world features, and near power-law scaling of the node degree. A putative model possibly underlying the latter phenomenon, involving the appearance of strange loops in the hierarchical categorization of words, has also been proposed. While extensive additional investigations are required in order to confirm such preliminary results, it is felt that the identified phenomena are likely to provide a reasonably formal scaffolding for further investigating and understanding word associations by humans and even more sophisticated brain dynamics Motter et al. (2002).

Several are the possibilities implied by the reported developments. First, it is important to note that the specific measurements extracted from digraphs obtained from different subjects can be possibly correlated to individual features or even for diagnosis. At the same time, it is likely that the obtained graphs will present a core shared by several subjects, corresponding to those more established and invariant collective concepts, while the graph difference residuals could provide interesting information about intrinsic individual features and preferences. Another interesting task would be to extend the reported approach in order to investigate associations in visual language, for instance by using eye-tracking systems. Several possibilities for further investigation can be defined by considering modified versions of the adopted psychophysical experiment. For instance, it would be interesting to study situations where the subject is allowed to enter a continuous flow of associated words, without any interference from the computer, except the presentation of the first word. Although more complex, given the additional degrees of freedom, such investigations could provide additional insights about long time memory effects, which are poised to reduce the number of word repetitions in respective Monte Carlo simulations. It would be interesting to compare how such extended context modifies the properties of the respectively obtained networks.

###### Acknowledgements.

The author is grateful to FAPESP (processes 99/12765-2 and 96/05497-3) and CNPq for financial support.