text:
# Calculation of the energy spectrum of a two-electron spherical quantum dot

## I Introduction

Progress in nanotechnology has allowed the development of small devices like quantum dots. The confinement potential can be safely approximated by a harmonic one , which has boosted the study of quantum dots with parabolic confinement during the last years (see e.g. Refs. 3 and 4, and references therein). The presence of many interacting electrons render the computation of the electronic states and properties a very complicated many-body problem. The first non-trivial exactly solvable problem of many-electrons is the one of two electrons confined in a parabolic potential, which made it a very attractive workbench for testing all kind of approximations (see e.g. ).

The Hamiltonian describing the pair of interacting electrons in a parabolic quantum dot in the effective masss approximation can be written as

$$\widehat{H}=\frac{\mathrm{}^2}{2m^{}}(_1^2+_2^2)+\frac{1}{2}m^{}\omega ^2(r_1^2+r_2^2)+\frac{e^2}{ϵ|\stackrel{}{r}_1\stackrel{}{r}_2|}$$
(1)
where $`m^{}`$ is the effective mass, $`\omega `$ the confinement strength, and $`ϵ`$ the dielectric constant of the host material and $`^2`$ is the Laplacian operator.

Equation (1) can be separated into a centre-of-mass and a relative motion component. Furthermore, due to the radial symmetry of the components, only those parts of the corresponding Schrödinger equations have to be solved. The centre-of-mass part will give the well-known harmonic oscillator problem. In three dimensions, the relative motion part may admit exact solutions for special choices of the parameters (see e.g. ). For the two dimensional case, similar separation and substitution can be made, and again there are no general solutions expressible in terms of special functions. Nevertheless, in Ref.  and more recently in Ref.  it was shown that there exist analytic solutions for special choices of the confinement constant.

Many-body effects are usually divided into exchange and correlation components . Exchange-only effects are considered in Hartree-Fock (HF) and differently, in Kohn-Sham (KS) approaches (although in KS, the correlation effects can be included), which typically amounts for around $`10\%`$ of the total energy. Correlation is about one order of magnitude smaller. Nevertheless, it has been shown that in two dimensions for two electrons in a harmonic field, correlation may play a bigger role specially for singlet states (see e.g. ).

The shifted-1/N method  has been applied to various condensed matter problems. Also, the two-dimensional relative motion Schrödinger problem have been solved using this technique in Refs. .

The article has been structured as follows: in Sect II we describe the Schrödinger, Hartree-Fock, and Kohn-Sham approaches we use. In Sect. III we solve the exact, HF, and KS-LDA equations using an accurate numerical technique and the shifted-1/N method. We also apply perturbation theory up to first order in both limits of confinement for the exact and mean field cases, and we propose some interpolation formulae for the energy. We discuss the accuracy of the mean field approaches, and of the shifted-1/N method for the present case.

## II Method

Throughout the paper the units of energy will be given in terms of the effective Rydberg constant $`^{}=\mathrm{}^2/(2m^{}a^2)`$ and the effective Bohr radius $`a^{}=\mathrm{}^2ϵ/m^{}e^2,`$ respectively.

In centre-of-mass and relative coordinates and measuring in reduced units the Hamiltonian reads

$$\widehat{H}=(\frac{1}{4}_\stackrel{}{R}^2+_\stackrel{}{r}^2)+\gamma ^2R^2+\frac{\gamma ^2}{4}r^2+\frac{1}{r}$$
(2)
where we have chosen the centre-of-mass $`\stackrel{}{R}=(\stackrel{}{r}_1+\stackrel{}{r}_2)/2`$ and the relative coordinate $`\stackrel{}{r}=(\stackrel{}{r}_1\stackrel{}{r}_2)`$.

The separation leads to the harmonic oscillator problem for the centre-of-mass coordinate, with energy

$$E_{NL}=\gamma (2N+L+3/2).$$
(3)
and eigenfunctions

$$\mathrm{\Psi }_{NLM}(\stackrel{}{R})=𝒩_{NL}\mathrm{exp}(\gamma R^2)(2\gamma )^{L/2}R^LL_N^{L+1/2}(2\gamma R^2)Y_{LM}(\theta _R,\varphi _R).$$
(4)
For the relative coordinate equation the wavefunction can be separated into radial and angular components

$$\mathrm{\Psi }(\stackrel{}{r})=\frac{u(r)}{r}Y_{lm}(\theta ,\varphi ),$$
(5)
where $`Y_{lm}(\theta ,\varphi )`$ are the spherical harmonics which are eigenfunctions of the angular momentum operator $`L_z`$, and $`L^2`$ with eigenvalues $`m`$ and $`l`$. This separation makes the corresponding radial Schrödinger equation the following second order ordinary differential equation

$$\left[\frac{d^2}{dr^2}+l(l+1)\frac{1}{r^2}+\frac{1}{4}\gamma ^2r^2+\frac{1}{r}E\right]u(r)=0.$$
(6)
It is well known that exact solutions of equation (6) cannot be expressed in a closed form in terms of special functions. There are analytic expressions for the energy for particular values of $`\gamma `$ and $`l`$ as it was pointed out in Refs. , among others.

The electrons should satisfy the Fermi-Dirac statistics which means in this case that for singlet states ($`s=0`$) the spatial part of the wavefunction should be antisymmetric and for triplet states ($`s=1`$) symmetric. As the centre-of-mass coordinate remains the same after exchanging to electrons, the antisymmetry requirement will be in the relative part. Because of the separation in radial and angular components of the relative-coordinate wavefunction it will mean that singlet states are associated with odd $`l`$ and triplet states with $`l`$ even, respectively.

It is interesting to compare the results of exact calculations with independent-electron models like Hartree-Fock (HF), and Kohn-Sham (KS)  in order to assess the relative importance of many-body effects like exchange and correlation, and also to evaluate the performance of the local-density approximation. For two paired electrons the electronic density is $`\rho _{HF}=2|\varphi _{HF}|^2`$, where $`\varphi _{HF}`$ is the orbital, and the exchange potential is equal to half of the Hartree one with opposite sign. The HF equation can be written as

$$\left[\frac{1}{2}^2+v(r)+\frac{1}{2}v_H[\rho _{HF}]\right]\varphi _{HF}=\epsilon _{HF}\varphi _{HF},$$
(7)
$`v(r)=\frac{1}{2}\gamma ^2r^2`$ and $`\epsilon _{HF}`$ is the HF orbital energy. The total HF energy is written as

$$E_{HF}=2\epsilon _{HF}\frac{1}{2}𝑑\stackrel{}{r}\rho _{HF}v_H[\rho _{HF}],$$
(8)
where $`v_H`$ is the Hartree potential given by

$$v_H[\rho ]=𝑑\stackrel{}{r}^{}\frac{\rho (\stackrel{}{r}^{})}{|\stackrel{}{r}\stackrel{}{r}^{}|}.$$
(9)
The KS equation can be written as

$$\left[\frac{1}{2}^2+v(r)+v_H[\rho _{KS}]+v_x(\rho _{KS})\right]\varphi _{KS}=\epsilon _{KS}\varphi _{KS}$$
(10)
$`\epsilon _{KS}`$ and $`\varphi _{KS}`$ are the KS orbital energy and eigenfunction, respectively, and again $`\rho _{KS}=2|\varphi _{KS}|^2`$. We take here the local-density approximation for which $`v_x(\rho )=\frac{4}{3}c_x\rho ^{1/3}`$, with $`c_x=\frac{3}{4}\left(\frac{3}{\pi }\right)^{1/3}`$ (see e.g. ). The total KS energy is given by

$$E_{KS}=2\epsilon _{KS}\frac{1}{2}𝑑\stackrel{}{r}\rho _{KS}v_H[\rho _{KS}]\frac{1}{3}𝑑\stackrel{}{r}\rho _{KS}v_x(\rho _{KS}).$$
(11)
Eqs. (7) and (10) have asymptotae controlled by the harmonic potential, so the asymptotic density looks like

$$\rho _a(\stackrel{}{r})\mathrm{exp}(\gamma r^2).$$
(12)
## III Results and Discussion

First we make an analysis of the solution of the exact case using the 1/N approximation. For details in the derivation of the formulae related to the shifted-1/N method, we refer to the literature (see e.g. ).

In order to establish the accuracy of the results obtained by the application of the shifted $`1/N`$ method we compare them with those obtained by using the Schwartz’s numeric method . The method is based in a numerical approximation of functions on a mesh and gives very accurate results . There are only empirical estimates of the error  which turns out to be exponentially decaying with the number of points given the mesh step. The interpolation function is chosen as

$$f(r)=\underset{m}{}f_m\frac{u(r)}{(rr_m)a_m},$$
(13)
where

$$u(r)=\mathrm{sin}[\pi (r/h)^{1/2}].$$
(14)
Here $`r_m`$ is a zero of $`u(r)`$, $`a_m`$ is a zero of its derivative, and $`h`$ is the step of the mesh which turns out to be quadratically spaced. The choice of the step $`h`$ was made after estimating the characteristic length of the effective potential, and then multiplying the obtained estimate by five and dividing it by the square of the number of points in the mesh, usually around 300. This guarantees that $`h`$ is minimal for a given $`\omega `$, and also that the function value at the last mesh point is practically zero.

In figure 1 we show the behavior of the error of the energy of few lowest eigenvalues calculated using the shifted-1/N method, compared with the accurate results obtained using the Schwartz’s method. We plotted the relative error defined as $`\delta E=E_{approx}/E_{exact}1`$, as function of the reduced variable $`\gamma ^{}=\gamma /(\gamma +1)`$. The $`0s,0p`$, and $`0d`$ states are the three lowest energy states of the relative motion with angular momentum $`l=0,1`$, and $`2`$, respectively. It can be seen that the error remains bounded in $`0.5\%`$ for the first eigenvalue, in $`0.075\%`$ for the second one and $`0.022\%`$ for the third one. It is noticeable that the method always underestimates the correct values of the energy, and that the error decreases with the increase of the angular quantum number $`l`$, as it should be expected. Furthermore, the error has some maximum between the two limiting cases, after which it decreases as expected, since the shifted-1/N method reproduces exactly the oscillator case.

Figure 2 shows the behavior of the total energy (the sum of the centre-of-mass and relative motion energies) as a function of $`\gamma ^{}`$, as calculated using the shifted-1/N method. The first two symbols are the indices of the centre-of-mass component of the energy, and the last two correspond to the relative motion. The lowest six states are depicted. The inset is a magnification for $`\gamma ^{}`$ from $`0.01`$ to $`0.15`$. Here, an apparent linear behavior can be observed, but a more careful analysis indicates that, for instance, for the ground-state first six-eight points, the effective power in terms of $`\gamma ^{}`$ is about $`0.8`$. In the last few points the curve can be fitted well with a function of the type $`\gamma ^{}/(1\gamma ^{})`$ which is the expected behavior for electrons in a harmonic field. The relative error is not shown, but it is estimated in roughly one half of the one shown in figure 1, since the centre-of-mass energy can be calculated exactly and it is typically of the size of the relative motion energy.

Although we are able to solve the Schrödinger equation very accurately for this system, perturbation analysis can give some more insight of the behavior of the electrons under weak and strong confinement. For weak confinement ($`\gamma 0`$) we have that the kinetic energy term can be neglected (see e.g. ), which corresponds to the strong interaction limit (Wigner crystallization). Then, the energy is approximately taken as the minimum of the effective potential (this is the zeroth order approximation to the energy),

$$V_l(r)=l(l+1)\frac{1}{r^2}+\frac{1}{4}\gamma ^2r^2+\frac{1}{r},$$
(15)
the minimum is reached for $`r_0^l`$ that satisfies the equation

$$\gamma ^2r_0^{l\mathrm{\hspace{0.33em}4}}2r_0^l4l(l+1)=0.$$
(16)
For $`l=0`$ $`r_0^0=(2/\gamma ^2)^{1/3}`$ and for $`l\gamma /2`$ $`r_0^l(2l/\gamma )^{1/2}`$. The minimum of the potential is then

$$U_0^l=l(l+1)\frac{1}{r_0^{l\mathrm{\hspace{0.33em}2}}}+\frac{1}{4}\gamma ^2r_0^{l\mathrm{\hspace{0.33em}2}}+\frac{1}{r_0^l},$$
(17)
which is $`2^{1/3}3/2\gamma ^{2/3}`$ for $`l=0`$ and $`\gamma (2l+1)/2+\gamma ^{1/2}/(2l)^{1/2}`$ for large enough values of $`l`$.

The next order in the approximation is to get the effective frequency $`\gamma _l^2=\frac{1}{2}d^2V_l(r)/dr^2|_{r=r_0^l}`$,

$$\gamma _l^2=\frac{1}{4}\gamma ^2+\frac{1}{r_0^{l\mathrm{\hspace{0.33em}3}}}+\frac{3l(l+1)}{r_0^{l\mathrm{\hspace{0.33em}4}}},$$
(18)
which is $`\gamma _0^2=3/4\gamma ^2`$ and $`\gamma _l^2=\gamma ^2+3/4\gamma ^2/l+2^{1/2}/4(\gamma /l)^{3/2}`$ for large $`l`$. Now we can estimate the energy levels of this effective harmonic field. So the total energy in the weak confinement limit is

$$E_{nl}l(l+1)\frac{1}{r_0^{l\mathrm{\hspace{0.33em}2}}}+\frac{1}{4}\gamma ^2r_0^{l\mathrm{\hspace{0.33em}2}}+\frac{1}{r_0^l}+2\gamma _l(n+1/2),$$
(19)
where $`n=0,1,..`$. For $`l=0`$

$$E_{n0}\frac{3}{2}2^{1/3}\gamma ^{2/3}+\frac{3^{1/2}}{2}\gamma (n+1/2).$$
(20)
For $`l`$ large enough

$$E_{nl}\frac{1}{2}\gamma (2l+1)+\left(\frac{\gamma }{2l}\right)^{1/2}+2\left[\gamma ^2+\frac{3}{4}\frac{\gamma ^2}{l}+\frac{2^{1/2}}{4}\left(\frac{\gamma }{l}\right)^{3/2}\right]^{1/2}(n+1/2).$$
(21)
Equation (20) gives the explanation why the effective power of the ground-state energy for small $`\gamma ^{}`$ is approximately $`0.8`$: it is between $`2/3`$ and $`1`$, the effective powers for the weak field limit for the relative motion and the centre-of-mass energies in terms of $`\gamma ^{}\gamma `$.

In the strong confinement regime, the zeroth order approximation amounts to neglect the electron-electron interaction, so it corresponds to the oscillator’s energy. Application of the first order of the perturbation theory  for strong confinement ($`\gamma \mathrm{}`$) for the relative coordinate equation leads to

$$E_{nl}E_{nl}^{osc}+\varphi _{nl}^{osc}|r^1|\varphi _{nl}^{osc},$$
(22)
where $`E_{nl}^{osc}=2\gamma (2n+l+3/2)`$ and substituting $`\varphi _{nl}^{osc}`$ (which is similar to equation (4) ) into the above equation gives

$$E_{nl}2\gamma (2n+l+3/2)+\frac{\gamma ^{1/2}}{2^{1/2}}\underset{k,k^{}=0}{\overset{n}{}}a_{nl,2k}a_{nl,2k^{}}\mathrm{\Gamma }(l+k+k^{}+1)𝒩_{nl}^2$$
(23)
where $`a_{nl,2k}`$ are the coefficients of the generalized Laguerre polynomials (of Eq.(4)) that satisfy the recursion

$$a_{nl,2k}=a_{nl,2(k1)}\frac{kn1}{k(k+l+1/2)}$$
(24)
with $`a_{nl,0}=1`$. Furthermore, the normalization constant is given by

$$𝒩_{nl}^2=\underset{k,k^{}=0}{\overset{n}{}}a_{nl,2k}a_{nl,2k^{}}\mathrm{\Gamma }(l+k+k^{}+3/2).$$
(25)
Based on the perturbative results of Eqs. (17), (19) and (23) we propose the following interpolation formula

$$E_{nl}^{(int)}(\gamma )=\frac{\gamma ^1E_{nl}^{(0)}+\gamma E_{nl}^{(\mathrm{})}}{\gamma ^1+\gamma },$$
(26)
where the superscripts $`(0)`$ and $`(\mathrm{})`$ correspond to the zero and infinite confinement limits. For $`\gamma 0`$ equation (26) will return approximately the weak confinement limit and for $`\gamma \mathrm{}`$ the strong confinement limit. For the zeroth order approximations equation (26) leads to

$$E_{nl}^{(int)}(\gamma )=\frac{\gamma ^1U_0^l+2\gamma ^2(2n+l+3/2)}{\gamma ^1+\gamma },$$
(27)
and for the first order pertubative results, it will be

$$E_{nl}^{(int)}(\gamma )=\frac{\gamma ^1[U_0^l+2\gamma _l(n+1/2)]+\gamma [2\gamma (2n+l+3/2)+2^{1/2}\gamma ^{1/2}\mathrm{\Delta }_{nl}]}{\gamma ^1+\gamma },$$
(28)
where $`\mathrm{\Delta }_{nl}`$ denotes the summation in the right hand side of equation (23). The interpolation scheme of equation (27) performs consistently bad, except for the very weak, and very strong fields. The error goes up to $`41\%`$ for the ground-state, underestimating the correct values. The introduction of first order corrections, corresponding to the interpolation scheme of equation (28) brings a dramatic improvement on the values: for the ground-state the relative error is never worse than $`3.3\%`$ for the energy of the relative motion, which means that the relative error of the total energy is around $`1.7\%`$.

In the case of the mean-field approximations like Hartree-Fock and Kohn-Sham, we have considered the paired-electron ground-state case, and since the confining potential is radially symmetric, the orbitals and the density are also radially symmetric. We implemented the shifted-1/N technique for the HF and KS equations. Here, a word about accuracy is needed: although the solution of the Schrödinger-like equations using Schwartz’s method is very accurate, the estimate of the Hartree potential is not as accurate anymore, nevertheless, five to six figures are always guaranteed. The results are shown in table 1 in the first and third columns, indicated as $`HF1/N`$ and $`KS1/N`$, respectively. Since the resulting wavefunction from the application of the 1/N has a complicated form which makes difficult a direct evaluation of the density or of the Hartree potential, we assumed that as initial guess the non-interacting density, which is also correct asymptotically. For comparison purposes, we have used again the Schwartz technique self-consistently to solve both the HF and KS equations. Numerical results are shown in the second and fourth columns of table 1, indicated as $`HFnum`$ and $`KSnum`$, respectively. Also for comparison purposes, we included the results of solving the full Schrödinger equation using the shifted-1/N method and the numeric solution, which are the last two columns indicated as $`Exact1/N`$ and $`Exactnum`$, respectively.

Table 1 shows very good agreement between the results from the accurate numerical result and the ones calculated with the 1/N method. We should not be too enthusiastic about the accuracy, since the remarkable agreement is probably a result of the compensation of errors from the calculation of the Hartree potential and energy (due to its simplicity) and the intrinsic error of the 1/N method, specially for the weak field case. It can also be seen that the relative accuracy improves from typically few times $`10^2`$ to $`10^410^5`$ with the increase of the strength of the field. Also the systematic difference between the HF and KS methods is reduced with stronger confinement. Both behaviors can be understood by taking into account that the confinement potential dominates over the decaying electron-electron interaction potential, and with the increase of the strength of confinement the problem becomes just a harmonic potential problem, for which the shifted-1/N is designed to give the exact energy, although the quality of the wavefunction is not too good.

In order to better understand the behavior of the pair of electrons in strong and weak confinement within these mean-field theories, we can recourse to pertubation analysis. In the strong confinement limit ($`\gamma \mathrm{}`$) the system will behave basically as a pair of non-interacting electrons, and the electron-electron interaction (Hartree and exchange potentials) can be considered as a perturbation. The HF and KS orbitals become

$$\varphi (\stackrel{}{r})=\left(\frac{\gamma }{\pi }\right)^{3/4}\mathrm{exp}(\frac{1}{2}\gamma r^2),$$
(29)
so, the Hartree energy will be $`2(2\gamma /\pi )^{1/2}`$ and the LDA exchange $`3c_x2^{1/3}(3\gamma /\pi )^{1/2}/4`$. Then the total HF energy will be

$$E_{HF}^{(\mathrm{})}=3\gamma +\left(\frac{2\gamma }{\pi }\right)^{1/2},$$
(30)
and the total KS energy

$$E_{KS}^{(\mathrm{})}=3\gamma +\left(2^{3/2}+\frac{3}{4}2^{1/3}3^{1/2}c_x\right)\left(\frac{\gamma }{\pi }\right)^{1/2}.$$
(31)
In the weak confinement limit ($`\gamma 0`$), due to scaling arguments, we can neglect the kinetic energy, and we can assume constant density, at least within a certain radius $`R`$ (this is asymptotically true for the HF approximation and arguably for the KS one). Now taking Eqs. (7) and (10), and using Poisson’s equation, we get that $`\rho _{HF}\gamma ^2/(2\pi )`$ and $`\rho _{KS}\gamma ^2/(4\pi )`$. Then from normalization condition $`R=[3/(2\pi \rho )]^{1/3}`$. The Hartree potential will take the form

$$v_H(r)=\frac{3}{R}\frac{r^2}{R^3},ifrR,$$
(32)
and $`v_H(r)=2/r`$ if $`rR`$. The Hartree energy will be equal to $`12/(5R)`$. Substituting the above result into equation (8) we get that the Hartree-Fock energy will be

$$E_{HF}^{(0)}(3\gamma )^{2/3}.$$
(33)
Substituting the result for the Hartree energy into equation (11) we find that in this limit the Kohn-Sham energy will be equal to

$$E_{KS}^{(0)}\gamma ^{2/3}(6^{3/2}+\frac{2c_x}{5(4\pi )^{1/3}}).$$
(34)
Again, we can use interpolation schemes like the one of equation (26), for the HF ground-state energy we will have

$$E_{HF}^{(int)}(\gamma )\frac{(3\gamma )^{2/3}+\gamma ^2[3\gamma +(2\gamma /\pi )^{1/2}]}{1+\gamma ^2}.$$
(35)
In the Kohn-Sham case we can write down

$$E_{KS}^{(int)}(\gamma )\frac{\gamma ^{2/3}[6^{3/2}+2c_x/5/(4\pi )^{1/3}]+\gamma ^2[3\gamma +(2^{3/2}+\mathrm{3\hspace{0.33em}2}^{1/3}3^{1/2}c_x/4)(\gamma /\pi )^{1/2}]}{1+\gamma ^2}.$$
(36)
Comparing the asymptotics of the exact equation for the ground-state in the strong confinement limit with the asymptotics of the HF equation, we observe that they coincide. In the Kohn-Sham case the zeroth order is the same, but the coefficient of the first order perturbation is slightly higher ($`1.62`$ compared to the exact $`\sqrt{2}1.41`$). In the low density limit (weak confinement), although the three asymptotes are proportional to $`\gamma ^{2/3}`$ in the zeroth order, the coefficients differ substantially, the exact is $`1.19`$, the HF $`2.08`$, and the KS $`3.17`$. This indicates a consistent overestimation of the energy by the independent-electron approximations, which is the expected behavior at least for the HF approximation. Our findings for this system for the weak and strong asymptotae, and the numerical results, suggest that KS-LDA energies are always higher than the HF ones, which is consistent with numerical experience on atoms and molecules . From the results shown in table 1 we can estimate the correlation energies for the ground-state for different confinement strengths: for $`\gamma ^{}=0.1`$, the correlation energy is about $`5.6\%`$ of the total energy, meanwhile for $`\gamma ^{}=0.9`$ (strong confinement) it is about $`0.17\%`$ (for the lightest many-electron atom, Helium, it is about $`1.4\%`$, and for Argon with $`Z=18`$, it is only $`0.14\%`$). This adds evidence to the suggestion that for harmonic fields correlation effects are more important, also in three dimensions, especially for weak confinement.

## IV Concluding remarks

In the present article we have calculated the energy spectrum of a two-electron spherical quantum dot for the few lowest states, using the shifted-1/N method and the very accurate numerical Schwartz’s method. From the comparison of the numerical results we could assess the quality of the shifted-1/N method, which consistently underestimates the correct values, although the error is rather small, and it decreases with the increase of the relative angular momentum. We have also applied perturbation theory up to first order in both limits of confinement of the electron, and we have proposed some interpolation formulae for the energy. Inclusion of first order perturbation allowed to construct an interpolation expression that performs reasonably well. We also solved the mean-field Hartree-Fock and local-density Kohn-Sham problems for the ground-state. Using the shifted-1/N method we got from reasonable to high accuracy already in the first iteration, compared to the self-consistent numeric solution using the Schwartz’s method. We made an analysis of the strong and weak confinement limits, and proposed interpolation formulae for both the Hartree-Fock and Kohn-Sham ground-state energies. It was shown that the correlation energy is relatively big for this systems, especially for weak confinement.

## V Acknowledgments

One of the authors (V.M.V.) would like to acknowledge support by CONICIT, Venezuela, under project 96000061.


format: text

text:
# Anomalous Response Variability in a Balanced Cortical Network Model

## 1 Introduction

Cortical neuronal responses often exhibit a puzzling variability (see, e.g., ): spike count distributions obtained for repeated presentations of a stimulus are frequently broader than Poisson distributions with the same mean counts. In this paper we study the statistics of neuronal responses in a simple model network with balanced excitation and inhibition and find that they also show large variability. We focus on the so-called Fano factor $`F`$: the ratio of the variance of the spike count to its mean. It is equal to 1 for any Poisson process, homogeneous or inhomgeneous. We find that in our model it depends sensitively on the neuronal firing threshold: $`F<1`$ for high thresholds and $`F>1`$ for low ones, while mean spike counts are hardly affected by changes in threshold.

To study the responses in our model, we utilize mean field theory, which is exact in the limit of a large network with homogeneous connection probabilities . It allows us to reduce the full network problem to one of a single neuron subject to a Gaussian synaptic current, the mean and covariance of which are self-consistently related to those characterizing the the neuron’s firing. It is Gaussian because of the central limit theorem: it is the sum of a large number of (what can be proved to be) independent contributions from the other neurons.

The mean field theory of balanced cortical networks has been studied by a number of authors , but generally in or near steady firing states and assuming that the self-consistent current input is uncorrelated. Here, like ref. , we consider time-dependent external drive (as in an experimental trial) and color the noise correctly.

## 2 The Model

For these exploratory investigations, we use a network of $`N`$ mutually inhibitory neurons with a 10% connection probability. This is probably the simplest model for which one can achieve balanced asynchronous activity. The nonzero connections are of equal strength, and the synaptic currents are assumed to act instantaneously, i.e., each presynaptic spike depresses the postsynaptic potential discontinously by a fixed amount. We do not include transmission delays. This and other features, such as the presence of two neuronal populations, inhibitory and excitatory, can easily be incorporated into the model. A modest extension of this form might serve as a serious model of a cortical (mini)column. However, we believe that the basic dynamical features of balanced asynchronous networks are already captured in our simple all-inhibitory model.

We find it convenient to scale variables in the way used by van Vreeswijk and Sompolinsky : When the mean number of neurons presynaptic to a given one is $`K1`$, each synaptic strength is given a value $`J/\sqrt{K}`$. Then the average synaptic current scales like $`K\times 1/\sqrt{K}=\sqrt{K}`$. It is inhibitory and is counterbalanced by an external excitatory current which is also proportional to $`\sqrt{K}`$. The fluctuations in the synaptic current are smaller than the mean by a factor of order $`1/\sqrt{K}`$, i.e., of order 1. In steady state, the two currents nearly cancel, leaving a net current of order 1. If the steady-firing state is stable, disturbances of the rates relax very rapidly, in a time of order $`1/\sqrt{K}`$. While this picture is something of a caricature, we believe it has some relevance to cortical dynamics.

Mean field theory can be used for any kind of neuron, but here we use leaky integrate-and-fire neurons. Thus, the (subthreshold) membrane potentials are described by

$$\dot{u}_i(t)=\frac{1}{\tau }u_i+\sqrt{K}I_i^0(t)\frac{J}{\sqrt{K}}\underset{k,s}{}c_{ij}\delta (tt_j^s),$$
(1)
where $`\tau `$ is their (common) membrane time constant, $`I_i^0(t)`$ is the excitatory external input current felt by neuron $`i`$, $`c_{ij}`$ is 1 or 0 according to whether there is a connection from neuron $`j`$ to neuron $`i`$, and $`t_j^s`$ is the $`s`$-th spike time of neuron $`j`$. The external input $`I_i^0(t)`$ represents input from elsewhere in the brain (e.g., the preceding stage in a sensory pathway), and as such is noisy itself. However, as the recurrent connections in the randomly-diluted network generate dynamical noise on their own, this extrinsic noise does not have a big qualitative effect, so we take $`I_i^0(t)`$ to be constant here. We also take it to be uniform across the population.

## 3 Mean Field Theory

As indicated above, in mean field theory, one studies a single neuron for which the recurrent synaptic current (the last term in (1)) is replaced by a self-consistent Gaussian current with self-consistent mean and variance. Explicitly, the subthreshold membrane potential of this single neuron obeys

$$\dot{u}=\frac{1}{\tau }u_i+\sqrt{K}[I^0+I^1(t)].$$
(2)
Whenever $`u`$ reaches a threshold $`\theta `$, it fires a spike and $`u`$ is reset to zero. The effective recurrent current $`I^1(t)`$ has a mean

$$I^1(t)=Jr(t)$$
(3)
proportional to the instantaneous firing rate $`r(t)`$ of the neuron and a covariance

$$\delta I^1(t)\delta I^1(t^{})=\frac{1K/N}{K}C(t,t^{}),$$
(4)
where $`C(t,t^{})=[S(t)r(t)][S(t^{})r(t^{})]`$ is the autocorrelation funtion of the neuronal firing $`S(t)=_s\delta (tt^s)`$. The $`K`$ in the denominator in (4) comes from the averaging over $`K`$ independent inputs, and the $`1K/N`$ in the numerator is a correction for finite connection concentration $`K/N`$. It can be derived using the methods of ref. .

This model can not be solved exactly analytically, but it is simple to solve numerically, using the method first introduced for spin glasses by Eisfeller and Opper . We start with a guess at the form of the mean and covariance function of the random current $`I^1(t)`$ and run a series of “trials”, in each of which we integrate (2) for an independent realization of $`I^1(t)`$. By averaging the output of the neuron over the trials, we get an estimate of $`r(t)`$ and $`C(t,t^{})`$, which is used to generate new examples of $`I^1(t)`$ for another set of trials. This loop is then iterated until the statistics converge to self-consistency. In our calculations here, we used 10000 trials per iteration and up to 30 iterations.

Each trial was 100 integration steps (which we call “milliseconds”) long. We chose parameters $`K=500`$, $`N=5000`$, and $`\tau =10`$ ms. The external excitatory input $`I^0(t)`$ was constant at a low value (evoking a background firing rate of 3 Hz) during the first and last 10 ms. In the middle 80 ms, an additional “stimulus” input, which peaked about 15 ms after onset, was added. It typically evoked about 4-5 spikes, with peak rates around 100 Hz. The spike count distributions, PSTHs and covariance functions were computed for a number of values of the firing threshold $`\theta `$.

## 4 Results

For all cases studied, the relaxation of the network to its state of balanced excitation and inhibition was very rapid, as expected, so the response tracked the time course of the excitatory external input closely (Fig. 1). The overall response strengths vary only weakly with threshold: a factor of 16 difference between the smallest and largest threshold values produced only a 20% difference in mean response. This is because the increased firing that would be produced by lowering the threshold is largely compensated by the comcomitant increase in inhibition.

However, varying the threshold had a strong effect on the irregularity of the firing. Fig. 2 shows the spike count distributions for three threshold values with ratio 1:4:16. While the intermediate value fits a Poisson distribution well, the low threshold leads to an anomalously broad distribution ($`F=2.4`$) and the high one to an anomalously narrow one ($`F=0.25`$).

These differences are also evident in the autocorrelation function $`C(t,t^{})`$. Fig. 3) shows $`C(t,t^{})`$ as a function of $`t^{}`$ for a fixed value of $`t`$ for the same three threshold values as in the preceding figures. For the lowest threshold, there is a “hill” centered around $`t=t^{}`$, while for the highest there is a valley. These are indicative of spike “bunching” and “antibunching”, respectively, leading naturally to the higher and lower spike count fluctuations seen in Fig. 2. The intermediate threshold value shows very little correlation (apart from the delta-function peak at $`t=t^{}`$), consistent with the nearly-Poisson count distribution found in this case.

## 5 Discussion

Measured Fano factors in visual and IT cortex  vary over a range at least as large as the one-order-of-magnitude difference between those for the smallest and largest thresholds described above. Of course, threshold differences are not the only possible source of such response variability. We have also explored the effects of varying synaptic strengths, with similar results, and it seems likely that diffferences in a wide variety of single-neuron properties can have the same kind of effect.

Neurons in a local cortical network can not all be expected to have the same threshold, and, furthermore, their thresholds (or other parameters) may fluctuate (uncontrollably) from trial to trial. We have also found large variations in the Fano factor in a model where these fluctuations are assumed independent for different neurons and in different trials.

All these results are only suggestive, and more systematic work, both experiments and modeling, is called for. However, they do point to the possibility that the observed response variability of cortical neurons may be accounted for in terms of natural variations in properties from neuron to neuron and trial to trial.


text:
# Equation of state in a small system: Violation of an assumption of Maxwell’s demon

## Abstract

An equation of state for an ideal gas with a small number of particles is studied. The resulting equation is found to differ from that expected in conventional thermodynamics, which is strikingly illustrated when considering the traditional thermodynamic problem of Maxwell’s demon. We clarify the mechanism of this different feature of thermodynamics arising in small systems.

The thought experiment by MaxwellMaxwell  known as Maxwell’s demon has attracted the attention of physicists for more than a centuryreview . Maxwell’s demon is a well-known paradox based on the fundamentals of thermodynamics, which seemingly permits a violation of the second law of thermodynamics. In the original thought experiment, we are asked to imagine a demon that controls a gate between two compartments. Being able to determine the speed of the particles in each compartment, the demon opens the gate in such a way as to collect more high-speed particles in one compartment, thus decreasing the entropy of the system, in seeming violation of the second law. About half a century after Maxwell, Szilard extended Maxwell’s demon, devising a model which allowed quantitative description on the extracted workSzilard .

The framework of Szilard’s formulation is depicted in Fig.1. One thermal molecule is in an isothermal system and the cycle of Maxwell’s demon is considered to be composed of two processes: 1) measurement (observation) of the molecule position and 2) isothermal expansion in the process from State B to C. By measurement, the demon determines whether the molecule is on the left or right side and inserts a wall to create a piston, which is accompanied by a mechanical load (weight) for extraction of work. In this way, work of an amount $`W=k_\mathrm{B}T\mathrm{ln}2`$ can be extracted from the system by isothermal expansion. Through cyclic operation, i.e. following the process from States A $``$ B $``$ C $``$ A, work is perpetually extracted out of the isothermal system. The operation transfers ambient thermal energy into work through the isothermal expansion of the “piston”, resulting in a decrease in entropy, $`k_\mathrm{B}\mathrm{ln}2`$. In order to save the second law, Szilard described how the entropy decrease would be compensated “… if the execution of such a measurement were, for instance, always accompanied by production of $`k\mathrm{ln}2`$ units of entropy”Szilard .

As has recently been pointed outMagnasco , Szilard did not prove his compensation theory, but assumed that his heat engine would always preserve the second law. Most studies of Maxwell’s demon have followed the theoretical framework of Szilardreview , and have focused on the energy cost of observation 11 . In such studies, it is generally assumed that, for ideal gases, one can extract work, $`W`$, even for a one-molecule isothermal piston, of an amount given by the relation $`W=_{V_i}^{V_f}P(V)𝑑V=nRT\mathrm{ln}(V_f/V_i)`$, where $`P`$ is the pressure, $`n`$ is the number of molecules in a mole, $`R`$ is the gas constant, $`T`$ is the temperature, and $`V_i`$, and $`V_f`$ are the initial and final volumes of the piston, respectively22 . Thus, the conventional equation of state for an ideal gas has been implicitly assumed: $`PV=nRT`$. In other words, studies of Maxwell’s demon have generally been performed within the framework of conventional macroscopic thermodynamics, regardless of the small number of moleculesHill .

On the other hand, recent experimental developments have enabled experimental observation of even a single molecule in some situations. In molecular motorsFiner ; Miyata ; Molloy  the value of the work extracted in an elementary process fluctuates over repeated experiments, which is in contrast to the thought experiment of Maxwell’s demon. This suggests that there may be an essential fault in the basic assumption of Maxwell’s demon on this scale, because an essential experimentally observed feature of a small-scale system, thermal fluctuation, has not been properly accounted for. In addition, recent studies that describe the energetics in thermally fluctuating systems have shown that there is another way of describing the thermodynamics of small systems kinetically without employing entropyKenJPS . From this background, a question was recently raised by Hatano and SasaHatano  as to the validity of the assumption of the work extracted in Szilard’s model. The study was an attempt at quantitative estimation of the work extracted in an operational manner, which is in contrast to conventional thought experiments. Although their study was a pioneering one, the model resulted in broken detailed balance even in the equilibrium state.

In this paper, we analyze the equation of state in an extremely small system, within the boundary of kinetics. We discuss an equilibrium process here; however, an interesting aspect of relaxation into equilibrium is discussed by Crosignani and PortoCrosignani . We show that the equation of state that Szilard and others have assumed for isothermal expansion is not appropriate. We clarify the theoretical origin of the difference between an extremely small system and conventional macroscopic systems towards the construction of a system of thermodynamics for small systemsHill .

We begin with an estimation of the equation of state of one molecule because typical features of small systems are expected at this limit. In order to derive the equation, we must confine ourselves to a concrete model that permits us to treat the system only within the framework of kinetics. As in the study by Hatano and SasaHatano , we consider a one-dimensional piston in thermal equilibrium (Fig. 2), where the dynamical variables are $`N`$ “thermal particles” and a movable piston on one side. By replacing the volume $`V`$ by the distance between the wall and the movable piston $`X`$, $`n`$ by $`N/n_\mathrm{a}`$, and $`R`$ by $`n_\mathrm{a}k_\mathrm{B}`$, we obtain the equation of state for the system, $`PX=Nk_\mathrm{B}T`$, where $`n_\mathrm{a}`$ and $`k_\mathrm{B}`$ are the Avogadro constant and the Boltzmann constant, respectively.

Hatano and Sasa claimed that for ideal gases the equation of state depends on the masses of the thermal particles and the piston and that the conventional equation of state is valid if the ratio of the mass of the piston to that of the particle is sufficiently large. However, we found, by the method of stochastic energeticsKenJPS , that in the study by Hatano and Sasa, there is a finite current between two thermal baths, one of which is in contact with a thermal wall and the other of which is in contact with a movable piston, resulting in the broken detailed balance. We can learn from this that we cannot naively utilize a “hybrid” model, for example, Langevin dynamics (for a piston) and Hamiltonian dynamics (for a thermal particle) simultaneously, because Langevin dynamics is a coarse-grained description of Hamiltonian dynamicsMori . In order to maintain self-consistency in a class of description, we should return to the framework of elementary mechanics.

In most literature on Maxwell’s demon, it is assumed that we can fix a movable piston by applying an external force of the same value as the conventional (internal) pressure. This assumption fails in a small system, as discussed later. Thus, we replace the coarse-grained quantity “pressure” $`P`$ by the mechanical force $`f`$ externally exerted on the piston. This yields a reference equation, $`fX=Nk_\mathrm{B}T`$, which we will use to discuss the validity of the conventional equation of state, and hence the validity of the frameworks of the previous studies. We assume that collisions between the variables (particle(s) and a piston) are perfectly elastic and that the movable piston is frictionless, since “friction” is also a coarse-grained quantity which must be excluded for consistency. Then, the total energy of the system of one thermal particle ($`N=1`$) is written as

$$E=\frac{p_\mathrm{t}^2}{2m}+\frac{p_\mathrm{p}^2}{2M}+fX,$$
(1)
where $`m(p_\mathrm{t})`$ and $`M(p_\mathrm{p})`$ are the masses (momenta) of the thermal particle and the piston, respectively, and $`f`$ is a force applied to the movable piston. To keep the system isothermal, a thermal wallBergmann  is introduced. The particle is reflected with a positive random velocity according to Maxwell’s velocity distribution at the boundary of the system, $`x=0`$, where $`x`$ is the position of a thermal particle. The thermal wall realizes thermal equilibrium of the system. Because the system is one-dimensional, the thermal particle is always lower than the piston, $`xX`$.

One naively expects that the equation of state for one molecule should be obtained by replacing $`N`$ by $`1`$ in the reference equation $`fX=Nk_\mathrm{B}T`$; namely $`fX=k_\mathrm{B}T`$. The extension of conventional thermodynamics into those of one molecule has been performed since SzilardSzilard . In the following, we will investigate the validity of this extension in order to uncover the nature of thermodynamics for a few molecules. Three methods of analyses are performed independently: 1) numerical simulation, 2) Master equation analysis, and 3) Gibbs’ statistical mechanics.

In the simulation, the system obeys Hamiltonian dynamics with the exception of the thermal wall, in which the thermal particle is reflected with a random velocity within the velocity distribution $`\frac{m|v|}{k_BT}e^{mv^2/2k_BT}(v>0)`$ when the particle arrives at the boundary $`x=0`$. The condition of the thermal wall is known to direct the target system into thermal equilibriumBergmann . The numerical simulation revealed that the equation of state for one molecule is in fact

$$fX=2k_\mathrm{B}T,$$
(2)
where the angular brackets indicate an average value. This result is obviously different from that which is conventionally assumed.

We obtained the same analytical result using the Master equation of the distribution function in a phase space, which consists of Liouville terms and collision terms. Let $`\rho (x,v,X,v_p)`$ be the probability density of the particles (a thermal particle and a movable piston), and $`v`$ and $`v_p`$ the velocities of the thermal particle and movable piston, respectively. Here we use velocities instead of momenta for simplicity. The density function obeys the following equation:

$`{\displaystyle \frac{\rho (x,v,X,v_p)}{t}}=`$ (3)
$`v{\displaystyle \frac{\rho }{x}}v_p{\displaystyle \frac{\rho }{X}}+{\displaystyle \frac{f}{M}}{\displaystyle \frac{\rho }{v_p}}`$
$`\theta \left(vv_p\right)\left(vv_p\right)\rho (x,v,X,v_p)\delta \left(xX\right)`$
$`+\theta \left(v_pv\right)\left(v_pv\right)`$
$`\times \rho (x,{\displaystyle \frac{mM}{M+m}}v+{\displaystyle \frac{2M}{M+m}}v_p,X,{\displaystyle \frac{2m}{M+m}}v+{\displaystyle \frac{Mm}{M+m}}v_p)`$
$`\delta \left(xX\right).`$
The first three terms of the R.H.S. of this equation come from the Liouville equation and the last two terms come from collision effects. A stochastic boundary condition is applied to the system at one end, that is, $`\rho (x=+0,v)=\frac{f}{T}\frac{\sqrt{m}}{\sqrt{2\pi T}}\mathrm{exp}[mv^2/2T]`$ for $`v>0`$. The stochastic condition corresponds to the thermal wall in the numerical simulation. From straightforward calculation, one obtains a stationary solution, $`\rho (x,p_\mathrm{t},X,p_\mathrm{p})=\frac{f^2}{2\pi (k_\mathrm{B}T)^3\sqrt{mM}}\mathrm{exp}\{(\frac{p_\mathrm{t}^2/2m+p_\mathrm{p}^2/2M+fX}{k_\mathrm{B}T})\}\theta (Xx),`$ where $`\theta `$ is a Heaviside step function. From this equation, we again obtain Eq. (2). Note that the result is independent of the masses of both the thermal particle and the movable piston, which is in contrast to the result by Hatano-Sasa Hatano .

Although we introduced a thermal wall and, correspondingly, a stochastic boundary condition, one may obtain the same generalized result without applying such boundary conditions, but instead using conventional Gibbs’ statistical mechanics. Since the present system has two constant parameters, temperature and force (applied to the piston), the system can be described by a pressure ensemble (P-T ensemble) of Gibb’s statistical mechanicsKubo . The quantity “pressure” in a conventional pressure ensemble is replaced here by a force. This replacement does not alter any general result of the pressure ensemble, because one can reprove any derivation in the pressure ensemble with this replacement. The Laplace transformation between a partition function of the canonical ensemble $`Z`$ and that of the pressure ensemble $`Y`$ is given as $`Y=_0^{\mathrm{}}𝑑V\mathrm{exp}[PV/(k_\mathrm{B}T)]Z`$Kubo , from which we obtain the partition function of the pressure ensemble of the system of $`N`$ particles, $`Y=F(\frac{k_\mathrm{B}T}{f})^{N+1}`$, where $`F`$ is a factor that will drop out later. The average position of the movable piston is obtained as $`X=\frac{1}{k_\mathrm{B}T}\frac{\mathrm{ln}Y}{f}=(N+1)k_\mathrm{B}T/f.`$ Thus, we have $`fX=(N+1)k_\mathrm{B}T.`$ The previous result, $`fX=2k_\mathrm{B}T`$, is again obtained by replacing $`N`$ by 1. One finds that in a sufficiently large system, $`N1`$ (thermodynamic limit), $`N+1`$ may be replaced by $`N`$. However, in the small system, the difference between $`N+1`$ and $`N`$ is crucial, as shown in Fig. 3.

This result is apparently in contrast to the basic assumption made by Szilard and others, in which the macroscopic thermodynamic relation is assumed to hold even in a small system, namely, $`fX=k_\mathrm{B}T`$ for one molecule. In such studies, we find that it has been assumed that one can “fix” or quasi-statistically “control” the position of the piston, where the existence of the working bath (or ”working reservoir”)55  behind the piston is implicitly assumed. However, for a small system, the question as to how one can fix or control the piston in an operational or experimental manner should be carefully examinedKenPRE . The validity of the concept of the working bath is not evident in small systemsKenJPS . If one fixes the position of the piston, that is, if one eliminates the variability of the piston, work cannot be extracted out of the system because the “fixed” piston cannot transfer thermal energy into work. Such a framework is in contradiction to, and thus is no longer physically valid in, studies of Maxwell’s demon. This is precisely the reason why we introduced an external load to the piston instead of using the concepts of “pressure” or a “working bath”. It should also be noted that the movable piston satisfies the thermal property (Maxwell’s velocity distribution) merely through kinetic collision with a thermal particle. The evidence supports the physical consistency of the present description.

In the following paragraphs, we will present the role of thermal fluctuation as it essentially relates to the equation of state0807 . We have the distribution function in phase space, $`\rho =F^{}e^{fX/k_\mathrm{B}T}\theta (Xx_N)\theta (x_Nx_{N1})\mathrm{}\theta (x_2x_1),`$ where $`F^{}`$ and $`x_i`$ are a factor and the position of the $`i`$-th thermal particle ($`Xx_N\mathrm{}x_2x_1`$), respectively. By integrating all the spatial coordinates, except for that of the movable piston, $`X`$, we obtain the distribution function of the movable piston:

$$\rho \left(X\right)=\frac{1}{N!}\left(\frac{k_\mathrm{B}T}{f}\right)^{N+1}\mathrm{exp}\left\{\frac{fXNk_\mathrm{B}T\mathrm{ln}X}{k_\mathrm{B}T}\right\}.$$
(4)
Now, we define the coordinate $`X=X^{}`$ where the distribution function $`\rho (X)`$ is maximum. From Eq. (4), we get the relation $`X^{}=Nk_\mathrm{B}T/f`$. Thus, the conventional thermodynamic relation $`fX=Nk_\mathrm{B}T`$ still remains for the most probable position of the piston, $`X^{}`$, although it does not hold for the average position. The reason why the average position of the piston deviates from the most probable position, $`X^{}`$, is that thermal fluctuation of the piston is substantial for small systems, as seen by the relation $`\sqrt{(XX)^2}/X=1/\sqrt{N+1}.`$ Figure 4 shows the distribution function of the movable piston for three different quantities of thermal particles ($`N=1`$, $`10`$, and $`100`$), where the horizontal axis is normalized by $`N`$. As the number of particles decreases, the fluctuation from the most probable position, $`X=X^{}`$, and simultaneously the asymmetry of the distribution function, increase, while the most probable position of the piston continues to obey the conventional equation of state. Thus, the average equation of state deviates from the conventional one as the number of particles decreases. An interesting question relevant to any technique that reduces the thermal fluctuation of one molecule, as applied in probe microscopyYanagida , may arise here: Which position does the frozen position of a molecule correspond to, the most probable position or the average position obtained before the reduction in thermal fluctuation?

We can learn here the features of the concept of a “ quasi-static process” in extremely small systems. Szilard and others implicitly assumed that the expansion process could be described in the same quasi-static framework as its macroscopic counterpart. However, the condition of a quasi-static process requires a completely different framework at the microscopic level from that at the macroscopic levelKenPRE . A quasi-static condition always requires that a system should itinerate all the possible phase space sufficiently under the current parameters of the system. For this requirement, the system parameters should vary slowly enough during a quasi-static thermodynamic process. In our case, what should vary slowly is not the position of the piston but the force, $`f`$. In the expansion process, the piston moves rapidly following a slow change in the force. Thus, the conventional picture of a gradually expanding piston does not correspond to quasi-static expansion in small systems. Szilard and others had not anticipated this fact.

In this paper, we showed in the framework of Szilard’s model that the equation of state for a few molecules is different from that previously assumed by Szilard and others. In the thermodynamic limit, the equation of state is found to coincide with the conventional equation of state0809 . The result reveals that the macroscopic framework of thermodynamics has been incorrectly extended into the one molecule world, which has been assumed in the literature of Maxwell’s demon.

This counterexample to the basic assumption of Maxwell’s demon brings to light the need to develop the thermodynamics of small systemsKenJPS . This is important not only for Maxwell’s demon but also for proper interpretation and modeling of one-molecule experiments, in which the experimental data vary among observations. One-molecule thermodynamics is not a rescaled version of its macroscopic counterpartEvans . Filling the gap between one-molecule thermodynamics and macroscopic thermodynamics should be an interesting challenge from both philosophical and practical viewpoints of physics.

The author would like to thank K. Sekimoto for making him aware of the problem and introducing him to the method of Master equations, J. Lebowitz for useful comments on the pressure ensemble and other fruitful discussions, and S. Sasa, M. Kikuchi, M. Sano, Y. Hayakawa, S. Takagi, F. Jülicher, J. Prost, T. Kawakatsu for helpful comments. This work was supported in part by the Inamori Foundation, and by a Grand-in-Aid from JSPS (Nos. 60261575 and 17654082).


text:
# FUSE Observations of Degree-Scale Variations in Galactic Halo O VI

## 1 Introduction

Observations of absorption from O VI in the interstellar medium (ISM) towards early-type stars and extragalactic objects provide a useful means of studying the hot gaseous medium of galaxies. Because 113.9 eV are required to ionize O V to O VI, photoionization is unlikely produce the observed quantities of O VI in galactic environments (see Savage et al., 2000; Sembach et al., 2000). This fact, coupled with the high intrinsic abundance of O and the strength of the O VI absorption doublet at 1031.926 and 1037.617 Å, makes this ion a sensitive probe of the conditions of the hot ISM in the Milky Way and other nearby galaxies.

The Far Ultraviolet Spectroscopic Explorer (FUSE), launched in June 1999, is producing a large database of high-resolution ($`R15,00020,000`$) spectra covering the far-ultraviolet (FUV) wavelength range 905–1187 Å, which includes the important O VI doublet. FUSE can easily observe stars in the Galactic halo, stars in the Magellanic Clouds, and distant extragalactic sources with sufficient signal to noise to produce reliable O VI column density and line profile measurements. Early reports of FUSE observations include studies of interstellar O VI absorption in the Galactic halo (Savage et al., 2000; Oegerle et al., 2000), in high-velocity clouds (Sembach et al., 2000; Murphy et al., 2000), and in the Galaxy and the Large Magellanic Cloud towards the star Sk–67$`^{}`$05 (Friedman et al., 2000).

In this paper we discuss FUSE observations of interstellar O VI in the halo of the Milky Way as observed towards 12 early-type stars in the Large Magellanic Cloud \[LMC; $`(l,b)=(280.5,32.9)`$; $`d=50`$ kpc\] and 11 early-type stars in the Small Magellanic Cloud \[SMC; $`(l,b)=(302.8,44.3)`$; $`d=65`$ kpc\]. These observations have previously been used to study the O VI content of these galaxies (Howk et al., 2002; Hoopes et al., 2002). Our emphasis in this paper is on the degree of variability in the O VI column densities over the small angular scales ($`0\stackrel{}{\mathrm{.}}05\mathrm{\Delta }\theta 5\stackrel{}{\mathrm{.}}1`$) probed by these Magellanic Cloud observations.

The remainder of this paper is laid out as follows. In §2 we discuss the FUSE observations and the properties of the sight lines towards each of the Magellanic Clouds. We present the small-scale O VI measurements in §3, and we briefly discuss the kinematics of the O VI profiles along these directions in §4. We discuss the implications of these measurements for the the highly-ionized Galactic halo in §5 and summarize our results in §6.

## 2 FUSE Observations of the Magellanic Cloud Sight Lines

### 2.1 The LMC Sight Line

The first high-resolution UV absorption line measurements along sight lines towards extragalactic objects were were International Ultraviolet Explorer (IUE) observations of early-type stars in the LMC (Savage & de Boer, 1979). This high-latitude sight line ($`b33^{}`$) probes the thin interstellar disk of the Milky Way only along the first few hundred parsecs ($`z200`$ pc at a distance of $`d400`$ pc from the Sun). The first $`100`$ pc in this direction are contained within the boundaries of the local bubble, where there is relatively little material (Sfeir et al., 1999). Most of the LMC sight line probes the “thick disk” or “halo” of the Milky Way, and many of the interstellar clouds along this sight line show abundances reminiscent of other Galactic halo clouds (Welty et al., 1999; Friedman et al., 2000). Gas participating in Galactic rotation (Clemens, 1985) within the first $`5`$ kpc of the Sun ($`z2.5`$ kpc) in this direction will be at low negative velocities ($`5v_{\mathrm{LSR}}0`$ km s-1). More distant gas will be found at positive velocities, reaching $`v_{\mathrm{LSR}}+100`$ km s-1 at 25 kpc ($`z12.5`$ kpc).

Several clouds or cloud complexes have been identified along the sight line to the LMC, mostly in tracers of neutral gas. McGee, Newton, & Morton (1983) summarize H I 21-cm observations of the cloud complexes in the general direction of the LMC. They find several components towards the LMC. The main emission complexes are seen at $`v_{\mathrm{LSR}}=32.8,+0.2,+65.5,\mathrm{and}+131.7`$ km s-1, where the velocities are a spatial average over the face of the LMC, with spatially-averaged H I column densities of $`7.1\times 10^{18}`$, $`4.3\times 10^{20}`$, $`4.3\times 10^{18}`$, and $`6.0\times 10^{18}`$ atoms cm-2. The component nearest the LSR can be associated with local gas in the Galactic disk. The $`v_{\mathrm{LSR}}33`$ and $`+65`$ km s-1 components are referred to as intermediate velocity clouds (IVCs), while the $`+130`$ km s-1 cloud is cataloged as a high velocity cloud (HVC). These components can also be seen in metal absorption lines (e.g. Ferlet, Dennefeld, & Marice, 1985; Savage & de Boer, 1981), including moderately- and highly-ionized species (Welty et al., 1999; Wakker et al., 1998), and Danforth et al. (2002) have discussed their distribution in FUSE observations of metal absorption lines. Each of these aforementioned “components” is likely made up of several individual clouds that are unresolved in many datasets (Welty et al., 1999), and the IVC and HVC complexes probably reside in the Galactic halo (Welty et al., 1999).

The highly-ionized gas along the sight line to the LMC at velocities below those associated with the LMC itself likely traces material in the vertically-extended Galactic halo. IUE observations of C IV and Si IV along LMC sight lines have been presented by Sembach & Savage (1992), Savage et al. (1989), Savage & de Boer (1981), and others, while Bomans et al. (1996) and Wakker et al. (1998) have presented HST observations of Galactic C IV along LMC sight lines. We also note that Widmann et al. (1998) have derived O VI column densities for the Galactic halo seen towards the LMC using ORFEUS data. However, high-quality measurements of the highly-ionized Galactic material towards the LMC have been sparse.

In this work we discuss the Galactic O VI absorption as observed with FUSE towards 12 early-type stars in the LMC. The observations and analysis of the data used in this work are described in detail by Howk et al. (2002, hereafter Paper I) who make use of the observations to study the hot gas content of the LMC. All of the FUSE observations used for this work have a resolution of $`20`$ km s-1 (FWHM). The LMC stars used as background continuum sources in Paper I are all either WR stars or O stars with spectral types earlier than O7 and were chosen to have well-defined continua in the region of interstellar O VI $`\lambda 1031.926`$ Å absorption. The (6-0) P(3) and (6-0) R(4) transitions of molecular hydrogen can overlap Galactic O VI absorption towards the LMC given the rest wavelengths of these lines and the velocity structure of these sight lines. Models for molecular hydrogen transitions that overlap the Galactic O VI absorption have been divided out in the processing of the spectra. The exposure times for the LMC stars ranged from 3.6 to 33 ksec per star, yielding signal-to-noise ratios in the range $`12`$ to 25 per resolution element. The continuum placement, molecular hydrogen correction, and other details are described in Paper I.

Table 1 gives the Milky Way O VI column densities derived by integrating the apparent column density, $`N_a(v)`$, profiles (see Savage & Sembach, 1991) of the 1031.926 Å transition. The $`N_a(v)`$ profiles themselves are displayed in Figure 1. Table 1 gives the O VI column densities integrated over two ranges (all velocities are approximately LSR; see Paper I): the first over the velocity range $`|v|50`$ km s-1, and the second over the entire velocity range of the Milky Way (the full range is indicated in the last column of Table 1). As discussed in Paper I the separation between the Milky Way and LMC absorption is not always clear. In these cases we adopt an average separation point of $`+175`$ km s-1. The quoted column density errors for the columns integrated over the entire range include the effects of changing the upper velocity limit by $`\pm 20`$ km s-1, which is sufficient to cover the range of separations observed in cases where there is a clear break between the Milky Way and LMC absorption.

The full range of Milky Way velocities likely includes O VI associated with the IVCs and HVCs present along the LMC sight lines. Examining the strongest transitions of Fe II and O I in the FUSE bandpass (at 1144.938 Å and 1039.230 Å, respectively), we find that 11 of the 12 sight lines studied here exhibit IVC absorption, while at least 10 of the 12 contain HVC absorption.111Danforth et al. (2002) present the absorption line profiles of several species towards almost 100 stars in the Magellanic Clouds, including the strong Fe II 1144.938 Å and O I 1039.230 Å transitions. We direct the reader to that work for better views of the sight line velocity structure along the lines of sight studied in this work. None of the sight lines is simultaneously free of absorption from IVC and HVC material (i.e., absorption from one or the other is seen along all of the sight lines). The column densities derived by integrating over $`|v|50`$ km s-1 should exclude much of the IVC and HVC absorption.

### 2.2 The SMC Sight Line

The higher-latitude ($`b44^{}`$) sight line to the SMC is more heavily biased towards low-density Galactic halo gas than that towards the LMC. Towards the SMC, gas within $`10`$ kpc ($`z7`$ kpc) that is participating in Galactic rotation (Clemens, 1985) would be found at negative velocities, and gas within $`20`$ kpc of the Sun ($`z14`$ kpc) should be found at velocities $`|v_{\mathrm{LSR}}|40`$ km s-1.

Similar to the case for the LMC, the sight line towards the SMC contains a number of low-ionization Galactic gas clouds or complexes spread over $`65`$ km s-1 (Welty et al., 1997; Mallouris et al., 2001). Material associated with the disk and low halo of the Galaxy is seen over the range $`25v_{\mathrm{LSR}}+40`$ km s-1. High-velocity Galactic material, such as that seen towards the LMC, would be blended with SMC material for velocities $`v_{\mathrm{LSR}}+60`$ km s-1.

Observations of the highly-ionized material in the direction of the SMC with IUE are discussed by Fitzpatrick & Savage (1983) and Fitzpatrick (1984). More recently, Koenigsberger et al. (2001) and Hoopes et al. (2001) have reported on HST and FUSE observations of the sight line towards HD 5980 in the SMC. These papers include comments on the highly-ionized species C IV, Si IV, and O VI in the Galactic halo.

In this work we make use of FUSE observations of Galactic O VI towards 11 early-type stars in the SMC. The observations and analysis of the data are described in Hoopes et al. (2002, hereafter Paper II), who use these observations to study the hot gas content of the SMC. The principal properties of the data are the same as those described above for the LMC observations.

Table 2 gives the Milky Way O VI column densities derived from an integration of the apparent column density profiles of the 1032.926 Å transition. The O VI $`N_a(v)`$ profiles for the SMC sight lines are shown in Figure 2. For the SMC sight lines, contamination of the O VI absorption profiles by H2 is not as significant as for the LMC sight lines (see Paper II, ). The H2 has not been removed before integrating the $`N_a(v)`$ profiles. The separation between Galactic and SMC O VI absorption is much cleaner along the SMC sight lines than the LMC sight lines. For this reason, the error budget includes the effects of changing the upper integration limit (last column of Table 2) by only $`\pm 5`$ km s-1.

## 3 Total O VI Column Densities and Degree-Scale Variations

The column densities quoted in Tables 1 and 2 can be seen to vary significantly from sight line to sight line. Towards the LMC the full range of O VI columns spans a factor of $`2.8`$ range when integrated over the full Milky Way velocity limits indicated in the table, or a factor of $`4.2`$ over $`|v|50`$ km s-1. The SMC sight lines show a factor of $`4.2`$ full range as well.

Figure 3 shows the absolute value of the logarithmic differences in column densities for each pairing of sight lines in our sample, i.e., $`|\mathrm{\Delta }\mathrm{log}N_{jk}(\text{}\text{VI})|=|\mathrm{log}N_j(\text{}\text{VI})\mathrm{log}N_k(\text{}\text{VI})|`$ for each pair of LMC sight lines $`j`$ and $`k`$, versus the angular separation of the pair, $`\mathrm{\Delta }\theta _{jk}`$. The top panel of this figure shows the difference in O VI column density integrated over the full Milky Way velocity range for each pairing of sight lines, while the bottom panel shows the column density differences when integrating the apparent column density profiles over $`|v|50`$ km s-1. Figure 4 shows a similar plot for the SMC sight lines.

Figures 3 and 4 demonstrate that there are significant column density variations on all of the angular scales probed by our observations. The smallest separations of LMC sight lines are of order $`\mathrm{\Delta }\theta _{jk}0\stackrel{}{\mathrm{.}}5`$, while the smallest angular scales probed by the SMC sight lines are $`\mathrm{\Delta }\theta _{jk}0\stackrel{}{\mathrm{.}}05`$.

Table 3 summarizes the statistics of the O VI column density measurements towards the LMC and SMC. For the LMC direction, the results are shown for integrations over the full Milky Way velocity range and over $`|v|50`$ km s-1, corresponding to the analogous entries in Table 1. Because of the presence of IVCs and HVCs along the LMC sight lines (§2.1), the variation in the upper panel of Figure 3 trace variations in both the low-velocity Milky Way O VI and in the O VI associated with these anomalous-velocity clouds. The IVC and HVC absorption along the LMC sight lines is to a large extent excluded in the integration over the velocity range $`|v|50`$ km s-1. The variations over this low-velocity material are on average larger than those integrated over the full velocity range (see Table 3). Indeed, it is interesting to note that the statistics of the LMC and SMC directions are almost indistinguishable if only the low-velocity gas towards the LMC is included.

The standard deviations of the individual O VI column density measurements are a significant fraction of the mean for the LMC and SMC sight lines. Table 3 shows that the standard deviation of each sample is $`27\%`$ to $`38\%`$ of the mean. For comparison, the standard deviation of the $`N(\text{}\text{VI})`$ determinations towards 11 bright AGNs in the study of Savage et al. (2000) is $`41\%`$ of the mean (with $`\mathrm{log}N(\text{}\text{VI})=14.50`$). Thus, the LMC/SMC measurements that each sample sight lines within a few degrees of one another show a degree of (fractional) variation that is similar to that seen along sight lines distributed across the entire celestial sphere. Indeed, compared with the uncertainties in the individual measurements, the variation between the LMC and SMC sight lines is more significant than that in the Savage et al. (2000) dataset (as judged by the values of $`\chi ^2`$ derived by comparing the data with the mean within each sample).

Figures 3 and 4 show that significant variations in the O VI column densities of the Galactic halo ($`|\mathrm{log}\mathrm{\Delta }N_{jk}(\text{}\text{VI})|>0`$ at greater than $`2\sigma `$ significance) can be found for sight lines over the entire range of angular scales probed by our observations. This includes pairs of sight lines separated by $`1^{}`$, many of which which show greater than 0.1 dex ($`25\%`$) variations.

To demonstrate the degree of variation of and any systematics in the distribution of column densities, Figures 5 and 6 show “maps” of the derived Galactic O VI column densities towards the LMC and SMC, respectively. The positions of the probe stars are marked with circles whose radii are linearly proportional to the column densities summarized in Tables 1 and 2. Scales showing the minimum and maximum columns are given. With only 12 and 11 sight lines, the maps are too sparsely populated to say whether or not coherent structures are visible in our data.

Figure 7 shows the angular autocorrelation function (ACF), $`w(\theta )`$, for the O VI column densities towards the LMC and SMC. The ACF was calculated in a manner analagous to the studies of the X-ray sky by Chen et al. (1994) and Kuntz et al. (2001). We define the ACF as:

$$w(\theta )\frac{\underset{i}{}\underset{ji}{}(N_iN)(N_jN)(\sigma _i\sigma _j)^{1/2}}{N^2\underset{i}{}\underset{ji}{}(\sigma _i\sigma _j)^{1/2}}.$$
(1)
The summation is calculated for each pair of sight lines $`i`$ and $`j`$ separated by an angle $`\theta `$ on the sky. In this expression $`N_i`$ and $`N_j`$ are the O VI column densities and $`\sigma _i`$ and $`\sigma _j`$ the column density uncertainties for the sight lines, while $`N`$ is the average column density for the entire sample of 12 or 11 stars. The $`w(\theta )`$ data displayed in Figure 7 have been summed over $`0\stackrel{}{\mathrm{.}}6`$ bins for the LMC dataset; the summation for the SMC dataset is over $`0\stackrel{}{\mathrm{.}}3`$ bins for separations smaller than $`2^{}`$ and $`0\stackrel{}{\mathrm{.}}6`$ bins for larger angular separations. The uncertainties in $`w(\theta )`$ represent $`95\%`$ confidence limits. These were determined through 1000 simulations wherein the column densities were randomly shuffled among the observed directions and $`w(\theta )`$ rederived, similar to the method described by Carrera et al. (1991). The confidence limits shown in Figure 7 include $`95\%`$ of the simulated $`w(\theta )`$ values. This approach to estimating uncertainties has the advantage that the simulated data are drawn from the same population as the real data, and hence no intrinsic column density distribution function need be assumed. To remove any systematic effects from the ACF, we have subtracted the mean simulated $`w(\theta )`$ from the real values, though the mean simulated values were always very small compared with the uncertainties.

Figure 7 shows that there is little evidence for significant power at any scale in the ACF for the O VI column densities towards the LMC and SMC. Indeed, those points that seem to be different from $`w(\theta )=0`$ are not robust to changes in the bin spacing. When binned over the same angular scale, the LMC and SMC $`w(\theta )`$ values are completely consistent over all scales probed.

There is no evidence for significantly greater power in the ACFs on small angular scales compared with large scales along these two sight lines. We caution that the small number of sight lines makes the uncertainties large over all scales. Furthermore, the ACFs are insensitive to structure on scales smaller than the adopted bin size for each sight line. Thus, both larger numbers of sight lines and finer sampling would be needed to probe the smallest structures. However, our data could be used with further measurements of O VI towards AGNs to constrain the largest-scale O VI structures in the sky.

## 4 Kinematics of Gas in the Galactic Halo

The kinematic profiles of the O VI absorption lines towards the LMC and SMC could, in principle, yield information on the distribution of gas in the Galactic halo assuming it corotates with the underlying thin disk and adopting an appropriate rotation curve (see, e.g., Savage & de Boer, 1981; Savage, Sembach, & Lu, 1997). Examining the distribution of apparent column density versus velocity for LMC and SMC sight lines in this manner can yield estimates of the turbulent velocity, midplane density, and scale height of the absorbing layer.

There are several caveats to this approach for the particular sight lines studied here, however. Figure 8 shows a comparison of the apparent column density, $`N_a(v)`$, profiles of O VI and Fe II (derived from the 1125.448 Å transition) for the sight lines towards Sk–67$`^{}`$211 and Sk–68$`^{}`$80 in the LMC. The Fe II absorption mostly traces gas associated with the neutral ISM (i.e., H I-bearing gas). Also shown is the Clemens (1985) rotation curve for the direction of the LMC.222The reader should be aware that the velocities presented in the O VI and Fe II profiles are expected to be accurate to $`\pm 10`$ to 15 km s-1 when comparing the observed profiles with the displayed rotation curve. These plots demonstrate the difficulty in analyzing the O VI kinematics in this direction. Although the sharp Fe II absorption in the main Milky Way component suggests the Fe II-bearing material is confined to relatively small distances from the midplane (consistent with the results of Edgar & Savage, 1989), determining the extent of the Milky Way halo O VI absorption is difficult given the overlap with IVC, HVC, and LMC material.

Furthermore, examining Figure 8 shows that the O VI $`N_a(v)`$ profiles for these two sight lines are significantly different at even low velocities, i.e., velocities that should be free from contamination by IVC and HVC contamination. The Sk–67$`^{}`$211 sight line shows much larger apparent column densities over all velocities $`v+50`$ km s-1 than the sight line towards Sk–68$`^{}`$80. Furthermore, the shape of the profiles is quite different over the lowest velocities, with the Sk–67$`^{}`$211 profile exhibiting a sharp peak at low positive velocities and significant absorption at negative velocities (neither of which are easily identified along the Sk–68$`^{}`$80 sight line). Given the differences seen between the two LMC sight lines examined above, which are separated by $`1\stackrel{}{\mathrm{.}}5`$ on the sky, models of smooth, plane-parallel halo material corotating with the underlying Galactic disk do not appropriately describe the complicated nature of the observed O VI profiles.

Figure 9 shows the O VI and Fe II profiles for the stars AV 232 and AV 321 in the SMC. These stars are separated by only $`0\stackrel{}{\mathrm{.}}3`$ on the sky and again significant variations between the sight lines are seen at low velocity in the O VI $`N_a(v)`$ profiles. The two $`N_a(v)`$ profiles exhibit differences not only in the peak and total O VI column densities, but also in the velocity distribution of the gas. While there is no clear evidence for IVC absorption in these directions (although Welty et al., 1997, have found weak IVC absorption in the direction of the SMC star Sk 108), absorption associated with the SMC begins to overlap the Milky Way absorption at velocities nearing +100 km s-1. The low-velocity tail of the broad O VI distribution in the SMC likely overlaps the Galactic O VI at velocities below this.

While the derivation of reliable estimates of scale height, midplane densities, and turbulent velocities using the O VI absorption in these directions may give ambiguous results, examination of the kinematic profiles of Milky Way O VI seen towards the Magellanic Clouds do lead to several important conclusions regarding the nature of highly-ionized Galactic halo gas. That the individual absorbing components are seen more distinctly in the Fe II profiles (Figure 8) implies the O VI absorption contains more and/or broader component complexes than the Fe II. This is consistent with our expectations given the temperature of the gas the two species are thought to trace.

In the direction of the SMC, the Milky Way O VI absorption is clearly broader than the Fe II absorption. This suggests that the O VI arises in a layer much more extended than that traced by Fe II, which has a vertical extension similar to that of H I (Edgar & Savage, 1989). Towards the LMC, the fraction of the gas arising at velocities associated with IVC and HVC absorption is quite different for the O VI and Fe II profiles. A larger fraction of the O VI absorption is found at high velocities. Given that the IVC and HVC material towards the LMC likely resides far from the Galactic plane, this implies a greater fraction of the O VI is found in the Galactic halo than Fe II.

This all suggests that while the O VI-bearing medium is quite inhomogeneous, the current observations certainly suggest the O VI layer of the Galaxy is significantly different than the Fe II-bearing layer. O VI is found at preferentially larger distances from the Galactic plane than is the Fe II. Stated another way, the effective scale height of the O VI-bearing gas is larger than that of the Fe II-bearing material.

## 5 Discussion

We have presented measurements of O VI in the Galactic halo observed towards early-type stars in the LMC and the SMC. These show significant O VI column density variations on degree scales (and smaller) over velocity ranges corresponding to gas in the Galactic disk and halo. The dispersion of halo O VI measurements is $`27\%`$ to $`38\%`$ of the averages in these directions. The ratio of the highest to lowest O VI column densities ranges from 2.8 to 4.2, with the latter applying to the low-velocity gas towards the LMC and the entire Milky Way range towards the SMC.

Several properties of O VI in the Galactic halo can be inferred from the properties of these sight lines towards the Magellanic Clouds. The most fundamental inference to be drawn from our observations is that the O VI-bearing medium of the Galactic halo is better described by spatially-patchy distribution than by a smooth layer (e.g., like that discussed by Spitzer, 1956). That the medium is patchy favors models in which the O VI resides in clouds of $`3\times 10^5`$ K material or in the interfaces between cooler ($`<10^5`$ K) clouds and a hotter ($`10^6`$ K) medium. The evidence for this conclusion comes both from the observed column density variations and the kinematic differences between sight lines (see Figures 1 and 2).

The greater velocity breadth of O VI than Fe II in the direction of the SMC and the much larger contribution of intermediate- and high-velocity halo clouds to the O VI profiles than to the Fe II profiles towards the LMC imply that the Galactic O VI in these directions is extended further from the plane than the Fe II. Effectively, these inferences imply the scale height of Galactic O VI is larger than that of Fe II. Although the models of smoothly-distributed, exponentially-stratified gas layers assumed for determining the properties of the O VI layer seem inappropriate (§4), our data do require that the Galactic O VI and Fe II be distributed in very different ways.

That the Galactic halo O VI-bearing medium is poorly described by a smoothly distributed layer of gas is a conclusion one must also draw from the Savage et al. (2000) observations of O VI absorption towards extragalactic FUV sources. However, the current observations show that the distribution of O VI within the halo must have significant variations over much smaller angular, and hence physical, scales. To what physical scales might the observed variations correspond? The distribution of O VI absorption seems to be more tightly concentrated towards the Galactic midplane than that of C IV and Si IV (Savage et al., 2000). The latter ions have distributions with best fit scale heights of $`4.4`$ and 5.1 kpc (Savage et al., 1997) – under the assumption that an exponentially-stratified layer with some degree of irregularities (patchiness) is a good description of the intrinsic distribution. If we assume that the observed O VI absorption in the Galactic halo towards the Magellanic Clouds resides at heights lower than $`z5`$ kpc (i.e., within one Si IV scale height), then the $`0\stackrel{}{\mathrm{.}}55\stackrel{}{\mathrm{.}}0`$ scales probed by our LMC observations correspond to linear distances of $`80800`$ pc while the $`0\stackrel{}{\mathrm{.}}053\stackrel{}{\mathrm{.}}3`$ scales probed by the SMC measurements correspond to distances of $`6400`$ pc. Much of the gas is likely to be significantly closer to the Sun than this limit, making these estimated size scales upper limits. In particular, if the effective scale height of the O VI distribution is $`h_z2.7`$ kpc (as derived by Savage et al., 2000), then $`85\%`$ of the Galactic O VI is expected to be at heights smaller than the adopted 5 kpc.

If we assume that the O VI variations are due to clouds at this $`5`$ kpc distance from the plane with angular sizes approximately equal to the transverse separation of the pairing of stars observed (giving physical sizes of $`80800`$ pc towards the LMC and $`5400`$ pc towards the SMC), we can estimate the required densities of O VI for the observed column density differences between the sight lines. The mean density derived under these assumptions is $`n_{\mathrm{OVI}}1.5\times 10^7`$ cm-3 for the LMC sight lines with significant ($`2\sigma `$) column density differences; the full range of densities is $`(0.4\mathrm{to}4.4)\times 10^7`$ cm-3. Towards the SMC the mean density derived in this manner is $`n_{\mathrm{OVI}}2.7\times 10^7`$ cm-3 with a full range of $`(0.6\mathrm{to}27)\times 10^7`$ cm-3. These values give estimated particle densities of $`n_\mathrm{H}1.4\times 10^3`$ cm-3 and $`2.5\times 10^3`$ cm-3 for the LMC and SMC directions, respectively, assuming the abundance of oxygen relative to hydrogen is $`5.5\times 10^4`$ (Holweger, 2001) and the fraction of oxygen in the form of O+5 is $`0.2`$ (Sutherland & Dopita, 1993).

The average density estimates can be used to estimate the pressure of the O VI-bearing gas. For $`P/k=nT=2.3n_\mathrm{H}T`$, where the factor of 2.3 allows for the total number of particles in a fully ionized gas containing $`10\%`$ helium, we obtain $`P/k10^3`$ to $`1.7\times 10^3`$ K cm-3 for the directions to the LMC and SMC, respectively. These estimates adopt the average densities derived in the previous paragraph, and we have assumed $`T=3\times 10^5`$ K, the temperature where O VI peaks in abundance in collisional ionization equilibrium. The densities and pressures derived in this way will be higher for clouds closer than the assumed 5 kpc height from the plane of the Galaxy.

This approach to estimating densities and pressures is overly simplistic, and we caution the reader not to put too much emphasis on these idealized calculations – neither the true distances to nor the morphology of the absorbing material is known. If the morphology of the O VI-bearing ISM is more sheet- or filament-like, then much of the variation in column density will be the result of the projection of the individual absorbing components along the line of sight (see, e.g., Heiles, 1997). Indeed, the interpretation of the absorption arising in individual clouds at a common distance is highly unsatisfactory given the apparent lack of any any preferred scales for the observed variability.

We suggest that the observations presented in this work, which show degree-scale column density and kinematic variations in the O VI-bearing medium of the Galactic halo, fit readily into pictures of the Galactic halo wherein the O VI is contained in complicated cloud- or sheet-like distributions of material. In particular, models that place the O VI gas in interfaces between more neutral clouds or sheets and a more pervasive very hot ($`10^6`$ K) ISM would seem to allow for the observed variations so long as there are not so many such interfaces along typical paths through the halo that the expected variations along random adjacent sight lines is expected to be small.

The morphological information contained in the absorption line datasets presented here is incomplete, probing only random pencil beams through the highly-ionized halo. With our sparsely-populated “maps” it is difficult to discern the structure of the halo O VI. Our observations do give the first hint of the complexity of the small-scale distribution of hot, highly-ionized material within the halo. The structures giving rise to the O VI column density variations demonstrated in this work will be accessible to future O VI emission-line imaging programs. Indeed, that we find significant variations on the smallest angular scales probed suggests there is a wealth of structure present in the highly-ionized sky; imaging data may reveal how this highly-ionized material is related to other phases of the ISM.

We also note that the degree of structure demonstrated by our O VI absorption measurements has important implications for understanding current observations of diffuse O VI emission. The comparison of O VI emission line intensities with the O VI column density along the same line of sight can yield estimates of the physical properties of the O VI-bearing gas, e.g., the density and pressure (Shelton et al., 2001; Dixon et al., 2001). However, the determination of the appropriate column density for comparison with an emission line measurement is difficult. Our observations (e.g., Figures 3 and 4) suggest that O VI absorption line measurements along sight lines even a few arcminutes from the direction of an emission line measurement can only be assumed to yield column densities within a factor of approximately two of the value appriopriate for the emitting region. The large degree of variation in the absorption line column densities along closely-spaced sight lines through the Galactic halo must be considered when assessing the systematic uncertainties in comparing emission and absorption line measurements of O VI. However, one may look at this in another way: our analysis suggests that sight lines several degrees apart give column densities that are expected, on average, to give results that are no less consistent than sight lines separated by less than a degree. That is, unless the opportunity exists to observe O VI emission and absorption along essentially the same sight line, the determination of an O VI column density several degrees away from an emission line measurement will provide no worse an indicator than an O VI column density determined for a sight line tens of arcminutes away.

## 6 Summary

We have presented observations of Galactic halo O VI absorption towards 12 and 11 early type stars in the Large and Small Magellanic Clouds, respectively. These observations reveal strong column density variations on degree scales. The principle results of this study are as follows.

1. The Galactic halo O VI probed by our observations shows significant structure on all of the scales probed by our observations ($`0\stackrel{}{\mathrm{.}}055\stackrel{}{\mathrm{.}}0`$). There seems to be no preferred angular scale for O VI variations over this range. The transverse size-scales probed by our observations are likely $`6`$ to 800 pc.
2. The Galactic O VI in the directions of the Magellanic Clouds is very patchy, with the standard deviations of the measurements being equivalent to $`26\%`$ to $`38\%`$ of the average values in these directions. The maximum variations correspond to factors of $`2.9`$ to 4.2.
3. Galactic O VI is distributed in a very different manner than Fe II. The O VI-bearing layer of the Galaxy is significantly more extended than the Fe II-bearing layer, which has a vertical extent similar to that of H I. However, the existence of large variations in the O VI column densities over very small angular scales implies that models of a smoothly-varying, exponentially-stratified halo are poor descriptions of the true distribution of highly-ionized gas. The meaning of detailed quantities derived by comparing such models with the observed O VI kinematics along a given sight line (or sight lines) is ambiguous at best.
4. The existence of large column density variations over small angular scales can be used to constrain the origin of Galactic halo O VI ions. Models in which the O VI-bearing medium is composed of complicated cloud- or sheet-like distributions of material (particularly narrow interfaces) should be considered more likely than smoothly-distributed layers of highly-ionized material in light of our observations.

This work is based on data obtained for the Guaranteed Time Team by the NASA-CNES-CSA FUSE mission operated by the Johns Hopkins University. Financial support to U. S. participants has been provided by NASA contract NAS5-32985. JCH and KRS recognize support from NASA Long Term Space Astrophysics grant NAG5-3485 through the Johns Hopkins University.