text:
# 1 Introduction

## 1 Introduction

Cohomology is defined by *cochain complex*

$$0C^0\stackrel{d^0}{}\mathrm{}\stackrel{d^{k2}}{}C^{k1}\stackrel{d^{k1}}{}C^k\stackrel{d^k}{}C^{k+1}\stackrel{d^{k+1}}{}\mathrm{}.$$
(1)
Here $`C^k`$ are linear spaces (more generally, abelian groups), graded by the integer number $`k,`$ called *dimension* or *degree* (depending on the context). The elements of the spaces $`C^k`$ are called *cochains*.

The linear mappings $`d^k`$ are called *differentials* (or *coboundary operators*). The main property of these mappings is “their squares are equal to zero”: $`d^kd^{k1}=0`$.

The elements of the space $`Z^k=\mathrm{Ker}d^k`$ are called *cocycles*. The elements of the space $`B^k=\mathrm{Im}d^{k1}`$ are called *coboundaries*. Note that $`B^kZ^k.`$

The $`k`$th *cohomology* is the quotient space

$$H^k=Z^k/B^k\mathrm{Ker}d^k/\mathrm{Im}d^{k1}.$$

There are many cohomological theories designed for investigation of different mathematical structures and the space $`H^k`$ carries important information about peculiarities in these structures. The only difference between cohomological theories lies in the constructions of the cochain spaces and coboundary operator. These constructions depend on the underlying mathematical structures.

The cohomology of the Lie (super)algebra $`A`$ in the module $`X`$ is defined via cochain complex (1) in which (see, e.g., ) the cochain spaces $`C^k=C^k(A;X)`$ consist of super skew-symmetric $`k`$-linear mappings $`A\times \mathrm{}\times AX,`$ $`C^0=X`$ by definition. Super skew-symmetry means symmetry with respect to swapping of two adjacent *odd* cochain arguments and antisymmetry for any other combination of parities for adjacent pair.

The differential $`d^k`$ takes the form111This version of formula for differential corresponds to the algorithm used in the program LieCohomology.

$`(d^kc)(a_0,\mathrm{},a_k)`$ $`=`$ $`{\displaystyle \underset{0i<jk}{}}(1)^{s(a_i)+s(a_j)+p(a_i)p(a_j)}c([a_i,a_j],a_o,\mathrm{},\widehat{a_i},\mathrm{},\widehat{a_j},\mathrm{},a_k)`$ (2)
$`{\displaystyle \underset{0ik}{}}(1)^{s(a_i)}a_ic(a_o,\mathrm{},\widehat{a_i},\mathrm{},a_k),`$
where the functions $`c(\mathrm{})`$ are elements of cochain spaces; $`a_iA`$; $`p(a_i)`$ is the parity of $`a_i`$; $`s(a_i)=i,`$ if $`a_i`$ is even element and $`s(a_i)`$ is equal to the number of even elements in the sequence $`a_0,\mathrm{},a_{i1}`$, if $`a_i`$ is odd element. In the case of trivial module (i.e., if $`ax=0`$ for all $`aA`$ and $`xX`$) one uses as a rule the notation $`H^k(A).`$

In papers  we presented an algorithm for computation of Lie (super)algebra cohomologies. These papers contain also description of its C implementation and some results obtained with the help of codes designed. This algorithm computes cohomology of Lie (super)algebra $`A`$ over module $`X`$ in a straightforward way, i.e., for cochain complex (1) the algorithm constructs the full set of basis super skew-symmetric monomials forming the space $`C^k`$, generates subsequently all basis monomials in the space $`C^{k+1}`$, computes the differentials corresponding to these monomials to obtain the set of linear equations determining the space of cocycles

$$Z^k=\mathrm{Ker}d^k=\{C^k|dC^k=0\},$$
(3)
constructs the space of coboundaries

$$B^k=\mathrm{Im}d^{k1}=\{C^k|C^k=dC^{k1}\}.$$
(4)
Finally, the algorithm constructs the basis elements of quotient space

$$H^k(A;X)=Z^k/B^k.$$
(5)
This last step is based on the Gauss elimination procedure.

The main difficulty in computing cohomology results from the very high dimensions of the spaces $`C^k`$: for $`n`$-dimensional ordinary Lie algebra and $`p`$-dimensional module

$$dimC^k=p(\begin{array}{c}n\\ k\end{array}),$$

and for $`(n|m)`$-dimensional Lie superalgebra

$$dimC^k=p\underset{i=0}{\overset{k}{}}(\begin{array}{c}n\\ ki\end{array})(\begin{array}{c}m+i1\\ i\end{array})p(\begin{array}{c}n\\ k\end{array})+p\underset{i=1}{\overset{k}{}}(\begin{array}{c}n\\ ki\end{array})(\begin{array}{c}m+i1\\ i\end{array}).$$

In many cases it is possibly to extract some easier to handle subcomplexes of the full cochain complex (1). The partition of cochain complex for a graded algebra and module into homogeneous components is a typical example. In many papers (see, e. g., ) more special subcomplexes were used successfully to obtain new results in the theory of cohomology of Lie (super)algebras.222The main trick consists in imposing some restrictions on the elements of $`C^k`$ and proving the invariance of these restrictions with respect to the differential.

The main idea of the new algorithm presented in  is to extract the minimal possible subcomplexes from complex (1) and to carry computations within these subcomplexes. There are two versions of the algorithm. One of them is applied when the cochain spaces under consideration are infinite-dimensional (or their dimensions are too large to fit the available memory), but the minimal subcomplexes contain finite-dimensional spaces of $`k`$-cochains. Another version of the algorithm is applied when it is possible to construct the full space $`C^k.`$ Below we present this version in the pseudocode form.

Here the subalgorithm GenerateMonomials generates the full set $`M_g^k`$ of super skew-symmetric monomials

$$c(\alpha _{i_1},\mathrm{},\alpha _{i_k};\xi _\iota )c(\alpha _{i_1})\mathrm{}c(\alpha _{i_k})\xi _\iota \alpha _{i_1}^{}\mathrm{}\alpha _{i_k}^{}\xi _\iota $$

forming basis of the cochain space $`C^k`$ in the grade $`g;`$ $`\alpha _{i_j}A`$ and $`\xi _\iota X`$ are basis elements of algebra and module; $`\alpha _i^{}`$ is the dual to $`\alpha _i`$ element. The subalgorithm ChooseMonomial takes some monomial $`m_g^kM_g^k.`$ This monomial is starting monomial for constructing the subcomplex $`s`$ by the subalgorithm ConstructSubcomplex. The subalgorithm ComputeCohomologyInSubcomplex computes basis cohomological classes $`BH_{g,s}^k`$ in the subcomplex $`s`$ by the straightforward algorithm described above.

## 2 Computation of $`H_g^k(\mathrm{H}(2|0))`$

In this section we present the results of computation of cohomology in the trivial module for Lie algebra $`\mathrm{H}(2|0)`$ of formal hamiltonian vector fields on the $`2`$-dimensional simplectic manifold. We describe also cohomological classes up to grade 8 for the Poisson algebra $`\mathrm{Po}(2|0)`$ which is a central extension of the algebra $`\mathrm{H}(2|0)`$.

The hamiltonian algebra $`\mathrm{H}(2\mathrm{n}|\mathrm{m})`$ is an algebra of vector fields (see, e.g., ) acting on the $`(2n|m)`$ supermanifold and preserving the following 2-form

$$\underset{i=1}{\overset{n}{}}dp_idq_i+\underset{j=1}{\overset{m}{}}d\theta _jd\theta _j,$$

where $`p_1,\mathrm{},p_n;q_1,\mathrm{},q_n`$ and $`\theta _1,\mathrm{}\theta _m`$ are even and odd local variables on the supermanifold, respectively. The elements of $`\mathrm{H}(2\mathrm{n}|\mathrm{m})`$ can be expressed in terms of *generating function* $`f(p_1,\mathrm{},p_n;q_1,\mathrm{},q_n;\theta _1,\mathrm{}\theta _m)`$ by the formula

$$\underset{i=1}{\overset{n}{}}\left(\frac{f}{p_i}\frac{}{q_i}\frac{f}{q_i}\frac{}{p_i}\right)(1)^{p(f)}\underset{j=1}{\overset{m}{}}\frac{f}{\theta _j}\frac{}{\theta _j},$$
(6)
where $`p(f)`$ is parity of the function $`f`$ (this function is called usually *hamiltonian*). Thus one can consider the formal hamiltonian vector fields as linear combinations of monomials in the variables $`p_i,`$ $`q_i`$ and $`\theta _j`$ (except of the monomial 1). Considering these monomials as basis elements of $`\mathrm{H}(2\mathrm{n}|\mathrm{m})`$ and using prescribed $`\mathrm{𝖹𝖹}`$-grading for the variables $`p_i,`$ $`q_i`$ and $`\theta _j`$ one can impose $`\mathrm{𝖹𝖹}`$-grading $`\mathrm{gr}()`$ on the algebra $`\mathrm{H}(2\mathrm{n}|\mathrm{m}).`$ The standard grading is $`\mathrm{gr}(p_i)=\mathrm{gr}(q_i)=\mathrm{gr}(\theta _j)=1.`$ For the standard grading the grade of algebra element corresponding to some monomial is equal to the grade of this monomial minus 2 (due to two differentiations in the terms of expression (6)).

Since the hamiltonian algebra is very important in both classical and quantum physics many efforts were applied to investigation of its cohomology. Most advanced results were obtained for the finite-dimensional algebras of the form $`\mathrm{H}(0|\mathrm{m})`$ . Nevertheless the hamiltonian algebras on supermanifolds with nonzero even dimension are more important in applications but computation of their cohomologies is much more difficult task. Some results about cohomologies of such algebras were obtained in . In the paper  some elements of $`H_g^k(\mathrm{H}(2|0))`$ were calculated by considering special subcomplexes (and using computer partially). We present here all cohomological classes (without discussing their meaning and interpretation) from $`H_g^k(\mathrm{H}(2|0))`$ up to grade 8.

The results of computation are summarized in Table 2. The boxes of this table corresponding to cochain degree $`k`$ and cochain grade $`g`$ contain the following information: $`dimC_g^k,`$ dimension of the full space of $`k`$-cochains in grade $`g`$; number of minimal subcomplexes $`C_{g,s}^{k1}\stackrel{d_{g,s}^{k1}}{}C_{g,s}^k\stackrel{d_{g,s}^k}{}C_{g,s}^{k+1}`$ extracted by the algorithm from the full complex; $`\mathrm{max}dimC_{g,s}^k,`$ maximum dimension of the subspace of $`(k,g)`$-cochains among all subcomplexes. The empty box means that $`dimC_g^k=0,`$ i.e. the space of $`(k,g)`$-cochains is empty. The boxes marked by the bullet $``$ contain nontrivial 1-dimensional cohomological classes. For example, the box corresponding to the pair $`(k,g)=(7,8)`$ tells that $`dimC_8^7=25488,`$ number of subcomplexes is 21, $`\mathrm{max}dimC_{8,s}^7=3148`$ and $`dimH_8^7=1.`$

More detailed information about computation in $`(k,g)=(7,8)`$ is given in Table 1.

In this table the columns $`dimZ_{g,s}^k,`$ $`dimB_{g,s}^k`$ and $`dimH_{g,s}^k`$ contain dimensions of cocycle, coboundary and cohomology spaces in subcomplexes, respectively. On can see that there are 10 pairs of subcomplexes with repeated structure and the only single subcomplex containing nontrivial cohomological class.

The full set of nontrivial (1-dimensional) cohomological classes is: $`H_2^2,H_2^5,H_0^7`$ (computed earlier) and $`H_8^7,H_6^{10}`$ (computed by the new program). As to the Poisson algebra $`\mathrm{Po}(2|0),`$ the part of its cohomological classes up to grade 8 coincides with those for $`\mathrm{H}(2|0)`$ except of $`H_2^2`$.333The cocycle $`H_2^2`$ describes the central extension of the algebra $`\mathrm{H}(2|0)`$ to $`\mathrm{Po}(2|0).`$ $`H_{g8}^k(\mathrm{Po}(2|0))`$ contains also four additional classes: $`H_4^6,H_2^8,H_6^8,H_6^{11}.`$ But all these classes are multiplicative consequences of the classes $`H_2^5,H_0^7,H_8^7,H_6^{10}.`$ These classes can be expressed in the form $`H_g^k=H_{g+2}^{k1}c(𝒵)`$ due to the general property  of cohomology of algebras containing central element $`𝒵.`$

## 3 Conclusion

Our new algorithm demonstrates substantially higher efficiency in comparison with the old one. For example, the program described in  computes the case $`(k,g)=(6,5)`$ in 35 min 45 sec = 2145 sec whereas the new program takes 54 sec for this task. For both runs we used PC Pentium III, 667MHz, 256MB RAM. The superiority of the new program grows with increasing of the task complexity. Nevertheless, due to rapidly increasing computational complexity the presented results are not sufficient to derive any general idea about the structure of cohomology ring $`H_{}^{}(\mathrm{H}(2|0))`$. Our computation was carried over the field of rational numbers $`\text{ }\mathrm{Q}.`$ As profiling shows, the most time consuming part of computation by the program LieCohomology is multiprecision arithmetic. This is common difficulty for almost all problems in computer algebra. It seems that carrying computation over the finite fields, say $`\mathrm{𝖹𝖹}_p,`$ we can go to the grade 40-50 for the problem considered here, but the results obtained in this way can be considered merely as hints.

## Acknowledgements

This work was partially supported by the grants RFBR 01-01-00708, RFBF 00-15-96691 and INTAS 99-1222.


format: text

text:
# Self-assembly of Nanometer-scale Magnetic Dots with Narrow Size Distributions on an Insulating Substrate

## Abstract

The self-assembly of iron dots on the insulating surface of NaCl(001) is investigated experimentally and theoretically. Under proper growth conditions, nanometer-scale magnetic iron dots with remarkably narrow size distributions can be achieved in the absence of a wetting layer. Furthermore, both the vertical and lateral sizes of the dots can be tuned with the iron dosage without introducing apparent size broadening, even though the clustering is clearly in the strong coarsening regime. These observations are interpreted using a phenomenological mean-field theory, in which a coverage-dependent optimal dot size is selected by strain-mediated dot-dot interactions.

Clustering on surfaces by nucleation and growth during atom deposition has been an important subject in basic and applied science for decades Venables ; Zinke . Recent efforts have been focused on searching for methods to obtain nanometer-scale clusters with narrow size distributions. Such clusters or quantum dots are potentially valuable for optical, electronic, and magnetic device applications, but mass production of such structures by lithography or etching-based fabrication has proved to be exceptionally challenging Wolf ; review . Alternatively, it has been realized that the strain energy associated with the lattice mismatch between the dot and the substrate materials can be exploited to induce self-assembled formation of quantum dots with narrow size distributions. This has generated much excitement, particularly in the area of semiconductor quantum dots Wolf ; review ; Mo ; Jesson ; Xie ; Xin . In such cases, the growth of the dots often proceeds in the Stranski-Krastanow (SK) mode, which is characterized by the presence of a wetting layer prior to three-dimensional (3D) clustering. To date, the precise mechanism for size selection in semiconductor quantum dot systems remains a subject of active debate review ; Mo ; Jesson ; Daruka : some attribute them to strain-induced thermodynamic equilibrium states, while others associate them with metastable configurations due to kinetic limitations.

Although improved size uniformity can be achieved in quantum dot formation via the SK growth mode, the presence of a wetting layer is often undesirable, particularly for electronic and magnetic device applications of metallic/magnetic quantum dots. For this reason, it is preferred to fabricate quantum dots in the Volmer-Weber (VW) growth mode, which is characterized by immediate 3D clustering on the substrate surface. Indeed, considerable recent efforts have been devoted to metallic/magnetic quantum dot formation on various substrates Fruchart ; huang ; Baumer , but no significant size uniformity has been achieved in such studies.

In this Letter we investigate the self-assembly of iron dots on NaCl(001), an insulating substrate, by thermal deposition and variable-temperature atomic force microscopy in ultrahigh vacuum. We show that, by properly choosing the growth conditions, nanometer-scale magnetic iron dots with remarkably narrow size distributions can be achieved in the absence of a wetting layer (VW growth). Moreover, by changing the dosage of iron, we can tune both the vertical and lateral sizes of the dots without introducing apparent size broadening, even though the clustering is already in the strong coarsening regime, signified by the decrease in dot density as a function of the iron dosage. The preserved narrowness in the island size distributions is in clear contradiction with the expectations of existing understanding of clustering on surfaces in the coarsening regime Zinke . We interpret these observations within a phenomenological mean-field theory, in which a coverage-dependent optimal dot size is selected by the competition between the self-energy of a dot, i.e. the total energy of an isolated Fe dot on NaCl(100), and the energy of the strain-induced dipolar interaction between the dots.

The experiments were performed in an ultrahigh vacuum (UHV) system with base pressure of 1$`\times `$1010 torr. The system is equipped with electron beam sources, laser molecular beam epitaxy, and an in-situ Omicron variable-temperature UHV beam deflection atomic force microscope (AFM)/ scanning tunneling microscope (STM) with cooling and heating facilities covering a temperature range of (13-1500) K. The noncontact mode AFM was used to study the surface morphology in this work. The NaCl single crystal substrates were cleaved in air, then were immediately loaded into the UHV chamber and were annealed to 530K for one hour to remove surface contamination prior to the experiments. AFM images of the cleaved NaCl(001) surface show monatomic height steps and large terraces with no detectable adsorbates. The iron was evaporated from an Fe wire (5N purity) heated by electron beam bombardment at a rate of 0.04 ML/min (1ML is equivalent to the nominal surface atomic density of bcc Fe(110), 1.7$`\times `$1015 atoms/cm2). The iron dots can be formed only within a finite temperature window, since it is known that the morphology of the substrate changes dramatically if the temperature is above 720K Szymonski , whereas at low temperatures, only random clusters or percolated iron films form Gai . The dots discussed in this Letter are grown at a substrate temperature of 530K.

Nanometer-scale dots can be formed directly on the NaCl(001) surface without a wetting layer, indicating that the growth proceeds in the Volmer-Weber mode. As shown in the AFM images in Figs. 1 and 2, deposited Fe atoms tend to nucleate and form dots even at submonolayer coverage. The preferred VW growth mode is not too surprising, given the large surface free energy difference (2.48 J/m2 for Fe versus 0.18 J/m2for NaCl)Boer  and the large lattice mismatch between Fe and NaClltt . The most eye-catching characteristic of the dots is their size uniformity. Fig. 1 (a) shows a representative AFM image of the dots with a nominal Fe thickness of 1.7ML. The very uniform dots are randomly distributed on the terraces. Fig. 1 (b) is a close up together with a line profile of the dots. Estimated from the line scan, the dots are 3.5 to 4 nm in height, and around 14 nm in diameter ldia . Fig. 1 (c) and (d) are the lateral size and height distributions of the dots shown in Fig.1 (a) along with the corresponding Gaussian fits. Consistent with the line profile, the centers of the Gaussian fits of the height and diameter distributions of the dots appear at 3.76 nm and 13.5 nm, respectively. Each distribution has a very narrow width. The dispersions, $`\mathrm{\Delta }h=(h^2h^2)^{1/2}`$, of the height and diameter distributions are only 0.36 nm and 1.1 nm, respectively, which are less than 10% of the average height and diameter values. Such narrow dispersions are comparable with the narrowest dispersions achieved in semiconductor quantum dot growth via the SK mode Xin , and are particularly striking because clusters grown via the VW mode normally have much broader distributions Wolf ; Fruchart ; huang ; Baumer .

Another striking observation lies in the coverage dependence of the dots shown in Fig. 2. As the Fe coverage increases, the dots grow larger and taller, but an optimal size of the dots is always selected while the narrowness of their lateral and vertical size distributions is preserved. Our quantitative analysis of the dispersions does not yield any significant additional (coverage-dependent) broadening of the dot size distributions, consistent with the impression drawn by naked eye from the AFM images shown in Fig. 2. In Fig. 3, the optimal height and lateral size of the dots are plotted as functions of the iron coverage, together with the density of the dots. The decrease of the dot density with increasing coverage indicates that the growth of the dots is clearly in the coarsening regime. This observation on the variation of the island density makes the persistence of the narrow size distributions even more dramatic, because in the coarsening regime, the absolute size distribution of the islands is expected to broaden with increasing coverage Zinke .

Preliminary surface magneto-optical Kerr effect (SMOKE) measurements show that these iron dots are superparamagnetic in nature. The details of the magnetic investigation, however, are beyond the scope of the present paper, which focuses on the formation of the highly uniform dots. In order to understand the underlying formation mechanism of these nanoclusters, it is highly desirable to know both their detailed atomic structures and their stability. Unfortunately, the iron dots formed on NaCl cannot be imaged by STM with higher (atomic) resolution. On the other hand, in-situ time (from immediately after formation up to 10 days) and temperature (from room temperature up to 550 K) dependent AFM studies reveal little change in the morphology of the Fe dots, indicating that the dot arrays on the surface are in an energetically favorable configuration that is at least metastable.

At the present, the precise underlying mechanism for the formation of these narrow-sized iron dots is not completely clear, but the following phenomenological mean-field theory, stressing on strain-mediated dot-dot interactions, is able to reproduce most of the salient features observed in the experiment. First, the fact that the preferred size of the Fe dots is hardly changed upon annealing clearly indicates that the energy per atom $`E`$, equivalently, the chemical potential, of the Fe dot assembly must have a local minimum at the observed optimal size for a given coverage. This suggests that the qualitative form of the energy $`E`$ follows the generic behavior shown in Fig. 4 for a given coverage: When the volume approaches zero, $`E`$ is close to the energy of one adatom; when the dot grows very large, $`E`$ is about the energy per atom in a bulk Fe crystal, which is much smaller than the energy of one adatom. Between the two limits there exists a local minimum, mainly induced by the elastic relaxation caused by the discontinuity of the intrinsic surface stress tensor at the dot edges review . Secondly, the observation of different optimal sizes at different coverages shown in Fig. 2 suggests that the local energy minimum is a function of the iron coverage, implying that dot-dot interactions should contribute substantially to the energy $`E`$. Based on these considerations, we can express $`E`$ as a function of both the iron coverage $`\theta `$ and the dot size $`V`$ within the mean field approximation as

$$E(\theta ,V)=E_0(V)+P(\theta ,V),$$
(1)
where $`E_0`$ is the self-energy of an isolated dot and $`P`$ is the dot-dot interaction energy.

The interaction energy per atom,$`P(\theta ,V)`$, is originated from misfit-induced strain within the system. By adopting the standard dipole-dipole interaction between two dots as $`P(V_1,V_2)V_1V_2/r^3`$review , we have

$$P(\theta ,V)=g\frac{\theta ^{3/2}}{V^{1/2}},$$
(2)
where $`g`$ is a constant related to the interaction strength and the effective number of the nearest neighbors of a given dot. We note that, on one hand, the long range dot-dot interaction is important enough in influencing the energy per atom of the dot assembly, thereby the existence of an optimal dot size; on the other hand, it is not strong enough to cause the center of mass of a given dot to change, thereby the lacking of the spatial order among the dotsLiu .

The self-energy $`E_0(V)`$ of an isolated dot of size $`V`$ includes all possible contributions to the dot, such as the strain energy, interface energy, step energy, and kink energy, but it does not depend on the total coverage $`\theta `$. In principle, an explicit, but rather complex, expression can be obtained from existing studies of strain-mediated dot formation in heteroepitaxy, such as Eq. (4.15) of Ref.review . But surprisingly, our detailed analysis based on Eq. (4.15)review  shows that no qualitatively correct fitting could be obtained for the range of iron coverages explored in the present experiment. In the following, we choose the simple expression

$$E_0(V)=bV^\alpha .$$
(3)
as an alternative phenomenological approach, and show that the corresponding energy per atom is capable of explaining many of the salient features observed in the experiment.

First, by setting the first derivative of the energy per atom $`E(\theta ,V)`$ with respect to the volume $`V`$ equal to zero, we obtain a coverage-dependent local energy minimum, at the volume $`b\alpha V^{\alpha +1/2}=\frac{g}{2}\theta ^{3/2}`$. This, in turn, yields the relation between the density of the dots $`\theta /V`$ and the coverage $`\theta `$

$$b\alpha \left(\frac{\theta }{V}\right)^{(\alpha +1/2)}=\frac{g}{2}\theta ^{1\alpha }.$$
(4)
With $`\alpha `$ = 0.46, our theory gives an excellent fit to the experimental data as shown in Fig. 3(a). Next, assuming that the equilibrium shape of the dots does not change with the volume Duport , we can also obtain the dot height (h) and the diameter (D) as functions of the coverage, given by $`h\theta ^{0.52}`$ and $`D\theta ^{0.52}`$, respectively. These two relationships also find very good fits to the experimental data as shown in Figs. 3(b) and 3(c). The fact that the exponent $`\alpha `$ as a single fitting parameter yields very good fittings to the slopes of three independent sets of experimental observations hints on the physical validity of the choice made in Eq. (3). Whereas the precise underlying physics for the success of Eq. (3) remains to be explored, we suspect that it could be related to the unusually large lattice mismatch in the present system ltt . Finally, we note that the dispersion of the distribution of the island heights (or the diameters) shows little change with the coverage, whereas our mean field theory predicts that the dispersion should increase with the coverage, $`\mathrm{\Delta }h\theta ^{0.4}`$. This discrepancy is not too surprising, because within the mean-field theory, the resulting optimal size of the islands can be reliably obtained, but the exponent of the scaling function around the optimal value is typically overestimated, a trend well recognized in critical phenomena Wilson .

In summary, we have shown that, under proper growth conditions, the self-assembly of iron dots on the insulating surface of NaCl(001) leads to the formation of nanometer-scale magnetic iron dots with narrow size distributions in the absence of a wetting layer. We have also demonstrated that the vertical and lateral sizes of the dots can both be changed by the iron coverage while the narrow size distributions are preserved, even though the clustering is clearly in the strong coarsening regime. These striking observations have been interpreted successfully with a phenomenological mean-field theory.

We are grateful to E.W. Plummer for invaluable discussions. This work was supported in part by the LDRD of ORNL, managed by UT-Battelle, LLC for the USDOE ( DE-AC05-00OR22725), and by the NSF (DMR-0071893 (BW) and DMR-0105232 (JP)).


text:
# Star formation in the warped outer pseudoring of the spiral galaxy NGC 3642

## 1 Introduction

NGC 3642 was classified as unbarred (SAbc from RC3 - De Vaucouleurs et al. 1991, Sb(r)I from the RSA - Sandage & Tammann 1981) and has three rings (De Vaucouleurs & Buta 1980, hereafter DVB80). This combination was our main motivation for the study of this galaxy, as part of our analysis of ringed non-barred galaxies (Verdes-Montenegro et al. 1995, 1997, 2000). However, based on optical images, we have found that the outer ring is just a not very well defined pseudoring lying in the faint outer parts of the optical disk. Furthermore, when compared with the statistical study performed by Athanassoula et al. (1982) we find that the size ratio of the outer and inner ring is far larger than the one needed to link them to the ratio of a ring at the outer Lindblad radius and a ring near corotation due to a single pattern speed.

New HI data of this galaxy show that its outer disk is warped, and that star formation occurs beyond R25 (at 1.4 $`\times `$ R25, where R25 = 161′′ $`\pm `$ 3′′ from RC3). Star formation in the faint outer parts of galactic disks is a relevant and little understood issue: very few examples have been so far studied (Ferguson et al. 1998 and references therein). These combined characteristics make NGC 3642 an interesting galaxy, and in this paper we present a detailed optical and HI (WSRT) study.

The systemic heliocentric radial velocity of this galaxy, obtained from our HI data (Sect. 2.3) is 1572 km s-1, and corrected to the centroid of the Local Group is 1682 km s-1. A Hubble constant of 75 km s-1 Mpc-1 gives a distance of 22.4 Mpc for NGC 3642.

## 2 Observations and data analysis

### 2.1 Photometry

We have analysed optical data of NGC 3642 in the B and R bands. Boroson (1981) has obtained photographic surface photometry of this galaxy in the B-band. We have supplemented Boroson’s data with CCD surface photometry extracted from an R-band image obtained from the ING archive. This image was obtained on February 3, 1995 at the INT telescope with a 1024 $`\times `$ 1024 pixel TEK1 CCD camera, each pixel being 24$`\mu `$m in size and 0.58′′. The measured seeing is 1$`\text{.}^{\prime \prime }`$7. Flatfielding was achieved by producing a superflatfield from the median of all obtained sky flats. Several fields of Landolt phometric standards were also retrieved from the ING archive, for the calibration of the R image. We find good agreement with the profile given by Boroson for the inner parts. For the outer parts some disagreement is found and we checked that this is not due to a bad sky subtraction of our R image.

We have also inspected the HST image published by Barth et al. (1998) with the F547M (V band) filter, which we retrieved from the archive. This image does not extend with sufficient signal to noise ratio to the inner ring but is well suited to explore the central part of NGC 3642, including the nuclear structure.

### 2.2 Radio observations

We have observed NGC 3642 in the 21 cm line of HI with the Westerbork Synthesis Radio Telescope (WSRT111The Westerbork Synthesis Radio Telescope is operated by the Netherlands Foundation for Radio Astronomy with financial support from the Netherlands Foundation for the Advancement of Pure Research (N. W. O.)) in 1983. We used 40 interferometers with spacings ranging from 36 m to 2844 m in steps of 72 m. This results in a synthesized beam of 13$`\text{.}^{\prime \prime }`$0 $`\times `$ 15$`\text{.}^{\prime \prime }`$1 ($`\alpha \times \delta `$) and first grating response at 10 $`\times `$ 12 ($`\alpha \times \delta `$). We used a digital backend (Bos et al. 1981) resulting in 31 channel maps at heliocentric velocities 1463.3 to 1710.7 km s-1. Hanning smoothing was applied on–line, giving a velocity resolution of twice the channel spacing of 8.2 km s-1.

The data were edited and calibrated as explained in Verdes-Montenegro et al. (1995). A rms noise level of $``$ 1.13 mJy/beam was achieved after 12 hours of integration. In order to get a higher signal–to–noise ratio in the integrated HI distribution and associated radial velocity field, we convolved the map data with a gaussian, leading to a beam size of 21$`\text{.}^{\prime \prime }`$4 $`\times `$ 18$`\text{.}^{\prime \prime }`$4 ($`\alpha \times \delta `$). Primary beam corrections have been applied to our maps.

The total HI emission in this galaxy was measured to be 66.16 $`\pm `$ 5.98 Jy km s-1 by Haynes & Giovanelli (1991) with the Green Bank radiotelescope and corrected for source to beam extent to 83.84 Jy km s-1. The HI spectrum obtained by integrating the flux density in each channel map over an area containing the line emission is shown in Fig. 1, and the derived total flux is 44 Jy km s-1. Hence we are missing emission, and comparison with the single dish spectrum indicates that it comes from those velocities where the emission is more extended (velocity range 1545 - 1590 km s-1, see Sect. 3.2). From our spectrum we derive a systemic velocity, obtained as the mean velocity at 20% of the peak, of 1572 km s-1. The spectrum has a FWHM at the 20% level of 74 km s-1.

## 3 Results and discussion

### 3.1 Optical emission

In Fig. 2a we show the R band image of NGC 3642. The central parts have been saturated in order to show better the outer spiral structure, that extends up to 220′′. We have marked with an ellipse the outer ring size given by DVB80, $``$ 5 $`\times `$ 4. We do not find a complete outer ring in our image, although part of the ellipse seems to coincide with the wide outer spiral structure, that extends from about 145 to 325. It is not obvious whether the thinner structure that extends from about 90 to 145 at smaller radius is part of the wider structure, or whether it corresponds to an inner spiral.

Boroson (1981) described the galaxy as having “at first glance a large bulge and a very faint disk”, but indicated that “after closer examination it turns that much of the bulge region shows spiral structure, and is actually a section of the disk with a shorter scalelength”. This (flocculent) spiral structure is shown in Fig. 2b, extending between radius 20′′ and 50′′, where the faint outer disk starts.

The spiral structure of NGC 3642 has a larger contrast in the outer parts than in the inner parts, and the difference between the two parts is more evident after substraction of a smoothed R image from the original one, that enhances the small scale structure (Fig. 2c). DVB80 reported the existence of an inner ring with a size of 25$`\text{.}^{\prime \prime }`$8 $`\times `$ 22$`\text{.}^{\prime \prime }`$8. This ring is faint in our R band image, and enhancement with unsharp masks is required in order to isolate it (Fig. 2d). We have measured its size as 28$`\text{.}^{\prime \prime }`$6 $`\times `$ 25$`\text{.}^{\prime \prime }`$4, in good agreement with DVB80, and its position angle as 0.

The third ring reported by DVB80, a nuclear one with a diameter of 6′′, cannot be resolved in our R band image, where the central part is saturated due to the bright LINER nucleus (Heckman 1980). The inner 20′′ of NGC 3642 have bright emission in the HST-V image, that shows sections (three quarters of a circle) of what could be the nuclear ring, with a diameter of 5′′. However a sharpening of the HST-V image (Fig. 2d) suggests that instead of a ring this emission enhancement is part of the nuclear spiral structure, since a single m=1 spiral is observed, reaching the active nucleus, and extending between 1′′ and 6′′ from the center (109 to 654 pc at 22.4 Mpc), with several dust lanes crossing the outer part of the minispiral. At larger radii the arm seems to enter the flocculent structure that we have indicated above.

Hence NGC 3642 exhibits three different types of spiral structure: a nuclear one-armed minispiral, an inner flocculent spiral and an outer spiral that winds into an outer pseudoring extending for about half a circle. The winding direction of all three is the same (clockwise).

The radial brightness profile shown in Fig. 3 has been calculated from the R band by integrating over circular annuli with thickness of 1$`\text{.}^{\prime \prime }`$5 on the deprojected image of the galaxy (see Table 1 for the deprojection parameters), after blotting of the stars that project on the galaxy. Two bumps are evident on the profile at radii of 110′′ and 170′′, respectively, both associated with bright features in the spiral structure, as indicated in the upper right inlay. It is not surprising that such features appear in the R band since it includes the H$`\alpha `$ line.

Ellipticities ($`ϵ`$) and position angles (p.a.) of the ellipses fitted to the R band isophotes are shown in Fig. 4a and b respectively as a function of radius, up to the larger radius at which we were able to perform such fitting (radius of 70′′). The results are very noisy since the isophotes are strongly influenced by the spiral arms. The ellipticity shows a mean value of 0.06 $`\pm `$ 0.02, corresponding to an inclination of $`i`$ = 20.4 $`\pm `$ 3.5, and the position angle does not show particularly abrupt changes with radius as expected for a bar embedded in a disk. There is also no obvious difference of the value of the inclination as function of radius, although the inner parts seem a bit more circular. Furthermore, there are no clearly marked minima in the light distribution corresponding to the Lagrangian points L4 and L5, a characteristic of the morphology of an oval distortion, nor is there any other sign of the presence of a bar. The low inclination makes the morphological determination of the position angle highly uncertain. Inside the inner pseudoring (r $`<`$ 14′′) the p.a. ranges between 125 and 160. From 14′′ to 20′′, i.e. up to the beginning of the spiral structure, the p.a. drops steeply from 180 down to 110. Then it jumps up to $``$ 158 $`\pm `$ 44, between 20′′and 70′′.

### 3.2 Atomic gas distribution and kinematics

In Fig. 5 we display the channel maps containing the HI emission at the indicated heliocentric velocities, produced with a beam of 21$`\text{.}^{\prime \prime }`$$`\times `$ 18$`\text{.}^{\prime \prime }`$4. Each channel has a width of 8.2 km s-1 and an rms noise of 1.9 K, implying a level of 3.4 $`\times `$ 1019 cm-2. The HI contours have been superimposed on the R-band image. There is a good correspondence between the brighter areas of HI emission and the outer spiral structure, as is also seen in the integrated HI map (Fig. 6). The brighest HI/optical features have column densities of 2.7 $`\times `$ 1021 cm-2. At a level of 1.5 $`\times `$ 1020 cm-2 the galaxy extends $``$ 240′′ in radius (1.5 $`\times `$ R25), offcentered by $``$ 25′′ to the west. This is plausibly a lower limit to the real extent since we are missing nearly half of the flux (see Sect. 2.3). There is a central HI depression with N(HI) $``$ 1.0 $`\times `$ 1019 cm-2. The azimuthally averaged radial distribution of the HI column density is shown in Fig. 4c. It has been obtained by averaging the two-dimensional HI distribution in elliptical rings using the geometrical parameters given in Table 1.

Numerous channel maps show emission at both receding and approaching sides of the galaxy. This results in the velocity field shown in Fig. 7a, where redshifted emission is observed in the blueshifted part of the velocity field and viceversa. Between 40′′ and 70′′ the isovelocity contours show kinks that correlate with the HI clumps associated to the spiral arms. But most of the perturbations for the largest radii are characteristic of a warped disk. We have derived the kinematical major axis based on the most symmetrical part of the velocity field (r $`<`$ 40′′). The major axis has been obtained as the direction perpendicular to the isovelocity contours, resulting in a value of p.a. = 122.6. This direction coincides well with the line joining the peak of the blueshifted and redshifted velocities in the inner parts (Fig. 7). It differs significantly from the optically derived one, although this is not surprising due to the low inclination and the effect of the spiral structure on the isophotal fitting. We will thus adopt the p.a. of NGC 3642 derived from the velocity field. The isocontour at the systemic velocity (1572 km s-1, Sect. 2.3) intersects the major axis direction at the optical center position, which we will also adopt as the kinematical center of the galaxy (Fig. 7b).

A cut along the major axis direction (Fig. 8, filled circles) shows clearly the radial velocity inversion observed in the velocity field. As we show in the next section, this can be modelled as a warped disk: due to the low inclination of the galaxy, any small difference in the inclination between the inner and outer parts can produce a large rotation of the major axis. We have overlapped the velocities obtained from a cut along the major axis of the velocity field. Declining velocities are observed for radii larger than 100′′ down to values even lower than the systemic velocity.

### 3.3 Description of the warp model

We have modelled the velocity field fitting a tilted ring model (see Begeman 1987; the ROTCUR and GALMOD tasks in Gipsy). We have tried different combinations of the inclination, position angle and rotational velocity, as we describe below. In all cases the galaxy has been divided into concentric rings, each of them with a width of 5′′ along the major axis and a central position fixed to the optical center (Table 1). Points within a sector of 30 from the minor axis were excluded from the fits. The expansion velocities were set to zero and the systemic velocity fixed to the central velocity of the HI spectrum (Sect. 2.3, Table 1).

We could not allow the inclination and position angle to vary freely, since the errors for the inclination and for the rotation velocity are large for all radii. The galaxy is seen so face-on that a determination of the rotation curve is almost impossible : the uncertainties of the inclination correction are simply too large. Nevertheless, we have attempted to model the galaxy, assuming a flat rotation curve, as a continuation of what we find in the inner parts. From this we can estimate an approximate amplitude of the warping. In Fig. 9 we show the best combination of position angle and inclination that we have found in order to reproduce the velocity field. The inversion in the velocity field is reproduced although we could not obtain its most extreme features (Fig. 10). The fit along the major axis is however very accurate as shown in Fig. 8, where we compare the rotation curve along the major axis for the model and original data. The rotation curve, inclination and position angle obtained from our modelling are plotted in Fig. 9. The inclination derived from the modelling of the HI emission fits well with the one derived on the base of R band isophotal fitting. A warping of the disk of $``$ 25 is needed in order to reproduce the velocity field of NGC 3642. The bumps in the position angle at $``$ 90′′ and 140′′ coincide with the location of the spiral arms.

An estimate of the rotational velocity of this galaxy can be obtained from the Tully-Fisher relation (Pierce & Tully 1992). Using the observed $`M_B`$ we get $`W_R`$ = 379 km s-1, which, combined with the observed extrema of the radial velocities, implies an inclination of 17, in good agreement with our adopted value given the scatter in the Tully-Fisher relation.

### 3.4 Comment on the z-velocities in warps

The galaxies most useful for rotation curves studies are inclined typically between 50° and 80° (e.g. Begeman et al. 1991), and if a warp is present in the outer parts, its effect is traditionally modelled using a tilted ring model (Rogstad et al. 1974, Bosma 1978, 1981, and Begeman 1987). Such a procedure, in particular the one implemented in ROTCUR, minimizes the importance of any peculiar motions, and maximizes the circular velocity in each ring independent of its spatial orientation. In reality, there ought to be some component of the motion of the inclined outer rings towards the main plane defined by the inner parts of the galaxy. We can get an upper limit to this motion by attributing to it all the peculiar motion seen in Fig. 8 (p-v-diagram along the major axis), which is 30 km s-1. This can be compared to the rotation velocity, which is certainly larger than 80 km s-1, and more likely of the order of 150 km s-1 or even 190 km s-1 (if the TF relation holds and the Sb(r)I classification of the RSA is retained). So at first order, z-motions of warps can be neglected.

### 3.5 Environment

We have studied the neighborhood of NGC 3642 within a radius of 500 kpc and 500 km s-1 and show the results in Fig. 11. From the NED222The NASA/IPAC extragalactic database (NED) is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration database we found that NGC 3642 has 10 neighbors with known redshift in the box under consideration. Among them, NGC 3610 and NGC 3613 have magnitudes similar to that of NGC 3642, while the rest are fainter by more than 2 magnitudes. The closest companion, NGC 3610, is located at an apparent separation of 34$`\text{.}^{}`$7, or 226 kpc, and has clear signs of merger: this is a quite perturbed elliptical galaxy, with a prominent, edge-on stellar disk along the major axis in the inner regions (Scorza & Bender 1990). It has several shells, plumes and boxy outer isophotes indicating a past merger involving two disk galaxies dated 4 Gyr ago (cf Whitmore et al. 1997) based on UBV colors. NGC 3613 is a disky elliptical with the disk fully embedded in a spheroidal component (Michard & Marchal 1994).

Garcia (1993) considered a larger scale, and finds that NGC 3642 belongs to a group composed by 5 galaxies. The brightest members of this group are NGC 3642 and 3610, followed by NGC 3619, a magnitude fainter, and by NGC 3674 and NGC 3683, which are two magnitudes fainter. The second nearby companion of NGC 3642 in mass, NGC 3619, is an early type disk galaxy with shells (cf. Schweizer & Seitzer 1988, Van Driel et al. 1989). The velocity dispersion of the group is 600 km s-1. NGC 3683a is close to the group, although with a larger redshift (see Fig. 11).

## 4 Discussion and conclusions

We have investigated the optical morphology as well as the gas distribution and kinematics of the nearly face-on spiral galaxy NGC 3642. As is typical for spiral galaxies, the HI gas extends a factor of 1.5 beyond the optical radius. What is more unusual is that the outer HI regions composing the pseudoring are warped. Furthermore the HI is more extended in the western side, where the more external star forming area is located, and higher column densities are associated to these regions. These perturbations affect only to the outer disk, since the kinematics within the main body conforms well to an ordinary differentially rotating disk.

The inner kpc of NGC 3642 shows a one-armed spiral structure. Spirals in the inner kpc of active galaxies have been considered as drivers of accretion into the central engine (Regan & Mulchaey 1999). While barred galaxies tend to show bisymmetrical patterns (Athanassoula 1992) non barred ones can have more irregular structures. The observed spiral could be driving nuclear accretion.

Both phenomena, nuclear activity and star forming warped outer parts, might be linked to the fact that NGC 3642 is located in a rich environment. Given the clear signs of merging in its closest bright companions, it is reasonable to assume that NGC 3642 also underwent recently an interaction, or, more likely, captured a low-mass, gas-rich dwarf, which disintegrated, and whose debris now form the outer HI and the outer patchy spiral structure in NGC 3642 itself. The differences in the kinematics between the approaching and receding sides of the galaxy are probably due to the accretion of a companion. The low surface brightness extension of the galaxy, as seen in the photometric structure, also points in this direction. Such low surface brightness extensions have been studied before, (e.g. Longmore et al. 1979, Bosma & Freeman 1993), and are not so uncommon. Longmore et al. (1979) indicated already that under favorable initial conditions, specially a low relative velocity encounter, a considerable fraction of the gas content of a galaxy can be drawn out in an interaction, and yet participate in subsequent star formation. We could therefore see star formation in the denser regions of redistributed material. i.e. in the outer pseudoring, which might be in the process of forming a continuous ring.

In conclusion, the observed outer and perhaps even the nuclear structure of the unbarred galaxy NGC 3642 can be attributed to accretion of an HI-rich dwarf galaxy that could be acquired from its rich environment. Star formation was triggered in its outskirts by incorporation of external gas that produced also a pronounced warp in the gaseous disk. Our modeling of the velocity field of this face-on galaxy allowed us to set a constrain on the amplitude of the z-motions in warps, which we find to be smaller than about 30 km s-1. It could well be that the formation of some of the faint outer extensions of galactic disks is due to this kind of destruction of a low-mass, gas-rich dwarf during capture.

###### Acknowledgements.

LV–M acknowledges the hospitality and financial support of the Observatoire de Marseille, where part of this work was carried out. LV–M is partially supported by AYA 2000-1564 and Junta de Andalucía (Spain).


text:
# Temperature dependence of density profiles for a cloud of non-interacting fermions moving inside a harmonic trap in one dimension

## I Introduction

The achievement of Bose-Einstein condensation in ultracold gases of bosonic atoms has given great impulse to the study of dilute quantum gases inside magnetic or optical traps . Similar cooling and trapping techniques are being used to drive gases of fermionic alkali atoms into the quantum degeneracy regime . It is also possible to experimentally generate and study strongly anisotropic atomic fluids effectively approaching dimensionality $`D=2`$ or $`D=1`$, at very low temperature and with very high purity, inside magnetic traps where the transverse and the longitudinal confinement are vastly different .

In the relevant conditions of temperature and dilution the atomic interactions become negligible in a gas of fermions placed in a single Zeeman sublevel inside a magnetic trap . One can obtain in this way a close laboratory realization of an inhomogeneous, non-interacting Fermi gas, which has special significance in regard to the foundations of density functional theory and to tests of the Thomas-Fermi approximation . In the quasi-onedimensional (1D) case the ground state and the excitation spectrum of such an ideal gas of spin-polarized (or ”spinless”) fermions can be mapped into those of a gas of hard-core impenetrable bosons . The latter is known as the Tonks gas, from the early work of Tonks  on the equation of state of hard-object adsorbates. Advances in atom waveguide technology, with potential applications to atom interferometry and integrated atom optics, especially motivate theoretical studies of dilute gases in a regime where the dynamics becomes essentially 1D .

The wave functions of a cloud of spinless fermions under harmonic confinement are, of course, Slater determinants of harmonic-oscillator single-particle orbitals. Their representation in terms of Hermite polynomials has, however, limited usefulness for numerical calculations on mesoscopic clouds . Brack and van Zyl  have developed a more powerful method for non-interacting fermions occupying a set of closed shells under isotropic harmonic confinement in D dimensions, leading to analytical expressions for the particle and kinetic energy densities at temperature $`T=0`$ in terms of Laguerre polynomials. These expressions are especially useful for numerical calculations on isotropic systems in $`D=2`$ and $`D=3`$. A Green’s function method, which altogether avoids the use of wave functions in favor of the matrix elements of the position and momentum operators, has been developed for similar purposes  and extended to calculations of the pair distribution function in the ground state for the 1D ideal Fermi gas . This method has also been extended to the ground state of ideal Fermi gases under harmonic confinement of arbitrary anisotropy in higher dimensionalities .

The purpose of the present paper is to extend the Green’s function method to an inhomogeneous 1D cloud of non-interacting fermions at finite temperature and to illustrate its usefulness by numerical calculations of the particle and kinetic energy density profiles as functions of temperature in the case of harmonic confinement. Analytical expressions have been derived by Wang for the same system at high temperature, where the chemical potential is lower than the single-particle ground-state energy. The case of non-interacting fermions under 3D harmonic confinement has been treated by Schneider and Wallis  through the use of Laguerre polynomials. The emphasis of our numerical calculations will be to illustrate how the characteristic quantum features of the fermion cloud, i.e. its shell structure and the spill-out of particles at the boundaries beyond the Thomas-Fermi radius, evolve with increasing temperature as a semiclassical regime is being approached.

The paper is organized as follows. Section II reports some essential definitions, starting from the one-body density matrix at $`T0`$, and gives the expressions of the density profiles in the semiclassical regime. Sections III and IV present the essential details of our method and our numerical results, respectively. A brief summary concludes the paper in Section V. Although we refer through the paper to the results as being appropriate to spinless fermions in 1D, they are equally valid for the Tonks gas of impenetrable bosons.

## II ESSENTIAL DEFINITIONS AND THE SEMICLASSICAL LIMIT

The generalized grand-canonical density matrix for fermions  can be written as

$$D(x_1,x;\beta ,\mu )=\underset{i=1}{\overset{\mathrm{}}{}}\frac{1}{1+\mathrm{exp}[\beta (E_i\mu )]}\psi _i^{}(x_1)e^{i\widehat{p}(xx_1)}\psi _i(x_1).$$
(1)
Here $`\beta =1/k_BT`$ and $`\mu `$ is the chemical potential, while $`\psi _i`$ and $`E_i`$ are the single-particle orbitals and the corresponding energy eigenvalues. The zero temperature limit of Eq. (1) leads to the Dirac density matrix for $`N`$ ideal spinless fermions,

$$\rho (x_1,x)=\underset{i=1}{\overset{N}{}}\psi _i^{}(x_1)e^{i\widehat{p}(xx_1)}\psi _i(x_1).$$
(2)
The particle density profile $`n(x)`$ of the gas at temperature $`T`$ and chemical potential $`\mu `$ is the zero-order moment of the matrix $`D(x,x_1)`$,

$`n(x)`$ $`=`$ $`D(x_1,x;\beta ,\mu )|_{x_1=x}`$ (3)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{1+\mathrm{exp}[\beta (E_i\mu )]}}\psi _i|\delta (x\widehat{x})|\psi _i.`$ (4)
The kinetic pressure $`P(x)`$ is given by a specific second-order moment of $`D(x,x_1)`$ ,

$`P(x)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{m}}{\displaystyle \frac{^2}{x_1^2}}D(x_1,x;\beta ,\mu )|_{x_1=x}`$ (5)
$`=`$ $`{\displaystyle \frac{1}{2m}}{\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{1+\mathrm{exp}[\beta (E_i\mu )]}}\psi _i|\widehat{p}^2\delta (x\widehat{x})+\delta (x\widehat{x})\widehat{p}^2|\psi _i.`$ (6)
The kinetic pressure $`P(x)`$ is twice the kinetic energy density of the fermion cloud. In these equations $`\widehat{p}`$ and $`\widehat{x}`$ are the momentum and position operators, and the chemical potential $`\mu `$ is determined from the condition

$$n(x;\beta ,\mu )𝑑x=N$$
(7)
where $`N`$ is the average number of particles.

In the semiclassical regime the particle density and the kinetic pressure of an ideal Fermi gas confined by a 1D potential $`V(x)`$ can be calculated by treating the energy levels as a continuum. The appropriate condition of validity is that the level spacing $`\mathrm{\Delta }E`$ is sufficiently smaller than the thermal energy $`k_BT`$. The appropriate expressions in the grand-canonical ensemble are

$`n_{sc}(x)`$ $`=`$ $`{\displaystyle 𝑑E\frac{dp}{2\pi \mathrm{}}\delta \left(E\frac{p^2}{2m}V(x)\right)\frac{1}{1+\mathrm{exp}[\beta (E\mu )]}}`$ (8)
$`=`$ $`{\displaystyle \frac{dp}{2\pi \mathrm{}}\left\{\mathrm{exp}\left[\beta \left(\frac{p^2}{2m}+V(x)\mu \right)\right]+1\right\}^1}`$ (9)
and

$$P_{sc}(x)=\frac{dp}{2\pi \mathrm{}}\frac{p^2}{m}\left\{\mathrm{exp}\left[\beta \left(\frac{p^2}{2m}+V(x)\mu \right)\right]+1\right\}^1.$$
(10)
The chemical potential is again determined by normalization to the average number of fermions.

## III The Green’s function method

Equations (4) and (6) can be rewritten in terms of the Green’s function $`\widehat{G}(x)=lim_{\epsilon 0^+}(x\widehat{x}+i\epsilon )^1`$ in coordinate space,

$`n(x)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}\underset{\epsilon 0^+}{lim}\mathrm{Im}{\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{1+\mathrm{exp}[\beta (E_i\mu )]}}\psi _i|\widehat{G}(x)|\psi _i`$ (11)
$`=`$ $`{\displaystyle \frac{1}{\pi }}\underset{\epsilon 0^+}{lim}\mathrm{ImTr}(𝒯\widehat{G}(x))`$ (12)
and

$`P(x)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}\underset{\epsilon 0^+}{lim}\mathrm{Im}{\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{1+\mathrm{exp}[\beta (E_i\mu )]}}\psi _i|{\displaystyle \frac{\widehat{p}^2}{m}}\widehat{G}(x)|\psi _i`$ (13)
$`=`$ $`{\displaystyle \frac{1}{\pi m}}\underset{\epsilon 0^+}{lim}\mathrm{ImTr}(𝒯\widehat{p}^2\widehat{G}(x)).`$ (14)
Here we have used a matrix formalism, by introducing the temperature matrix $`𝒯`$ whose diagonal elements $`[𝒯]_{i,i}=1/\{1+\mathrm{exp}[\beta (E_i\mu )]\}^1`$ contain the statistical Fermi factors, while the off-diagonal elements $`[𝒯]_{i,j}`$ with $`ij`$ are null.

The evaluation of the expressions in Eqs. (12) and (14) is carried out by an immediate extension of the procedure developed in  for the case of zero temperature. The trace of a generic matrix $`A`$ is connected to the elements of the inverse matrix $`A^1`$ by the relation

$$\mathrm{Tr}A=\frac{}{\lambda }[\mathrm{ln}det(A^1+\lambda 𝕀)]|_{\lambda =0}.$$
(15)
We get from Eqs. (12) and (14)

$$n(x)=\frac{1}{\pi }\underset{\epsilon 0^+}{lim}\mathrm{Im}\frac{}{\lambda }[\mathrm{ln}det(x\widehat{x}+\lambda 𝒯+i\epsilon )]|_{\lambda =0}$$
(16)
and

$$P(x)=\frac{1}{\pi m}\underset{\epsilon 0^+}{lim}\mathrm{Im}\frac{}{\lambda }[\mathrm{ln}det(x\widehat{x}+i\epsilon +\lambda 𝒯\widehat{p}^2)]|_{\lambda =0}.$$
(17)
In the specific case of harmonic confinement we make use of the representation of the position and the momentum operators in the basis of the eigenstates of the harmonic oscillator with energy $`E_n=(n1/2)\mathrm{}\omega `$: that is, $`\widehat{x}=(a+a^{})/\sqrt{2}`$ and $`\widehat{p}=i(a^{}a)/\sqrt{2}`$ with $`a|\psi _n=\sqrt{n1}|\psi _{n1}`$ and $`a^{}|\psi _n=\sqrt{n}|\psi _{n+1}`$. The tridiagonal form of the matrices representing $`\widehat{x}`$ and $`\widehat{p}`$ allows us to express the determinants in Eqs. (16) and (17) by a recursive algorithm  as products of an infinite number of matrices having dimension $`1\times 1`$ for the particle density and $`2\times 2`$ for the kinetic pressure. If we write the matrix $`𝒦^\nu =\widehat{x}\lambda 𝒯\widehat{p}^\nu `$ in the tridiagonal form

$$𝒦^\nu =\left(\begin{array}{cccc}𝒜_1& _{1,2}& & \\ _{2,1}& 𝒜_2& _{2,3}& \\ & \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right)$$
(18)
where the indices refer to blocs of dimension $`1\times 1`$ if $`\nu =0`$ and $`2\times 2`$ if $`\nu =2`$, we obtain

$$det(x𝒦^\nu +i\epsilon )=\underset{j=1}{\overset{\mathrm{}}{}}det(x\stackrel{~}{𝒜}_j+i\epsilon ).$$
(19)
Here $`\stackrel{~}{𝒜}_1=𝒜_1`$ while the renormalized blocs $`\stackrel{~}{𝒜}_j`$ for $`j>1`$ are given by the recursive formula

$$\stackrel{~}{𝒜}_j=𝒜_j+_{j,j1}(x\stackrel{~}{𝒜}_{j1}+i\epsilon )^1_{j1,j}.$$
(20)
As in the calculation of density profiles at $`T=0`$  the evaluation of the determinant in Eq. (19) converges quite rapidly, yielding accurate results for mesoscopic fermion clouds with a moderate amount of numerical effort.

## IV Numerical results

We report in this section some numerical results for the particle density and the kinetic pressure that we have obtained by the Green’s function method. Our main purpose is to illustrate the approach of the profiles to their semiclassical expressions and we consider for clarity clouds containing rather small average numbers of particles, i.e. $`N=4`$ and $`N=20`$. Some sharpening of the structures in the profiles should be expected in a canonical ensemble viewpoint, where fluctuations in the particle number are suppressed.

In the calculations reported below we have approximated the determinant in Eq. (19) by the product of its first $`M`$ terms with $`M=10^7`$ as an upper bound. This choice allows for up to $`10^7`$ thermally excited states and the corresponding value of the spectral resolution parameter is chosen as $`\epsilon =10^3`$. These two choices evidently limit to some minor extent the accuracy of our numerical results relatively to the formally exact theory given in Sect. III. Of course, for such small numbers of particles one may also approach the same calculations by using directly the expression of the harmonic-oscillator orbitals in Eq. (1). However, the routines that are available to evaluate Hermite polynomials also involve numerical approximations which are increasingly severe for polynomials of high degree, as needed to describe the gas at finite temperature.

The first step in our calculations is to evaluate the chemical potential at fixed $`N`$ by means of a self-consistent solution of Eqs. (7) and (9. The resulting semiclassical value of $`\mu `$ is used in Eqs. (16) and (17) for the evaluation of the particle density and kinetic pressure profiles. The procedure is justified a posteriori by the fact that integration of the spatial density $`n(x)`$ obtained at this stage reproduces the correct average number $`N`$ of fermions.

In Figures 1 and 2 we report the particle density profiles at various temperatures for $`N=4`$ and $`N=20`$, respectively. The corresponding profiles of kinetic pressure are given in Figures 3 and 4. In each Figure the left-hand panel shows how the quantum effects, which consist of oscillations from shell structure in the profiles and in particle spill-out and negative kinetic pressure at the edges, are washed away by thermal excitations. The right-hand panels show instead how the profiles approach the semiclassical regime: the latter already holds to a good approximation when $`k_BT\mathrm{}\omega `$, although some trace of a negative kinetic pressure in the spill-out region remains in the “exact” profile at this temperature.

## V CONCLUSIONS

In conclusion, in this paper we have extended to finite temperature the Green’s function method that was previously developed for the formally exact evaluation of ground-state properties of an inhomogeneous ideal Fermi gas in 1D. The method dispenses with knowledge of the single-particle orbitals in favour of the single-particle energy levels and of the matrix elements of the position and momentum operators. Our numerical applications to clouds containing small numbers of particles under harmonic confinement have shown that the commonly used semiclassical approach to a thermally excited fermion cloud yields (i) accurate values of the chemical potential except at the lowest temperatures, and (ii) a reliable account of particle and kinetic energy density profiles once their main quantum features become essentialy washed out at $`k_BT\mathrm{}\omega `$.

Two further comments are in order. Firstly, the profiles shown in Figures 1 and 2 also describe (aside from an immediate rescaling of the units used on the axes) the momentum distribution of the fermion cloud under harmonic confinement. Secondly, the same profiles also refer to the Tonks gas of impenetrable bosons in 1D, as was explicitly shown by Yang and Yang .

###### Acknowledgements.

This work was partially supported by INFM through the PRA-Photonmatter program.


text:
# Steering effect on the shape of islands for homoepitaxial growth of Cu on Cu(001)

## Abstract

The steering effect on the growth of islands is investigated by combining molecular dynamics (MD) and kinetic Monte Carlo (KMC) simulations. Dynamics of depositing atoms and kinetics of atoms on a substrate are realized by MD and KMC, respectively. The reported experimental results on the asymmetric island growth \[van Dijken et al., Phys. Rev. Lett. 82, 4038 (1999).\] is well reproduced. A salient phenomenon, the reversal of the asymmetry, is found as the island size increases, and attributed to the asymmetric flux on the lower terrace of island.

The growth of thin film is an essential step for many modern technologies and scientific investigation, and a lot of efforts have been made to understand and tailor the growth process. There have been many studies on the effect of energeticBauer  and kineticReviews  variables to the growth of thin film, while minor attention has been paid to the role of dynamic variables such as deposition conditions. Recently, van Dijken et al.Dijken  reported the growth of rectangular Cu islands on square-symmetric Cu(001) when the deposition was made at grazing incidence angle. Furthermore, they controlled the shape of Co island on Cu(001) by varying the deposition angle and successfully manipulated its magnetic anisotropyDijken2 . Their work clearly reveals the importance of dynamic variable and expand the adjustable parameters for the growth of thin film.

An atomistic picture for the asymmetric island growth was also proposed by van Dijken et al.Dijken ; Incident atoms see a modified potential by pre-existing island and as a result the incoming flux is focused on the upper terrace of the island near the front edge, while depleted on the lower terrace near the rear edge of the island. Such steering effect is argued to result in the shortened edge length along the deposition direction, $`x`$-axis, and leaves the edge perpendicular to the deposition direction, $`y`$-axis, relatively longer. The model, however, is based on a qualitative argument without any detailed kinetic description on island growth from the inhomogeneous flux distribution, and the steering effect is considered only in the incident plane rather than in the three dimensional space as it should be. Independently, Zhong et al.Zhong  proposed that the enhanced corner crossing diffusion at the front edge due to the transient mobility of the deposited atom increased the growth speed along the y-direction relative to that along the $`x`$-direction. However, the very existence and the role of the transient mobility still remains controversialSanders . Steering effect has also been studied in the roughening in amorphous thin film growth by MDLandman  and non-linear stochastic equationRaible . For the growth of a crystalline film, only recently Montalenti and VoterMontalenti  reported the steering induced instability of Ag film on Ag(001), but the study was restricted to the roughening of the film on a substrate at 0 K and to the deposition at normal incidence.

The purpose of the present study is to investigate the steering effect on the growth of asymmetric islands by performing a realistic simulation combining KMC and MD simulation; MD simulation is executed to calculate the trajectory of depositing atoms in three dimensional space, only when a deposition event is selected in KMC. Before the next deposition event is selected, various diffusion events of adatoms are realized by KMC simulationClancy . In the present study, the asymmetric island shape observed in the experimental studyDijken  is well reproduced. Moreover, a salient phenomenon, the reversal of the asymmetry of island from elongation in the $`y`$-direction to the $`x`$-direction, is found as the coverage or the substrate temperature increases. The asymmetric shape of the island is attributed mainly to the asymmetric deposition flux on the lower terrace of island that in turn depends on its size.

In MD simulation, Lennard-Jones potential in the form of $`U(r)=4D[(\sigma /r)^{12}(\sigma /r)^6]`$ with $`D=`$0.4093 eV and $`\sigma =`$2.338 $`\AA `$ is used for the interaction between the incident atom and surface atomsSanders , and Verlet algorithm is adapted. The simulation box is composed of 6 layer high empty space in $`z`$-direction on an fcc(001) surface of 400$`\times `$400 atoms with periodic boundary condition in the plane. To reduce computational efforts, during deposition, the atoms on the substrate are frozen at each lattice site. (The surface lattice constant, $`a_0=`$2.56 $`\AA `$ and the interlayer spacing, $`d=`$1.805 $`\AA `$) The incident atom along  direction has initial kinetic energy of 0.15 eV determined from the melting temperature of Cu and follows the trajectory determined from the interaction potential starting at the initial height of $`10a_0`$ until the atom experiences the repulsive force by pre-existing atoms on the substrate. In the present simulation, no transient mobility is taken into account.

For the diffusion processes, 12 different ones are taken into account including those suggested by Furman et al.Furman , and their diffusion barriers are adopted. Some of the most influential diffusion barriers are listed in Table 1.

From the simulated morphology, the mound radius, which can directly be compared with the diffraction resultsDijken , is determined as the first zero of the height-height correlation function $`<h(𝐫)h(0)><h>^2`$. To increase the statistical reliability, simulations are performed 60 times under the same growth condition, and the mound radius is calculated over the whole set of simulated morphologies.

In Fig. 1 we summarize the simulation results for the growth of Cu atoms on Cu(001) under the identical growth conditions to those of the previous experimentDijken . (We denote the mound radius along the $`x`$-axis by $`L`$ and that along the $`y`$-axis by $`W`$.) For the normal deposition, $`W`$ and $`L`$ are identical at all coverages, while for the deposition at an angle, $`80^o`$ off the surface normal, the symmetry is broken and the difference between $`W`$ and $`L`$ becomes larger as the coverage increases. At 0.5 monolayer (ML), $`W`$ is larger than L by 5% reproducing the previously reported resultDijken . It is noteworthy that we reproduce the experimental resultDijken  without taking any transient mobility into account. The inset of Fig. 1 shows that for constant coverage of 0.5 ML, the aspect ratio, $`W/L`$, increases, as does the deposition angle. This result demonstrates that the asymmetry of the island shape is closely related to the reduced symmetry of the deposition geometry.

The top figure in Fig. 2(a) reproduces the flux distribution along the $`x`$-axis through the center of a pre-existing island as reported in Ref. 3, which herein is called the front and rear edge(FR-edge) steering effect. Also found in the present study, as one may see in the following figures of Fig. 2(a), is the strong dependence of the flux distribution on the island size; as the island size becomes small, the portion of the island with enhanced flux is enlarged. Besides, in Fig. 2(b), an additional steering effect called herein the side edge(S-edge) steering is observed; the trajectories of incident atoms are curved toward both side edges and enhanced flux is found along both side edges. As will be discussed later, both the dependence of deposition flux on the island size and the relative strength of FR-edge steering to S-edge steering have great impact on shaping the island.

To have a comprehensive picture on the asymmetric island growth, a series of simulation is performed with various substrate temperatures. Fig. 3 shows the dependence of $`L`$ and $`W`$ on coverage at various growth temperatures and the inset summarizes the temperature dependence of the aspect ratio at 0.5 ML. A stunning new phenomenon, the reversal of the asymmetry of island around 230K, is found; $`W>L`$ above 230K as observed in Ref. 3, but $`W<L`$ below 230K. Such a reversal of the island shape has never been predicted before.

We find that aforementioned, apparently diverse island growth can be understood by exploring the nature of the inhomogeneous deposition flux. The inhomogeneous flux is expected to manifest itself in two different ways; 1) asymmetric downhill current from the top of the island and 2) asymmetric deposition flux around the lower terrace of the island. As for the first one, despite the inhomogeneous deposition flux on the upper terrace of the island as shown in Fig. 2, the ratio of the downhill current per unit edge length over $`x`$-directional edge to that over $`y`$-directional edges is very close to 1.0 when we count the atoms coming over the ES barrier during growth simulation. This direction-independence of the downhill diffusion current occurs because the terrace diffusion on top of the island is more frequent by $`10^2`$ times than the downhill diffusion over ES barrier, and effectively homogenizes the inhomogeneous deposition flux on the island before atoms come down over the barrier.

As for the second one, we investigate the flux distribution on the lower terrace of the island. In Fig. 4(a), shown is the average deposition flux on terraces within 3 $`a_0`$ from the side edges (S-terrace) and that from the front and rear edges (FR-terrace), respectively, as a function of the island size. The flux on the FR-terrace decreases monotonically, but that on the S-terrace shows a maximum as the island becomes larger. For the island with its edge length smaller than 10 $`a_0`$, the flux on FR-terrace is larger than that on S-terrace, but the other way is true for the larger island.

Aforementioned dependence of the deposition flux on the island size can be understood from the flux distribution shown in Figs. 2 and 5. For small islands(Fig. 2(a) middle and bottom), the perturbed potentials at front and rear edges interfere and result in the increased flux on the FR-terrace. As the island size increases, the flux near the rear edge becomes relatively small (Fig. 2(a) top) due to the uninterfered rear edge potential. Hence, the average flux on the FR-terrace monotonically decreases as the island size increases as found in Fig. 4(a).

On the other hand, the average flux on the S-terrace is influenced by two counter-acting factors; firstly, the flux on the S-terrace increases with the island size, because the longer the edges the more atoms moving parallel to the side edge are attracted to the edge. However, the flux on S-terrace is not homogeneous along the side edge and another intricate variable is involved (Fig. 5). For small islands (Fig. 5 top and middle), the deposition flux is enhanced over the whole edge. As the island size increases, the relative portion of the side edge with pronounced flux enhancement gradually decreases(Fig. 5 bottom). This suggests that the average flux over the side edge will gradually decrease as the edge length increases. Hence, the competition between these two factors results in the maximum of the average flux at the island size of 17-20 $`a_0`$, explaining the dependence of the average flux on the island size observed in Fig. 4(a). The origin of the inhomogeneous enhancement of flux along the side edges on S-terraces can be found by the comparison of Figs. 2(a) and 5. The flux profile along the side edge(Fig. 5) shows quite a resemblance with that along the deposition direction through the center of the island (Fig. 2(a)). It suggests that the side edge flux is determined not solely by the S-edge steering but cooperatively by the FR-edge steering.

The island-size dependence of the flux ratio as shown in Fig. 4(b) is quite similar to the temperature dependence of the aspect ratio of the islands as shown in the inset of Fig. 3. Keeping in mind that the mean island size monotonically increases with the substrate temperature for the same coverage, it suggests that the temperature dependence of the aspect ratio of the island may originate from the dependence of the flux ratio on the island size; At low substrate temperature where small islands form, the average flux on FR-terrace is larger than that on S-terrace and results in asymmetric island having $`L>W`$, and vice versa. We examine this idea further by plotting the aspect ratio of the islands as a function of the mean size of islands formed with varying growth conditions. As shown in Fig. 4(b), even though the data are collected from uncorrelated simulations, there exists a strong correlation between the aspect ratio and the island size. Furthermore, it reproduces all the features found in the dependence of the flux ratio on the island size, especially the crossover (aspect ratio$`=`$1), maximum, and thereafter monotonic decrease of the ratio. (The island sizes where the crossover and the maximum in the aspect ratio occur are larger than the corresponding ones in the flux ratio. This is understood from the fact that the island size is determined by the integration of deposition flux.) It suggests that the most important factor determining the shape of the island is the deposition flux on the lower terrace of the island, which in turn is determined by the size of the island. If that conjecture is correct, the asymmetric island shape ($`L>W`$) should also be observed in the initial deposition stage and so does the reversal of the asymmetry with further growth unless the temperature is too low to form islands larger than the crossover size, because the island size is always small at the early stage of the island growth, irrespective of the substrate temperature. As indicated with $`+`$ signs in Fig. 3, the reversal of the asymmetry is indeed observed at similar island sizes, regardless of the substrate temperatures above 230K, reassuring that the island-size dependent flux ratio at the lower terrace of the island is the critical factor for shaping the island.

In summary, the present simulation combining MD and KMC properly reproduces the previously reported rectangular island growth of Cu/Cu(001) under deposition at a grazing incidence angle, and shows that the transient mobility is not a necessary condition for the asymmetric island growth. Instead, the asymmetry of the island shape is attributed mainly to the asymmetry of the deposition flux on the lower terrace of islands. Also found is the existence of the reversal of the aspect ratio of asymmetric islands, depending on the substrate temperature or the island size. This, so far, unanticipated phenomenon remains to be proven by further experimental study. It is also found important that, for a proper interpretation of the experimental results, the steering effect should be treated in full three dimension.

###### Acknowledgements.

This work is support by KISTEP(00-B-WB-07-A-01) (JK).