text:
# Critical thermodynamics of the two-dimensional ±𝐽 Ising spin glass

## The model and our measurements —

The Hamiltonian of our two-dimensional ($`2d`$) spin glass is

$$H_J(\{\sigma _i\})\underset{ij}{}J_{ij}\sigma _i\sigma _j$$
(1)
where the sum runs over all nearest neighbor pairs of Ising spins ($`\sigma _i=\pm 1`$) on a square lattice of volume $`V=L\times L`$ with periodic boundary conditions. The quenched random couplings $`J_{ij}`$ take the value $`\pm J`$ with probability $`1/2`$. The partition function at inverse temperature $`\beta T^1`$ is $`Z_J=_{\{\sigma _i\}}e^{\beta H_J(\{\sigma _i\})}`$ and can be written as

$$Z_J(\beta )=e^{2L^2\beta J}P_J(X=e^{2\beta J}).$$
(2)
Here $`P_J(X)`$ is the polynomial whose coefficient of $`X^p`$ is the number of spin configurations of energy $`E=(2L^2+2p)J`$. Saul and Kardar Saul and Kardar (1993, 1994) showed that determining $`P_J`$ can be reduced to computing determinants which they did using exact arithmetic of arbitrarily large integers. More recently a more powerful approach has been developed Galluccio et al. (2000, 2001), based on the use of modular arithmetic to compute pfaffians. With this algorithm, one first finds the coefficients modulo a prime number, thereby avoiding costly arbitrary precision arithmetic. Then the computation is repeated for enough different primes to allow the reconstruction of the actual (huge) integer coefficients using the Chinese remainder theorem.

The algorithm proposed and implemented in Galluccio et al. (2000, 2001) is powerful enough to solve samples with $`L100`$; the total CPU time needed to compute $`Z_J`$ grows approximately as $`L^{5.5}`$. In our study we have determined $`Z_J`$ for a large number of disorder samples at different lattice sizes: for instance we have $`400000`$ samples at $`L=6`$, $`100000`$ at $`L=10`$, $`10000`$ at $`L=30`$, $`1000`$ at $`L=40`$ and $`300`$ at $`L=50`$. The total computation time used is equivalent to about 40 years of a 1.2 GHz Pentium processor. For each sample we derive from $`Z_J`$ various thermodynamic quantities such as the free energy $`F_J(\beta )=\beta ^1\mathrm{ln}Z_J`$, the internal energy $`U_J(\beta )=H_J`$, and the specific heat $`dU_J/dT`$. We also study in detail the number of ground states and of excited states. Note that flipping any spin changes the energy by $`0`$, $`\pm 4J`$ or $`\pm 8J`$; the gap between the ground state and the first excited state is thus $`4J`$.

## Low temperature behavior of $`c_V`$

The study of $`2d`$ Ising spin glasses has a long history. We will only discuss here results about the $`\pm J`$ model. It is generally agreed that this model is paramagnetic for $`T>0`$, spin-glass ordering arising only as $`T0`$. The critical region thus corresponds to $`TT_c=0`$. Since there is an energy gap $`4J`$, the free energy should have a singularity of the form $`\mathrm{exp}(4J/T)`$. This is difficult to check, in particular via Monte Carlo where the free energy is not directly measurable. Instead, it is better to concentrate on the specific heat density $`c_V`$. For that observable, the difference between the models with bimodal ($`J_{ij}=\pm J`$) and continuous couplings is striking: in the first case $`c_V`$ goes to zero rapidly as $`T0`$ while in the second there is a clear linear behavior.

Even though our computations provide us with the free energy, we also prefer to work with $`c_V`$. Note that $`c_V`$ is related to a second derivative of the free energy so the corresponding singularities are directly related. Also, $`c_V`$ should provide a cleaner signal as irrelevant “constants” such as the ground state energy that fluctuate from sample to sample have been subtracted out. Consider now any given sample. As $`T0`$, we have the scaling

$$c_V\frac{\left[HH\right]^2}{L^2T^2}\frac{16J^2e^{S_1S_0}e^{4J/T}}{L^2T^2}$$
(3)
where $`S_0`$ and $`S_1`$ are the logarithms of the degeneracies of the ground state and first excited state energy levels for the given sample. ($`S_0`$ and $`S_1`$ are microcanonical entropies; we have dropped the index $`J`$ denoting a sample dependence.) Note that $`4J`$ appears because it is the energy gap in our system. It thus seems unavoidable that $`c_V`$ will have an $`\mathrm{exp}(4J/T)`$ singularity. Surprisingly, in 1988, Wang and Swendsen Wang and Swendsen (1988) postulated that instead

$$c_VT^p\mathrm{exp}\left(AJ/T\right)$$
(4)
with $`A=2`$. They performed a Monte Carlo study in which $`A3`$ for most of the temperatures they could access, but their effective $`A`$ drifted and their final prediction was $`A=2`$ from an analogy with a one dimensional model (we shall come back to this later). This issue was taken up a few years later by Saul and Kardar Saul and Kardar (1993, 1994) who claimed $`A=4`$; their work is based on exact computations of partition functions and thus does not suffer at low $`T`$ from the thermalization problems of the Monte Carlo approach. We are aware of no specific heat measurements in this model since. How could $`A`$ not be $`4`$? The subtlety is that we must take $`L\mathrm{}`$ at fixed $`T`$, and only *after* can we take $`T0`$; indeed Eq. (4) assumes $`L=\mathrm{}`$ whereas Eq. (3) assumes $`T0`$ at fixed $`L`$.

Using the algorithm in Galluccio et al. (2000, 2001), together with the availability of cheap and powerful computers, we have extended significantly the study of Saul and Kardar. For the sake of comparison, they had $`80`$ samples at $`L=20`$, 22, and 24, and 4 samples at $`L=32`$ and $`L=36`$. (They also had samples for $`L18`$.) We go much beyond that, both in lattice sizes and in the number of samples we consider. In the left part of Fig. 1 we show our first analysis of $`c_V`$ as follows.

Let’s set $`J=1`$. When $`T0`$, if naïve scaling ($`A=4`$) holds, $`\mathrm{ln}\left(T^2c_V\right)+4/T\text{const}`$, while $`\mathrm{ln}\left(T^2c_V\right)+4/T\frac{4A}{T}`$ if $`AA_{\text{naïve}}`$ and $`p=p_{\text{naïve}}=2`$. (The $`c_V`$ resulting from our exact partition function computations has been averaged over disorder samples.) In the plot we see that for any given lattice size, when $`T`$ becomes small enough there is a saturation toward the naïve scaling behavior, i.e., the points go to a constant value. The physically relevant regime is the thermodynamic limit, given by the envelope of these curves; this envelope does appear and seems to be linear in $`1/T`$. Note that the envelope emerges only on quite large lattices ($`L30`$); because of this, the true scaling escaped detection by Saul and Kardar. The straight line in the left part of Fig. 1 is our best linear fit to the $`L=50`$ data when $`\beta [2.5,5.5]`$. It is a very satisfying fit and gives $`A=2.02\pm 0.03`$, close to the integer value $`A=2`$.

We can also present the data in a slightly different fashion. In the right part of Fig. 1 we plot $`T\mathrm{ln}\left(T^2c_V\right)`$ versus $`T`$. Here the coefficient $`A`$ is given by the intercept of the envelope’s extrapolation to $`T=0`$, the left axis of the picture. We can distinguish three regions. The first region is for very low $`T`$ values. Here the naïve (non-thermodynamic scaling) with $`A=A_{\text{naïve}}=4`$ is very clear. This region, where the intercept at $`T=0`$ is $`4`$, shrinks to zero with increasing lattice size. In a second region we have the physical scaling; for the large lattice sizes we have, the value $`A2`$ emerges. Notice that this is the same region where in the left part of Fig. 1 the $`L=50`$ data lie on a straight line. The third and last region corresponds to “high” $`T`$ ($`T0.4`$) where one is far from the critical point and no scaling is apparent.

Our conclusion here is that thanks to the larger sizes available to us and to a technique that does not suffer from low temperature critical slowing down, the thermodynamic scaling of $`c_V`$ is now finally clarified.

## Ground state properties —

Our computations also give the ground state energies and degeneracies. Theoretical arguments Bouchaud et al. (2003) suggest that the mean ground state energy density $`e_0`$ has power corrections in $`1/L`$:

$$e_0(L)=e_0^{}+aL^{2+\mathrm{\Theta }^{(e)}}.$$
(5)
We have $`e_0^{}=1.4017(3)`$ which agrees well with previous work. We also find $`\mathrm{\Theta }^{(e)}=0.08(7)`$; note that the prediction in Bouchaud et al. (2003) is that $`\mathrm{\Theta }^{(e)}=\theta _{DW}`$, the exponent associated with domain wall energies. Following the work of Hartmann and Young Hartmann and Young (2001), there is general agreement that $`\theta _{DW}=0`$ in the $`2d`$ $`\pm J`$ model. Thus our estimate for $`\mathrm{\Theta }^{(e)}`$ is in excellent agreement with the conjecture in Bouchaud et al. (2003).

We have performed a similar study for the mean ground state entropy density $`s_0(L)`$. We find $`s_0^{}=0.0714(2)`$ which compares well with the recent work of Blackman et al. (1998) in which $`s_0^{}=0.0709(4)`$. The fit also gives $`\mathrm{\Theta }^{(s)}=0.42(2)`$, though if we take into account systematic effects we cannot rule out $`\mathrm{\Theta }^{(s)}=1/2`$. We believe that this large value, unrelated to $`\theta _{DW}`$, denotes the presence of a subtle organization of the ground states.

## Anomalous density of excitations —

The microcanonical entropy $`S(E)`$ of an energy level $`E`$ is defined as the logarithm of the number of spin configurations having exactly that energy. Clearly, $`S(E)`$ is obtained from the knowledge of $`P_J`$ as computed in Eq. (2). Of major interest is $`S_1S_0`$, the increase of entropy when going from the ground state energy $`E_0`$ to the lowest excitation energy. In the pure ferromagnetic model, the lowest excitation corresponds to taking the ground state (all spins parallel) and flipping a single spin. This gives $`S_1S_0=\mathrm{ln}\left(V\right)`$. One says that the excitations are “elementary”, and the system at low temperature is accurately described as a gas of independent excitations.

The situation changes dramatically when a large enough fraction of $`J_{ij}`$ are negative, taking the system from a ferromagnetic to a spin-glass phase. In such a phase, the large $`V`$ law of $`S_1S_0`$ is modified. This is illustrated in Fig. 2 where we plot our numerical estimate of $`\overline{S_1S_0}`$ as a function of $`\mathrm{ln}\left(V\right)`$. The dotted line is $`\mathrm{ln}\left(V\right)`$ while the dash-dotted one is $`2\mathrm{ln}\left(V\right)`$. We see that the true scaling behavior emerges only for large lattices and that the large $`V`$ behavior is compatible with a $`2\mathrm{ln}\left(V\right)`$ growth. How can one interpret this anomalous growth? Imagine classifying all excitations in terms of the size of the cluster of spins flipped when comparing to a given ground state configuration. (Naturally, one may ignore all clusters of zero excitation energy, and it is enough to focus on connected clusters.) Just as in the ferromagnetic case, some of the excitations correspond to single spin flips; there are $`O(V)`$ such objects. For any bounded-size cluster, the number of objects is $`O(V)`$, leading to $`S_1S_0\mathrm{ln}\left(V\right)`$. Since one has instead $`\overline{S_1S_0}2\mathrm{ln}\left(V\right)`$, finite-size clusters are irrelevant: necessarily large scale excitations dominate the set of excitations of lowest energy.

It is important to understand the nature of these large scale excitations, but unfortunately our computational approach does not give us configurations, it merely counts their number. There are other ways to gain insight into this problem. To begin, we consider as in Wang and Swendsen (1988) an analogy with the $`1d`$ pure Ising model.

In that system, when using periodic boundary conditions, the lowest excitation is composite, corresponding to a pair of kinks with a total energy $`4J`$; however the “true” elementary excitations are single kinks, necessarily absent when using periodic boundary conditions. It is easy to see that for this $`1d`$ model the quantity $`S_1S_0`$ grows as $`2\mathrm{ln}\left(V\right)`$, i.e., as in our $`2d`$ system.

How may objects of energy $`2J`$ appear in our 2d lattices with periodic boundary conditions for which the gap is $`4J`$? To answer this question, consider in a ground state configuration any connected cluster of spins and associate to its surface the corresponding closed path $`𝒫`$ on the dual lattice Schuster (1979). (The cluster’s surface is the set of edges connecting the cluster to its complement.) When flipping the cluster, the change in the configuration’s energy comes only from those bonds crossing $`𝒫`$; in fact, for each such bond that is satisfied ($`J_{ij}S_iS_j=1`$), the energy increases by $`2J`$, and otherwise it is decreased by $`2J`$. It is easy to see that all clusters lead to $`𝒫`$ with an even number of bonds, and thus excitation energies are quantized in units of $`4J`$. However there are closed paths that are *not* associated with clusters: an example is a path that winds around one of the directions of the lattice! Such topologically non-trivial paths are called domain walls; when comparing periodic and anti-periodic boundary conditions, the set of bonds in the ground state that are changed from satisfied to unsatisfied or vice-versa form exactly such a path. When $`L`$ is odd, domain walls have energies $`\pm 2J`$, $`\pm 6J`$, …and the the quantum $`2J`$ appears. Of course, to have a physical excitation, one needs to introduce domain walls in *pairs*; then the flipped cluster of spins is topologically a strip with a surface in two pieces, one for each domain wall, while its energy is a multiple of $`4J`$. Note that this is exactly what happened in the one dimensional case, the domain walls there being simply kinks which also arise in pairs.

To justify the anomalous scaling of $`S_1S_0`$, we appeal to the fact realized a few years ago by Hartmann and Young Hartmann and Young (2001) that low energy domain walls proliferate in the $`\pm J`$ spin glass. Let $`\delta S`$ be the typical entropy of a single domain wall; if we focus on those excitations of energy $`4J`$ associated with two domain walls of energy $`2J`$, we have an excess entropy $`\mathrm{\Delta }S=2\delta S`$. Our data thus suggest $`\delta S=\mathrm{ln}\left(L^2\right)`$; this law can be interpreted by saying that $`\delta S`$ is the sum of a $`\mathrm{ln}\left(L\right)`$ term coming from the $`L`$ possible mean transverse positions of the domain wall and of an additional $`\mathrm{ln}\left(L\right)`$ term coming from the degeneracy (proliferation) at a given position. To extend this reasoning to the case of $`L`$ even, we first remark that the domain walls there have energies $`0`$, $`4J`$, …To have a “strip” excitation, we need one domain wall of energy $`4J`$ and one of $`0`$ energy. Undoubtedly, the entropy of these domain walls increases with their energy; a simple pattern is obtained if we conjecture that the excess entropy increases by $`\mathrm{ln}\left(L\right)`$ every time the energy increases by the quantum $`2J`$. If this is so, the first domain wall will contribute $`3\mathrm{ln}\left(L\right)`$ to the excess entropy and the second $`\mathrm{ln}\left(L\right)`$, leading again to the desired $`2\mathrm{ln}\left(V\right)`$ result. Such a conjecture is quite elegant and should be amenable to testing using a recent Monte Carlo method Houdayer (2001).

## Finite-size scaling —

Given the result for $`S_1S_0`$, we go back to Eq. (3) to understand the finite-size scaling of $`c_V`$. When $`T0`$ and $`L\mathrm{}`$ simultaneously, standard arguments lead to

$$T^2c_V(L,T)e^{2\beta J}\left[L/\mathrm{\Lambda }(T)\right].$$
(6)
Here $`\mathrm{\Lambda }(T)`$ is a temperature dependent length that determines the cross-over between the thermodynamic scaling of $`c_V`$ (going as $`\mathrm{exp}\left[2\beta J\right]`$ when $`L=\mathrm{}`$) and the “naïve” scaling as in Eq. (3). $``$ is a finite-size scaling function; when its argument is large, $`L\mathrm{\Lambda }(T)`$, we recover the thermodynamic limit and thus necessarily $``$ must tend toward a constant. (Since $`c_V`$ is intensive, the $`L`$ dependence must drop out.) On the contrary, when $`L\mathrm{\Lambda }(T)`$, we recover the behavior of Eq. (3) where $`c_V`$ goes as $`\mathrm{exp}(4\beta J)`$ *but* diverges as $`L\mathrm{}`$. Interestingly, in this regime $`c_V`$ is not self-averaging and so one should apply finite-size scaling for the whole probability distribution of $`c_V`$. Just as before where we computed $`\overline{S_1S_0}`$ and not $`\mathrm{log}(\overline{\mathrm{exp}(S_1S_0)})`$:, we focus on the typical behavior and so we consider the median rather than the average of $`c_V`$ (this distinction is relevant only in the very low $`T`$, unphysical region, while it is irrelevant in our scaling region for $`T`$, say, close to $`0.3`$). When using this data we have a very reasonable data collapse, consistent with Eq. (6) as shown in Fig. 3. We find that this median scales as $`L^2\mathrm{exp}(2\beta J)`$ at low $`T`$ and thus $`(x)x^2`$ as $`x0`$. This then gives

$$\mathrm{\Lambda }(T)\mathrm{exp}(\beta J).$$
(7)
## Summary and discussion —

We have investigated the critical thermodynamics of the $`2d`$ Ising spin glass with binary couplings. Our main conclusion is that the specific heat density scales as $`c_V\mathrm{exp}(2J\beta )`$. This scaling is “anomalous” in the sense that it does not follow from the size of the energy gap (which is $`4J`$). To find this scaling law, it is necessary to go to rather large systems, $`L30`$. We also found that the typical degeneracy of the first excited level grows about $`L^4`$ times faster than that of the ground state level. We believe this high degeneracy has its roots in the proliferation of domain walls, two domain walls enabling one to define a composite excitation. Such a picture justifies the analogy with kink pairs proposed many years ago by Wang and Swendsen Wang and Swendsen (1988): each domain wall may indeed play the role of a kink, albeit with an additional entropy contribution. Finally, using finite-size scaling, we found a cross-over length scale $`\mathrm{\Lambda }(T)`$ that grows as $`\mathrm{exp}(J/T)`$; this divergence is exactly as expected from hyperscaling arguments.

This work was supported in part by the European Community’s Human Potential Programme under contracts HPRN-CT-2002-00307 (DYGLAGEMEM) and HPRN-CT-2002-00319 (STIPCO).


format: text

text:
# Equation of state and opacities for hydrogen atmospheres of magnetars

## 1. Introduction

Neutron star atmospheres differ from the atmospheres of ordinary stars because of the high gravity and magnetic fields. The majority of the known pulsars possess magnetic fields $`B10^{11}`$$`10^{13.5}`$ G (Manchester et al., 2003). Most popular models of the soft gamma repeaters and anomalous X-ray pulsars assume they are magnetars, neutron stars with superstrong magnetic fields $`B10^{14}`$$`10^{15}`$ G (e.g., Thompson 2003).

The properties of matter are strongly modified if the electron cyclotron energy $`\mathrm{}\omega _{\mathrm{ce}}=\mathrm{}eB/m_\mathrm{e}c`$ exceeds 1 a.u. – i.e., the field strength $`B`$ is higher than $`B_0=m_\mathrm{e}^2ce^3/\mathrm{}^3=2.3505\times 10^9`$ G (for a review, see Lai 2001). Shibanov et al. (1992) presented the first detailed model of hydrogen atmospheres for neutron stars with strong magnetic fields. Later it was developed and used by many authors (e.g., Zavlin et al., 1995; Pavlov et al., 1995; Page, Shibanov, & Zavlin, 1996; Zane et al., 2001; Ho & Lai, 2003; Özel, 2003; and references therein). These works have played a valuable role in assessing the observed spectra of the neutron stars. Despite the recognition of the possibility that the absorption by atoms can significantly contribute to the opacities, the atoms were not included in the atmosphere model in a thermodynamically consistent manner in the listed works.

Theoretical equation of state for partially ionized hydrogen plasma with strong magnetic fields was developed by Potekhin, Chabrier, & Shibanov (1999). The nonperturbative effects of atomic thermal motion across the magnetic field were taken into account and shown to be important. Based on this theory, in a previous paper (Potekhin & Chabrier, 2003, hereafter Paper I) we performed extensive calculations of the thermodynamic functions and calculated monochromatic and Rosseland mean opacities, taking into account the bound-bound and bound-free transitions. These calculations were done at $`11.9\mathrm{log}B[G]13.5`$. In addition, we reevaluated the free-free opacities, taking into account the motion of both interacting particles, electron and proton, in arbitrary magnetic fields.

Here we report an extension of the previous calculations of the equation of state and opacities of a partially ionized hydrogen plasma to higher field strengths, $`13.5\mathrm{log}B[G]15`$.

## 2. Theoretical model

A major complication when incorporating bound species in the models of strongly magnetized neutron-star atmospheres arises from the strong coupling between the center-of-mass motion of the atom and the internal atomic structure. If an atom rests without motion in a strong magnetic field, there are two distinct classes of its quantum states: at every value of the Landau quantum number $`n`$ and the magnetic quantum number $`s`$ ($`n0`$, $`sn`$), there is one tightly bound state, with binding energy growing asymptotically (at $`B\mathrm{}`$) as $`(\mathrm{ln}B)^2`$, and an infinite series of loosely bound states (numbered by $`\nu =0,1,\mathrm{}`$) with binding energies approaching the energies of a field-free H atom (e.g., Canuto & Ventura, 1977). The Landau number $`n`$ is not a “good” quantum number, in the sense that the Coulomb interaction mixes different Landau orbitals of a free electron, but this numbering is unambiguous and convenient at $`BB_0`$. The atom is elongated: its size along the magnetic field $`𝑩`$ either decreases logarithmically (for the tightly bound states) or remains nearly constant (for the loosely bound states), while the transverse radius decreases as $`B^{1/2}`$. The radiative transition rates are different for the three basic polarizations: the linear polarization along the field and the two circular polarizations in the transverse plane.

This simplicity is destroyed when atomic motion is taken into account. The electric field, induced in the comoving frame of reference, breaks down the cylindrical symmetry. Separating the center-of-mass motion and choosing an appropriate gauge of the vector potential (Gor’kov & Dzyaloshinskii, 1967; Vincke, Le Dourneuf, & Baye, 1992; Potekhin, 1994), one comes to an effective one-particle Schrödinger equation with an effective Hamiltonian depending on the pseudomomentum $`𝑲`$, which is a conserved quantity related to the center-of-mass motion in a magnetic field. If the component of $`𝑲`$ perpendicular to $`𝑩`$, $`K_{}`$, is finite, then not only $`n`$, but also $`s`$ is not a “good” quantum number; nevertheless, it is convenient to retain the quantum level numbering by $`n`$, $`s`$, and $`\nu `$, specified above. In the adiabatic approximation, widely used in the past (Gor’kov & Dzyaloshinskii, 1967; Ipatova, Maslov, & Subashiev, 1984), the transverse part of the wave function is postulated to be the same as for an electron without Coulomb interaction. We perform calculations without this approximation, using the technique described by Potekhin (1994) for the discrete atomic spectrum and by Potekhin & Pavlov (1997) for the continuum. This technique uses the basis of Landau functions, which represent the transverse parts of the wave functions of a free electron in a magnetic field and are the same for the solutions of the Schrödinger and Dirac equations. Therefore this technique can be used for the ordinary ($`\mathrm{}\omega _{\mathrm{ce}}<m_\mathrm{e}c^2`$) as well as superstrong fields ($`\mathrm{}\omega _{\mathrm{ce}}>m_\mathrm{e}c^2`$), provided the atomic binding energy $`m_\mathrm{e}c^2`$, which is always the case at $`B10^{15}`$ G.

An atom moving across the strong magnetic field acquires a constant dipole moment perpendicular to $`𝑩`$ and $`𝑲`$. Those radiative transitions, which were dipole-forbidden for an atom at rest because of conservation of the $`z`$-projection of the angular momentum, become allowed and should be taken into account in the atmosphere models. The binding energies decrease with increasing $`K_{}`$. Asymptotically, at large $`K_{}`$, the binding energies tend to $`e^3B/(cK_{})`$.

The equation of state for partially ionized hydrogen in strong magnetic fields was constructed and discussed by Potekhin et al. (1999). The treatment is based on the free energy minimization. The free energy model is in essence a generalization of the nonmagnetic model of Saumon & Chabrier (1991, 1992) to the case of a strong magnetic field. We consider a plasma composed of $`N_\mathrm{p}`$ protons, $`N_\mathrm{e}`$ electrons, $`N_\mathrm{H}`$ hydrogen atoms, and $`N_{\mathrm{mol}}`$ molecules in a volume $`V`$, the number densities being $`n_jN_j/V`$. The Helmholtz free energy is written as the sum

$$F=F_{\mathrm{id}}^\mathrm{e}+F_{\mathrm{id}}^\mathrm{p}+F_{\mathrm{id}}^{\mathrm{neu}}+F_{\mathrm{ex}}^\mathrm{C}+F_{\mathrm{ex}}^{\mathrm{neu}},$$
(1)
where $`F_{\mathrm{id}}^\mathrm{e}`$, $`F_{\mathrm{id}}^\mathrm{p}`$, and $`F_{\mathrm{id}}^{\mathrm{neu}}`$ are the free energies of ideal gases of the electrons, protons, and neutral species, respectively, $`F_{\mathrm{ex}}^\mathrm{C}`$ takes into account the Coulomb plasma nonideality, and $`F_{\mathrm{ex}}^{\mathrm{neu}}`$ is the nonideal contribution which arises from interactions of bound species with each other and with the electrons and protons. Ionization equilibrium is given by minimization of $`F`$ with respect to particle numbers under the stoichiometric constraints, keeping $`V`$ and the total number $`N_0`$ of protons (free and bound) constant. The latter number is determined by the total mass density $`\rho `$: $`n_0N_0/V\rho /m_\mathrm{H},`$ where $`m_\mathrm{H}=m_\mathrm{e}+m_\mathrm{p}`$. The formulae for each term in Eq. (1) are given in Potekhin et al. (1999) and in Paper I. The employed minimization technique is similar to the technique presented by Potekhin (1996) for the nonmagnetic case.

Once the free energy is obtained, its derivatives over $`\rho `$ and $`T`$ and their combinations provide the other thermodynamic functions.

The atomic number fraction $`x_\mathrm{H}=n_\mathrm{H}/n_0`$, evaluated in the course of the free energy minimization, can be used in calculations of atmospheric opacities. One should take into account that the strong magnetic field affects the polarization properties of radiation (e.g., Ginzburg, 1970). At photon energies $`\mathrm{}\omega `$ much higher than $`(4\pi \mathrm{}^2e^2n_\mathrm{e}/m_\mathrm{e})^{1/2}28.7\rho _0^{1/2}\mathrm{eV},`$ where $`\rho _0\rho /(1\text{ g cm}\text{-3})`$, radiation propagates in the form of two so-called *normal modes*. These modes have different polarization vectors $`𝒆_j`$ and different absorption and scattering coefficients, which depend on the angle $`\theta _B`$ between the propagation direction and $`𝑩`$ (e.g., Kaminker, Pavlov, & Shibanov, 1982). The two modes interact with each other via scattering. Gnedin & Pavlov (1973) formulated the radiative transfer problem in terms of these modes.

Let the magnetic field be directed along the $`z`$-axis. At a fixed photon frequency $`\omega `$, the absorption opacity $`\kappa _j^\mathrm{a}(\theta _B)`$ in each mode $`j`$ and scattering opacities $`\kappa _{jj^{}}^\mathrm{s}(\theta _B)`$ from mode $`j`$ into mode $`j^{}`$ can be presented as (e.g., Kaminker et al., 1982)

$`\kappa _j^\mathrm{a}(\theta _B)=m_\mathrm{H}^1{\displaystyle \underset{\alpha =1}{\overset{1}{}}}|e_{j,\alpha }(\theta _B)|^2\sigma _\alpha ^\mathrm{a},`$ (2)
$`\kappa _{jj^{}}^\mathrm{s}(\theta _B)={\displaystyle \frac{3}{4}}{\displaystyle \underset{\alpha =1}{\overset{1}{}}}|e_{j,\alpha }(\theta _B)|^2`$
$`\times {\displaystyle \frac{\sigma _\alpha ^\mathrm{s}}{m_\mathrm{H}}}{\displaystyle _0^\pi }|e_{j^{},\alpha }(\theta _B^{})|^2\mathrm{sin}\theta _B^{}\mathrm{d}\theta _B^{},`$ (3)
where $`\alpha =0,\pm 1`$, $`e_{j,0}=e_{j,z}`$ is the $`z`$-component of $`𝒆_j`$, and $`e_{j,\pm 1}=(e_{j,x}\pm \mathrm{i}e_{j,y})/\sqrt{2}`$ are the circular components. The cross sections $`\sigma _\alpha `$ depend on $`\omega `$, but not on $`j`$ or $`\theta _B`$.

The total scattering opacity from mode $`j`$ is $`\kappa _j^\mathrm{s}=\kappa _{j1}^\mathrm{s}+\kappa _{j2}^\mathrm{s}`$, and the total extinction opacity is $`\kappa _j=\kappa _j^\mathrm{a}+\kappa _j^\mathrm{s}`$. According to equations (2), (3), we can write

$$\kappa _j(\theta _B)=\underset{\alpha =1}{\overset{1}{}}|e_{j,\alpha }(\theta _B)|^2\widehat{\kappa }_\alpha ,$$
(4)
where $`\widehat{\kappa }_\alpha `$ ($`\alpha =1,0,1`$) do not depend on $`\theta _B`$. With a good accuracy, $`\widehat{\kappa }_\alpha (\sigma _\alpha ^\mathrm{a}+\sigma _\alpha ^\mathrm{s})/m_\mathrm{H}`$.

In a partially ionized atmosphere, the opacity is contributed by electrons, ions, and bound species. The scattering cross section includes contributions from the electrons and protons: $`\sigma _\alpha ^\mathrm{s}=\sigma _\alpha ^{\mathrm{s},\mathrm{e}}+\sigma _\alpha ^{\mathrm{s},\mathrm{p}}`$. With a good accuracy, $`\sigma _\alpha ^{\mathrm{s},\mathrm{e}}`$ and $`\sigma _\alpha ^{\mathrm{s},\mathrm{p}}`$ are described by simple analytical formulae (e.g., Pavlov et al. 1995). The absorption cross section $`\sigma _\alpha ^\mathrm{a}`$ includes contributions from absorption by plasma electrons and protons (mainly by free-free transitions due to the electron-proton collisions, $`\sigma _\alpha ^{\mathrm{ff}}`$), transitions between discrete states of an atom (bound-bound absorption, $`\sigma _\alpha ^{\mathrm{bb}}`$) and photoionization (bound-free, $`\sigma _\alpha ^{\mathrm{bf}}`$). Thus, for the hydrogen atmosphere, we can write

$$\sigma _\alpha ^\mathrm{a}x_\mathrm{H}(\sigma _\alpha ^{\mathrm{bb}}+\sigma _\alpha ^{\mathrm{bf}})+(1x_\mathrm{H})\sigma _\alpha ^{\mathrm{ff}}.$$
(5)
For the hydrogen atoms, electrons, and protons moving in a strong magnetic field, the cross sections $`\sigma _\alpha ^{\mathrm{bb}}`$ were studied by Pavlov & Potekhin (1995), $`\sigma _\alpha ^{\mathrm{bf}}`$ by Potekhin & Pavlov (1997), and $`\sigma _\alpha ^{\mathrm{ff}}`$ in Paper I.

Ho & Lai (2003) presented convenient formulae for $`|e_{j,\alpha }(\theta _B)|^2`$ for a fully ionized, strongly magnetized electron-ion plasma (taking into account vacuum polarization). However, the bound species affect the dielectric tensor of the medium and hence the polarization properties of the normal modes. For a completely neutral hydrogen gas in a strong magnetic field, polarization modes were studied by Bulik & Pavlov (1996), but for a partially ionized gas the problem is not yet solved. For the time being, we use the formulae of Ho & Lai (2003). Since the neutral fraction is typically small, we expect that the resulting error should not be large for the effective Rosseland mean opacities presented in the tables.

## 3. Results

### 3.1. Calculation of tables

Since our model of the ionization equilibrium and equation of state is computationally expensive, it is not possible to use our computer code “on line” in any practical application. The alternative is to tabulate thermodynamic quantities covering the density, temperature, and magnetic field domain of interest and to rely on an interpolation of the tabular values. In Paper I, we presented the tables which cover a range of $`\rho `$, $`T`$, and $`B`$ appropriate for most typical neutron stars, such as isolated radio pulsars. We used the fitting formulae for the binding energies, quantum-mechanical sizes, oscillator strengths, and electron collision widths derived by Potekhin (1998). Calculations at $`B10^{13.5}`$ G were hampered by the absence of the relevant fits.

In order to overcome this difficulty, in the present work we have calculated the tables of the required atomic quantities as functions of the transverse pseudomomentum $`K_{}`$, using the computer code described by Potekhin (1994). A representative grid of $`50`$ equally spaced values of $`\mathrm{log}K_{}`$ has been used for every relevant set of atomic quantum numbers and every considered field strength. Handling of these tables is easier at superstrong than at weaker fields, because the number of involved discrete states is smaller. For example, there remain no bound states with $`s=2`$ and $`s=1`$ at $`B>2\times 10^{13}`$ G and $`B>6\times 10^{13}`$ G, respectively. The latter point is illustrated by Fig. 1, which shows the energies of the most important atomic transitions at $`K_{}=0`$ as functions of $`B`$. With increasing $`B`$, the tightly bound levels $`s=3,2,`$ and 1 consecutively merge with the continuum and cease to contribute to the bound-bound spectrum; instead, they appear as autoionization resonances. There still remain excited loosely bound states (with $`s=0`$ and any $`\nu `$) at any $`B`$, but, at the typical neutron-star atmosphere densities, most of them are destroyed by interactions with surrounding particles.

We use a cubic interpolation across these precalculated tables in our computer code for obtaining the number fractions of atoms and molecules, thermodynamic functions, and cross sections $`\sigma _\alpha ^{\mathrm{bb}}`$. As previously, we employ the approximate expression for $`\sigma _\alpha ^{\mathrm{ff}}`$ \[Eq. (51) in Paper I\] and numerical tables of $`\sigma _\alpha ^{\mathrm{bf}}`$ as functions of $`K_{}`$ and $`\omega `$ (§4.2 of Paper I).

In this way we have calculated and tabulated thermodynamic functions and opacities of the hydrogen plasma at $`13.5\mathrm{log}B\text{ [G]}15.0`$ (preliminary results for $`B=10^{14}`$ G were presented by Chabrier, Douchin, & Potekhin, 2002). We have also updated the Fortran program for the cubic-polynomial interpolation of the tabulated data in the 3-parameter space of $`B`$, $`T`$, and $`\rho `$.111The tables and the program are available at http://www.ioffe.ru/astro/NSG/Hmagnet/.

As discussed by Potekhin et al. (1999), our model becomes less reliable at relatively low $`T`$ and high $`\rho `$, particularly because of formation of molecules and chains Hn, which are treated in an approximate manner. In this domain, the partial number fractions and thermodynamic quantities are strongly model-dependent. The stronger magnetic field, the higher the temperature below which the molecules and chains dominate (Lai, 2001). Trying to avoid the domain of uncertainty, we shifted the lower bound on the temperature in the tables, as compared to the case of weaker fields (Paper I), from $`2\times 10^5`$ K to $`5\times 10^5`$ K. On the other hand, the temperature at which the atoms become thermally ionized increases with increasing $`B`$, therefore we have shifted the upper temperature bound in the tables from $`10^7`$ K to $`4\times 10^7`$ K. It would be meaningless to consider still higher $`T`$, because hydrogen undergoes efficient thermonuclear burning at such temperatures (Ergma, 1986). The third input parameter in our tables, in addition to $`B`$ and $`T`$, is the “astrophysical density parameter” $`R=\rho _0/T_6^3`$, where $`T_6=T/10^6`$ K. In order to avoid the rather uncertain region of possible phase separation (see below), we restricted the public tables by $`\mathrm{log}R3.4`$, although higher $`R`$ values were also considered. This range of $`R`$ should be sufficient for modeling the atmospheres of most neutron stars, except the unusually cold ones mentioned below.

Along with the pressure, our tables contain internal energy, entropy, specific heat, and the logarithmic derivatives of pressure with respect to $`V`$ and $`T`$. Other second-order thermodynamic quantities can be calculated using the Maxwell relations (Landau & Lifshitz, 1993).

The format of the tables is the same as in Paper I. In addition to the thermodynamic quantities, the tables contain the number fractions of different chemical species, and also the effective Rosseland mean opacities for diffusion of radiation along and across the magnetic field.

### 3.2. Ionization equilibrium and equation of state

A strong magnetic field generally increases the fraction of bound species at a given $`T`$. Figure 2 illustrates this point. In this figure, ionization equilibrium curves are plotted as functions of density for different values of $`B`$ and $`T`$. At the lowest $`T`$ and highest $`B`$ shown in Fig. 2, the molecular H2 fraction (dot-dashed line on the left panel) becomes significant at $`\rho 1\text{ g cm}\text{-3}`$. At the same $`T`$ and $`B`$ and higher density $`\rho 100\text{ g cm}\text{-3}`$, the H2 molecules disappear, but the Hn chains and clusters (not shown in the figure) become abundant. As mentioned above, our results are not well justified under such conditions.

Figure 3 shows the domains of partial ionization in the $`T`$$`\rho `$ plane at three values of $`B`$. With increasing $`B`$, the domains where the atomic fraction is above a specified level expand significantly towards higher temperatures.

Figure 4 shows pressure $`P`$ as function of $`\rho `$ at $`T=10^{6.5}`$ K and at different field strengths. At very low $`\rho `$, we have nearly fully ionized, almost ideal nondegenerate gas of electrons and protons. The magnetic field does not affect the equation of state in this case. With increasing density, at $`10^2\text{ g cm}\text{-3}\rho 10^3\text{ g cm}\text{-3}`$, the electrons become degenerate at $`B=0`$, but remain nondegenerate at $`B10^{12}`$ G. This explains why the dashed line in the figure (pressure at $`B=0`$) is higher than the other lines in this regime. Meanwhile, because of the partial recombination of atoms, the solid curves ($`P`$ of the nonideal magnetized plasma) are slightly lower than the dotted ones ($`P`$ of the ideal fully ionized plasma). At still higher density, however, the atoms become pressure ionized: in this region the increase of $`P`$ due to the increased number of free electrons and protons competes with a negative nonideal contribution, which is mainly due to the Coulomb term $`F_{\mathrm{ex}}^\mathrm{C}`$ in the free energy.

At sufficiently high $`B`$ or low $`T`$, the Coulomb interaction leads to a discontinuity, exemplified in Fig. 4 by the line $`B=10^{15}`$ G. It occurs at $`T`$ below the critical temperature $`T_\mathrm{c}`$ and at $`\rho `$ around the critical density $`\rho _\mathrm{c}143B_{12}^{1.18}\text{ g cm}\text{-3}`$ (Potekhin et al., 1999), where $`B_{12}=B/10^{12}`$ G. This phase transition is directly related to the magnetic field assisted separation between the condensed hydrogen and vapor, discussed by Lai & Salpeter (1997). Note that the density of the condensed phase at zero pressure and zero temperature, estimated by Lai & Salpeter (1997) as $`\rho _\mathrm{s}560B_{12}^{6/5}\text{ g cm}\text{-3}`$, scales with $`B`$ similarly to $`\rho _\mathrm{c}`$. At $`B=10^{15}`$ G, this estimate gives $`\rho _\mathrm{s}2.2\times 10^6\text{ g cm}\text{-3}`$, whereas we have numerically $`\rho _\mathrm{s}1.7\times 10^6\text{ g cm}\text{-3}`$ in Fig. 4. However, in general, the validity of the free-energy models in the framework of the chemical picture of plasmas is questionable near the plasma phase transition. In order to study the matter properties near the phase transition, one needs to consider accurately the interaction of Hn molecules in the gas phase and their interaction with the surface of the condensed phase. This problem goes beyond the scope of the present paper, although it certainly deserves a future study. If this phase separation is real, it can be important for interpretation of the emission of neutron stars with unusually low $`T`$ and/or high $`B`$, which in this case do not possess an optically thick atmosphere above the condensed surface. This hypothesis was suggested by Lai & Salpeter (1997); a possible example of such “naked neutron star” has been recently considered by Turolla et al. (2003).

### 3.3. Opacities

Figure 5 shows the monochromatic extinction coefficients $`k_\alpha =\rho \widehat{\kappa }_\alpha `$ for the three basic polarizations $`\alpha =0,\pm 1`$ at $`B=10^{13.5}`$ G, $`T_6=1`$ or 4, and $`\rho =0.1\text{ g cm}\text{-3}`$ or 50 g cm-3. Figure 6 demonstrates the same functions for 10 times higher magnetic field and densities. These densities correspond to the outer and inner layers of the atmosphere, where the optical depth is of the order of 1 for ordinary (high-opacity) and extraordinary (low-opacity) modes of radiation, respectively (see, e.g., Ho & Lai 2003). At the higher densities, the opacities are featureless, except the proton cyclotron resonance at $`\mathrm{}\omega _{\mathrm{cp}}=\mathrm{}\omega _{\mathrm{ce}}m_\mathrm{e}/m_\mathrm{p}=6.305B_{12}`$ eV, because the plasma at such density is pressure-ionized. On the contrary, absorption features are clearly visible at the lower density.

At the lower field strength (Fig. 5), one sees the broad absorption feature for $`\alpha =\pm 1`$ (the left and central panels) due to the $`s=01`$ bound-bound transition at $`\mathrm{}\omega 0.3`$ keV. For $`\alpha =+1`$ (the left panel), this absorption merges with the strong proton cyclotron free-free resonance at lower energies and with the photoionization continuum at higher energies. This situation is rather special, because at this field strength $`\mathrm{}\omega _{\mathrm{cp}}`$ is rather close to the photoionization threshold, whereas the limiting energy for the bound-bound transition falls in between (see Fig. 1). Both the photoionization edge and the absorption line are considerably broadened because of the energy shifts due to the atomic motion (the magnetic broadening; see Pavlov & Potekhin 1995 and Potekhin & Pavlov, 1997), therefore they effectively merge together. On the other hand, the free-free absorption coefficient at the lower $`T`$ (the solid curve) reveals the bumps at $`\mathrm{}\omega 1`$ keV, which arise from the opening of ionization channels with higher $`s`$ values at higher energies, and from autoionization resonances just below these partial ionization thresholds (Potekhin & Pavlov, 1997). On the right panel ($`\alpha =0`$, the polarization vector along $`𝑩`$), one notices only the photoionization edge at relatively small $`\rho `$ and $`T`$ (the lower solid line).

At the higher field strength (Fig. 6), all transitions changing $`s`$ belong to the continuum (see Fig. 1) and form a series of resonances on the photoionization profile at the energies above $`\mathrm{}\omega _{\mathrm{cp}}`$, particularly visible for $`\alpha =+1`$ at low $`\rho `$ and $`T`$ (the lower solid line on the left panel). Near the bottom of the photosphere ($`\rho =500\text{ g cm}\text{-3}`$), the bound states are destroyed and do not contribute to the spectrum.

## 4. Conclusions and outlook

We have calculated the equation of state of fully and partially ionized hydrogen plasmas in a wide range of densities, temperatures, and magnetic fields expected for photospheres of the magnetars. We have also calculated monochromatic radiative opacities for three basic polarizations. The first- and second-order thermodynamic functions, non-ionized fractions, and effective Rosseland mean opacities are published in the electronic form.

These results can be useful for construction of the atmosphere models and calculation of the spectra of outgoing radiation from warm neutron stars with superstrong magnetic fields. This work is under way; some results have been published by Ho et al. (2003a, b). However, there remain unsolved problems which hamper an immediate application of our present results to modeling accurate and reliable outgoing spectra of magnetars with hydrogen atmospheres. First, vacuum polarization alters the dielectric tensor of the medium and polarization of photon modes (see Pavlov & Gnedin 1984), which is important for magnetar atmospheres. In particular, the “vacuum resonance” phenomenon may lead, under certain conditions, to an adiabatic mode conversion and “mode collapse” (Ho & Lai, 2003, and references therein). In the latter case, the description in terms of the two modes fails, and one has to solve more general equations for the evolution of electromagnetic waves (Lai & Ho, 2003). Second, as noted in §2, the presence of the bound species also changes the dielectric properties of the medium and the normal mode polarization, which necessitates derivation of a polarization tensor for partially ionized, strongly magnetized plasmas. Third, for relatively cold neutron stars with strong magnetic fields, the phase separation mentioned in §3.2 should be studied in a thermodynamically consistent way. We are planning to study these problems in near future.

We thank Dong Lai for useful discussions, and Yura Shibanov and Wynn Ho for valuable remarks. A.P. gratefully acknowledges hospitality of the theoretical astrophysics group at the Ecole Normale Supérieure de Lyon and the Astronomy Department of Cornell University. The work of A.P. is supported in part by RFBR grants 02-02-17668 and 03-07-90200.


text:
# Contribution of the VLBA Network to Geodynamics

## 1. Introduction

Method of VLBI first proposed by Matveenko, Kardashev and Sholomitsky (1965) is used for two major applications: making images of compact radio sources and measuring antenna’s motion with respect to extragalactic source. Fringe amplitudes and phases measured with precision $`10^1`$$`10^3`$ are used as primary observables for imaging; group delays measured with relative precision $`10^8`$$`10^{10}`$ are used as observables for geodynamics applications. Group delay is sensitive to the projection of the baseline vector to the source direction. Having observed many sources in different directions and different baselines we are able to solve for source positions, site positions, Earth orientation parameters (EOP) and the rate of their changes in a least square solution. The paper of Overs, Fanselow and Jacobs (1998) provides an overview of geodetic VLBI technique.

## 2. Precision and accuracy

There are several measures of how precise results are. Formal uncertainties, which are calculated according to a law of propagation of errors of an individual observation, are the first measure. This estimate gives us an upper limit of precision, since an estimation model usually does not account correlations correctly. Another measure of precision is repeatability — the root mean square deviation of the series of estimates of the targeted parameter with respect to the average. And the last measure is accuracy: the measure of how far our estimate of the targeted parameter is from the true value. Table 1 summarizes the current level of precision and accuracy of geodetic VLBI.

## 3. Scientific objectives

Although site positions and source coordinates are determined in processing VLBI observations, and they can be useful for particular applications, such as source mapping, determination of positions is not the objective per se. The objective of observations is to measure changes in site positions and source coordinates. Investigation of changes in positions of the sites allows us to make judgments about the forces which caused the motion and the response the Earth to external forcing.

When the whole VLBI dataset is analyzed, the following parameters are estimated:

* Secular site velocities. Causes: 1) plate tectonic; 2) po-seismic slips; 3) post-glacial isostatical adjustment; 4) environmental changes 5) ??
* amplitudes of harmonic site position variations. Causes: 1) solid Earth tides; 2) ocean tides; 3) ??
* anharmonic site position variations. Causes: 1) atmospheric pressure loading; 2) hydrology loading; 3) non-tidal ocean loading; 4) ??
* variations in the Earth orientation parameters. Causes: 1) torques exerted by the Sun and planets; 2) exchange of angular momentum with the atmosphere; 3) exchange of angular momentum with the ocean; 4) earthquakes? 5) ??

One of the examples of achieving these objectives are results of analysis of the amplitudes of the Earth’s nutation fitted to VLBI data.

In the case if the Earth responses to external torques as an absolutely rigid body, the complex transfer function would be (1,0). The presence of the liquid core causes the resonance at the retrograde free core nutation frequency. Dissipation in the ocean and in the liquid body causes appearance of the non-zero imaginary part of the transfer function. Fitting empirical transfer function to a model allowed Buffett, Mathews, & Herring, (2002) to estimate so-called basic Earth paramenters which includes the core mantle boundary flattening and the magnitude of the magnetic field at the core mantle boundary (0.7mT).

VLBI observations allowed us to detect a secondary tidal effect: site displacements caused by water mass redistribution due to ocean tides. Mass loading results in crust deformation even thousands kilometers away from the mass load. It was found that observed amplitudes are in a rather good agreement with models (fig. 2. Comparison of the observed amplitudes with theoretical amplitudes at many VLBI sites, including VLBA antennas, allowed Petrov and Ma (2003) to discriminate between competing models of ocean tides. It is relatively easy to validate ocean tides models close to the shore where tide gauges are located, but it is a challenge to validate ocean tides models at open ocean.

Similarly to mass loading, changes of atmospheric mass above the surface cause crust deformation due to pressure loading. Unlike tides, changes in atmospheric pressure are not predictable. However, global numerical weather models based on past observations of the world meteorological network allow us to compute displacements due to atmospheric pressure loading with latency several days. Examples of the atmospheric pressure loading are presented in fig. 3. Comparison of the admittance factors of the theoretical time series for all sites, including VLBA sites, showed surprisingly well agreement between the theoretical model and observations, even for horizontal loading which does not exceed 1 mm for most of the sites. This allowed Petrov and Boy (2003) to set the upper limit (4%) of possible errors of loading Love numbers which describe elastic properties of the Earth.

## 4. Contribution of VLBA

VLBA stations started to participate in geodetic program even before construction of the VLBA network has been completed. By June 2003 the total number of observations used in analysis was 4 032 371. The number of observations with at least one VLBA site was $`\mathit{894\hspace{0.17em}690}`$. Thus, approximately 1/5 of geodetic observations are done at the VLBA network.

Results of estimation of secular station velocity are presented in table 2. The station of MK-VLBA belongs to the Pacific tectonic plate, SC-VLBA belongs to the Caribbean plate. They move rapidly (in geophysical meaning of this word) with respect to the North American plate where 8 other stations are located. The station of OV-VLBA is located rather close to the plate boundary and this station, and in a lesser degree BR-VLBA, are affected by the plate deformation near the edge of the plate. Peculiar horizontal velocity of Pietown is not yet completely understood. As Gordon (2003) indicated, it may be related to antenna’s tilting. Six other stations do not exhibit relative horizontal motion with respect to each other at the level exceeding 1 mm during entire history of the VLBA.

Exceptionally high precision of the Earth orientation parameters makes 20 station experiments correlated at the VLBA correlator especially attractive for investigation of free core nutation – a poorly understood phenomenon of variable amplitude with the period around 430 days. Unlike forced nutation caused by torques exerted by Moon, Sun and planets, the mechanism of excitation of the free core nutation is not known. Currently it cannot be predicted and should be monitored like polar motion.

## References

Buffett, B.A., Mathews P.M., & Herring T.A., 2002, J. Geophys. Res., 107(B4), 10.1029/2000JB000056.

Gordon, D. 2003, this volume.

Matveenko, L.I., Kardashev N.S., & Sholomitsky G.B., 1965, Izvestija Vysshih Uchebnyh Zavedenij, Radiofizika, 8, 651.

Petrov, L., & Ma, C., 2003, J. Geophys. Res., 108(B4), 10.1029/2002JB001801.

Petrov, L., & Boy, J.-P., 2003, Submitted to J. Geophys. Res., 10.1029/2003JB002500.

Sovers, O.J., Fanselow, J.J., & Jacobs, C.S., 1998, Phys.Rev.D, 70(4), 1393.


text:
# A NEW METHOD TO CONSTRAIN THE IRON ABUNDANCE FROM COOLING DELAYS IN CORONAL LOOPS

## 1  Introduction

Coronal EUV emission is mainly produced by radiative decay of collisionally excited, highly-ionized iron ions, i.e. by Fe8+ and Fe9+ in the 171 Å TRACE passband ($`T1.0`$ MK) and by Fe11+ at 195 Å ($`T1.5`$ MK). Coronal loops undergo various phases of heating and cooling. When steady heating operates, the loops evolve into a steady-state, where heating input is balanced by thermal conduction losses to the chromosphere and radiative losses into space, as described by the energy equation derived by Rosner, Tucker, & Vaiana (1978), generalized for gravity and non-uniform heating by Serio et al. (1981). When heating stops, coronal loops cool off by thermal conduction and radiative losses. One consequence of this cooling process is that the EUV emission peaks first in a high-temperature filter 1, e.g. in TRACE 195 Å, and later in a lower-temperature filter 2, e.g. in TRACE 171 Å , with a time delay that we call cooling delay $`\mathrm{\Delta }t_{12}`$. This cooling delay was found to scale proportionally to the loop length $`L`$ in a recent study (Winebarger et al. 2003). Under the assumption that the cooling time is dominated by radiative cooling, we can relate this observable time delay to the radiative loss function $`\mathrm{\Lambda }(T)`$, which allows us to constrain the absolute abundances of iron ions, which dominates the radiative loss function in this coronal temperature range of $`T12`$ MK. The new method of determining iron abundances provides an important diagnostic for coronal heating mechanisms that involve preferential ion heating. This study demonstrates also that much more physical information can be inferred from the temporal evolution of the EUV intensity than from the intensity measured in a single image.

## 2  Theoretical Model

Plasma cooling through a narrow temperature range can be approximated by an exponential function,

$$T_e(t)=T_1\mathrm{exp}\left(\frac{t}{\tau _{cool}}\right),$$
(1)
where $`T_1`$ is the initial temperature at time $`t=0`$ and $`\tau _{cool}`$ is the cooling time. Cooling over large temperature ranges would require the full consideration of the hydrodynamic equations, but the exponential approximation is fully justified for the narrow temperature range of $`T_e1.01.4`$ MK we are considering here, during the cooling of coronal flare loops through the two TRACE 171 and 195 Å filters (see also measurements in Warren et al. 2003). So, when a cooling plasma reaches the temperature $`T_1`$ of the peak response of the hotter filter ($`T_1=1.4`$ MK for TRACE 195 Å), the time delay $`\mathrm{\Delta }t_{12}`$ to cool down to the cooler filter $`T_2`$ (e.g. $`T_2=0.96`$ MK for 171 TRACE Å), can be expressed with Eq.1 as

$$\mathrm{\Delta }t_{12}=t_2t_1=\tau _{cool}\mathrm{ln}\left(\frac{T_1}{T_2}\right).$$
(2)
The cooling time scale could be dominated by thermal conduction losses in the initial phase, but is always dominated by radiative losses in the later phase (e.g., Antiochos & Sturrock 1982). Here we make the a priori assumption that the cooling of EUV loops in the $`T_e1`$ MK temperature range is dominated by radiative losses. This working hypothesis is particularly justified near the almost isothermal loop tops and is also corroborated by the observational result found in Aschwanden et al. (2000), where a median value of $`\tau _{cool}/\tau _{rad}=1.02`$ was obtained from the statistics of 12 nanoflare loops observed with TRACE 171 and 195 Å. Moreover we show later that the radiative cooling time is always significantly shorter than the conductive cooling time for the cases analyzed here. Thus we set the cooling time $`\tau _{cool}`$ equal to the radiative cooling time $`\tau _{rad}`$,

$$\tau _{cool}\tau _{rad}=\frac{E_{th}}{dE_R/dt}=\frac{3n_ek_BT_e}{n_en_H\alpha _{FIP}\mathrm{\Lambda }(T_e)}$$
(3)
where $`n_e`$ is the electron density, $`n_H`$ the hydrogen density, $`T_e`$ the electron temperature, $`k_B=1.38\times 10^{16}`$ erg K-1 the Boltzmann constant, $`\alpha _{FIP}`$ the abundance enhancement factor for low first-ionization potential elements (at $`<10`$ keV), and $`\mathrm{\Lambda }(T_e)`$ the radiative loss function, which can be approximated with a constant in the limited temperature range of $`T_e0.52.0`$ MK, according to the piece-wise powerlaw approximation of Rosner et al. (1978), $`\mathrm{\Lambda }(T)\mathrm{\Lambda }_0=10^{21.94}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^3`$, for T=0.5-2.0 MK. The computation of the radiative loss function at a given temperature depends on the elemental atomic abundances, and thus we define a reference value $`\mathrm{\Lambda }_{0,Ph}`$ for photospheric abundances. Coronal abundances generally show a density enhancement for low first-ionization potential (FIP) elements, which we express with an enhancement factor $`\alpha _{FIP}`$. Since iron (Fe) contribution strongly at these temperatures, the radiative loss function in the corona scales as $`\mathrm{\Lambda }_{0,cor}=\mathrm{\Lambda }_{0,ch}\times \alpha _{Fe}`$. The cooling delay $`\mathrm{\Delta }t_{12}`$ as function of the coronal iron abundance $`\alpha _{Fe}`$ is thus (with Eqs. 1-3), assuming full ionization in the corona ($`n_H=n_e`$),

$$\mathrm{\Delta }t_{12}=\frac{3k_BT_e}{n_e\alpha _{Fe}\mathrm{\Lambda }_{0,Ph}}\mathrm{ln}\left(\frac{T_1}{T_2}\right).$$
(4)
When a loop cools through a passband, the maximum of the flux $`F(t)`$ is detected at the time when the loop temperature matches the peak of the response function, so the peak flux $`F_2`$ of the light curve in the lower filter corresponds to the emission measure $`EM`$ at the filter temperature $`T_2`$,

$$F_2=EM\times R_2=(n_e^2wq_w)\times R_2$$
(5)
with the flux $`F_2`$ in units of DN/(pixel s), $`w`$ is the loop width or diameter, $`q_w`$ is the linear filling factor in case of unresolved substructures, and $`R_2`$ is the response function, which is $`R_2=0.37\times 10^{26}`$ cm5 DN/(pixel s) for 171 Å and photospheric abundances (see Appendix A in Aschwanden et al. 2000). The value $`R_2`$ of the response function refers to the time at the beginning of the mission, while the degradation decreased this value by a factor of 0.78 in August 1999, the latest date of analyzed observations used here. The corresponding correction by a factor of $`\sqrt{0.78}=0.88`$ is neglected in the numerical values given in Table 1. Inserting the density from Eq.(5) into Eq.(4) we find the following expression for the iron abundance $`\alpha _{Fe}`$,

$$\alpha _{Fe}=\frac{3k_BT_2}{\mathrm{\Lambda }_{0,Ph}\mathrm{\Delta }t_{12}}\sqrt{\frac{R_2wq_w}{F_2}}\mathrm{ln}\left(\frac{T_1}{T_2}\right)=4.17\left(\frac{wq_w}{1\mathrm{Mm}}\right)^{1/2}\left(\frac{F_2}{10\mathrm{DN}/\mathrm{s}}\right)^{1/2}\left(\frac{\mathrm{\Delta }t_{12}}{1\mathrm{min}}\right)^1.$$
(6)
For a filling factor of unity ($`q_w=1`$), the iron enhancement factor can be determined with an accurcay of about $`\stackrel{<}{}20`$%, because the observables $`w`$, $`F_2`$, and $`\mathrm{\Delta }t_{12}`$ can each be measured better than $`\stackrel{<}{}10`$%. In case of unresolved fine structure, i.e. filling factors of $`q_w<1`$, we obtain with Eq.(6) an upper limit for the iron enhancement.

Recent EUV observations with TRACE have shown that the cooling delay $`\mathrm{\Delta }t_{12}`$ is roughly proportional to the loop length $`L`$. In order to understand such a correlation we use the energy balance equation, which is valid in a steady-state, e.g., before the cooling process, at the turning point from dominant heating to dominant cooling, or at the turning point from dominant conductive cooling to radiative cooling. The resulting scaling law is according to Rosner et al. (1978),

$$T_{max}1400(p_0L)^{1/3}\times q_{Serio}$$
(7)
with $`p_0`$ the pressure and $`L`$ the loop half length. This scaling law has been generalized for gravity and non-uniform heating by Serio et al. (1981), modified by the correction factor

$$q_{Serio}=\mathrm{exp}\left(0.08\frac{L}{s_H}0.04\frac{L}{s_p}\right),$$
(8)
where $`s_H`$ is the heating scale length and $`s_p=47,000\times T_{MK}`$ km the pressure scale height. With the ideal gas law ($`p=2n_ek_BT_{max}=p_0q_p`$), corrected for the pressure at the loop top \[$`q_p=\mathrm{exp}(h/s_p)=\mathrm{exp}(2L/\pi s_p)`$\], we can eliminate the pressure $`p`$ in the RTV scaling law and find the following expression for the density $`n_e`$,

$$n_e=\frac{T_{max}^2}{2k_BL(1400q_{Serio})^3q_p}.$$
(9)
Inserting this density into the relation for the cooling delay (Eq. 4) we find indeed a proportional relation $`\mathrm{\Delta }t_{12}L`$,

$$\mathrm{\Delta }t_{12}=L\times \left[\frac{6(1400q_{Serio})^3q_pk_B^2}{T_{max}\mathrm{\Lambda }_{0,Ph}\alpha _{Fe}}\mathrm{ln}\left(\frac{T_1}{T_2}\right)\right].$$
(10)
which should show up for cooling loops with similar maximum temperatures $`T_{max}`$. The relation $`\tau _{cool}L`$ was also derived in Cargill et al. (1995; Eq. 14E therein) and Serio et al. (1991; Eq. 14 therein).

## 3  Data Analysis

We are using three data sets for which the cooling delay between the TRACE 195 Å and 171 Å filters has been measured: 11 nanoflare loops analyzed in Aschwanden et al. (2000), 4 medium-sized EUV loops analyzed in Schrijver (2002), and 5 large EUV loops analyzed in Winebarger at al. (2003).

In the first study (Aschwanden et al. 2000), nanoflare loops were measured on 1999 Febr 17, 02:15-03:00 UT, with a cadence of $`2`$ min in both 171 Å and 195 Å with TRACE. The time delay between the appearance in the two filters was measured by cross-correlation of the two time profiles, shown in Fig.6 and listed in Table 1 of Aschwanden et al. (2000). The time delay was found to be positive in 11 out of 12 cases. We use only the 11 cases with positive time delay, sorted according to the loop length in Table 1. The loop length was determined from a geometric model of a projected semi-circular cylindrical loop. The electron density $`n_e`$ is measured from the flux $`F_2`$ and loop width $`w`$ according to Eq.(5).

In the second study (Schrijver 2001), time delays between the maximum intensity in the 171 and 195 Å TRACE passbands were measured in 4 active region loops observed above the limb on 2000 May 26. The loop half length was estimated with $`L=(\pi /2)h+h_{limb}`$ based on the height $`h`$ of the loop top above the limb and a height correction $`h_{limb}5`$ Mm for the offset between optical limb and the portion of the plage inside the limb. The fluxes $`F_{171}`$ and loop widths $`w`$ could not be reliably measured for this subset, because of confusion problems in the crowded limb regions.

In the third study (Winebarger et al. 2002), 5 loops were measured on 1999 Aug 18, 1998 Jul 04, 1998 Jul 25, 1998 Aug 17, and 1998 Jul 25, in both the 171 Å and 195 Å filters with TRACE. The time delay between the two filters was measured from the peak times of asymmetric gaussian curves fitted to the light curves in the two filters. The loop lengths were measured from the best fit of elliptical and dipolar geometric models to the projected loop shapes.

Inserting the measured values $`F_2`$, $`w`$, and $`\mathrm{\Delta }t_{12}`$ into Eq.(6) yields the iron abundance enhancement factors $`\alpha _{Fe}`$ listed in Table 1, which have a mean and standard deviation of $`\alpha _{Fe}4.9\pm 1.7`$ for the dataset of Aschwanden et al. (2000), and $`\alpha _{Fe}1.4\pm 0.4`$ for the dataset of Winebarger et al. (2003), respectively. Note that the mean iron enhancement is significantly higher for the small-scale nanoflares analyzed in Aschwanden et al. (2000) than for the large-scale loops of Winebarger et al. (2003). The variable degree of iron enhancement could be related to different physical conditions in freshly-filled small-scale loops compared with longer-lived large-scale loops. The effect of gravitational settling has been observed in observations of coronal streamers with SUMER (Feldman et al. 1999) and with UVCS (Raymond et al. 1997).

The values of $`\alpha _{Fe}`$ have to be considered as upper limits if the filling factor $`q_w1`$ is lower than unity. We can turn the argument around and assume that the iron abundance enhancement has to be larger than or equal unity, which would then constitute lower limits for the filling factors: $`q_w0.23\pm 0.08`$ (Aschwanden et al. (2000) and $`q_w0.78\pm 0.22`$ (Winebarger et al. (2003).

Figure 1 shows the correlation plot between the cooling delay $`\mathrm{\Delta }t_{12}`$ and the loop half length $`L`$. A linear regression fit between the logarithic values yields the power-law relation $`\mathrm{\Delta }t_{12}L^{1.08\pm 0.16}`$ which is fully consistent with the theoretical prediction of a linear relationship $`\mathrm{\Delta }t_{12}L^1`$ (Eq. 10).

In order to verify our initial assumption of dominant radiative loss, we estimate also the conductive cooling time,

$$\tau _{cond}=\frac{3n_ek_BT_{max}}{F_C}=1.1\times 10^9n_eT_{max}^{5/2}L^2,$$
(11)
where we assign a mean value of $`T_{max}1.2`$ MK for the looptop temperature, when the loop cools through the two TRACE passbands. The last two columns in Table 1 show that the conductive cooling time is always much larger than the radiative cooling time, $`\tau _{cond}\tau _{rad}`$, which corroborates our a priori assumption of dominant radiative cooling.

## 4  Discussion

We have developed a simple method to constrain the iron abundance or filling factor in a coronal loop, based on the cooling delay measurement between two EUV filters. This model predicts a linear relationship between the cooling delay and the loop length, i.e. $`\mathrm{\Delta }t_{12}L`$, which is consistent with the observed relation, i.e. $`\mathrm{\Delta }t_{12}L^{1.08\pm 0.16}`$.

The only underlying assumption is that the cooling is dominated by radiative loss in the temperature range of EUV loops ($`T1.01.4`$ MK here), rather than by conductive loss. Theoretical models of cooling in flare loops predict that conductive cooling is only dominant in the initial phase of very hot plasma seen in soft X-rays, say at $`T\stackrel{>}{}10`$ MK, while the later cooling phase seen in EUV is dominanted by radiative cooling (Antiochos 1980; Antiochos & Sturrock 1982; Cargill et al. 1995). Our measurements of the density allows us to estimate upper limits for the radiative cooling time. If there would be a filling factor $`q_w<1`$, the density would be higher and the radiative cooling time shorter. But even for a filling factor of unity, we find that the radiative cooling time in the EUV temperature range is significantly shorter than the conductive cooling time, which corroborates our a priori assumption. This is also consistent with other observations of EUV loops, where the ratio of the cooling time to the radiative cooling time was found to be $`\tau _{cool}/\tau _{rad}=1.02`$ (Aschwanden et al. 2000).

The main result of this study is the estimation of the iron abundance. We show a compilation of radiative loss functions in Fig. 2, which has been calculated for photospheric abundances (Meyer 1985), with an absolute iron abundance of $`log(A_{Fe})=7.59`$ relative to hydrogen $`log(A_H)=12.0`$ (i.e. $`A_{Fe}/A_H=3.9\times 10^5`$), as well as for coronal abundances (Feldman 1992), which have an iron enhancement by a factor of $`\alpha _{Fe}=3`$. This enhancement factor of $`\alpha _{Fe}=3`$ in density produces a change of the radiative loss rate that can be seen between the two curves calculated by Martens et al. (2000) in Fig. 2, at a temperature of $`T1.0`$ MK. Recent measurements of the absolute abundance of iron based on comparisons or EUV and radio data yielded a value of $`A_{Fe}/A_H=1.56\times 10^4`$, or a coronal iron enhancement by a factor of $`\alpha _{Fe}=4.0`$ (White et al. 2000). Our measurements from 16 different loops in many different active regions yield a median value of $`\alpha _{Fe}=4.0`$, or a mean and standard deviation of $`\alpha _{Fe}=3.7\pm 2.2`$. Because of this reasonable agreement of iron enhancements with radio methods (White et al. 2000) and spectroscopic measurements (Feldman 1992), our result corroborates the notion that low-FIP elements such as Fe are enhanced in the corona relative to photospheric values (Feldman 1992). By the same token we can argue the filling factor is close to unity in coronal EUV loops, otherwise we would have disagreement with spectroscopic and radio iron abundance measurements.

Acknowledgements: Part of this work was supported by NASA contracts NAS5-38099 (TRACE) and NAS8-00119 (SXT).

## References

Antiochos,S.K. 1980, ApJ 241, 385

Antiochos,S.K. & Sturrock,P.A. 1982, ApJ 254, 343.

Aschwanden,M.J., Tarbell,T.T., Nightingale,W., Schrijver,C.J., Title,A., Kankelborg,C.C., Martens,P., and Warren,H.P. 2000, ApJ 535, 1047

Aschwanden,M.J. and Schrijver,C.J. 2002, ApJS, 142, 269

Cargill,P.J., Mariska,J.T., and Antiochos,S.K. 1995, ApJ 439, 1034

Cook,J.W., Cheng,C.C., Jacobs,V.L., and Antiochos,S.K. 1989, ApJ 338, 1176

Feldman,U. 1992, Physica Scripta 46, 202

Feldman,U., Doschek,G.A., Schühle,U., and Wilhelm,K. 1999, ApJ 518, 500.

Golub,L. and Pasachoff,J.M. 1997, The Solar Corona, Cambridge: Cambridge University Press

Martens,P.C.H., Kankelborg,C.C., and Berger,T.E. 2000, ApJ 537, 471

Meyer,J.-P. 1985, ApJS 57, 173

Raymond,J.C. and 26 co-authors 1997, Solar Phys. 170, 105,

Rosner,R., Tucker,W.H., and Vaiana,G.S. 1978, ApJ 220, 643

Schrijver,C.J. 2001, Solar Phys. 198, 325

Schrijver,C.J. and Zwaan,C. 2000, Solar and stellar magnetic activity, 2000, Cambridge: Cambridge University Press

Serio,S., Peres,G., Vaiana,G.S., Golub,L., and Rosner,R. 1978, ApJ 243, 288

Serio,S., Reale,F., Jakimiec,J., Sylwester,B., and Sylwester,J. 1991, AA 241, 197.

Warren,H.P., Winebarger,A.R., and Mariska,J.T. 2003, ApJ, (subm.)

White,S.M., Thomas,R.J., Brosius,J.W., and Kundu,M.R. 2000, ApJ 534, L203

Winebarger,A.R., Warren,H.P., and Seaton,D.B. 2003, ApJ, (subm.)


text:
# The partial-fractions method for counting solutions to integral linear systems

## 1. Euler’s generating function

We are interested in computing the number of integer solutions of the linear system

$$𝐱_0^d,𝐀𝐱=𝐛,$$

where $`𝐀`$ is a nonnegative $`(m\times d)`$-integral matrix and $`𝐛^m`$. We think of $`𝐀`$ as fixed and study the number of solutions $`\varphi _𝐀(𝐛)`$ as a function of $`𝐛`$. (Strictly speaking, this function should only be defined for those $`𝐛`$ which lie in the nonnegative linear span of the columns of $`𝐀`$.)

The function $`\varphi _𝐀(𝐛)`$, often called a *vector partition function*, appears in a wealth of mathematical areas and beyond: Number Theory (partitions), Discrete Geometry (polyhedra), Commutative Algebra (Hilbert series), Algebraic Geometry (toric varieties), Representation Theory (tensor product multiplicities), Optimization (integer programming), as well as applications to Chemistry, Biology, Physics, Computer Science, and Economics.

Denote the columns of $`𝐀`$ by $`𝐜_1,\mathrm{},𝐜_d`$. The following lemma goes back to at least Euler :

###### Lemma 1 (Euler).

$`\varphi _𝐀(𝐛)`$ equals the coefficient of $`𝐳^𝐛:=z_1^{b_1}\mathrm{}z_m^{b_m}`$ of the function

$$f(𝐳)=\frac{1}{\left(1𝐳^{𝐜_1}\right)\mathrm{}\left(1𝐳^{𝐜_d}\right)}$$

expanded as a power series centered at $`𝐳=0`$.

###### Proof.

Expand each factor of the right-hand side into a geometric series. ∎

Equivalently, the coefficient of $`𝐳^𝐛`$ in $`f(𝐳)`$ equals the constant term in $`\frac{f(𝐳)}{𝐳^𝐛}`$, denoted by $`\mathrm{const}\frac{f(𝐳)}{𝐳^𝐛}`$. So Euler’s Lemma can be conveniently stated as

$$\varphi _𝐀(𝐛)=\mathrm{const}\frac{1}{\left(1𝐳^{𝐜_1}\right)\mathrm{}\left(1𝐳^{𝐜_d}\right)𝐳^𝐛}.$$

In a series of articles , we used complex integration of $`\frac{f(𝐳)}{𝐳^𝐛}`$ to compute $`\varphi _𝐀(𝐛)`$ for special cases of $`𝐀`$. Similar techniques were applied in . Here we expand $`\frac{f(𝐳)}{𝐳^𝐛}`$ into partial fractions to compute its constant term, and hence $`\varphi _𝐀(𝐛)`$. This work constitutes in a sense a refinement of the complex-integration methods, with the advantage that it is more flexible and—more importantly—applicable for *any* integral linear system.

## 2. Vector partition functions

The nature of the counting function $`\varphi _𝐀(𝐛)`$ is given by the following theorem. A *quasi-polynomial* is a finite sum $`Q(𝐛)=_𝐧c_𝐧(𝐛)𝐛^𝐧`$ with coefficients $`c_𝐧`$ that are functions of $`𝐛`$ which are periodic in every component of $`𝐛`$. The *degree* of $`Q`$ is the degree of the largest power $`𝐛^𝐧`$ appearing in $`Q`$. A matrix is *unimodular* if every square submatrix has determinant $`\pm 1`$.

###### Theorem 2 (Sturmfels ).

The function $`\varphi _𝐀(𝐛)`$ is a piecewise-defined quasi-polynomial in $`𝐛`$ of degree $`d\mathrm{rank}(𝐀)`$. The regions of $`^m`$ in which $`\varphi _𝐀(𝐛)`$ is a single quasi-polynomial are polyhedral, that is, they are defined by linear constraints. If $`𝐀`$ is unimodular then $`\varphi _𝐀`$ is a piecewise-defined polynomial.

The unimodular case of this theorem is due to Dahmen and Micchelli \[11, Corollary 3.1\].

The computation of both $`\varphi _𝐀(𝐛)`$ and the chamber complex consisting of the regions of quasi-polynomiality give rise to challenging problems. The most powerful technique for computing $`\varphi _𝐀(𝐛)`$ which we are aware of is due to Brion and Vergne . (The methods described here are much more elementary.) The chamber complex is still much of a mystery. A promising approach can be found in .

We finish this section with a *reciprocity theorem*. Let $`\varphi _𝐀^{}(𝐛)`$ count the integer solutions of

$$𝐱>0,𝐀𝐱=𝐛.$$

Both $`\varphi _𝐀`$ and $`\varphi _𝐀^{}`$ are quasi-polynomials, and can hence be algebraically defined for arguments which are not integer vectors in the positive span of $`𝐀`$. The following identity shows the close relationship between the two functions.

###### Theorem 3 ().

The quasi-polynomials $`\varphi _𝐀`$ and $`\varphi _𝐀^{}`$ satisfy

$$\varphi _𝐀(𝐛)=(1)^{d\mathrm{rank}𝐀}\varphi _𝐀^{}(𝐛).$$

This identity gives rise to a symmetry property of $`\varphi _𝐀(𝐛)`$. Let $`r_k`$ denote the sum of the entries in the $`k^{\text{th}}`$ row of $`𝐀`$, and let $`𝐫=(r_1,\mathrm{},r_m)`$. Then the integer solutions of

$$𝐱>0,𝐀𝐱=𝐛$$

are in bijection (via $`x_kx_k1`$) with the integer solutions of

$$𝐱0,𝐀𝐱=𝐛𝐫,$$

and hence $`\varphi _𝐀^{}(𝐛)=\varphi _𝐀(𝐛𝐫)`$. This yields:

###### Corollary 4.

The quasi-polynomial $`\varphi _𝐀`$ satisfies

$$\varphi _𝐀(𝐛)=(1)^{d\mathrm{rank}𝐀}\varphi _𝐀(𝐛𝐫).$$

## 3. The partial-fractions method

This section describes the idea behind our computations. Recall that our goal is to derive

$$\varphi _𝐀(𝐛)=\mathrm{const}\frac{1}{\left(1𝐳^{𝐜_1}\right)\mathrm{}\left(1𝐳^{𝐜_d}\right)𝐳^𝐛}.$$

We start by expanding

$$\frac{1}{\left(1𝐳^{𝐜_1}\right)\mathrm{}\left(1𝐳^{𝐜_d}\right)𝐳^𝐛}.$$

into partial fractions in one of the components of $`𝐳`$, say $`z_1`$:

$$\frac{1}{\left(1𝐳^{𝐜_1}\right)\mathrm{}\left(1𝐳^{𝐜_d}\right)𝐳^𝐛}=\frac{1}{z_2^{b_2}\mathrm{}z_m^{b_m}}\underset{k=1}{\overset{d}{}}\frac{A_k(𝐳,b_1)}{1𝐳^{𝐜_k}}+\underset{j=1}{\overset{b_1}{}}\frac{B_j(𝐳)}{z_1^j}.$$

Here $`A_k`$ and $`B_j`$ are polynomials in $`z_1`$, rational functions in $`z_2,\mathrm{},z_m`$, and exponential in $`b_1`$. (We are tacitly assuming that there are no multiple poles (besides $`𝐳=0`$), which generally holds unless $`m=1`$, a case which can be handled easily.) The two sums on the right-hand side correspond to the analytic and the meromorphic part with respect to $`z_1=0`$. The latter does not contribute to the $`z_1`$-constant term, whence

$`\varphi _𝐀(𝐛)`$ $`=\mathrm{const}_{z_2,\mathrm{},z_m}\left({\displaystyle \frac{1}{z_2^{b_2}\mathrm{}z_m^{b_m}}}\mathrm{const}_{z_1}\left({\displaystyle \underset{k=1}{\overset{d}{}}}{\displaystyle \frac{A_k(𝐳,b_1)}{1𝐳^{𝐜_k}}}\right)\right)`$
(1) $`=\mathrm{const}\left({\displaystyle \frac{1}{z_2^{b_2}\mathrm{}z_m^{b_m}}}{\displaystyle \underset{k=1}{\overset{d}{}}}A_k(0,z_2,\mathrm{},z_m,b_1)\right).`$
The effect of one partial fraction expansion is to eliminate one of the variables of the generating function, at the cost of replacing one rational function by a sum of such. The following idea is now evident.

###### Algorithm.

Apply (1) repeatedly to eliminate $`z_1`$, then $`z_2`$, etc., up to $`z_{m1}`$.

The constant term of the remaining rational functions in the one variable $`z_m`$ can be computed with the methods introduced in . This implies, in particular, that any vector partition function is a quasi-polynomial whose nontrivial ingredients are *Fourier-Dedekind sums* . The constant-term computation for the last variable also explains the quasi-polynomial character of $`\varphi _𝐀`$, since in this last step roots of unity will appear, with the components of $`𝐛`$ as exponents.

We note that this algorithm computes $`\varphi _𝐀(𝐛)`$ *as a function of $`𝐛`$*, that is, it allows symbolic computation. Secondly, it is not very hard to deduce Sturmfels’s Theorem 2 from this algorithm. What might be more important, however, is the fact that the constraints which define the regions of quasi-polynomiality of $`\varphi _𝐀`$ are obtained “on the go” as one computes $`\varphi _𝐀`$: When expanding into partial fractions, one has to check where the poles of a rational function are. The components of $`𝐛`$ will appear (linearly) in the exponents of these rational functions, and hence one will automatically have to split the computation into cases which give rise to different regions in the chamber complex of quasi-polynomiality.

While we hope that our algorithm gives a computational tool, in particular, for obtaining the chamber complex defining the regions of quasi-polynomiality of a vector partition function, it is not clear to us why the chambers obtained from the partial fractions analysis coincide with Sturmfels’s predicted chambers.

Our algorithm is best illustrated by going through an actual example.

## 4. An illustrating example

Let $`𝐀=\left(\begin{array}{cccc}1& 2& 1& 0\\ 1& 1& 0& 1\end{array}\right)`$, and write $`𝐛=(a,b)`$, so $`\varphi _𝐀(a,b)`$ counts the integer solutions of

$$x_1,x_2,x_3,x_40,\begin{array}{ccc}\hfill x_1+2x_2+x_3& =& a\hfill \\ \hfill x_1+x_2+x_4& =& b.\hfill \end{array}$$

By Euler’s Lemma 1,

$$\varphi _𝐀(a,b)=\mathrm{const}\frac{1}{(1zw)(1z^2w)(1z)(1w)z^aw^b}.$$

We first expand into partial fractions with respect to $`w`$:

$$\frac{1}{(1zw)(1z^2w)(1w)w^b}=\frac{\frac{z^{b+1}}{(1z)^2}}{1zw}+\frac{\frac{z^{2b+3}}{(1z)(1z^2)}}{1z^2w}+\frac{\frac{1}{(1z)(1z^2)}}{1w}+\underset{k=1}{\overset{b}{}}\frac{\mathrm{}}{w^k}.$$

Taking constant terms gives

$`\varphi _𝐀(a,b)`$ $`=\mathrm{const}_z\left({\displaystyle \frac{1}{(1z)z^a}}\mathrm{const}_w\left({\displaystyle \frac{1}{(1zw)(1z^2w)(1w)w^b}}\right)\right)`$
$`=\mathrm{const}\left({\displaystyle \frac{1}{(1z)z^a}}\left({\displaystyle \frac{z^{b+1}}{(1z)^2}}+{\displaystyle \frac{z^{2b+3}}{(1z)(1z^2)}}+{\displaystyle \frac{1}{(1z)(1z^2)}}\right)\right)`$
(2) $`=\mathrm{const}\left({\displaystyle \frac{z^{ba+1}}{(1z)^3}}+{\displaystyle \frac{z^{2ba+3}}{(1z)^2(1z^2)}}+{\displaystyle \frac{1}{(1z)^2(1z^2)z^a}}\right)`$
(At this point, we could interpret each of the three constant terms as counting integer solutions to new linear systems. For example, the last term gives the number of integer points $`(x,y)0`$ satisfying $`2x+ya`$. To keep a general flavor, we continue with our general algorithm.)

We compute the constant term of each of the three terms separately. For starters,

$$\mathrm{const}\frac{z^{ba+1}}{(1z)^3}=0$$

if $`ba+1>0`$, equivalently (since $`a`$ and $`b`$ are integers) $`ba`$. If $`b<a`$, we use

$$\frac{1}{(1z)^3}=\underset{k0}{}\left(\genfrac{}{}{0pt}{}{k+2}{2}\right)z^k,$$

which gives

$$\mathrm{const}\frac{1}{(1z)^3z^{ab1}}=\left(\genfrac{}{}{0pt}{}{ab+1}{2}\right)=\frac{(ab)^2}{2}+\frac{ab}{2}.$$

Incidentally, this identity is true not only for $`b<a`$, but also for $`b=a`$ and $`b=a+1`$, because the right-hand side vanishes then. This suggests that the regions of quasi-polynomiality *overlap*, as was proved by Szenes and Vergne . Hence for the first constant term, we obtain

$$\mathrm{const}\frac{1}{(1z)^3z^{ab1}}=\{\begin{array}{cc}0\hfill & \text{ if }ba,\hfill \\ \frac{(ab)^2}{2}+\frac{ab}{2}\hfill & \text{ if }ba+1.\hfill \end{array}$$

For the second term in (2),

$$\mathrm{const}\frac{z^{2ba+3}}{(1z)^2(1z^2)}=0$$

if $`2ba+20`$. If $`a2b+3`$ we expand into partial fractions again:

$`\mathrm{const}{\displaystyle \frac{z^{2ba+3}}{(1z)^2(1z^2)}}`$ $`=\mathrm{const}\left({\displaystyle \frac{1/2}{(1z)^3}}+{\displaystyle \frac{\frac{a2b3}{2}+\frac{1}{4}}{(1z)^2}}+{\displaystyle \frac{\frac{(a2b3)^2}{4}+\frac{a2b3}{2}+\frac{1}{8}}{1z}}+{\displaystyle \frac{(1)^{a+1}/8}{1+z}}\right)`$
$`={\displaystyle \frac{(a2b)^2}{4}}+{\displaystyle \frac{2ba}{2}}+{\displaystyle \frac{1+(1)^{a+1}}{8}}.`$
Similar to the first constant-term computation, this identity is also valid for $`a=2b+2,2b+1,2b`$, whence

$$\mathrm{const}\frac{z^{2ba+3}}{(1z)^2(1z^2)}=\{\begin{array}{cc}0\hfill & \text{ if }a2b+2,\hfill \\ \frac{(a2b)^2}{4}+\frac{2ba}{2}+\frac{1+(1)^{a+1}}{8}\hfill & \text{ if }a2b.\hfill \end{array}$$

The computation for the last term in (2) is almost identical. (Note that this term always contributes.)

$`\mathrm{const}{\displaystyle \frac{1}{(1z)^2(1z^2)z^a}}`$ $`=\mathrm{const}\left({\displaystyle \frac{1/2}{(1z)^3}}+{\displaystyle \frac{\frac{a}{2}+\frac{1}{4}}{(1z)^2}}+{\displaystyle \frac{\frac{a^2}{4}+\frac{a}{2}+\frac{1}{8}}{1z}}+{\displaystyle \frac{(1)^a/8}{1+z}}\right)`$
$`={\displaystyle \frac{a^2}{4}}+a+{\displaystyle \frac{7+(1)^a}{8}}.`$
Summing up all terms in (2) gives finally

$`\varphi _𝐀(a,b)`$ $`=\{\begin{array}{cc}\frac{a^2}{4}+a+\frac{7+(1)^a}{8}\hfill & \text{ if }ab,\hfill \\ \frac{(ab)^2}{2}\frac{ab}{2}+\frac{a^2}{4}+a+\frac{7+(1)^a}{8}\hfill & \text{ if }\frac{a}{2}1ba+1,\hfill \\ \frac{(ab)^2}{2}\frac{ab}{2}+\frac{(a2b)^2}{4}+\frac{2ba}{2}+\frac{1+(1)^{a+1}}{8}+\frac{a^2}{4}+a+\frac{7+(1)^a}{8}\hfill & \text{ if }b\frac{a}{2}\hfill \end{array}`$
$`=\{\begin{array}{cc}\frac{a^2}{4}+a+\frac{7+(1)^a}{8}\hfill & \text{ if }ab,\hfill \\ ab\frac{a^2}{4}\frac{b^2}{2}+\frac{a+b}{2}+\frac{7+(1)^a}{8}\hfill & \text{ if }\frac{a}{2}1ba+1,\hfill \\ \frac{b^2}{2}+\frac{3b}{2}+1\hfill & \text{ if }b\frac{a}{2}.\hfill \end{array}`$
It is a fun exercise to show that $`\varphi _𝐀(a,b)=\varphi _𝐀(a4,b3)`$, as promised by Corollary 4.

We pause for a moment to introduce the geometry behind our linear system. By thinking of $`x_3`$ and $`x_4`$ as slack variables, we can see that

$$x_1,x_2,x_3,x_40,\begin{array}{ccc}\hfill x_1+2x_2+x_3& =& a\hfill \\ \hfill x_1+x_2+x_4& =& b.\hfill \end{array}$$

is equivalent to

$$x_1,x_20,\begin{array}{ccc}\hfill x_1+2x_2& & a\hfill \\ \hfill x_1+x_2& & b.\hfill \end{array}$$

Depending on the relationship between $`a`$ and $`b`$, the geometric figure described by this linear system is a quadrilateral or triangle.

The inequalities which define the different cases for $`\varphi _𝐀(a,b)`$ determine the chamber complex in the parameter space, and in this example also the combinatorial type of the polytope: Except for boundary cases, in the first ($`ab`$) and last ($`b\frac{a}{2}`$) case, it is a triangle, and in the second case ($`a>b>\frac{a}{2}`$) a quadrilateral. In fact, the geometry of the triangle in the first case only depends on $`a`$, which is reflected in $`\varphi _𝐀(a,b)`$, and similarly the last case shows only dependency on $`b`$. In this last case, the triangle actually has integer vertices (the linear system is unimodular), whence we get a polynomial. In the first two cases the vertices of the polytope are half-integral, which is reflected in the period-2 quasi-polynomials. The following picture illustrates the three overlapping chambers.

## 5. Ehrhart quasi-polynomials

A *convex polytope* $`𝒫`$ in $`^d`$ is the convex hull of finitely many points in $`^d`$. Alternatively (and this correspondence is nontrivial ), one can define $`𝒫`$ as the bounded intersection of affine halfspaces. A polytope is *rational* if all of its vertices have rational coordinates. (A *vertex* of $`𝒫`$ is a point $`𝐯𝒫`$ for which there is a hyperplane $`H`$ such that $`\{v\}=𝒫H`$.) We denote by $`𝒫^{}`$ the relative interior of $`𝒫`$. For a positive integer $`t`$, let $`L_𝒫(t)`$ denote the number of integer points (“lattice points”) in the dilated polytope $`t𝒫=\{tx:x𝒫\}`$. The fundamental result about the structure of $`L_𝒫`$ is as follows.

###### Theorem 5 (Ehrhart ).

If $`𝒫`$ is a convex rational polytope, then the functions $`L_𝒫(t)`$ and $`L_𝒫^{}(t)`$ are quasi-polynomials in $`t`$ whose degree is the dimension of $`𝒫`$. If $`𝒫`$ has integer vertices, then $`L_𝒫`$ and $`L_𝒫^{}`$ are polynomials.

Ehrhart conjectured and partially proved the following *reciprocity law*, which was proved by Macdonald .

###### Theorem 6 (Ehrhart-Macdonald).

The quasi-polynomials $`L_𝒫`$ and $`L_𝒫^{}`$ satisfy

$$L_𝒫(t)=(1)^{dim𝒫}L_𝒫^{}(t).$$

The computation of Ehrhart quasi-polynomials is only slightly easier than that of $`\varphi _𝐀`$. Recent work includes .

Suppose the convex polytope $`𝒫^d`$ is given by an intersection of halfspaces, that is,

$$𝒫=\{𝐱^d:𝐀𝐱𝐛\},$$

for some $`(m\times d)`$-matrix $`𝐀`$ and $`m`$-dimensional vector $`𝐛`$. We may convert these inequalities into equalities by introducing slack variables. If $`𝒫`$ has rational vertices, we can choose $`𝐀`$ and $`𝐛`$ in such a way that all their entries are integers, without loss of generality nonnegative ones. In summary, we may assume that a convex rational polytope $`𝒫`$ is given by

$$𝒫=\{𝐱_0^d:𝐀𝐱=𝐛\},$$

where $`𝐀M_{m\times d}(_0)`$ and $`𝐛_0^m`$. (If we are interested in counting the integer points in $`𝒫`$, we may assume that $`𝒫`$ is in the nonnegative orthant, i.e., the points in $`𝒫`$ have nonnegative coordinates, as translation by an integer vector does not change the lattice-point count.)

The connection to vector partition functions is now evident. Since $`t𝒫=\{𝐱_0^d:𝐀𝐱=t𝐛\}`$, we obtain $`L_𝒫(t)=\varphi _𝐀(t𝐛)`$ as a special evaluation of $`\varphi _𝐀`$. Note that $`t𝐛,t=1,2,\mathrm{}`$ lie in the same chamber of quasi-polynomiality of $`\varphi _𝐀`$. Ehrhart’s Theorem 5 is therefore a special case of Sturmfels’s Theorem 2, and Theorem 6 is a special case of Theorem 3.

As an example, the quadrilateral $`𝒬`$ described by

$$x,y0,\begin{array}{ccc}\hfill x+2y& & 5\hfill \\ \hfill x+y& & 4\hfill \end{array}$$

(a special case of the polygons appearing in Section 4) with vertices $`(0,0),(4,0),(3,1),(0,5/2)`$ has the Ehrhart quasi-polynomial

$$L_𝒬(t)=\varphi _𝐀(5t,4t)=\frac{23}{4}t^2+\frac{9}{2}t+\frac{7+(1)^t}{8}.$$

Ehrhart quasi-polynomials are easier to compute, since one does not need to derive/know the chamber complex of quasi-polynomiality of $`\varphi _𝐀`$. Our algorithm, naturally, works just as well for $`L_𝒫`$.

## 6. Concluding remarks

Many open problems and questions remain untouched. Most importantly, we hope that our ideas will give rise to practical implementations of computing vector partition functions (including the chamber complex defining the regions of quasi-polynomiality) in their various disguises in form of a computer software. Again, the partial-fraction method is more general and more flexible than algorithms using complex residues, for example in .

Not unrelated is the question of computational complexity. It is known that the rational generating function of the Ehrhart quasi-polynomial of a $`d`$-polytope can be computed in polynomial time if $`d`$ is fixed . We have not analyzed the complexity of our algorithm. This seems interesting in the light that our algorithm depends on the number of contraints, geometrically corresponding to the facets (codimension-1 faces) of the polytope. Barvinok’s algorithm, in contrast, depends on the number of vertices of the polytope.

Another venue which should be explored is the following: After each step in our algorithm, one could reinterpret the constant term of each summand as counting integer solutions of a linear system. It might be interesting to simplify each of these linear systems (as far as this is possible without changing the number of integer solutions), and then proceed with the algorithm.


text:
# Reply to Lindgren and Salomonson

## Abstract

In the accompanying Comment \[Phys. Rev. A 67, 056501 (2003)\], I. Lindgren and S. Salomonson claim to prove for the Kohn-Sham kinetic energy functional of ground state electron density that a Frëchet functional derivative exists, equivalent to a multiplicative local potential function. If true, this result would imply an exact Thomas-Fermi theory for ground states of noninteracting electrons. However, such a theory is not consistent with the exclusion principle for more than one electron of each spin. The simplest counterexample is the lowest triplet state of a noninteracting two-electron atom. If only the total electron density were normalized, as in Thomas-Fermi theory, the lowest state would collapse into a doubly-occupied $`1s`$ spin-orbital. Two independent parameters $`ϵ_{1s}`$ and $`ϵ_{2s}`$ are required to maintain independent subshell normalization. The argument presented by these authors is discussed in the light of this unphysical implication.

For Phys.Rev.A

The underlying issue here is whether Hohenberg-Kohn theory implies an exact ab initio Thomas-Fermi theory (TFT) for ground states. For noninteracting electrons, the question reduces to the existence of a Fréchet functional derivative of the Kohn-Sham kinetic energy functional, equivalent to a multiplicative local potential function. If such a functional derivative did exist, it could be expressed as

$`{\displaystyle \frac{\delta T}{\delta \rho }}=\mu v(𝐫),`$ (1)
which simply states the TFT Euler-Lagrange equation. Constant $`\mu `$ is determined by the normalization condition $`\rho =N`$ for N electrons. In the context of density functional theory (DFT), $`\rho =_in_i\rho _i=_in_i\varphi _i^{}\varphi _i`$, a sum of orbital subshell densities. Density variations generated by unconstrained independent variations in the orbital Hilbert space are necessary and sufficient to determine the Euler-Lagrange equations. Spin indices and sums are assumed but suppressed in the notation here. For noninteracting electrons, the correct Euler-Lagrange equations for Fermi-Dirac electrons, as used in Kohn-Sham theory, are a set of Schrödinger equations,

$`\widehat{t}\varphi _i(𝐫)=\{ϵ_iv(𝐫)\}\varphi _i(𝐫),`$ (2)
where $`\widehat{t}=\frac{1}{2}^2`$, for orthogonal but independently normalized orbital functions $`\varphi _i(𝐫)`$. Independent orbital normalization requires independent parameters $`ϵ_i`$, and orbital orthogonality implies the exclusion principle. The restriction to one free parameter in Eq.(1) means that the exclusion principle cannot be enforced unless all energy eigenvalues $`ϵ_i`$ are equal. Hence, in a compact system, these equations cannot be reconciled for more than one electron of each spin. The simplest specific example is the lowest $`1s2s^3S`$ state of a two-electron system. TFT implies bosonic condensation with two electrons in the $`1s`$ spin-orbital. The argument of Lindgren and Salomonson (LS) must fail, because it implies that TFT is an exact theory for noninteracting electrons. The proposed argument also disagrees with DFT, in which the Kohn-Sham equations for noninteracting electrons are just Eqs.(2), valid for atomic shell structure and consistent with the exclusion principle. The TFT and DFT Euler-Lagrange equations are inconsistent in the $`1s2s^3S`$ example, because two independent constants cannot be deduced from one. This conclusion has been disputed and defended in recent publications.

LS put forward two variational derivations that set out to determine the density functional derivative of the Kohn-Sham kinetic energy functional $`T_s[\rho ]`$, constructed as the ground-state limit of the Schrödinger orbital functional $`T[\{\varphi _i\}]=_in_i(i|\widehat{t}|i)`$. In Section IIA, they consider an N-electron wave function $`\mathrm{\Psi }`$. For noninteracting electrons, the Schrödinger equation is separable, and the wave function can be expressed as a serial product of orbital functions, $`\mathrm{\Psi }=\mathrm{\Pi }_i\varphi _i(𝐫_𝐢)`$. Fermi-Dirac statistics are not implied, but must be imposed by explicit orthogonalization and by independent normalization of the partial densities $`\rho _i=\varphi _i^{}\varphi _i`$. In terms of partial densities, LS Eq.(18) should be $`\delta T=_in_i(ϵ_iv)\delta \rho _i`$. Instead, LS extend the functional derivative of Eq.(22) to unnormalized densities by making the implicit assumption that the undetermined constant has no orbital subshell structure, which is inconsistent with the orbital Schrödinger equations. This point can be examined in detail by considering specific wave functions, as in LS Section IIB.

In Section IIB, LS propose a variant of standard Hartree-Fock theory in which the product of orbital overlap integrals is retained, in contrast to the usual simplified argument based on the orbital functional $`T`$. If done correctly, this variant leads to standard Hartree-Fock equations, equivalent for noninteracting electrons to Kohn-Sham equations, in which the kinetic energy operator $`\widehat{t}`$ takes its usual form. The problem here is that all terms in the total energy expression for a Slater determinant must be considered together. To enforce orbital orthonormality, as required by the definition of a Slater determinant, a matrix of Lagrange multipliers $`ϵ_{ij}`$ is required. It is customary (for closed-shell systems) to assume a diagonalized canonical representation, leaving the diagonal elements $`ϵ_i`$ to be determined as eigenvalues of the canonical Hartree-Fock equations. These are Lagrange multipliers for independent normalization constraints $`(i|i)=1`$. Derivations using the simplified orbital functional $`T`$ are standard in applications of DFT and are consistent with Hartree-Fock theory. Products of the overlap factors are unity for variational solutions and can consistently be omitted. If retained in the derivation, overlap variations are multiplied by a factor that vanishes. An appropriate variational functional for noninteracting electrons is

$`{\displaystyle \underset{i}{}}n_i(i|\widehat{t}+v(𝐫)ϵ_i|i)\mathrm{\Pi }_{ji}(j|j).`$ (3)
This is derived directly for a serial product wave function, with Lagrange multipliers for independent orbital normalization. Orbital orthogonality is implied by the Euler-Lagrange equations or by symmetry. The first factor in each term here vanishes in stationary states, so all variations of the overlap product drop out, and variations of the first factor are multiplied by unity. Written out in full, the correct differential for variations about a stationary state is

$`{\displaystyle \underset{i}{}}n_i`$ $`[{\displaystyle d^3𝐫\delta \varphi _i^{}(𝐫)\{\widehat{t}+v(𝐫)ϵ_i\}\varphi _i(𝐫)\mathrm{\Pi }_{ji}(j|j)}+cc]`$ (4)
$`+{\displaystyle \underset{i}{}}n_i`$ $`(i|\widehat{t}+v(𝐫)ϵ_i|i)\delta \mathrm{\Pi }_{ji}(j|j).`$ (5)
This vanishes for a stationary state, but cannot do so unless all Lagrange multipliers $`ϵ_i`$ are included, precluding an equivalent Thomas-Fermi theory.

The argument presented by LS fails to consider the full variational expression including $`(i|v(𝐫)ϵ_i|i)`$ for each occupied orbital. LS Eq.(33) should be written as

$`\delta T_s`$ $`=`$ $`\delta \varphi _1|ϵ_1v(𝐫)|\varphi _12|2`$ (6)
$`+`$ $`\delta \varphi _2|ϵ_2v(𝐫)|\varphi _21|1+cc,`$ (7)
plus terms involving $`\delta (i|i)`$ that drop out of Eqs.(4). A correct derivation must take into account that the normalization product modifies all terms of Eq.(3), and modifies the usual DFT definition $`\rho =_in_i\rho _i`$. The functional derivative must be consistent with the Hartree-Fock equations implied for Fermi-Dirac electrons. Failure to include a separate Lagrange multiplier for normalization of each occupied orbital is inconsistent with the constraints needed to specify a Slater determinant as the variational trial function.

The chain rule connecting orbital and density functional derivatives of the functional $`T_s[\rho ]`$ is a consistency condition implied by the Kohn-Sham construction of $`T_s`$ using the Schrödinger orbital functional $`T`$. In their discussion of the chain rule appropriate to DFT, LS invoke an expression, LS Eq.(36), that is undefined unless a Fréchet functional derivative exists. The integral in LS Eq.(36) is inappropriate, since $`\rho `$ is an explicit function of $`\varphi _i^{}`$. LS Eq.(38) is inconsistent with the orbital Euler-Lagrange equations. Such logic can be avoided because the ground-state functional, extended to include density variations unrestricted by normalization, defines Gâteaux derivatives such that

$`{\displaystyle \frac{\delta T_s}{n_i\delta \rho _i}}\varphi _i={\displaystyle \frac{\delta T_s}{n_i\delta \varphi _i^{}}}=\widehat{t}\varphi _i=\{ϵ_iv\}\varphi _i.`$ (8)
A Fréchet derivative is implied only if all $`ϵ_i`$ are equal. Alternatively, since these Gâteaux derivatives are operationally equivalent to the Schrödinger operator $`\widehat{t}`$, one can define a generalized functional derivative as such a linear operator. Using this generalized definition, the Kohn-Sham equations can be derived entirely in terms of density functional derivatives.

LS misstate the mathematical implications of the fact that the noninteracting Kohn-Sham or Schrödinger equations imply by direct construction the existence of Gâteaux functional derivatives of the extended functional $`T_s`$. Unless all $`ϵ_i`$ are equal, this implies that a Frechét functional derivative does not exist. If a Fréchet derivative is simply assumed, in order to define an extended functional of total density, the theory is inconsistent with Fermi-Dirac statistics. As a consequence of the exclusion principle, there is no physically correct definition of such an extended functional. However, because Gâteaux derivatives exist, a unique extended functional of the subshell densities $`\rho _i`$ is defined. This suffices to determine the Euler-Lagrange equations of the theory. For practical applications of DFT, which must assume a parametrized functional $`E_{xc}[\rho ]`$, the implied universal functional $`F_s`$ has the same properties as $`T_s`$. An extended functional $`F_s[\rho ]`$ and the Fréchet derivative required by TFT do not exist in general for more than two electrons, but a unique extended functional of the partial densities $`\rho _i`$ does exist, as do the Gâteaux derivatives required to determine Euler-Lagrange equations. This generalizes TFT to a theory consistent with Fermi-Dirac statistics and with electronic shell structure.