text:
# Quantum size correction to the work function and the centroid of excess charge in positively ionized simple metal clusters

## I Introduction

In the SJM, the total energy of a cluster with jellium radius $`R=\overline{r}_s^{1/3}N^{}`$ is given by

$`E_R[n]=T_s[n]+E_{xc}[n]`$ $`+`$ $`{\displaystyle 𝑑𝐫\{\frac{1}{2}\varphi _R[n;𝐫]+\delta v_{\mathrm{WS}}\mathrm{\Theta }(𝐫)\}[n(𝐫)n_+(𝐫)]}`$ (1)
$`+`$ $`(\epsilon _\mathrm{M}+\overline{w}_R){\displaystyle 𝑑𝐫n_+(𝐫)}.`$ (2)
The first and second terms in the right hand side of Eq.(1) are the non-interacting kinetic and the exchange-correlation energies, respectively and

$$\varphi _R([n];𝐫)=𝑑𝐫^{}\frac{n(𝐫^{})n_+(𝐫^{})}{\left|𝐫𝐫^{}\right|}.$$
(3)
All equations throughout this paper are in atomic units ( $`\mathrm{}=e^2=m=a_0=1`$) unless otherwise explicitly expressed. $`n_+(𝐫)`$ is the positive jellium background density, $`n_+(𝐫)=\overline{n}\theta (Rr)`$ where $`\overline{n}=3/4\pi \overline{r}_s^3`$ is the bulk density. The quantity $`\delta v_{\mathrm{WS}}`$ is the average of the difference potential over the Wigner-Seitz cell and the difference potential, $`\delta v`$, is defined as the difference between the pseudopotential of a lattice of ions and the electrostatic potential of the jellium positive background. $`\mathrm{\Theta }(𝐫)`$ takes the value of unity inside the jellium background and zero, outside. $`\epsilon _\mathrm{M}`$ and $`\overline{w}_R`$ are the Madelung energy and the repulsive part of the pseudo-potential, respectively.

In the continuum approximation, the ground-state density, $`n_R(𝐫)`$, of the neutral cluster is given by the solution of the Euler equation

$$\mu (R)=\frac{\delta E_R[n]}{\delta n(𝐫)}|_{n=n_R}=\varphi _R([n_R];𝐫)+\delta v_{\mathrm{WS}}\mathrm{\Theta }(𝐫)+\frac{\delta E_{kxc}[n]}{\delta n(𝐫)}|_{n=n_R}.$$
(4)
In Eq (4), $`\mu (R)`$ is the chemical potential of electrons which is constant throughout the cluster for the exact ground-state density $`n_R`$, and $`E_{kxc}=T_s+E_{xc}`$. The ground-state density satisfies the costraint $`𝑑𝐫n_R(𝐫)=N^{}`$. Removing an electron from the neutral cluster gives rise to a new ground-state density $`n_R^{}(𝐫)=n_R(𝐫)+\delta n(𝐫)`$ with $`𝑑𝐫\delta n(𝐫)=1`$.

The ionization energy is defined as the difference in the ground-state energies of the two systems with $`N^{}`$ and $`(N^{}1)`$ electrons which may be represented by the expansion

$`I(R)`$ $`=`$ $`E_R[n_R^{}]E_R[n_R]`$ (5)
$`=`$ $`{\displaystyle 𝑑𝐫\frac{\delta E_R[n]}{\delta n(𝐫)}}|_{n=n_R}\delta n(𝐫)+{\displaystyle \frac{1}{2}}{\displaystyle 𝑑𝐫𝑑𝐫^{}\frac{\delta ^2E_R[n]}{\delta n(𝐫)\delta n(𝐫^{})}}|_{n=n_R}\delta n(𝐫)\delta n(𝐫^{})+\mathrm{}`$ (6)
The first term in the right hand side of Eq. (5) is $`\mu (R)`$ and in local density approximation for the kinetic and exchange-correlation energies we have

$$\frac{\delta ^2E_R[n]}{\delta n(𝐫)\delta n(𝐫^{})}|_{n=n_R}=\frac{1}{|𝐫𝐫^{}|}+\delta (𝐫𝐫^{})\frac{^2}{n^2}\{n\epsilon _{kxc}(n)\}|_{n=n_R(𝐫)}.$$
(7)
Here, $`\epsilon _{kxc}(n)`$ is the sum of the bulk kinetic and exchange-correlation energies per electron at the density $`n`$. Choosing a spherical distribution for $`\delta n(𝐫)`$

$$\delta n(𝐫)=\frac{\xi (r)}{4\pi (R+a)^2},$$
(8)
and interpretting $`(R+a)`$ as the centroid of the excess charge, $`\xi (r)`$ satisfies the following two equations:

$$_0^{\mathrm{}}\frac{\xi (r)}{4\pi (R+a)^2}4\pi r^2𝑑r=1$$
(9)
$$_0^{\mathrm{}}r\frac{\xi (r)}{4\pi (R+a)^2}4\pi r^2𝑑r=(R+a).$$
(10)
The simplest choice for $`\xi (r)`$ which stisfies the above two constraints, is

$$\xi (r)=\delta (rRa).$$
(11)
Using the Eqs. (7)-(11) in the second term of the right hand side of the Eq. (5) one obtains

$$I(R)=\mu (R)+\frac{1}{2(R+a)}+O(R^2).$$
(12)
Since $`\mu (R)`$ is a constant, independent of $`𝐫`$, for the exact ground-state density of the neutral cluster, $`n_R(𝐫)`$, it is possible to evaluate it at a point deep inside the cluster where $`n_R(𝐫)=\overline{n}`$. The density profile $`n_R(r)`$ and the potential $`\varphi _R([n_R],r)`$ for large clusters can be expressed as corresponding quantities for the planar surface (i.e. $`R=\mathrm{}`$) plus a size correction

$$n_R(x)=n(x)+\frac{f(x)}{R}+O(R^2)$$
(13)
$$\varphi _R(x)=\varphi (x)+\frac{h(x)}{R}+O(R^2)$$
(14)
where $`x=rR`$ gives the distance from the jellium edge. At large distances from the jellium edge outside the cluster, we have $`f(+\mathrm{})=\varphi (+\mathrm{})=h(+\mathrm{})=0`$, and for a point deep inside the cluster ($`x=\mathrm{}`$) we have $`n(\mathrm{})=\overline{n}`$, $`f(\mathrm{})=0`$. Hence Eq. (4) reads

$$\mu (R)=\varphi (\mathrm{})+\frac{h(\mathrm{})}{R}+\delta v_{\mathrm{WS}}+\frac{}{n}\{n\epsilon _{kxc}(n)\}|_{n=\overline{n}}.$$
(15)
Upon inserting Eq. (15) into Eq. (12) and using the definition of the work function of a planar surface in the SJM

$$W=\mathrm{\Delta }\varphi \delta v_{\mathrm{WS}}\frac{}{n}\{n\epsilon _{kxc}(n)\}|_{n=\overline{n}},$$
(16)
where $`\mathrm{\Delta }\varphi =\varphi (+\mathrm{})\varphi (\mathrm{})`$ is the electrostatic dipole barrier, one obtains

$$I(R)=\left(W+\frac{c}{R}\right)+\frac{1}{2(R+a)}+O(R^2),$$
(17)
where $`c=h(\mathrm{})`$.

Equation (17) is also valid in the JM for which $`W`$ is obtained by putting $`\delta v_{\mathrm{WS}}=0`$ in Eq. (16) and recalculating $`\mathrm{\Delta }\varphi `$ with the JM density profile. There are two methods for calculating $`W`$ and $`c`$ called the ”Koopmans’ method” and the ”change in self-consistent field ($`\mathrm{\Delta }SCF`$) method” which give the same results for the correct $`n(x)`$, but the $`\mathrm{\Delta }SCF`$ method has the advantage that its results are less sensitive to the exact density profile $`n(x)`$. However, in our calculations, all the quantities $`W`$, $`c`$, and $`a`$ are obtained from fitting of our self-consistent ionization energies of the clusters to Eq. (17).

## II Calculational Scheme and Results

In this work, using the JM as well as the SJM with local spin-density approximation, we have solved the self-consistent Kohn-Sham equations and obtained the energies of neutral and singly-ionized $`N`$-atomic clusters of different sizes ($`2N100`$). The calculations have been performed for different $`\overline{r}_s`$ values of the jellium sphere ($`2\overline{r}_s7`$) in steps 0.5 bohr. For each $`\overline{r}_s`$ value, we have calculated the ionization energies

$$I(R)=E_R(N1)E_R(N),$$
(18)
for $`2N100`$ with $`R=N^{1/3}\overline{r}_s`$. Then, by fitting the $`I(R)`$ values to Eq. (12) we have obtained the values of $`W`$, $`c`$, and $`a`$ for that $`\overline{r}_s`$ value.

Figure 1(a) shows the work function obtained for our SJM data and is labeled as SJM(1) because, in this paper we have called the Eq. (17) as scheme 1. To compare our results with the bulk calculation results , SJM(bulk), we have also plotted the data from table IV of Ref. REFERENCES. As is seen, our results show a good agreement with the bulk calculations results.

In Fig. 1(b) we have compared our JM(1) results with the bulk JM results of Ref. REFERENCES. Here also the good agreement is obvious.

Figure 1(c) compares the JM(1) and the SJM(1) results. These two curves coincide at $`\overline{r}_s4.18`$ where the JM is mechanically stable. The difference $`W_{\mathrm{SJM}}W_{\mathrm{JM}}`$ is positive for $`\overline{r}_s<4.18`$ and changes sign for $`\overline{r}_s>4.18`$. This behavior can easily be predicted if one resorts to the corresponding expressions for the work functions: For $`\overline{r}_s4.18`$, we have $`\delta v_{\mathrm{WS}}0`$ and therefore, the density profiles are the same which results in $`\mathrm{\Delta }\varphi [n_{\mathrm{SJM}}]=\mathrm{\Delta }\varphi [n_{\mathrm{JM}}]`$ and in turn, by Eq. (16) leads to the same workfunctions. However, since the term $`\delta v_{\mathrm{WS}}`$ in Eq. (16) changes sign at $`\overline{r}_s=4.18`$, it gives rise to the change in the sign of $`W_{\mathrm{SJM}}W_{\mathrm{JM}}`$.

In Fig. 1(d), the values of the quantum size correction $`c`$ for the JM(1) and the SJM(1) are compared. These two plots have two intersections. The one at $`\overline{r}_s4.18`$ is obvious if we consider the dependence of the quantity $`c`$ on the density profile $`n(x)`$

$$c=h(\mathrm{})=4\pi _{\mathrm{}}^+\mathrm{}𝑑x\{x^2[n(x)\overline{n}\theta (x)]+xf(x)\},$$
(19)
$`{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑xf(x)`$ $`=`$ $`2{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑xx[n(x)\overline{n}\theta (x)]`$ (20)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}\mathrm{\Delta }\varphi .`$ (21)
The intersection at $`\overline{r}_s2.5`$ is because of the improper behavior of the properties at high densities in the JM ( See for example Fig. 6(a) of Ref. REFERENCES for the surface energy which becomes negative for $`\overline{r}_s<2.5`$).

In Fig. 2(a), we have shown the position of the centroid of excess charge, $`a`$, for the singly-ionized clusters as a function of $`\overline{r}_s`$ both for the JM(1) and the SJM(1). As is seen, both plots predict negative values which means that the centroid for a positively charged cluster lies inside the jellium \[ the jellium surface is taken as the origin, see Eq. (8)\]. This is consistent with our understandings that the excess charge in a metal resides on the surface. On the other hand, the JM(1) and the SJM(1) have opposite behaviors: the former is decreasing ( increasing in absolute value) and the latter is increasing ( decreasing in absolute value) and intersect each other at $`\overline{r}_s4.18`$. The explanation of these opposite behaviors is straightforward: In simple JM, if we take the thickness of the jellium shell at the edge as $`\lambda `$, then the volume of this shell contains a unit charge (for singly-ionized clusters)

$$\lambda (4\pi R^2)\overline{n}=1,$$
(22)
which gives

$$a\frac{1}{2}\lambda =\frac{1}{6R^2}\overline{r}_s^3.$$
(23)
By Eq. (23), for a fixed $`R`$, in JM as $`r_s`$ increases, the position of the centroid moves inside from the edge towards the center of the jellium sphere. However, in the SJM, there exist two competing effects: On the one hand, as in the simple JM, because the jellium background is the same, increasing $`\overline{r}_s`$ has the tendency to increase $`\lambda `$ but, on the other hand, the SJM Kohn-Sham effective potential contains an extra term $`\delta v_{\mathrm{WS}}`$ which is an increasing function of $`\overline{r}_s`$ ( See Table IV of Ref. REFERENCES ). This term is negative for $`\overline{r}_s<4.18`$ and hence deepens the effective potential ( compared to the JM ) which reduces the spill-out of the electrons, so that $`a_{\mathrm{SJM}}<a_{\mathrm{JM}}`$. On the other hand, for $`\overline{r}_s>4.18`$ the extra term is positive which shallowen the effective potential ( compared to the JM ) which in turn increases the range of the Kohn-sham orbitals from which the electron density is calculated. The result is that in the SJM, the overlap of the negative and positive charges increases (the shell thickness decreases) with increasing $`\overline{r}_s`$, so that $`a_{\mathrm{SJM}}>a_{\mathrm{JM}}`$.

The quantity $`a`$ calculated above is, of course, the centroid position for singly-ionized cluster. To calculate the centroid position for doubly-ionized clusters, we fitted the Kohn-Sham energies of doubly-ionized clusters to the Eq. (20) of Ref. REFERENCES which is obtained from Taylor expansion of the energy of a $`z`$-ply charged cluster about $`z=0`$. The resulted values for $`a`$ was +0.34, -0.18, -0.29 for $`\overline{r}_s`$ values 2, 4, 6, respectively. Obviousely, the result for $`\overline{r}_s=2`$ is incorrect because of its wrong sign. The predicted value for $`\overline{r}_s=6`$ is quite acceptable and consistent with the above arguments. One reason for the incorrect results may be in the use of Taylor expansion about $`z=0`$ for points $`z>>0`$, as in Eq. (20) of Ref. REFERENCES. A better approximation is obtained by using the $`z`$-th ionization energy

$`I_z(R)`$ $`=`$ $`E_R(N^{}z)E_R(N^{}z+1)`$ (24)
$`=`$ $`\epsilon _R^{\mathrm{HO}}(N^{}z+1)+{\displaystyle \frac{1}{2(R+a_z)}},`$ (25)
where we have used the Taylor expansion of the ground-state energy $`E_R(N^{}z)`$ at $`N^{}z+1`$. The quantity $`\epsilon _R^{\mathrm{HO}}(N^{}z+1)`$ is the highest occupied Kohn-Sham orbital energy for the $`(z1)`$-ply ionized cluster with radius $`R=\overline{r}_s(N^{})^{1/3}`$.

It is instructive to compare the behaviors of the position of the centroid of excess charge with the position of the image plane of a planar surface. In Fig. 2(b) we have compared these quantities in the JM and the SJM. As is seen, both of them are positive and intersect at $`\overline{r}_s4.`$ In the JM plot we have used the Eq. (8) of Ref. REFERENCES and for the SJM we have used the data from Fig. 1 of Ref. REFERENCES. Here, although the argument used for the finite cluster to calculate $`\lambda `$ is useless but, the effect of $`\delta v`$ is the same as in the finite cluster case which leads to similar behaviors.

Now, we show the results of fitting our self-consistent Kohn-Sham energies to Eq. (17) but neglecting the quantum size correction, as in Eq. (16) of Ref. REFERENCES. Putting $`c=0`$ in Eq. (17) we call the resulting equation as scheme 2.

In Fig. 3(a) we have compared the work functions obtained in schemes 1 and 2. As is seen, the scheme 2 predicts lower values than the scheme 1 otherwise the shapes are the same.

Figure 3(b) shows the plots of $`a`$ for JM and SJM in scheme 2. As we see, both of them are positive and have increasing behaviors. This means that the position of the centroid of the excess positive charge (for the singly-ionized clusters) lies outside the jellium edge which is incorrect because there exists no positive charge. Therefore, inclusion of the quantum size effect is vital for predicting correct values of the physical quantities.

Finally, we introduce the scheme 3 which is obtained by putting $`a=0`$ in Eq. (17) as Eqs. (1) of Refs. REFERENCES and REFERENCES. In Fig. 4(a) we have compared the SJM work functions in the schemes 1 and 3 with the bulk values. It is seen that scheme 3 improves the workfunction and reduces the distance from the bulk values at higher densities.

The SJM fitted values of the size correction for schemes 1 and 3 are compared in Fig. 4(b). It is seen that in scheme 3, the values are nearly constant for the metallic densities $`c3.7eV`$ in good agreement with the result of Ref. REFERENCES.

## III summary and conclusions

We have calculated the self-consistent Kohn-Sham energies for neutral and singly-ionized clusters with different sizes and repeated the calculations for different $`\overline{r}_s`$ values in the metallic range. The ionization energies for each $`\overline{r}_s`$ has been separately fitted to Eq. (17), we called scheme 1. The results obtained for the position of the centroid of excess charge is consistent with the understanding that excess charges reside at the outer surface of a metal particle. However, for small doubly-ionized clusters, using the Taylor expansion about $`z=0`$ gives incorrect $`a`$ values for clusters of high electron-density metals. A better result would be obtained by using a Taylor expansion of the energy of $`z`$-ply ionized cluster about the point $`(z1)`$.

Next, in scheme 2, we have fitted the same ionization energies to Eq. (17) taking $`c=0`$. The results show incorrect values for the centroid position.

Finally, in scheme 3, we put $`a=0`$ in Eq. (17) and fit the ionization energies, which results in better values for the work function at higher densities.

To conclude, if one wishes to calculate the position of the centroid of excess charge of $`z`$-ply charged cluster, one should use the Taylor expansion about point $`(z1)`$ as in Eq. (25). On the other hand, if one wishes choose one of the versions of Eq. (17) to calculate the ionization energies ( or electron affinities) of simple metal clusters, the best choice would be the Eq. (17) with the $`\mathrm{\Delta }SCF`$ values for $`W`$ and $`c`$, and for sufficiently large clusters ($`R>>a`$) the scheme 3 with the $`\mathrm{\Delta }SCF`$ values for $`W`$ and $`c`$ gives the same results as scheme 2.

ACKNOWLEDGEMENT

The author would like to thank Bahram Payami for providing computer facilities.


format: text

text:
# Lensing effects of misaligned disks in dark matter halos

## 1 Introduction

Lensing by individual galaxies provides a plethora of observational signatures ranging from multiply imaged, highly magnified background sources to a weakly sheared background (see Schneider, Ehlers & Falco 1992, Blandford & Narayan 1992). Several authors have attempted to use these observed image geometries, positions and magnitudes to infer and reconstruct the mass profile of the lens using ground-based optical and radio surveys (Browne et al. 2001 \[JVAS/CLASS\] and references therein; Lehar et al. 2000 \[CASTLES\]). Attempts have been made to constrain both the density profiles of the mass distribution within the inner regions of galaxies (Keeton & Madau 2001; Cohn et al. 2001) as well as obtain constraints on the Hubble parameter $`H_0`$ (Kundić et al. 1997; Fassnacht et al. 2002). Both of these enterprises are complicated by the presence and precise configuration of structure within and exterior to the lensing galaxy in question. The degeneracies in such modeling attempts arise due to the presence of structural complexity in the inner regions of galaxies, the prevalence of disky components, bulges, central black holes (Möller & Blain 1998; Mao, Witt & Koopmans 2001), as well as the presence of external mass perturbations: nearby galaxies, groups or clusters in the vicinity of the lens (Keeton, Kochanek & Seljak 1997; Möller, Natarajan, Kneib & Blain 2002). In this work, we explore a further source of uncertainty - the presence of misaligned disks within dark halos.

The motivation for studying galaxy scale systems that have their baryonic disks misaligned slightly with their dark matter halos is two-fold: first, given what is currently accepted regarding the heirarchical assembly of structure in cold dark matter dominated structure formation models (Frenk & White 1991) i.e. that merging of galaxies is frequent - then dynamically induced small misalignments between the disk and halo are expected during the merger process. Since the lensing properties of a system (an individual galaxy in this case) depend on the projected surface mass density, we demonstrate in the following sections that these small misalignments are likely to have very interesting and important observable consequences. Secondly, it has been recently claimed in the literature (Schneider & Mao 1998; Metcalf & Madau 2000; Dalal & Kochanek 2002) that the observed magnification ratios of several if not most multiply imaged systems (notably B1422+231 and PG1115+08) can only be explained as arising due to the presence of small scale substructure within galaxy halos. Here, we illustrate that there are in fact other physical scenarios such as an anisotropic mass distribution within the lensing galaxy: non-zero quadrupole and octopole components induced in the lensing potential due to a misaligned disky component can also give rise to magnification ratios and image characteristics that differ from those predicted by the standard isothermal/pseudo isothermal mass models currently in use. While the evidence for substructure is compelling it is useful to explore the degeneracy with primary lens properties.

## 2 When are halos and disks misaligned?

In the standard paradigm for structure formation in cold dark matter dominated models, galaxies are assembled over time via the merger of smaller sub-clumps (Frenk & White 1991). The merger rate for halos of a given mass can be computed semi-analytically using the extended Press-Schechter formalism (see Cole et al. 1995 and subsequent refinements by Sheth, Mo & Tormen 2001), and can be directly compared with cosmological N-body simulations (Raig et al. 2001). When dark matter halos merge their baryonic components also merge; detailed numerical simulations of this process have been carried out by Mihos & Hernquist (1994) and Barnes et al. (1998). In these simulations it is found that the baryonic gas flows to the center of the more massive halo and the inner regions of the halo become baryon-dominated. The tidal torquing during the merger process is significant and can easily cause the embedded disk and halo to be misaligned. Misalignments of less than an arcsecond are highly probable. It is also likely that during the ‘relaxation’ process the disk wanders transiently within the halo (Nelson & Tremaine 1995; Dubinski & Kuijken 1995). Finally, interactions that do not result in mergers can also perturb the disk-halo configuration and induce transient misalignments. These are expected to be damped well within a Hubble time (Dubinski & Kuijken 1995). Since massive galaxies are assembled via mergers, and because such galaxies are particularly efficient as lenses, it is likely that even transitory misalignments could be important for a number of lens systems.

### 2.1 Estimating the number of misaligned disks in halos

Lens galaxies tend to sample the high mass end of the galaxy mass function. As claimed above, such systems are expected to be built up by mergers. Misalignments are highly probable during this stage and are likely to persist post-merger. We use the observationally determined galaxy merger rate to estimate the expected frequency of misaligned systems. Patton et al. (2002) report a merger rate of $`\frac{dN}{dz}(1+z)^{2.3\pm 0.7}`$ from a study of a sample of bright galaxies, $`18<M_B<21`$, in the CNOC2 redshift survey. We assume here that all mergers lead to misalignments due to tidal torquing during the merger process and that these misalignments persist for the duration of the merger. From detailed N-body simulations of interacting galaxies, the merger time-scale is found to be of order 1 Gyr (see Fig. 4 of Mihos & Hernquist 1994). Since typical lenses lie within a redshift range $`z=0.10.5`$, the misaligned fraction is roughly,

$`f_{\mathrm{misalign}}({\displaystyle \frac{t_{\mathrm{merger}}}{\delta t}}){\displaystyle \frac{V(\mathit{0.1},\mathit{0.5})}{V(\mathit{0},\mathit{5})}}{\displaystyle \frac{_{z=0.1}^{z=0.5}(\frac{dN}{dz})𝑑z}{_{z=0.1}^{z=0.5}𝑑N}}10\%,`$ (1)
where $`V(z_\mathit{1},z_\mathit{2})`$ is the comoving volume between redshifts $`z_1`$ and $`z_2`$; $`t_{\mathrm{merger}}1\mathrm{Gyr}`$ is the merger timescale, and $`\delta t`$ the total elapsed time from $`z=0.1\mathrm{\hspace{0.17em}0.5}`$. With the assumptions given above, $`f_{\mathrm{misalign}}`$ is found to be about 10% in a $`h=0.6`$, $`\mathrm{\Omega }_\lambda =0.7`$, $`\mathrm{\Omega }_\mathrm{m}=0.3`$ cosmology. Note that the merger rate found by Patton et al. (2002) is consistent with semi-analytic estimates from models that are calibrated to observational data, as shown in the top panels of Fig. 6 in Somerville et al. (2000).

At present, there are approximately 50 confirmed cases of individual galaxies lensing high redshift background quasars or galaxies. We expect that at least a few of these are likely to contain misaligned disks. Note that while the majority of lenses detected thus far are typically early-type galaxies, observationally it is becoming clear that early-type galaxies also do often contain embedded compact disky components (Rest et al. 2001). In fact, we demonstrate below that even a highly compact misaligned disk that lies interior to the Einstein radius can produce the anomalous image geometries. Therefore, our analysis is generic and applicable to all lensing galaxies that contain a disky mass component.

## 3 Methodology: The Ray-tracing method

In order to study the strong lensing properties of a misaligned halo and disk we used the multiple plane ray-tracing code developed by Möller & Blain (1998). Computing the lensing properties of multi-component lens systems is complicated and can be very time-consuming if this is not done with an optimized algorithm. The code we employ here is based on the ‘grid search method’ (Schneider, Ehlers & Flaco 1992) and uses an adaptive grid that is optimized to calculate deflection angles for parametric multiple component lens models quickly and accurately. A regular grid of triangle pairs on the image plane is mapped to the source plane. The magnification at each point on the source plane is calculated as the sum of the ratios of the area of the unmapped triangles on the image plane to that of the mapped triangles on the source plane. Similarly, shears are calculated from the distortions of mapped (source plane) and unmapped (image plane) triangles. In order to reduce the computation time the mapping process is iterated adaptively. The numerical inaccuracies on magnifications and shears are very small, usually less than 0.1%. The numerical uncertainty when calculating the number of images for a given point source is greater; faint images near a caustic line may be missed by this method (Schneider, Ehlers & Falco 1992). However, increasing the resolution of the grid enough can overcome this problem, and it does not affect any of the results discussed in this paper.

The illustrative single lens systems studied here consists of a dark halo and an embedded disk. The halo model is a cored isothermal sphere (CIS). We have have also investigated pseudo-isothermal elliptical mass distributions (PIEMDs) with an ellipticity of $`ϵ=0.3`$ (Kassiola & Kovner 1993) and find that finite ellipticity in the halo does not significantly alter our results. The two parameters required to specify the CIS are the core-radius and the mass within a radius $`R`$. The former is not well-determined by observations of real galaxies. In order to test the effects of our choice of $`r_\mathrm{c}`$ on the results we explore a wide range of values, including $`r_\mathrm{c}=0`$. The total halo mass contained within a 20 $`\mathrm{kpc}`$ radius was held constant for all the models considered here to enable a sensible comparison. The primary effect of varying $`r_\mathrm{c}`$ is a change in the relative significance of the halo in the lensing system, which affects the size of the elliptical caustic (see § 4.1). This is due to the fact that a small core radius implies that the convergence $`\kappa `$ in the central region is very high. In this case, the mass profile is centrally concentrated and thus the cross section for strong lensing is significantly increased.

## 4 Lens properties of simulated spiral systems

We take the Milky Way as the template for our spiral lens, and hence model the halo as a CIS with a small core radius $`r_\mathrm{c}=0.5\mathrm{kpc}`$ and mass $`M_\mathrm{h}=3.64\times 10^{11}M_{\mathrm{}}`$ within a radius of 20 $`\mathrm{kpc}`$. This corresponds roughly to a velocity dispersion of 160 km/s and the Einstein radius of the halo is 0.4”. The disk has mass $`M_\mathrm{d}=6.4\times 10^{10}M_{\mathrm{}}`$ and scale length $`r_\mathrm{d}=3.2\mathrm{kpc}`$. Note that models with less massive and more compact disks are also studied and the results are presented in a subsequent section. An inclination angle of 70 degrees was assumed for the disk which is roughly the minimal inclination needed to significantly increase the cross-section for multiple imaging (see Fig. 2 in Möller & Blain 1998). The lens is located at $`z=0.2`$ and the source plane at $`z=1.5`$. The values assumed for cosmological parameters are $`H_0=60\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, $`\mathrm{\Omega }_\mathrm{m}=0.3`$, and $`\mathrm{\Omega }_\lambda =0.7`$. We have also experimented with alternate cosmologies, such as $`\mathrm{\Omega }_\mathrm{m}=1.0,\mathrm{\Omega }_\lambda =0.0`$ (EdS) and $`\mathrm{\Omega }_\mathrm{m}=0.3,\mathrm{\Omega }_\lambda =0.0`$ (OCDM), and find that our results are not significantly affected by the particular choice of these parameters.

In this paper we study simple spatial offsets of the disk relative to the halo. We assume that the disk retains an exponential profile and that the halo retains an isothermal profile. Offsets in the plane of the disk are denoted by $`\mathrm{\Delta }x`$ and offsets outside the plane of the disk are denoted $`\mathrm{\Delta }y`$. The following sections assume offsets of $`\mathrm{\Delta }x=0.6^{\prime \prime }`$, $`\mathrm{\Delta }y=0.6^{\prime \prime }`$, and $`\mathrm{\Delta }x=\mathrm{\Delta }y=0.6^{\prime \prime }`$. We consider this to be a reasonable approximation to the disturbed and asymmetric mass distribution that will result from mergers or strong interactions.

### 4.1 Magnification maps, cross sections and image geometries

In this section we present magnification maps in the source plane for several lens configurations, varying the position of the disk relative to the halo in each case, i.e. the extent of spatial misalignment. We also generate the corresponding magnification maps in the image plane, where the value of each pixel is the magnification of the observed image at that location. The caustics (points of infinite magnification) in the source plane map onto critical curves in the image plane. Magnification cross sections and representative image geometries are also presented.

Fig. 1 shows the magnification maps in the source plane for various representative values of the disk offset. It is evident that a realistic shift in the disk position relative to the halo does indeed distort the caustic shapes severely. A shift in the x-direction, or in the plane of the disk, yields a shift of the diamond-shaped curves that form the *tangential* caustic, which is caused by the presence of the disk. The elliptical (*radial*) caustic, which is due to the peak in the central mass concentration of the combined halo and disk, reduces in size and is distorted in the direction of the offset. Offsets outside the plane of the disk lead to more complicated caustic geometries. The tangential caustic shifts and the upper cusp of the caustic appears to fold over into the diamond itself. The radial caustic changes its shape drastically, and is pinched inward near the lower cusp. But a comparison of Panel (c) and Panel (d) of Fig. 1 shows that an offset in the x-direction can partially reduce the effect of a shift in the y-direction.

Fig. 2 shows the magnification maps in the image plane. The small inner critical curve corresponds to the radial caustic and the outer one corresponds to the tangential, diamond-shaped caustic. It is clear from these images that the locations of highly magnified images can be substantially affected by misalignments. As with the magnification maps in the source plane, offsets out of the plane of the disk also create complicated geometries. In Panel (c), Fig. 2, the upper, roughly triangular shaped critical curve has split from the inner critical curve. In Panel (d), Fig. 2, the region within this new critical curve has nearly collapsed. If the disk were shifted further in the x-direction then this new critical curve would be absorbed into the outer critical curve. The evolution of caustics and critical curves is continuous, i.e. the curves gradually distort, pinch off and absorb each other.

To investigate the effects of misalignment on total magnification, we calculate the magnification cross section ratio of the misaligned cases to the aligned case. The magnification cross section for a lens is derived by adding up the total number of pixels with a value greater than some threshold value $`A`$ in a magnification map. The results are shown in Fig. 3. An offset in the plane of the disk is seen to have little effect. But the cross section for high magnification $`(A>\mathrm{\hspace{0.17em}30})`$ events is significantly increased for offsets outside the plane of the disk. This is a manifestation of the large surface area very near the caustic curves in Panel (c) of Fig. 2. The offset in the y-direction of 0.6” corresponds approximately to the maximum possible enhancement in high magnification cross section for our lens model. The observed trend that offsets in the x-direction tend to compensate for the more dramatic effects of shifts in the y-direction is apparent here also.

This effect on the cross section for high magnifications is likely to have a measurable effect on the statistics of multiple image systems. Due to the preferential selection of high luminosity images in a survey, any flux limited sample of lens systems is biased towards systems with a large cross section for high magnifications (magnification bias, e.g. Borgeest, 1991). Therefore the fraction of lens systems with misaligned disks as estimated in § 2.1 is conservative, and the actual number of systems in which a misaligned disk is important might be larger.

The number of images of a source will increase or decrease by two if the source crosses a caustic (see Schneider, Ehlers & Falco 1992 for a mathematical illustration of the same). An aligned halo and disk are capable of producing at most five separate images, in the configuration of a four-image, possibly asymmetric, Einstein cross with a demagnified central image. However, our distorted caustic geometries for the misaligned cases have regions where even higher image numbers are possible. For instance, the lens configuration of Panel (c), Fig. 1 has a region where seven images are formed and other lens configurations can produce even more images. Such a configuration with a large multiplicity of images has been observed; B1359+154 is a 6-image system (Rusin et al. 2001) where the lens is a compact group of galaxies at $`z=1`$. These asymmetries in caustic geometries are expected to yield complicated, asymmetric image geometries. Fig. 4 illustrates examples of the kind of image geometries that can be produced with the increased parameter space provided by disk offsets. The inset within each panel shows the magnification maps in the source plane and location of the point source. A representative range of image geometries with a wide range of magnification ratios is shown. Current observations of such systems might not reveal their complexity entirely; for example the configuration shown in Fig. 4, Panel (c) might be easily confused for an Einstein cross because the separate demagnified image is likely to go undetected.

## 5 Lens properties of early-type systems

In the previous section we have concentrated on misalignment between the disk and halo component in late-type lensing galaxies with Milky Way-type spiral disks. The majority of lens systems known are, in fact, early type galaxies, as these tend to be more massive and hence have a larger cross section for lensing. Recent observations show that disks or disk-like structures might well be present in early-type lenses (Kelson et al. 2000). This suggests that early-type galaxies may contain disks that may also be misaligned with respect to the bulge/halo. Whether there is a possible misalignment between a disk and the (baryonic) bulge depends on the exact cause of the misalignment, but a scenario as described in § 2 may cause a misalignment between disk and the bulge as well as with the dark matter halo.

### 5.1 Magnification maps, cross-sections and image geometries

We investigate here the lensing properties of a misaligned disk component that is similar to that observed in ellipticals. In particular, we choose disk parameters that are comparable to those found by fitting the light profile of ellipticals in CL 1358+62 with a bulge+disk model as in Kelson et al. (2000). For our early-type lens model, we choose a disk scale length of $`r_\mathrm{d}=0.5\mathrm{kpc}`$ which corresponds to the average disk scale length of the ellipitcals in CL 1358+62 at the cluster redshift of $`z=0.33`$. The disk mass is taken to be $`M_\mathrm{h}=2\times 10^{10}M_{\mathrm{}}`$. We model the bulge/halo as a CIS with $`r_\mathrm{c}=0.5\mathrm{kpc}`$ and increase the total mass within $`20\mathrm{kpc}`$ to $`M_\mathrm{h}=7\times 10^{11}M_{\mathrm{}}`$, to reflect the fact that early-type lens galaxies are usually more massive. This choice of parameters gives a velocity dispersion of 220 km/s and yields an Einstein radius of 1” for the bulge/halo.

Fig. 5 shows the magnification maps in the source plane for the elliptical lens. The values for disk offsets are the same as in the previous section. Although this disk is less massive than the disk in our late-type lens, it is also significantly more compact, meaning that the central $`\kappa `$ value is in fact high enough for the disk to produce an additional elliptical component to the overall caustic structure. So the magnification maps in Fig. 5 represent a superposition of two elliptical components due to the disk and the halo, as well as the tangential caustic due to the disk. Fig. 6 shows the corresponding magnification map in the image plane. As seen in Panel (a) of Figs. 5 and 6, our model of a compact disk has virtually no effect on the structure of the caustics or critical curves when the disk and halo are aligned. The circular critical curve in Fig. 6 is nearly a circular Einstein ring. The other panels of Fig. 6 show critical curves that are introduced by a misalignment of the compact disk.

Even though the disk components in the elliptical lens galaxy model contain a smaller fraction of the total mass than the disks in the late-type lens models discussed in previous sections, the effect of the misalignment between disk and halo is still very pronounced. In particular for misalignments in both the x and y direction, the effect on the structure of the caustics and critical curves is drastic. Since the majority of observed 4-image lens systems are expected to be produced by lensing of a point source that lies close to the high-magnification caustic lines, this indicates that the image properties of these systems are affected strongly by the presence of a misaligned disk, even if it is contributes only a small fratction of the total mass and has a scale length of less than a few hundred parsecs.

It is therefore likely that the image geometries will also change significantly. In Fig. 7 we show some typical image geometries of early-type lens systems with a misaligned disk for two different lens models. The boxes are centered on the image positions and have sizes that are proportional to the logarithm of the magnification. The numbers near the images give the magnifications of the image. Depending on the source position, both common and unusual image geometries may be produced. The bottom-left panel shows a fairly typical 3-1 geometry produced by a source near a cusp. In contrast, the top-left panel shows a case with non-standard magnification ratios; the central image is not the brightest image. These ‘inverted’ magnification ratios are a unique feature of complex multi-component lens models; as discussed below, single component mass models predict that the central image of the group of three images in 3-1 lens geometries is always the brightest. The top-right panel of Fig. 7 shows a lensed source near a cusp, but the 3-1 image geometry is not present. Instead there are four clustered images, with the two brightest ones nearly overlapping, and a fifth image lying at a distance of a few arcseconds. Finally, the source in the bottom-right panel does not lie on a cusp and the 3-1 image geometry is broken. All images have roughly similar magnification and there is no clear group of three closely grouped images.

## 6 Interacting halos

We point out in this section that even grazing interactions between neighboring halos are likely to produce transient misalignments due to tidal torquing. While these offsets might not persist for significant periods of time, they are worth exploring as a special case of the generic misaligned case. We use the disk$`+`$halo model described in § 4 but add a secondary halo with about 30% of the mass of the primary at a distance of one to several arc-seconds away from the primary lens. Simulations of such interacting halos suggest that the disk misalignment will tend to be along the direction of the secondary. We have explored a variety of halo separations and choose a distance of 1.7 arc-seconds to illustrate the effect. This separation in fact corresponds to a spatial distance of 2$`r_\mathrm{d}`$, where $`r_\mathrm{d}=3.2`$ kpc is the disk scale length for the proto-type Milky Way type spiral. The magnification maps in the source plane are shown in Fig. 8. The main effects are an increase in the area enclosed by the caustics, an overall shift of the caustic structure, and induced distortions in the direction of the secondary halo. Comparing Panel (c) of Fig. 1 to Panel (c) of Fig. 8 also shows that the radial caustic is less pinched at the lower edge. Close pairs of galaxies are in fact observed as lenses; B1608+656 (Koopmans & Fassnacht 1999) is a clear example of this. If one or both of the lens galaxies in this system contain discs, a scenario as discussed above is likely to arise.

## 7 Magnification ratios

There has been much interest recently in explaining some of the observed magnification ratios with lens models that include compact substructure with masses that are approximately $`10^6M_{}`$. Such substructure is ubiquitous in numerical simulations of structure formation in standard cold dark matter cosmologies, but has not yet been observed. An argument in favor for the possible ubiquity of dark clumps is that there are many viable mechanisms that could prevent star formation in such compact sub-halos (Somerville 2002) so that they could indeed exist but be unobservable directly.

Gravitational lensing does in principle provide a method to detect dark matter substructure on these small scales, and the recent arguments presented by Metcalf & Madau (2001), Metcalf & Zhao (2002), Dalal & Kochanek (2002) and Bradać et al. (2002) are compelling. The current claim is that in systems like B1422+231, the observed image magnification ratios cannot be obtained without substructure. The most general and convincing case was originally presented by Schneider & Mao (1998) for B1422+231 where there are 3 highly magnified images on one side of the lens centre (images A,B and C), forming a line. On the other side of the lens lies a single, fourth image D, that has a radio flux that is a factor of 50 lower than that of the other three images. This means that either image D is demagnified or the other 3 images are magnified very strongly. This, however would imply that the source has to lie close to a cusp. There exists a theorem that holds strictly only for sources on a cusp: the sum of the flux of the outer two images ought to be equal to the flux of the innermost image, or, written in terms of the magnifications:

$$\mu _\mathrm{A}+\mu _\mathrm{C}=\mu _\mathrm{B}\mathrm{for}\mu _{\mathrm{tot}}\mathrm{},$$
(2)
where B is the central, brightest, image. However, the images in B1422+231 do not obey this relation. Moving the source away from the cusp would mean the theorem need not hold but implies that the total magnification is small and hence that image D needs to be strongly demagnified. Mao & Schneider (1998) argue that this is unlikely as the image is far from the lens center and hence in a low $`\kappa `$ region, leaving only the possibility that the source does lie close to a cusp but that the magnification of one or more of the images is effected by substructure in the vicinity.

Recently, Keeton et al. (2002) investigated the lensing properties of sources lying in cusps in more detail. His analysis is based on a generic one-component elliptical power law lens model, including external shear. He finds that the image magnification ratios of a number of observed lens systems, B2045+265, RX J0911+0551 and B0712+472 can not be fit by such ”cusp lenses”. No detailed modeling is done in that paper and the argument is based on the deviant magnification ratios alone. The magnification ratio $`(\mu _\mathrm{A}+\mu _\mathrm{C})/\mu _\mathrm{B}`$ of observed systems lies typically between 1.5 and 2 for total magnifications of $`40`$. One-component lens models of these systems predict a magnification ratio $`(\mu _\mathrm{A}+\mu _\mathrm{C})/\mu _\mathrm{B}`$ that is within a few per-cent of unity. The question arises whether the mismatch between observed and simulated image magnification ratios is not due to the specific one-component lens model used; a more complicated lens mass profile may explain the observed magnification ratios.

As seen from the maps in Figs. 1-8 mergers and misalignments between a halo and disk can cause a strong change in the magnification pattern on the image plane. It is thus not unreasonable to propose that the magnification ratios might be explained by a more complex structure in the primary lens – for example the presence of a misaligned disk, as discussed in this paper.

In Figs. 9 and 10 we show the magnification ratio $`1/|(\mu _\mathrm{A}+\mu _\mathrm{C})/\mu _\mathrm{B}1|`$ for our Milky Way-type model and for the early-type lens respectively for different offsets and disk parameters. In all cases, only regions in the source plane that would produce 4 detectable images are shown. In the case of the Fig. 9 this is only true for regions in the source plane that form a total of 5 images. For the early-type lens model, regions that produce a total of 7 images in general also produce 4 bright images, the remaining three images being strongly demagnified (cf. Fig 7, top-left and bottom-left panels). To determine the magnification ratios for these systems we proceed as follows: for each source position as defined by a grid on the source plane consisting of $`400\times 400`$ pixels of size $`2.5\times 10^3`$arcsec, we identify the 4 most strongly magnified images. Within this group, the image with the smallest magnification, image D, is identified. Within the remaining group of three magnified images we identify the central image as image B. The remaining two images are assigned labels A and C randomly. Having thus identified the images A,B,C and D for each source position, we determine the magnification ratios from the individual magnifications as described in § 3.

The gray-scale in Figs. 9 and 10 is such that lightly-shaded regions indicate source positions for which the magnification ratios do not obey the cusp relation. The top left panel in Fig. 8 shows the case for no halo-disk misalignment for our Milky Way-type lens model; the disk in the remaining three panels has a moderate misalignment in the x, the y and the x-y direction. The caustic shape changes significantly once the disk is misaligned, leading to strong changes in the function of magnification ratio vs. source position. Especially for a misalignment in the y directory (top-left panel) and for misalignments of $`\mathrm{\Delta }x=\mathrm{\Delta }y=0.6\mathrm{"}`$, the caustic shape becomes very complicated and in several cusp regions the expected magnification differs drastically from that expected in ”normal” cusps. For the early-type lens, we show the magnifcation ratio only for offsets of $`\mathrm{\Delta }x=\mathrm{\Delta }y=0.6\mathrm{"}`$, but for different disk parameters. The top left panel is for the standard early-type lens model described in the previous section. The three other panels are for three different sets of disk masses and sizes. The bottom-right panel shows the results for a small, light disk. In all cases sources lying in the upper-right cusp would produce 3-1 image geometries with $`(\mu _\mathrm{A}+\mu _\mathrm{C})/\mu _\mathrm{B}`$ magnification ratios that differ significantly from the expected ratio of 1.

Note that for the Milky Way-type lens there is a strong effect due the presence of the disk even when it is not misaligned. We show the corresponding result for our standard early-type lens model in Fig. 11. The scale is drastically different than that of Figs. 9 and 10; the caustic structure is very small. Since we assume a spherical halo, the asymmetry in the potential that causes the inner caustic is solely due to the small disk embedded in the halo. However, there is not much qualitative change compared to the MW-type lens: a disk affects the magnification ratios of 4-image lens systems strongly also in E-type lens galaxies. This is in concordance with the findings of Möller, Hewett & Blain who studied the effects of such aligned disks in early-type lensing galaxies (MNRAS, submitted).

## 8 Conclusions

As has been illustrated in the plots in this paper, the presence of a misaligned disky component embedded in a dark matter halo produces a complex range of strong lensing effects - in the shapes and areas of the caustics, magnification ratios, and image geometries.

Most importantly, we show that for models of high-magnification ($`\mu _{\mathrm{tot}}50`$) systems the presence of a disky component can substantially alter the expected magnification ratios and image geometries. This is true for both Milky Way-type disks and small disks as may be found in early-type galaxies. But misalignments of the disk relative to the halo, such as might be caused by mergers or strong interactions, can compound this effect. The shape of the caustic structure changes dramatically, and depending on the ellipticity of the dark matter halo the enclosed area is enhanced by up to a factor of 10. Thus, misaligned disks increase the overall cross-section for the formation of four magnified images, while at the same time causing strong deviation from the magnification ratios expected from simple one-component lens models. Single component, elliptical mass models are unable to reproduce such features, which may be the reason why models based on simple elliptical mass models fail to reproduce the observed image magnification ratios as for B1422+231.

Due to the complexity of the lens model investigated here we have been concerned mainly with investigating the qualitative effect of misaligned disks. Even though detailed modeling of existing systems using a halo in combination with a misaligned disk lens is in principle possible, this will be a complicated and time-consuming task. Due to the large number of parameters and the relatively few observational lensing constraints, it is very unlikely that a definite model can be found; instead, it is to be expected that there will be a considerable degree of degeneracy. However, our results do show clearly that the simple one-component lens models currently used are likely to be too simplistic and that more complicated models including a possibly misaligned disk component can reproduce the observed properties of many lens systems without introducing the need for halo sub-clumps.

PN acknowledges a Research Fellowship from Trinity College, Cambridge. OM acknowledges support from the Marie Curie Fellowship programme.


text:
# 1 Introduction

## 1 Introduction

First-principles calculations of alloy thermodynamic properties have been successfully employed in a variety of contexts for metallic, semi-conductor and ceramic systems, including the computation of: composition-temperature phase diagrams, thermodynamic properties of stable and metastable phases, short-range order in solid solutions, thermodynamic properties of planar defects (including surfaces or antiphase and interphase boundaries), and the morphology of precipitate microstructures .

Although the formalism that allows the calculation of thermodynamic properties from first-principles has been known for decades , its practical implementation remains tedious. These practical issues limit the accuracy researchers are able to obtain without spending an unreasonable amount of their time writing input files for various computer codes, monitoring their execution and processing their output. These practical difficulties also limit the community of researchers that use these methods solely to those that possess the necessary expertise to carry out such calculations.

The Alloy Theoretic Automated Toolkit (ATAT)  drastically simplifies the practical use of these methods by implementing decision rules based on formal statistical analysis that free the researchers from a constant monitoring during the calculation process and automatically “glues” together the input and the output of various codes, in order to provide a high-level interface to the calculation of thermodynamic properties from first-principles. In order the make this powerful toolkit available to the wide community of researchers who could benefit from it, this article present a concise user guide to this toolkit.

## 2 Theoretical Background

While there exist numerous methodologies that enable the calculation of thermodynamic properties from first-principles, we will focus on the following two-step approach (see Figure 1). First, a compact representation of the energetics of an alloy, known as the cluster expansion , is constructed using first-principles calculations of the formation energies of various atomic arrangements. Second, the cluster expansion is used as a Hamiltonian for Monte Carlo simulations  that can provide the thermodynamic properties of interest, such as the free energy of a phase or short-range-order parameters as a function of temperature and concentration. This two-step approach is essential, because the calculation of thermodynamic quantities through Monte Carlo involves averaging the property of interest over many different atomic configurations and it would be infeasible to calculate the energy of each of these configurations from first-principles. The cluster expansion enables the prediction of the energy of any configuration from the knowledge of the energies of a small number of configurations (typically between 30 and 50), thus making the procedure amenable to the use of first-principles methods.

Formally, the cluster expansion is defined by first assigning occupation variables $`\sigma _i`$ to each site $`i`$ of the *parent lattice*, which is defined as the set of all the atomic sites that can be occupied by one of a few possible atomic species. In the common case of a binary alloy system, the occupation variables $`\sigma _i`$ take the value $`1`$ or $`+1`$ depending on the type of atom occupying the site. A particular arrangement of these “spins” on the parent lattice is called a *configuration* and can be represented by a vector $`\sigma `$ containing the value of the occupation variable for each site in the parent lattice. Although we focus here on the case of binary alloys, this framework can be extended to arbitrary multicomponent alloys (the appropriate formalism is presented in ).

The cluster expansion then parametrizes the energy (per atom) of the alloy as a polynomial in the occupation variables:

$$E(\sigma )=\underset{\alpha }{}m_\alpha J_\alpha \underset{i\alpha ^{}}{}\sigma _i$$
(1)
where $`\alpha `$ is a cluster (a set of sites $`i`$). The sum is taken over all clusters $`\alpha `$ that are not equivalent by a symmetry operation of the space group of the parent lattice, while the average is taken over all clusters $`\alpha ^{}`$ that are equivalent to $`\alpha `$ by symmetry. The coefficients $`J_\alpha `$ in this expansion embody the information regarding the energetics of the alloy and are called the effective cluster interaction (ECI). The *multiplicities* $`m_\alpha `$ indicate the number of clusters that are equivalent by symmetry to $`\alpha `$ (divided by the number of lattice sites).

It can be shown that when *all* clusters $`\alpha `$ are considered in the sum, the cluster expansion is able to represent any function $`E\left(\sigma \right)`$ of configuration $`\sigma `$ by an appropriate selection of the values of $`J_\alpha `$. However, the real advantage of the cluster expansion is that, in practice, it is found to converge rapidly. An accuracy that is sufficient for phase diagram calculations can be achieved by keeping only clusters $`\alpha `$ that are relatively compact (e.g. short-range pairs or small triplets). The unknown parameters of the cluster expansion (the ECI) can then be determined by fitting them to the energy of a relatively small number of configurations obtained through first-principles computations. This approach is known as the Structure Inversion Method (SIM) or the Collony-Williams  method.

The cluster expansion thus presents an extremely concise and practical way to model the configurational dependence of an alloy’s energy. A typical well-converged cluster expansion of the energy of an alloy consists of about 10 to 20 ECI and necessitates the calculation of the energy of around 30 to 50 ordered structures (see, for instance, ). Once the cluster expansion has been constructed, the energy of any configuration can be calculated using Equation 1 at a very small computational cost. This enables the use of various statistical mechanical techniques such as Monte Carlo simulations , the low-temperature expansion (LTE) , the high-temperature expansion (HTE) , or the cluster variation method (CVM)  to calculate thermodynamic properties and phase diagrams. The ATAT software implements Monte Carlo simulations, the LTE and the HTE.

Paralleling the two-step approach described in the previous section, ATAT consists of two main computer programs (see Figure 1). The cluster expansion construction is performed by the MIT Ab initio Phase Stability (MAPS) code , while the Monte Carlo simulations are driven by the Easy Monte Carlo Code (EMC2), developed at Northwestern University . Each of these codes will be discussed in turn.

While the present user guide describes how the ATAT software can be used to carry out all the steps necessary for the calculation of thermodynamic properties from first-principles, it must be emphasized that each part of the toolkit can be used as a stand-alone code. For instance, many users may have access to an existing cluster expansion obtained through the SIM or other popular methods, such as concentration-wave-based methods (see, for instance, ). It is then staightforward to setup the appropriate input files to run the emc2 Monte Carlo code. Alternatively, after obtaining a cluster expansion using the maps code, users could choose to calculate thermodynamic properties with the cluster variation method (CVM) , as implemented in the IMR-CVM code . The modularity of the toolkit actually extends below the level of the maps and emc2 codes — many of the subroutines underlying these codes can be accessed through stand-alone utilities .

## 3 Cluster expansion construction using the MAPS code

The maps code implements the so-called Structure Inversion Method (SIM), also known as the Connolly-Williams method . While the algorithms underlying the maps code are described in , the present section focuses on its practical use.

### 3.1 Input files

The maps code needs two input files: one that specifies the geometry of the parent lattice (lat.in) and one that provides the parameters of the first-principles calculations (xxxx.wrap, where xxxx is the name of the first-principles code used). The clear separation between the thermodynamic and first-principles calculations is a distinguishing feature of ATAT that enables the package to be easily interfaced with any first-principles code. Table 1 gives two annotated examples of a lattice geometry input file. The package includes ready-made lattice files for the common lattice types (e.g. bcc, fcc, hcp). It also includes an utility that automatically constructs multiple lattice geometry input files for common lattices. For instance,

* makelat Al,Ti fcc,bcc,hcp

creates 3 subdirectories containing the appropriate input files for each specified lattice.

The first-principles input file is usually less than 10 lines long, thanks to the dramatic improvements in the user-friendliness of most modern first-principles codes. For instance, in the case of the widely used VASP code , a typical input file is given in Table 2. Examples of such input files are provided with the package. Note that ATAT contains a utility that enables the automatic construction of $`k`$-point meshes from a single parameter defining the desired target $`k`$-point density, the number of $`k`$-point per reciprocal atom (KPPRA).

### 3.2 Running the code

The maps code is started using the command

* maps -d &

where the option -d indicates that all default values of the input parameters should be used, which is what most users will ever need. (The optional parameters can be displayed by typing maps by itself and further help is available via the command maps -h.) The trailing & character cause the command to execute in “background” mode. In this fashion, maps can continuously be on the lookout, responding to various “signals”, while the user performs other tasks. (The ongoing discussion assumes that the code is run under a UNIX environment within a shell such as sh, csh, tcsh or bash.)

The process of constructing a cluster expansion from first-principles calculations can be summarized as follows.

1. Determine the parameters of the first-principles code that provide the desired accuracy.
2. Let maps refine the cluster expansion.
3. Decide when the cluster expansion is sufficiently accurate.

Typically, one calibrates the accuracy of the first-principles calculations using the “pure”111In the case of pseudobinary alloys with spectator ions (e.g. the MgO-CaO system), the “pure” structures would correspond to the structures where the sublattice of interest is entirely filled with a single type of atom. structures of the alloy system of interest. To generate the two “pure” structures, type

* touch ready

This creates a file called ready which tells maps that you are ready to calculate the energy of a structure. Within 10 seconds, maps replies with

* Finding best structure...
* done!

maps has just created a directory called 0 and, within it, a file called str.out that contains the geometry of one of the two “pure” structures. If you type touch ready once more, the other “pure” structure is written to 1/str.out. You now need to launch the first-principles code to calculate the energy of each structure. Type

* cd 0

to go into the directory of the first structure. Assuming that your first-principles code is called xxxx, type

* runstruct\_xxxx &

After this command has successfully terminated, display the energy of that structure and go back to the initial directory

* cat energy
* cd ..

and edit the file defining the first-principles code parameters

* emacs xxxx.wrap &

so that the precision of the calculation is increased (e.g. increase the $`k`$-point density or the cut-off of the plane-wave energy). Then you rerun the calculations to check by how much the calculated energy has changed:

* cd 0
* runstruct\_xxxx &
* cat energy (After the calculations are completed)
* cd ..

This process is repeated until the user is satisfied with the precision of the calculation (that is, if the energy has become insensitive to changes in the input parameters within the desired accuracy).222A key number to keep in mind is that an error of 25meV corresponds to 300K on a temperature scale. A similar study should also be performed for the other “pure” structure (labeled structure 1) and, if one is really concerned with precision, for a few structures with intermediate concentrations.

Once the appropriate ab initio code parameters have been determined, the fully automated process can begin. From within the directory where maps was started, type

* pollmach runstruct\_xxxx

to start the job manager that will monitor the load on your local network of workstations and ask maps to generate new structures (i.e. atomic arrangements) whenever a processor becomes available. Note that the first time the command is run, instructions will appear on screen that explain how to configure the job dispatching system in accordance to your local computing environment. Once this configuration is complete, the above command should be invoked in the background by appending a “&” to it.

### 3.3 Output of MAPS

While the calculations are running, you can check on the status of the best cluster expansion obtained so far. The file log.out contains a brief description of the status of the calculations, such as the accuracy of the cluster expansion and various warning messages. Most of the messages pertains to the accurate prediction of the so-called ground states of the alloy system. The ground states, which are the structures that have the lowest energy for each given concentration, are extremely important to predict accurately because they determine which phases will appear on the phase diagram. The four possible messages are described below.

* Not enough known energies to fit CE. Before displaying any results, maps waits until enough structural energies are known to fit a minimal cluster expansion containing only nearest-neighbor pair interactions and test its accuracy. Thus, the first cluster expansion is typically constructed after at least 4 structural energies have been computed (this number may vary as a function of the symmetry of the lattice).
* Among structures of known energy, true ground states differ from fitted ground
states. The current cluster expansion incorrectly predicts which structures have the lowest energy for given concentrations, among structures whose first-principles energy is known. The code has built-in checks to avoid this. However, in rare instances, it may be mathematically impossible to satisfy all the constraints that the code imposes for a cluster expansion to be acceptable. This problem becomes increasingly unlikely as the number of calculated structural energies increases, so the user should just wait until the problem fixes itself.
* Among structures of known energy, true and predicted ground states agree. Opposite of the previous message. When this message is displayed, maps also displays either one of the following two messages.
* New ground states of volume less or equal to $`n`$ predicted, see predstr.out. This indicates that the cluster expansion predicts that, at some concentration, there exist other structures that should have an energy even lower than the one of the structures whose energy has been calculated from first-principles. In this case, maps will investigate the matter by generating those structures and requesting that their energy be calculated. Once again, the user should just wait for the problem to fix itself. The predicted ground states are flagged by a g in the predstr.out file, so that you can display their energy by typing

+ grep g predstr.out
* No other ground states of $`n`$ atoms/unit cell or less exist. The energies of all ground states predicted by the cluster expansion have been confirmed by first-principles calculations. Because it can be computationally intensive to perform a full ground state search when interactions extend beyond the nearest-neighbor shell , maps uses a search algorithm that merely enumerates every possible structures having $`n`$ atoms or less per unit cell and uses the cluster expansion to predict their energies. The upper limit $`n`$ increases automatically as calculations progress.

The log.out file also contains two other pieces of information:

* Concentration range used for ground state checking: \[$`a`$,$`b`$\] This displays the user-selected range of concentration over which ground state checking is performed (which can be specified as a command-line option of the maps command: -c0=$`a`$ -c1=$`b`$). It may be useful to relax the constraints that ground states be correctly reproduced over the whole concentration range when it is known that other parent lattices are stable in some concentration range. In this fashion, the code can focus on providing a higher accuracy in the concentration range where the user needs it.
* Crossvalidation score: $`s`$. This provides the predictive power of the cluster expansion. It is analogous to the root mean square error, except that it is specifically designed to estimate the error made in predicting the energy for structures not included in the least-squares fit . It is defined as

$$CV=\frac{1}{n}\underset{i=1}{\overset{n}{}}\left(E_i\widehat{E}_{\left(i\right)}\right)^2$$

where $`E_i`$ is the calculated energy of structure $`i`$, while $`\widehat{E}_{\left(i\right)}`$ is the predicted value of the energy of structure $`i`$ obtained from a least-squares fit to the $`\left(n1\right)`$ other structural energies.

The maps code also outputs quantitative data in various output files. The simplest way to analyze this data is by typing

* mapsrep

As illustrated in Figure 2, this command displays, in turn

* The log.out file described earlier.
* The formation energy of all structures whose energy is known from first-principles calculations, as well as the predicted energy of all structures maps has in memory. The convex hull of the ground states among structures of known energy is overlaid while the new predicted ground states (if any) are marked by an “$`\times `$”. (Note that this ground state line is only meaningful if the log.out file contains “Among structures of known energy, true and predicted ground states agree.”)
* The formation energy of all structures calculated from first-principles and associated ground state line.
* A plot of the magnitude of the Effective Cluster Interactions (ECI) as a function of the diameter of their associated cluster (defined as the maximum distance between any two sites in the cluster). Pairs, triplets, etc. are plotted consecutively. This plot is useful to assess the convergence of the cluster expansion. When the magnitude of the ECI for the larger clusters has clearly decayed to a negligible value (relative to the nearest-neighbor pair ECI), this is indicative of a well-converged cluster expansion.
* A plot of the residuals of the fit (i.e. the fitting error) for each structure. This information is useful to locate potential problems in the first-principles calculations. Indeed, when first-principles calculations exhibit numerical problems, this typically results in calculated energies that are poorly reproduced by the cluster expansion.

When the user is satisfied with the results (which are constantly updated), maps can be stopped by creating a file called stop in the current directory using the command:

* touch stop

while the job dispatching system can be stopped by typing:

* touch stoppoll

A cluster expansion can be considered satisfactory when

1. All ground states are correctly reproduced and no new ground states are predicted. (The log.out file would then indicate that Among structures of known energy, true and predicted ground states agree. No other ground states of $`n`$ atoms/unit cell or less exist.)
2. The crossvalidation score, as given in the log.out file, is small (typically less than 0.025 eV).
3. Optionally, it is instructive to verify that the magnitude of the ECI decays as a function of the diameter of the corresponding cluster and as a function of the number of sites it contains.

## 4 Monte Carlo simulations

The emc2 code implements semi-grand canonical Monte Carlo simulations , where the total number of atoms is kept fixed, while the concentration is allowed to adapt to an externally imposed difference in the chemical potential of the two types of atoms. The chemical potential difference will be simply referred to as the “chemical potential” in what follows. This ensemble offers the advantage that, for any imposed chemical potential, the equilibrium state of the system is a single phase equilibrium, free of interfaces.333If the simulation cell is commensurate with the unit cell of the phase under study, a requirement that the code automatically ensures. It also simplifies the process of calculating free energies through thermodynamic integration. While a detailed description of the algorithm underlying this code can be found in , the current section focuses on the practical usage of the code.

### 4.1 General input parameters

The Monte Carlo code needs the following files as an input

1. A lattice geometry file (lat.in), which is the same as the input for maps (see Table 1).
2. Files providing the cluster expansion (the clusters used are listed in the clusters.out file while the corresponding ECI are listed in the eci.out file.). These files are automatically generated by maps, although users can supply their own cluster expansion, if desired. A description of the format of these files is available by typing maps -h.
3. A list of ground states (gs\_str.out), which merely provide convenient starting configurations for the simulations. maps also automatically creates this file.

The parameters controlling the simulation are specified as command-line options. The first input parameter(s) needed by the code are the phase(s) whose thermodynamic properties are to be determined. There are two ways to invoke the Monte Carlo simulation code. When the command emc2 is used, a single Monte Carlo simulation is run to allow the calculation of thermodynamic properties of a single phase for the whole region of chemical potential and temperature where that phase is stable. The phase of interest is specified by a command-line option of the form

* -gs=$`n`$,

where $`n`$ is a number between $`1`$ and $`G1`$ (inclusive), where $`G`$ is the number of ground states. The value $`1`$ indicates the disordered phase while values ranging from $`0`$ to $`G1`$ indicate the phases associated with each ground states ($`0`$ denoting the ground state with the smallest composition). When the command phb is used, two Monte Carlo simulations are run simultaneously to enable the determination of the temperature-composition phase boundary associated with a given two-phase equilibrium. The two phases are specified by

* -gs1=$`n_1`$ -gs2=$`n_2`$.

It is possible to compute a two-phase equilibrium between phases defined on a different parent lattice. In this case, the user needs to specify the directories where the cluster expansions of each lattice resides using the options of the form

* -d1=directory 1 -d2=directory 2

The accuracy of the thermodynamic properties obtained from Monte-Carlo simulations is determined by two parameters: The size of the simulation cell and the duration of the simulation.

The size of the simulation cell is specified by providing the radius $`r`$ of a sphere through the command-line option

* -er=$`r`$.

As illustrated in Figure 3a, the simulation cell size will be the smallest supercell that both contains that sphere and that is commensurate with the unit cell of the ground state of interest. This way of specifying the simulation cell size ensures that the system size is comparable along every direction, regardless of the crystal structure of the ground state of interest. It also frees the user from manually checking the complicated requirement of commensurability. It is important that the user check that the simulation cell size is sufficiently large for the thermodynamic properties of interest to be close to their infinite-system-size limiting value. This can be done by gradually increasing the system size until the calculated quantities become insensitive to the further increases in system size, within the desired accuracy.

The duration of the simulations is automatically determined by the code from a user-specified target precision on the atomic composition of the phase, indicated by a command-line option of the form

* -dx=$`\mathrm{\Delta }x`$.

Alternatively, the user may also manually set the number $`n_{eq}`$ of Monte Carlo steps the system is allowed to equilibrate before thermodynamic averages are computed over a certain number $`n_{avg}`$ of Monte Carlo steps using the options

* -eq=$`n_{eq}`$ -n=$`n_{avg}`$.

The Monte Carlo code also needs additional parameters that specify which portion of a phase’s free energy surface needs to be computed. With emc2, the range of temperatures to be scanned are specified in either one of the following two ways:

* -T0=$`T_0`$ -T1=$`T_1`$ -dT=$`\mathrm{\Delta }T`$ (for steps in direct temperature)
* or -T0=$`T_0`$ -T1=$`T_1`$ -db=$`\mathrm{\Delta }(1/T)`$ (for steps in reciprocal temperature).

The temperature steps in reciprocal temperature ($`\mathrm{\Delta }(1/T)`$) can be useful when calculations are started from infinite temperatures down to a finite temperature. The -T1 and -dT (or -db) options can be omitted if calculations at a single temperature are desired. Since the program automatically stops when a phase transition is detected, it is not necessary to know in advance the temperature range of stability of the phase. The user only needs to ensure that the initial temperature lies within the region of stability of the phase of interest. An obvious starting point is $`T_0=0`$, since the ground state is then stable, by definition. With the phb code, the syntax is

* -T=$`T`$ -dT=$`\mathrm{\Delta }T`$

If the -T option is omitted, calculations start at absolute zero.444Temperature steps in reciprocal temperature are not needed, because a two-phase equilibrium never extents up to infinite temperature. The energy and temperature units used are set by specifying the Boltzman’s constant with the command-line option

* -k=$`k_B`$.

A value of $`8.617\times 10^5`$ corresponds to energies in eV and temperatures in Kelvin.

With emc2, the range of chemical potentials to be scanned needs to be specified. Once again, only the starting point really matters, because the code will stop when a phase transition is reached. By default, chemical potentials are given in a dimensionless form, so as to facilitate the link between the value of the chemical potential and the phase that it stabilizes. For instance, a chemical potential equal to $`3.0`$ is such that it would stabilize a two phase equilibrium between phase number $`2`$ and phase number $`3`$ at absolute zero (see Figure 3b). A chemical potential between $`3.0`$ and $`4.0`$ stabilizes phase number $`3`$ at absolute zero. While these ranges of stability are no longer exact at finite temperature, this dimensionless chemical potential still provides easy-to-interpret input parameters. The syntax is

* -mu0=$`\mu _0`$ -mu1=$`\mu _1`$ -dmu=$`\mathrm{\Delta }\mu `$

where $`\mathrm{\Delta }\mu `$ is the chemical potential step between each new simulation. Chemical potentials can also be entered in absolute value (say in eV, if the energies are in eV) by specifying the -abs option. Note that the output files always give the absolute chemical potentials, so that thermodynamic quantities can be computed from them. With phb, the initial chemical potential is optional when starting from absolute zero because the code can determine the required value from the ground state energies. It can nevertheless be specified (in absolute value) with the -mu=$`\mu `$ option, if a finite temperature starting point is desired.

A list of the command line options of either the emc2 or phb codes can be displayed by simply typing either command by itself. More detailed help is displayed using the -h option.

### 4.2 Examples

We now give simple examples of the usage of these commands. Consider the calculation of the free energy of the phase associated with ground state number $`1`$ as a function of concentration and temperature. Then, the required commands could, for instance, be

* emc2 -gs=1 -mu0=1.5 -mu1=0.5 -dmu=0.04 -T0=300 -T1=5000 -dT=50 -k=8.617e-5 -dx=1e-3 -er=50 -innerT -o=mc10.out
* emc2 -gs=1 -mu0=1.5 -mu1=2.5 -dmu=0.04 -T0=300 -T1=5000 -dT=50 -k=8.617e-5 -dx=1e-3 -er=50 -innerT -o=mc12.out

(The only difference in the two command lines is the value of -mu1 and the output file name, specified by the -o option.) These commands separately compute the two “halfs” of the free energy surface, corresponding to the values of the chemical potential below and above the “middle” value of $`1.5`$ which stabilizes ground state $`1`$ at absolute zero. This natural separation allows you to run each half calculation on a separate processor and obtain the results in half the time. The values of -dmu, -dT, -dx and -er given here are typical values. The user should ensure that these values are such that the results are converged. Note that, thanks to the way these precisions parameters are input, if satisfying values have been found for one simulation, the same values will provide a comparable accuracy for other simulations of the same system. The option -innerT indicates that the inner loop of the sequence of simulations scans the temperature axis while the outer loop scans the chemical potential. In this fashion, the point of highest temperature in the region of stability of the phase will be known early during the calculations. If the user is more interested in obtaining solubility limits early on, this option can be omitted and the inner loop with scan the chemical potential axis. In either cases, the code exits the inner loop (and the outer loop, if appropriate) when it encounters a phase transition.

The emc2 code thus enables the automated calculation of the whole free energy surface of a given phase, as illustrated in Figure 4a. Such free energy surfaces can be used as an input to construct thermodynamic databases or supplement existing ones. To facilitate this process, a utility that converts the output of emc2 into input files for the fitting module of ThermoCalc is provided.

While the above examples focus on the calculation of a phase’s thermodynamic properties over its whole region of stability, one may be interested in directly computing the temperature-composition phase boundary without first constructing a full free energy surface. To accomplish this task, a typical command-line for the phb program would be

* phb -gs1=0 -gs2=1 -dT=25 -dx=1e-3 -er=50 -k=8.617e-5 -o=ph01.out

This command computes the two phase equilibrium between phase 0 and phase 1, starting at absolute zero and incrementing temperature in steps of $`25`$ K. The output file ph01.out contains the temperature-composition phase boundary of interest, as well as the chemical potential stabilizing the two-phase equilibrium as a function of temperature. This output can be used to generate phase diagrams, as illustrated in Figure 4b.

The program automatically terminates when the “end” of the two-phase equilibrium has been reached. If the two-phase equilibrium disappears because of the appearance of a third phase, two new two-phase equilibria have to be separately calculated. To do so, one uses the final temperature $`T`$ and chemical potential $`\mu `$ given in the output file as a starting point for two new phb runs:

* phb -T=$`T`$ -mu=$`\mu `$ -gs1=0 -gs2=-1 -dT=25 -dx=1e-3 -er=50 -k=8.617e-5 -o=ph0d.out
* phb -T=$`T`$ -mu=$`\mu `$ -gs1=-1 -gs2=1 -dT=25 -dx=1e-3 -er=50 -k=8.617e-5 -o=phd1.out

In the above example, it is assumed that the new phase appearing is the disordered phase (indicated by the number $`1`$), which will usually be the case. Of course, it is also possible that a given two-phase equilibrium terminates because one of the two phases disappears. In this case, only one new calculation needs to be started, as in the following example:

* phb -T=$`T`$-mu=$`\mu `$ -gs1=0 -gs2=2 -dT=25 -dx=1e-3 -er=50 -k=8.617e-5 -o=phd1.out

Note that phase $`1`$ has been replaced by phase $`2`$. Finally, it is also possible that the two-phase equilibrium terminates because the concentration of each phase converges to the same value, a situation which requires no further calculations. The user can easily distinguish these three cases by merely comparing the final composition of each phase.

### 4.3 Interpreting the output files

The output file of emc2 reports the value of all calculated thermodynamic functions for each value of temperature and chemical potential scanned. The quantities reported include

* Statistical averages over Monte Carlo steps, such as energy, concentration, short-range and long-range order parameters.555For efficiency reasons, the long range order parameters are only calculated when starting from an ordered phase.
* Integrated statistical averages, such as the Gibbs free energy $`G`$ or the semi-grand-canonical potential $`\varphi =G\mu x`$.
* The result of common approximations, namely, the low temperature expansion (LTE) , the mean-field (MF) approximation and the high temperature expansion (HTE) (see, for instance, ).

While quantities obtained from statistical averages over Monte Carlo steps are valid for all temperatures and chemical potentials, caution must be exercised when interpreting the result of the various approximations or when looking at the integrated quantities. The LTE, MF and HTE approximations are only accurate in a limited range of temperature and it is the responsibility of the user to assess this range of validity. Also, the free energy or the semi-grand-canonical potential are obtained from thermodynamic integration and are thus only valid if the starting point of the integration is chosen appropriately. By default, the low temperature expansion value is used as a starting point whenever the phase of interest is a ground state, while the high temperature expansion is used when the phase of interest is the disordered state. Hence, to obtain absolute values of the semi-grand-canonical potential, one must ensure that the calculations are started at a sufficiently low temperature (or sufficiently high temperature, in the case of the disordered phase). This can be checked by comparing the Monte Carlo estimates with the LTE (or HTE) estimates and verifying that they agree for the first few steps of the thermodynamic integration. A user-specified starting point for $`\varphi `$ (e.g. obtained from an earlier Monte Carlo simulation) can be indicated using the option

* -phi0=$`\varphi `$

Note that, unlike emc2, the phb code automatically makes use of the low temperature expansion whenever it is sufficiently accurate in order to save a considerable amount of computational time.

By default, the code reports the thermodynamic quantities associated with the semi-grand-canonical ensemble, such as the semi-grand-canonical potential $`\varphi `$. The command-line option -can, instructs the code to add $`\mu x`$ to all appropriate thermodynamic quantities, so that the code outputs the more commonly used canonical quantities, such as the Gibbs free energy $`G`$ and the internal energy $`E`$.

## 5 Future developments

At the present time, ATAT only handles binary systems. However, a multicomponent version of the toolkit is under development. About 75% of the code has already been written with multicomponent systems in mind.

Although the present tutorial does not discuss the topic, ATAT also implements reciprocal space cluster expansions and, in particular, the constituent strain formalism . A tutorial on the use of this feature will be available shortly on our web site .

In its current stage of development, the ATAT package only focuses on configurational sources of entropy. To obtain quantitative agreement with experimental measurements (such as calorimetry), it may be necessary to explicitly add the vibrational and electronic contributions to the free energy . A simple correction that makes the calculated free energy comparable to experiments would be to add a concentration-weighted average of the known experimental free energies of the appropriate elements. Of course, such a shift does not affect the calculated phase diagram, unless it exhibits phases that are superstructures of distinct parent lattices. In the near future, a module that computes electronic and vibrational free energy contributions from first-principles will be added. (In fact, the Monte Carlo code already includes all the features required to accommodate nonconfigurational entropy, such as allowing for temperature-dependent ECI.)

## 6 Conclusion

The Alloy Theoretic Automated Toolkit (ATAT) drastically simplifies the practical implementation of the Connolly-Williams method, in combination with semi-grand-canonical Monte Carlo simulations, thus providing a high-level interface to the calculation of thermodynamic properties from first-principles. This toolkit enables researcher to focus on higher-level aspects of first-principles thermodynamic calculations by encapsulating the intricate details of the calculations under an easy-to-use interface. It also makes these powerful methodologies readily available to the wide community of researchers who could benefit from it.

## Acknowledgements

This work was supported by the NSF under program DMR-0080766 and by the U.S. Department of Energy, Office of Basic Energy Sciences, under contract no. DE-F502-96ER 45571.