text:
# Simulation of the Sedimentation of a Falling Oblate

## I Introduction

The way in which objects fall to the ground has been studied since antiquity. Objects were thought to return to “their natural” places by the ancient Greeks. Newton showed that the bodies fall on earth driven by a constant acceleration. But despite gravity’s undeniable attraction, not all falling objects travel downwards on straight trajectories. The tree leaves flutter to the ground in the autumn, exhibiting a complex motion and refusing to follow the shortest path.

A deep understanding of the motion of falling objects in a fluid is of great technical importance, and has been investigated in a variety of contexts, including meteorology , aircraft stability , power generation , chemical engineering , etc. Also Newton observed the complex motion of objects falling in both air and water . This phenomenon was also studied by Maxwell, who discussed the motion of a falling paper strip .

In the nineties, Aref and Jones  found through numerical solutions of Kirchhoff’s equations that the trajectory for an object moving through an incompressible inviscid and irrotational fluid, is chaotic. Tanabe and Kaneko , using a phenomenological model for the falling of a one-dimensional (1D) piece of paper, including lift and dynamical viscosity, but neglecting the inertia of the fluid, describe five falling regimes. Two of them are chaotic. In 1997 Stuart B. Field et al. , investigated experimentally the behavior of falling disks in a fluid, and identified different dynamical regimes as function of the moment of inertia and the Reynold’s number. They obtained experimental evidence for chaotic intermittency . In 1998 Andrew Belmonte, et al. , in an experiment with thin flat strips falling through a fluid, observed only two motions: side to side oscillation (flutter) and end-over-end rotation (tumbling). They proposed a phenomenological model including inertial drag and lift which reproduces this motion, and yields the Froude similarity, which describes the transition from flutter to tumble regime. Mahadevan et al , in 1999 made an experiment of dropping horizontal cards of thickness $`d`$ and width $`w`$, showing that the tumbling frecuency $`\mathrm{\Omega }`$ scales as $`\mathrm{\Omega }d^{1/2}w^1`$, consistent with a dimensional argument that balances the drag against gravity.

Given the difficulties to study this problem theoretically and experimentally, we took a computational approach simulating the falling of one oblate ellipsoid in a viscous fluid in a three dimensional container. An oblate is an ellipsoid for which the two largest principal radia are equal. We organize the paper in the following manner. In section 2 we give an review over the model that we use. In Sec. 3A we describe the main features of the falling oblate. In Sec. 3B-D we present the results for the change in the initial height, the oblate’s aspect-ratio and the dynamical viscosity, for the steady-falling regime. In sec. 3E we show the periodic behavior and compare to the results of reference . In Sec. 3F the sensitivity to tiny variations in the oblate’s initial orientation is presented in the chaotic regime. In sec. 3G the parameter phase space is sketched and compared to the results of reference . Section 4 summarizes present results and discusses possible further applications.

## II  Model

The general idea, proposed by Fogelson and Peskin , is to work with a simple grid for the resolution of the fluid motion at all times and represent the particles not as boundary conditions to the fluid, but by a volume force term or Lagrange multipliers in the Navier Stokes equations.

This technique was developed in the work of Kuusela, et al , Wachmann, et al  and Hoefler et al. , , and employs a numerical solver for the dynamical simulation of three-dimensional rigid particles in a Newtonian fluid, bounded by a rectangular container.

The motion of the fluid is described by the dimensionless Navier-Stokes equations:

$$\frac{\stackrel{}{v}}{t}+(\stackrel{}{v})\stackrel{}{v}=p+\frac{1}{Re}^2\stackrel{}{v}+\stackrel{}{f}$$
(1)
Here $`p`$ and $`\stackrel{}{v}`$ are the pressure and the velocity of the fluid, respectively, and $`\stackrel{}{f}`$ is an external force.

The Reynolds number $`Re`$ is defined as

$$Re:=\frac{UD\rho _f}{\nu }$$
(2)
where $`U`$ is the mean vertical oblate velocity, $`D`$ a characteristic length that in our case is the largest oblate’s diameter, $`\rho _f`$ the density and $`\nu `$ the dynamical viscosity of the fluid.

For an incompressible fluid, the continuity equation:

$$\frac{\rho _f}{t}+(\rho _f\stackrel{}{v})=0$$
(3)
reduces to

$$\stackrel{}{v}=0$$
(4)
Equation (1) is discretized on a regular, marker-and-cell mesh to second order precision in space. For the time stepping, we employ an operator-splitting-technique which is explicit and accurate to first order. The suspending fluid is subjected to no-slip boundary conditions at the surface of the suspended particles. More details of the solution procedure are presented in , , .

The geometry of the oblate ellipsoid is characterized by $`\mathrm{\Delta }r`$ its aspect-ratio defined as the ratio of the smallest radius over the largest one.

An oblate is represented by a rigid template connected to fluid tracer particles, which are moving on the trajectories of the adjacent fluid. The connection is made by using the body force term, in the Navier-Stokes equations, as constraints on the fluid such to describe the oblate.

The force density $`\stackrel{}{f}^c`$, is chosen elastic with a spring constant that, and guarantees that the elongation remains small against the grid spacing at all times , and it is zero in the exterior of the region outside the oblate. We can define this force density $`\stackrel{}{f}^c`$ as:

$$\stackrel{}{f}^c=f^c(\stackrel{}{x_{ij}}+\stackrel{}{ϵ}(\stackrel{}{x_{ij}}))=k\stackrel{}{ϵ}(\stackrel{}{x_{ij}})$$
(5)
where $`\stackrel{}{x_{ij}}`$ is the displacement field of the separation between the markers $`i`$ and their corresponding reference point $`j`$. The stiffness constant $`k`$, must be chosen large enough so that $`\stackrel{}{ϵ}(\stackrel{}{x_{ij}})h`$, ($`h`$ size grid), holds for all iterations.

In general the displacement field $`\stackrel{}{ϵ}(\stackrel{}{x_{ij}})`$ is defined as:

$$\stackrel{}{ϵ_i}(\stackrel{}{x_{ij}}^m)=\stackrel{}{x_{ij}}^m\stackrel{}{x_{ij}}^r$$
(6)
The vector $`\stackrel{}{x_{ij}}^m`$ is the position of a fluid tracer, whose motion is determined by the fluid local velocity, i.e.,

$$\dot{\stackrel{}{x_{ij}}^m}=\stackrel{}{u}(\stackrel{}{x}_{ij}^m)$$
(7)
The $`\stackrel{}{x_{ij}}^r`$ are the reference points associated to a template having the shape of the physical particle:

$$\stackrel{}{x_{ij}}^r=\stackrel{}{x}_i+O_i(t).\stackrel{}{r}_{ij}$$
(8)
Here $`\stackrel{}{x}_i`$ is the center of mass of the template, $`O_i(t)`$ is the rotation matrix that describes the present orientation of the oblate and $`\stackrel{}{r_{ij}}`$ denote the initial position of the reference points with respect to the center of mass. For the quaternion formulation of the rotation, we use the techniques described in ref. .

A velocity-Verlet integrator  serves to integrate the equations of motion for the translation and a Gear-predictor integrator  for the rotation on the template:

$$\stackrel{}{F}=Mg\widehat{j}+\rho _fVg\widehat{j}+\underset{i}{}\stackrel{}{f}_i^c$$
(9)
where $`\widehat{j}`$ is the unit vector along the vertical.

$$\stackrel{}{T}=\underset{i}{}(\stackrel{}{x_i}\stackrel{}{x_{cm}})\times \stackrel{}{f}_i^c$$
(10)
with respect to the template’s center of mass $`\stackrel{}{r}_{cm}`$.

The equations of motion of the particle template are:

$$\dot{\stackrel{}{U}}=\frac{\stackrel{}{F}}{M}$$
(11)
and

$$\dot{\stackrel{}{\mathrm{\Omega }}}=\frac{\stackrel{}{T}}{I}$$
(12)
where $`M`$ is the mass of the template particle; $`\stackrel{}{U}`$ and $`\stackrel{}{\mathrm{\Omega }}`$ are the linear and angular velocities of the template particle, respectively; $`I`$ is the moment of inertia; and $`\stackrel{}{T}`$ is the torque, , .

The boundary conditions at the container wall are zero for the normal velocity component of the fluid and no-slip condition for the tangential component, because the walls are assumed impenetrable, , . The interaction between the oblate and the walls is defined through a contact force, , where the walls are treated as a particle with infinite mass and infinite radius.

## III Results

For all our simulations, and in order to reduce the parameter space we use $`\rho _{fluid}=1\frac{g}{cm^3}`$ and $`\rho _{oblate}=3.5\frac{g}{cm^3}`$.

We found three basically different motions in our simulations: steady-falling, side-to-side periodic-oscillation known as ’flutter’ ref., and a chaotic motion as shown in fig. 2. The above phenomenology can be compared to the work of ref., for the case of dropping disks. In general, the trajectories depend strongly on the initial conditions and the properties of the system (oblate’s orientation $`\mathrm{\Theta }_o`$, dynamical viscosity $`\mu `$ and the oblate’s aspect-ratio $`\mathrm{\Delta }r`$, etc). We don’t find the tumbling motion described in the above references.

### A Phenomenology of the Steady-Falling Oblate

For the velocity of the center of mass, the vertical component decreases when the oblate approaches the container bottom and shows a damped wavering, with an amplitude that decreases with time (fig. 3 top).

The vertical trajectory is composed of succesive turning points, that correspond to the points where the trajectory changes the sign of the rate of change of the vertical velocity component marked in fig. 3(bottom) by the horizontal lines L2, … L5. We also see, that from $`L_1`$ to $`L_6`$ the amplitude for the trajectory and the velocity decrease.

The horizontal components of the velocity obey a different behavior as the vertical ones, and in general they have non-regular oscillations as seen in fig.4.

The oblate orientation is described through the three rotational degrees of freedom, around the center of mass. We present the time evolution of the angle between the oblate’s normal and the vertical direction that we will call vertical orientation fig. 5 (top). $`\mathrm{\Theta }=0`$ implies that the oblate’s principal axis will be parallel to the container’s horizontal fig. 5 (bottom). At the beginning of the movement, there is a larger angular change of the oblate’s normal $`\mathrm{\Delta }\mathrm{\Theta }`$ fig. 5(bottom), which is characterized by a larger peak-to-peak amplitude $`\mathrm{\Theta }_{pp}`$. In fig. 5 the definitions we illustrate the definition of the peak-to-peak amplitude as the distance between succesive turning points, which decreases while the oblate sinks. In the steady-falling regime, both the peak-to-peak amplitude $`\mathrm{\Theta }_{pp}`$ and the change $`\mathrm{\Delta }\mathrm{\Theta }=\mathrm{\Theta }_f\mathrm{\Theta }_i`$ , become smaller from the top towards the container’s bottom. The oblate tends to align its major axis along the vertical , presenting the lowest resistance to its descent in the fluid, and acquiring a limit vertical velocity fig. 3 (top).

The oblate generates shedding vortices in the fluid along its falling trajectory, an example is shown in figure 6 (right), that shows the velocity fluid field around the oblate, and the vortex is localized just in the top region above the oblate, and where the oblate has experienced the larger angular change $`\mathrm{\Delta }\mathrm{\Theta }`$ as shown in fig. 5 (top). The Reynolds number calculated from the oblate’s diameter and terminal velocity is $`Re=128`$. We point out that the vortex structure is obtained also in the work of Belmonte et al,  where a shedding vortex created by the zigzag motion of the falling strip is seen in the left figure 6.

### B Steady-Falling Oblate: Change in the Initial Height.

We choose three different initial heights fig. 7I, (fixing the rest of parameters). If we superpose the trajectories in fig. 7I, it can be shown that there is no variation in the wavelength or in the peak-to-peak amplitude. For all the heights, the trajectories generated are in good agreement with a damped harmonic oscillation.

For the range of falling heights used, the oblate attains the same terminal vertical velocity fig. 7II, ($`3\frac{cm}{s}`$), its magnitude is independent on the falling height. The vertical velocity finally, reaches a stable state (uniform and linear motion), after the same time ($`20s`$ fig 7II). We also see that the vertical velocity suddenly becomes zero when the oblate touches the bottom.

We see in fig. 7 III(c-a), an increase in the final angle with respect to the initial height. For $`h_o=56`$ the smaller height, we obtain $`85^o`$ and for larger one $`h_o=96`$, we have $`95^o`$. For the larger height the oblate is still in a transitory state before it arrives to its final angle fig. 7III.

### C  Steady-falling oblate: Change in the dynamical viscosity.

The oblate starts its swinging motion with a given initial orientation $`\theta _0=26.6^o`$. It glides downwards and to the side acquiring some amplitude, while the dynamical viscosity $`\mu `$ acts reducing the subsequent amplitudes of oscillation (fig. 8I). As we increase the dynamical viscosity from $`(a)\mu =0.025`$ to $`(c)\mu =0.100`$ the attenuation in the oscillatory trajectory becomes stronger. In fig. 8I a, there is a weak damping producing a long oscillatory behavior and in fig. 8I c, the trajectory is quickly attenuated in the first half of the falling height, and hereafter it follows a vertical trajectory. Similar decreasing oscillations, are observed when one small sphere suspended from a fixed point by a string or a rod, oscillates against the air and generates a damped harmonic oscillator.

The behavior of the vertical velocity is shown in fig. 8II. There is a clear damping in the peak-to-peak velocity amplitude, in all the cases. The attenuation in the vertical velocity amplitude and the time between two consecutive turning points are not very different from each other, as the fluid dynamical viscosity is changed. For the three viscosities the oblate adquires approximately the same final mean velocity $`v_y=3\frac{cm}{s}`$.

A larger variation in the velocity is observed in its horizontal component fig. 8II. There is a strong attenuation in the curve for larger viscosities fig. 8IIc. One explanation is that the interaction between the walls and the oblate is smaller when the viscosity is decreased . This is also supported by the larger angular variation $`\mathrm{\Delta }\mathrm{\Theta }`$ in fig. 8IIc, for smaller viscosities.

In fig. 8 III we see that the first peak-to-peak amplitude are nearly the same for all considered viscosities. For $`\mu =0.025`$ the subsequent peak-to-peak amplitudes attain a constant value of $`\mathrm{\Theta }_{pp}15^o`$ while for higher viscosities periodic oscillations are replaced by erratic motion of low amplitude.

### D Steady-Falling Oblate: Change in the Aspect-Ratio of the Oblate.

When the oblate’s falling motion begins, the oblate gains a larger oscillation amplitude in its oscillatory falling trajectory. And for all the three trajectories presented in fig. 9, this first amplitude decreases $`A_o=3.0,2.8,1.8cm`$ as the oblate’s aspect-ratio is increased $`\mathrm{\Delta }r=0.18,0.22,0.29`$ respectively. This first large oscillation is a common characteristic for these trajectories. As the aspect-ratio is incremented, the number of cycles and the amplitude of the trayectories change. For $`\mathrm{\Delta }r=0.22`$ the trajectory presents a well defined steady-falling behavior, fig. 9I(b), but, if the aspect-ratio increases, ($`\mathrm{\Delta }r=0,29`$), the trajectory varies, the peak-to-peak amplitude is quickly damped in the first half of the vertical trajectory, with the interesting observation that in the second half the trajectory doesn’t have a steady-falling behavior fig. 9 Ia. When $`\mathrm{\Delta }r=0.18`$, the smaller aspect-ratio, the trajectory has a oscillatory behavior fig. 9I(c), with a constant peak-to-peak amplitude of $`3cm`$.

Since the minor oblate radius is fixed in our simulations, when the aspect-ratio is increased $`\mathrm{\Delta }r=0.18,0.22,0.29`$, the oblate’s area gets smaller, and the final vertical velocity increases $`Vy=3.5,3.7,3.9\frac{cm}{seg}`$, respectively, fig. 9II. The final vertical velocity decreases with the decrement in the aspect-ratio, since the smaller aspect-ratio presents more area against the fluid. As the aspect-ratio is increased the peak-to-peak amplitude in the vertical velocity diminishes. For the smaller aspect-ratio $`\mathrm{\Delta }r=0.18`$, we have the larger amplitude $`1\frac{cm}{s}`$, and as the aspect-ratio is increased $`\mathrm{\Delta }r=0.22`$ and $`\mathrm{\Delta }r=0.29`$ the peak-to-peak amplitude tends to be much smaller.

The peak-to-peak amplitude for the vertical orientation $`\mathrm{\Theta }_{pp}`$ increases when the aspect-ratio decreases or the oblate’s area increases. For $`\mathrm{\Delta }r=0.29`$ the peak-to-peak amplitude is $`\mathrm{\Theta }_{pp}=15^o`$, (fig. 9a), and much smaller compared to $`\mathrm{\Theta }_{pp}=70^o`$, (fig. 8c), for $`\mathrm{\Delta }r=0.18`$. In all cases the oblate at the end orients vertically (fig. 9 III).

### E Periodic Behavior of a Falling Oblate.

We have found periodic behavior for smaller dynamical viscosity ($`Re=480`$) and smaller aspect-ratio ($`\mathrm{\Delta }r=0.133`$). The dynamics of the falling oblate is governed by inertial effects. In figure 10 (I), we show the transition from a quasi-periodic, or a long steady-falling trajectory($`\mu _1=0.100`$), to a periodic behavior fig. 10 (I,$`\mu _2=0.025`$), when the dynamical viscosity is varied from $`\mu _1=0.1`$ to $`\mu _2=0.025`$. The trajectory presented in fig. 10I, with dynamical viscosity $`\mu _2=0.025`$ has a wave length of 20 cm.

The vertical velocity shown in fig. 10 (II), has the same transition from a long steady-falling regime with a final average velocity of $`3.0\frac{cm}{s}`$ to the periodic regime where the velocity has a oscillation period of 3.3 s.

The vertical orientation presented in fig. 10 (III) has also the same transition from a long steady-falling regime to periodic behavior with a period of 6.6 s, and the angular values oscillate around $`\mathrm{\Theta }_0=90^o`$ with angular peak-to-peak amplitude $`\mathrm{\Theta }_{pp}=60^o`$.

We perform three simulations in the periodic regime with very different initial orientation angles and the corresponding trajectories are shown in fig. 11I. In the case of $`\mathrm{\Theta }_0=26^o`$ the peak-to-peak amplitude is 2.3 cm and for ($`\mathrm{\Theta }_0=90^o`$) it is 0.4cm, and the oscillatory behavior is observed for the three cases. The peak-to-peak amplitude of the oscillation in the trajectory fig.10I, decreases as $`\mathrm{\Theta }_o`$ is increased. In the case of the vertical velocity and orientation fig. 11II-III, the initial orientation angle also plays the same role, reducing the peak-to-peak amplitude of the curves, and for $`\mathrm{\Theta }_o=90^o`$ the amplitude of oscillation is the smallest of the three.

The final vertical velocity and orientation for the three trajectories are $`3.0\frac{cm}{s}`$ and $`85^o`$ respectively. We can say that the average final values for the vertical velocity and orientation are not modified by the variation of the initial orientation $`\mathrm{\Theta }_o`$.

We have found the largest peak-to-peak amplitude $`1.0\frac{cm}{s}`$ for the vertical velocity oscillation for an initial orientation of $`\mathrm{\Theta }_0=26^o`$ and it becomes smaller as the oblate’s initial orientation tends to $`\mathrm{\Theta }_0=90^o`$. For the peak-to-peak amplitude of the vertical orientation the oblate shows a similar behavior: we have the larger value for the amplitude $`\mathrm{\Theta }_{pp}=70^o`$, and it is obtained for an initial orientation of $`\mathrm{\Theta }_0=26^o`$. The smallest peak-to-peak amplitude $`\mathrm{\Theta }_{pp}4^o`$ is obtained for $`\mathrm{\Theta }_0=90^o`$.

In figure 12a (simulation) we show the time dependence of the vertical orientation with $`\mu _2=0.025`$. The value of the angular peak-to-peak amplitude is $`\mathrm{\Theta }_o=60^o`$. The vertical velocity fig. 12b (simulation) reaches its maximal value $`3.6\frac{cm}{seg}`$ as $`\mathrm{\Theta }`$ approaches $`\mathrm{\Theta }_{max}`$, presenting a minimal drag.

The smaller vertical velocity $`V_y=2.5\frac{cm}{s}`$ in the oscillation is reached at $`\mathrm{\Theta }_{min}105^o`$. The butterfly shape of fig. 12b was also measured in the experimental work of Belmonte et al.  (fig. 12b), exhibiting a vertical orientation $`\mathrm{\Theta }`$ that oscillates at double the period of $`V_y`$.

The horizontal velocity oscillates around zero with the same period as $`\mathrm{\Theta }`$, presenting its maximum value at $`Vx_{max}=1.5\frac{cm}{s}`$ and the minimum at $`Vx_{min}=1.5\frac{cm}{s}`$ at $`\mathrm{\Theta }90^o`$ as seen in fig. 12c. When the horizontal velocity is zero the oblate takes its maximum ($`120^o`$) and minimum ($`60^o`$) values in $`\mathrm{\Theta }`$ corresponding to a non-zero value of the vertical velocity $`V_y=3.2\frac{cm}{s}`$.

### F Sensitivity to a Change in the Initial Orientation In the chaotic Regime.

We can ask what is the sensitivity of the oblate to tiny changes in the initial orientation. Therefore we have simulated three trajectories shown in fig. 13, which have slightly different initial orientation. A tiny variation in the relative orientation ($`\mathrm{\Delta }\theta _o=10^3`$) produces, however, a significant variation in the shape of all the curves. These can be appreciated in the lower part of the trajectories.

In order to get better sensitivity, we have incremented the falling height to $`h_o=166cm`$. The resulting trajectories for four slightly different initial orientations in the vertical plane are presented in fig. 14. In this regime the system presents a high sensitivity to the initial orientation condition. For the four trajectories the relative angular variation is $`\mathrm{\Delta }\theta _o=10^3`$.

Due to this dependence on small changes in the initial orientation, we proceed to use as a tool of diagnosis, the Fourier power spectrum of the time series of the horizontal coordinate $`x(t),x(t+\delta t),x(t+2\delta t)\mathrm{}`$, and in our case $`\delta t=0.053566`$. A broad spectrum of frecuencies appears, as shown in fig. 15II, indicating chaotic motion.

The autocorrelation function, for the same time series (see fig. 15I), does not fall quickly to zero, it decreases linearly with time. The points are not independent of each other and a self similarity is present in the data.

In the figures 14III, we present a slices or Poincaré sections $`(p_x,x)`$, corresponding to the trajectories in fig. 14a, and which are quite irregular.

The orbits are quasi-periodic in the sense that they pass repeatedly and irregularly through the whole domain without ever closing on themselves, and without any particular time period associated with succesive passages.

The sensitivity to initial conditions is clear in these four figures. A small change in the initial orientation results in large changes in position and velocity.

We can now investigate quantitatively this sensitivity by studying the increment in the Euclidean distance $`d_{p_1p_2}=\sqrt{(x_1x_2)^2+(y_1y_2)^2}`$, between the curves presented in fig. 14 (a) and (c). Fig. 16, shows that the distance between nearby points has an overall exponential time dependence $`d(t)\mathrm{exp}(\lambda t)`$ and the fit gives an estimate for the Lyapunov exponent $`\lambda =0.052\pm 0.005`$. The positivity of the Lyapunov exponent is a clear indication for chaos.

### G Parameter Phase Diagram.

We explore the phase space in the dimensionless moment of inertia $`I^{}`$, which is the ratio of the moment of inertia of an oblate around its principal axis to the moment of inertia of a sphere of liquid with the same diameter and the Reynolds number $`Re`$. We do a similar analysis for our results as in the work of Field, et al. . It is important to remark that the mentioned experiment was for a falling disk, with small aspect-ratio, and we expect that the dynamics of the system will be close to that of an oblate ellipsoid.

The definitions of the dimensionless variables for our system are:

$$I^{}=\frac{I_{oblate}}{I_{sphere}}=\frac{5}{4}\frac{r_m}{r_M}\frac{\rho _{oblate}}{\rho _{fluid}}=\frac{5}{4}\frac{\rho _{oblate}}{\rho _{fluid}}\mathrm{\Delta }r$$
(13)
$$Re=\frac{U(2r_M)\rho _{fluido}}{\mu }$$
(14)
Figure 17 bottom, shows our results. At low values of $`I^{}`$ and small Reynolds number (high dynamical viscosity), the left-down corner of the diagram, the motion is overdamped and the oblate drops to the bottom container without any oscillation. If the Reynolds number increases $`(Re100)`$, fixing the moment of inertia, the trajectory is composed of successive oscillations that will decrease in amplitude until the oblate finally comes to stop at the lowest part of the container part. For all these cases this type of motion was called in our results a steady-falling, and was studied in section IIIB-B,C,D.

For small values of $`I^{}\mathrm{\Delta }r1`$, we have a flattened ellipsoid, and Reynolds number $`(Re400)`$, the trajectory, velocity and orientation are characterized by oscillations that repeat at equal intervals of time and space. This case was called in our simulations the oscillatory regime, and was studied in section IIIB-E.

As we increase $`I^{}`$, the object will become a sphere slightly flattened in the poles, and its dynamic becomes sensitive to smaller variations in the initial orientation, exihibiting a chaotic trajectory, which is explained in section IIIF.

If we compare our diagram with the experimental results obtained by S. Field et al, ref (fig. 17 top), we can see that in both pictures, the steady-falling region is located in the lower-left part, the oscillatory region is in the lower-right, and the chaotic regime is located along the top of the picture. We can say that the two diagrams are similar, but, with the difference, that the tumbling regime in the Field’s diagram is not present in our results because by the comparison of the two diagrams in fig. 17. The tumbling regime would be found for larger dimensionsless moment of inertia, in our case smaller aspect-ratio, demanding a larger three dimensional container and impliying a very expensive computational study, than in the other regimes.

If we take a different initial orientation angle, the new phase diagram exhibits the same dynamical behavior of figure 17 bottom. The coexistence of the dynamical phases, explained above, is independent of the initial orientation of the oblate.

## IV Conclusions and Outlook

The motion of a single oblate settling in a fluid in a three dimensional container has been studied. We found three basic regimes for the dynamics of the system (steady-falling, oscillatory, and chaotic). The steady-falling and the periodic motion exhibit a similar physical behavior as observed for flattened bodies , . With the exception that the tumbling motion is missing in our simulations.

We have characterized the dynamics of the steady-falling regime when the dynamical viscosity, dropping height, and oblate’s aspect-ratio are changed. Several conclusions can be drawn from this part of the work.

(a) The spatial trajectories $`(x,y)`$ are composed of oscillations that correspond to a damped harmonic motion. This regime is present for small values of $`I^{}0.51`$, $`Re100`$ and is shown in fig.6-7. There is no variation in the trajectories when we increase the initial height. When the aspect-ratio is varied, the trajectories change very much, fig. 8. The aspect-ratio dominates strongly the type of trajectory that is present in the system.

(b) The final vertical velocity $`V_y`$ does not depend on the initial falling height and the dynamical viscosity.

(c) The vertical orientation $`\mathrm{\Theta }`$ of the oblate, undergoes a rotational motion till its mayor axis is aligned with the direction of gravity. This tendency of finding the minimal resistance against the fluid, is present for all Reynolds numbers in the range $`Re30100`$ used in our simulations.

The periodic behavior in our simulations is found for ($`Re500`$), and small ($`I^{}0.5`$). The most important characteristic on this regime is that the vertical orientation $`\mathrm{\Theta }`$, oscillates at the double of the period of the vertical velocity $`V_y`$, and at the same period of the horizontal velocity $`V_x`$. This periodic motion is also present in the work of Belmonte et al. . In this regime the initial orientation determines the value of the oscillation amplitude in the spatial trajectory $`(x,y)`$, velocity $`V_y`$ and orientation $`\mathrm{\Theta }`$. For $`\mathrm{\Theta }_o=90^o`$ the amplitudes of the above quantities approach a smaller value.

The chaotic behavior is present for $`I^{}1`$ and in the entire range of Reynolds numbers used in the simulation. The separation between the spatial trajectories of the falling oblate will diverge for small variations in the initial orientation $`\mathrm{\Theta }_o`$, and grows exponentially in time. The value found for the Lyapunov exponent is around $`\lambda =0.052\pm 0.005`$.

The construction of the phase diagram shows three well-defined dynamical regions as in the case of ref. . But with the difference that the chaotic behavior in the above reference is associated with a transition to chaos through intermittency for which we have no indication in our simulations. The phase diagram is independent on the initial orientation.

More work to better understand the role of the fluid pressure and velocity fields as well as a more systematic study of the phase transition in the phase diagram seem necessary

## V Acknowledgment

This research is part of the SFB-404-Project A7. Helpful discussions with S. Schwarzer, R. Mück, and E. Kuusela are acknowledged and appreciated.


format: text

text:
# HST/STIS Spectroscopy of the Intermediate Polar EX Hydrae1footnote 11footnote 1Based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with proposal 8807.

## 1 INTRODUCTION

Intermediate polars (IPs) are a class of magnetic cataclysmic variables in which a white dwarf primary star sustains a magnetic field that is sufficient enough to divert the gravitationally controlled mass transfer. IPs may have a truncated accretion disk (the inner edge being wiped out by the magnetic field) or may be diskless accretors (accreting directly from the ballistic accretion stream), but the defining characteristic that separates IPs from the higher magnetic field polars is an asynchronously rotating white dwarf. As the magnetic field lines of the white dwarf sweep around the inner edge of the accretion disk, curtains of magnetically controlled material form and follow along the field lines until the material impacts the white dwarf surface at or near the magnetic poles, with both upper and lower poles in IPs accreting equally.

EX Hydrae is an IP with a spin period of 67 minutes and an orbital period of 98 minutes, which causes an unusual scenario in which the corotation radius is close to the $`L_1`$ point (King & Wynn, 1999). Magnetically controlled accretion occurs from the inner edge of the disk, and possibly from the outer edge of the disk, as evidenced by dips seen during the bulge dip in EUV light curves of EX Hya (Belle et al., 2002). Radiation from EX Hya is seen to be modulated on both the binary orbital and white dwarf spin periods. Binary modulation is seen in the form of a large bulge dip at phases $`\varphi _{98}=0.51.1`$ ($`\varphi _{98}`$ refers to the binary phase) in the EUV light curve, which is caused by absorption along the line of sight to the emitting poles due to the extended hot spot on the edge of the accretion disk. Continuum emission is modulated sinusoidally over spin phase for all wavelength regimes (e.g., Hellier et al., 1987; Rosen et al., 1991; Belle et al., 2002). Maximum spin phase occurs when the upper accretion pole points away from the observer, allowing a direct view of upper and lower accretion regions; minimum spin flux occurs when the upper pole points towards the observer and the lower pole is self-eclipsed by the white dwarf.

Previous UV observations of EX Hya consist of HUT spectra (Greeley et al., 1997), ORFEUS II and IUE spectra (Mauche, 1999), and HST/FOS spectra (Eisenbart et al., 2002). These observations show that the broad emission lines and the continuum flux are formed in different regions; the emission lines from an optically thick accretion curtain and the continuum from the inner most part of the accretion curtain and the white dwarf itself. Blue-shifted absorption was also seen in the C IV $`\lambda 1549`$Å emission line. The UV data of EX Hya have also been modeled using a variety of continuum emission sources. The continuum has been fit by a $`T=25000`$ K blackbody (Greeley et al., 1997), but this same work also showed that a bremsstrahlung model, with $`T=74000`$ K, or a white dwarf photosphere of $`T=26200`$ K also fit the data well. Eisenbart et al. (2002) found that a white dwarf photosphere of $`T_{\mathrm{eff}}=25000\pm 3000`$ K (determined from line ratios) fit their FOS spectra well.

## 2 OBSERVATIONS

Our HST/STIS observations, obtained through DDT time, comprised six orbits and used the E140M grating in Echelle mode. Each orbit produced a time-tagged spectrum consisting of 40 Echelle orders. The E140M grating spans a wavelength range of $`11501735`$Å and has a resolving power of $`R=45,800`$, which corresponds to a wavelength resolution of $`0.0130.019`$Å/pixel throughout the spectrum. The six orbit observation provided complete coverage of EX Hya’s spin orbit, but only $`2/3`$ of its binary orbit. Table 1 gives the specific phase coverage of each data set. Each spectrum was processed with IRAF/STSDAS software and the time-tagged data was utilized to bin the spectra into phase resolved components.

Spectra were extracted for a variety of orbital and spin phase bins. Spin phase bins were determined according to the spin modulation seen in the EUV light curve of EX Hya (Belle et al., 2002) and we have defined spin maximum phases as $`\varphi _{67}=0.8650.365`$ and spin minimum phases as $`\varphi _{67}=0.3650.865`$, where $`\varphi _{67}`$ refers to the 67 minute spin ephemeris, $`T=2437699.8914(5)+0.046546504(9)\mathrm{E}7.9(4)\times 10^{13}\mathrm{E}^2`$ (Hellier & Sproats, 1992). Spectra extracted throughout the spin period were binned into $`\mathrm{\Delta }\varphi _{67}=0.1=402.2`$ s phase bins; bins were chosen as such to preserve a moderate signal-to-noise ratio. Three types of binary phased spectra were created: binary phased, spin maximum binary phased, and spin minimum binary phased. For the binary ephemeris, we used $`T=2437699.94179+0.068233846(4)\mathrm{E}`$ (Hellier & Sproats, 1992). Spin maximum data existed for binary phases $`\varphi _{98}=0.00.7`$, but spin minimum data only existed for binary phases $`\varphi _{98}=0.00.6`$. Each binary phase bin has a resolution of $`\mathrm{\Delta }\varphi _{98}=0.1=589.5`$ s, where $`\varphi _{98}`$ refers to the 98 minute orbital period. We used IRAF/SPLOT to analyze the spectra and determine the $`1\sigma `$ errors we report here. Figure 1 presents our complete spectral data set divided into two phase bins: spin maximum and spin minimum. The spin maximum spectrum has on overall higher continuum level and stronger emission lines than the spin minimum spectrum. Each spectrum exhibits broad emission lines of C III $`\lambda 1175`$Å, Si III $`\lambda 1206`$Å, $`\lambda 1298`$Å, Ly$`\alpha `$ $`\lambda 1216`$Å, N V $`\lambda 1240`$Å, C II $`\lambda 1335`$Å, O V $`\lambda 1371`$Å, Si IV $`\lambda 1396`$Å, C IV $`\lambda 1549`$Å, and He II $`\lambda 1640`$Å, along with narrow emission lines N V $`\lambda \lambda 1238,1242`$Å and O V $`\lambda 1371`$Å, which appear in all phases except for the early binary phases ($`\varphi _{98}=0.00.1`$).

## 3 DISCUSSION AND ANALYSIS

### 3.1 Continuum

Continuum flux measurements were made at three emission line-free wavelength regions along each spectrum: $`\lambda \lambda 13521362`$Å, $`\lambda \lambda 14401520`$Å, and $`\lambda \lambda 15651620`$Å. The results of the $`\lambda \lambda 13521362`$Å measurements are shown in Figure 2 (spin phased). The other two wavelength regions exhibit similar behavior. Error bars on each measurement are RMS values as reported by IRAF/SPLOT. The spin phased continuum flux (Figure 2) is fit well by the sinusoidal function, $`f(10^{13}\mathrm{ergs}\mathrm{s}^1\mathrm{cm}^2\mathrm{\AA }^1)=A+B\mathrm{sin}2\pi (\varphi _{67}\varphi _0)`$, where $`A=1.68`$, $`B=0.38`$, and $`\varphi _0=0.76\pm 0.05`$, which places the peak of the light curve at $`\varphi _{67}=0.01\pm 0.05`$. The error here is dominated by the fact that each phase bin measurement was confined within the $`\mathrm{\Delta }\varphi _{67}=0.1`$ phase bin resolution of the data. The continuum flux is modulated by $`B/A=23\%`$. The peak in UV continuum flux is earlier than that of the EUV emission, which occurs at $`\varphi _{67}=0.115\pm 0.001`$ (Belle et al., 2002). As both the HST and EUVE observations of EX Hya occurred simultaneously, an error in the ephemeris is ruled out as a cause for phase difference in peak continuum fluxes; rather, the phase disparity is likely a result of distinct EUV and UV continuum emission sites. The sinusoidal fits to all of the spin phased UV continuum wavelength regions are given in Table 2.

The spectra were also measured over binary phase, and, in agreement with previous studies (e.g., Hellier et al., 1987), we found that the continuum flux did not exhibit coherent modulation over the binary orbital period.

### 3.2 Emission Lines

#### 3.2.1 Broad Emission Lines

The broad emission lines were fit with single Gaussian functions and a linear continuum, using IRAF/SPLOT, to determine parameters such as central wavelength, flux, FWHM, and equivalent width. The Ly$`\alpha `$ and Si III $`\lambda 1206`$Å lines were not measured because they are blended during the spin maximum phases.

Over spin phase, the emission line fluxes exhibit a sinusoidal shape (C IV is shown as an example in Figure 3), and can be fit well with a sinusoidal function of the form $`f(10^{12}\mathrm{ergs}\mathrm{s}^1\mathrm{cm}^2)=A+B\mathrm{sin}2\pi (\varphi _{67}\varphi _0)`$. The solutions for each line are given in Table 3. The emission peaks near $`\varphi _{67}0.050.1`$ for each line. Considering the error of $`0.05\varphi _{67}`$ on each phase bin, the emission line fluxes peak at roughly the same phase as the continuum flux ($`\varphi _{67}0`$). C III $`\lambda 1175`$Å has the greatest modulation, at $`73\%`$, but all of the lines are strongly modulated by around $`50\%`$. Compare this with the continuum with a flux modulation of only $`2025\%`$. While both the continuum and emission lines are modulated sinusoidally over spin phase, the difference in the modulation amplitude suggests that the continuum and emission lines have different formation sites. Over binary phase the emission line flux is not modulated coherently. Comparison of the emission line fluxes over binary and spin phase demonstrates that the flux modulations are clearly a function of spin phase, rather than binary phase.

Equivalent widths (EWs) of strong UV emission lines from optically thin plasmas are typically on order of 100s of Ångstroms. The EWs of the emission lines appearing in the EX Hya spectra are smaller than expected for optically thin plasma emission. C IV, which has the largest EW of all the broad emission lines, is shown as an example in Figure 4 (solutions for the sinusoidal fits to all of the lines are given in Table 4). There are two possible scenarios that may explain the observed EWs. The first is that optically thick material is responsible for the emission lines. The second is that the optically thin material is only weakly emitting and another form of radiation, such as cyclotron emission, may play an important role in the radiative cooling of the accretion curtains. If the latter is true, evidence should be seen in the form of cyclotron humps; for a magnetic field strength of $`B5`$ MG, the first cyclotron harmonic should appear around 10$`\mu `$. The almost triangular shapes of the lines implies that they are not from a purely optically thin material, thus it is likely that both optically thin and optically thick emission contributes to each emission line.

Figures 5 and 6 present two examples of the broad emission line evolution over binary and spin orbit. The N V $`\lambda 1240`$Å line shown in Figure 5 has been separated into spin maximum (thick line) and spin minimum (thin line) and phased over the binary orbit. The line retains roughly the same form throughout the binary orbit for both spin phases. Aside from binary phase bin $`\mathrm{\Delta }\varphi _{98}=0.00.1`$, the flux of the spin maximum line is constantly greater than the flux of the spin minimum line. This suggests that the lines are formed in a region that would be partially absorbed or occulted during spin minimum phases, i.e. the inner part of the accretion curtain. Also of interest is the lack of narrow emission lines during binary phases $`\mathrm{\Delta }\varphi _{98}=0.00.1`$, which correspond to the egress of the bulge dip; absorption features appear in their place. The narrow lines will be discussed in further detail in §3.2.2. The fact that the lines have roughly the same amount of absorption during $`\mathrm{\Delta }\varphi _{98}=0.00.1`$ implies that the lines are formed interior to the bulge.

If the change in flux of the broad emission lines between spin phases maximum and minimum is due to absorption through the accretion curtain, we can estimate the density of the curtains from the amount of absorption. Using the simple relation of $`I=I_0e^\tau `$, we may first deduce $`\tau `$ and then assume a Kramer’s opacity law to determine the density. The physical size of the accretion curtains must first be estimated. At most, the radial extent of the footprint of the curtain on the disk could be the entire truncated accretion disk and at the least, it would only be a small fraction of the disk size. The disk size was determined by measuring optical emission lines (using the optical data obtained by us during this campaign) to determine the inner and outer accretion disk radius, assuming a Keplerian accretion disk. We considered two white dwarf masses (see §3.3), $`M_{WD}=0.90M_{}`$ and $`M_{WD}=1.33M_{}`$, and calculated the disk radii corresponding to each white dwarf mass ($`R_{\mathrm{inner}}=1.62.3\times 10^9\mathrm{cm}`$ and $`R_{\mathrm{outer}}=2.55.9\times 10^9\mathrm{cm}`$). Several different physical sizes, $`S`$, of the accretion curtain were used for our density calculation. We take $`S=aR_{\mathrm{disk}}`$, where $`a=0.01,0.10,0.50`$ and $`1.0`$. In this way we were able to sample a range of physical depths for the accretion curtain. The emission line fluxes were measured during spin maximum and spin minimum, excluding binary phases of the bulge dip. This ensures that material associated with the bulge, which is also responsible for absorption at certain spin phases, will not be mistakenly identified with accretion curtain material. Five emission lines were used for this calculation; C III, N V, Si IV, C IV, and He II. Si III $`\lambda 1296`$Å and C II were not used due to their ambiguous line profiles during spin minimum. For each white dwarf mass and corresponding $`S`$ values, we find that $`\tau =0.71.1`$ and the density ranges between $`10^{10}10^9\mathrm{g}\mathrm{cm}^3`$ for all values of $`S`$. If we assume that the absorbing material is pure hydrogen, we obtain a number density of $`N_\mathrm{H}=7\times 10^{13}8\times 10^{14}\mathrm{cm}^3`$ within the accretion curtains.

The C IV $`\lambda 1549`$Å emission line has a blue-shifted absorption component that appears during spin phases $`\mathrm{\Delta }\varphi _{67}0.20.8`$ as evidenced in Figure 6. A similar feature was also observed in the C IV line in HUT data during spin phases $`\mathrm{\Delta }\varphi _{67}=0.360.48`$ (their data only extended through $`\varphi _{67}=0.5`$, Greeley et al. 1997). Although their C IV line is phased into finer orbital bins, our C IV line clearly shows an absorption component that begins prior to $`\varphi _{67}=0.3`$ and evolves into an absorption that eats away the blue wing of the emission line. At phases $`\varphi _{67}=0.20.3`$, the blue absorption has a FWHM of $`275\mathrm{km}\mathrm{s}^1`$ and is blue-shifted by $`615\mathrm{km}\mathrm{s}^1`$ from the center of the C IV doublet. This is slightly less than the value of $`1000\mathrm{km}\mathrm{s}^1`$ reported by Greeley et al. (1997) for their blue C IV absorption feature.

C III, N V, and C II also show a blue absorption component during spin phases $`\varphi _{67}0.20.7`$. The N V absorption component has a FWHM of $`400\mathrm{km}\mathrm{s}^1`$ and is blue-shifted by $`700\mathrm{km}\mathrm{s}^1`$ from the center of the line. The C III absorption component is really too shallow to accurately measure, though it is visible during spin phases $`\varphi _{67}0.20.7`$. The C II line seems to have a blue-shifted absorption component throughout the entire spin phase. Its FWHM ranges from $`655\mathrm{km}\mathrm{s}^1`$ at $`\varphi _{67}=0.40.5`$ to $`200\mathrm{km}\mathrm{s}^1`$ at $`\varphi _{67}=0.00.1`$ while the blue-shift from the center of the emission line remains fairly constant at $`1400\mathrm{km}\mathrm{s}^1`$. However, the C II absorption feature is associated with a white dwarf absorption line (see §3.4), which may influence the motion of the line.

Radial velocity curves for the broad emission lines are shown in Figure 7 phased over the spin period. The spin phased radial velocity curves present surprising results. In previous EX Hya observations (e.g., Hellier et al., 1987), the maximum blue-shift of the optical emission lines observed around $`\varphi _{67}=0.0`$ has been used to place the origin of these lines in the accretion curtain; while the upper pole points away from us, material flowing into the pole will be at its maximum blue-shift. We see that three lines exhibit this behavior: C III, N V, and C II. Si III, Si IV, and He II, on the other hand, have maximum blue-shifts closer to $`\varphi _{67}0.5`$. A double peaked structure is visible in all but the Si III radial velocity curve.

#### 3.2.2 Narrow Emission Lines

Three narrow emission lines appear in the HST spectra: N V $`\lambda \lambda 1238.8,1242.8`$Å and O V $`\lambda 1371.3`$Å, and represent the two highest ionization energy lines present in the spectrum, with $`\chi _{\mathrm{NeV}}=98`$ eV and $`\chi _{\mathrm{OV}}=114`$ eV. These lines have FWHM$`0.2`$Å and are visible over both spin and binary orbit, and only disappear during binary phases $`\varphi _{98}=0.00.1`$, which correspond to the egress of the bulge dip. A radial velocity curve over spin phase is shown for the narrow emission lines in Figure 8. Each radial velocity curve has a double-peaked structure, although the O V and N V $`\lambda 1238`$Å lines are less modulated than the N V $`\lambda 1242`$Å line. For each line the velocities tend to be blue for two reasons: one, only two-thirds of the binary orbit, $`\varphi _{98}=0.00.7`$, are (unevenly) sampled in each spin phase bin and two, a $`\gamma `$velocity has not been subtracted. Previous measurements of the $`\gamma `$velocity (derived from optical spectroscopy) range from $`9\pm 35\mathrm{km}\mathrm{s}^1`$ (Gilliland, 1982) to $`180\mathrm{km}\mathrm{s}^1`$ (Hellier et al., 1987), and we will show below that we obtain a value slightly different from the Gilliland (1982) and Hellier et al. (1987) results.

The radial velocities of the narrow emission lines over the binary phase are shown in Figure 9 and can be fit with a sinusoidal function of the form $`v=\gamma +K\mathrm{sin}2\pi (\varphi _{98}\varphi _0)`$, with $`\gamma =9.5\pm 3\mathrm{km}\mathrm{s}^1`$, $`K=59.6\pm 2.6\mathrm{km}\mathrm{s}^1`$ and $`\varphi _0=0.98\pm 0.05`$. The $`K`$ amplitude derived here is comparable with a previously determined value of $`K_1=69\pm 9\mathrm{km}\mathrm{s}^1`$ (Hellier et al., 1987). The small velocities, the clear sinusoidal modulation of the velocities, and the $`K`$ amplitude suggest that these lines are formed near or on the white dwarf surface. A high temperature emission region for these lines increases the likelihood that the lines are not contaminated by emission from the accretion disk. This narrow line radial velocity curve is therefore likely to be a true representation of the motion of the white dwarf, as compared to optical (disk) emission lines that have been used previously to determine the $`K_1`$ amplitude and the $`\gamma `$velocity (e.g., Hellier et al., 1987). We would like to note that the proximity of the narrow line emission region to the white dwarf surface should cause these lines to be gravitationally redshifted. Unfortunately, without an accurate $`\gamma `$-velocity measurement, it is impossible to separate the gravitationally redshift from the $`\gamma `$-velocity of the narrow emission line radial velocity curve.

The narrow emission lines appear to be double-peaked during certain spin phases. The N V doublet shows a strong double-peaked profile during spin phases $`\varphi _{67}=0.10.3`$, $`\varphi _{67}=0.50.6`$, and $`\varphi _{67}=0.91.0`$. O V is double-peaked during spin phases $`\varphi _{67}=0.00.3`$, $`\varphi _{67}=0.40.6`$ and $`\varphi _{67}=0.70.8`$. These phases roughly correspond to spin maximum phases. Close inspection of Figure 5 also shows the double-peaked nature of the N V line during spin maximum and that during spin minimum, it is only the blue component of the double-peaked line that is still visible. The radial velocities of these lines over binary phase tell us that they originate close to the white dwarf surface. The double-peaked nature of the lines suggest a two velocity component as their origin. We propose that these lines are formed in a wind or material outflowing from the poles. During spin maximum, we have a direct view of both poles (possibly not a direct view of the northern pole as it points away from us, but a wind would clearly be seen) giving us both red and blue-shifted components. At spin minimum, when the upper pole points towards the observer, the lower pole and its (red-shifted) emission is blocked from view; only the blue-shifted narrow line component would be seen at these phases. As discussed in §3.2.1, evidence of a wind is also seen in the form of blue-shifted absorption in several of the broad emission lines.

### 3.3 Mass Determinations

In order to calculate a white dwarf mass, a value for the secondary mass must be estimated. Previous studies have used an $`M_2`$ value determined from a zero age main sequence (ZAMS) mass-radius relation applied to the Roche lobe radius for a given binary orbital period. Hellier et al. (1987) found $`M_2=0.128M_{}`$ using the ZAMS mass-radius relation given in Patterson (1984) and an earlier study used the relation presented in Warner (1976) to find $`M_2=0.18M_{}`$. We use a value based on an up-to-date mass-radius relation (Howell et al., 2001) expressed in the form of

$$M_2=0.08f^{1.95}P_{\mathrm{orb}}^{1.3}$$
(1)
where $`P_{\mathrm{orb}}`$ is the binary orbital period in hours and $`f`$ is a factor that describes the bloating of the secondary star in systems above the period gap. For the EX Hya orbital period of 1.6 hr and $`f=1`$, we calculate $`M_2=0.152M_{}`$, a value slightly higher than the $`0.128M_{}`$ that has been used previously to determine the white dwarf mass.

Using the $`K_1`$ amplitude determined from our UV narrow emission line radial velocity curve, we expect to obtain an unambiguous determination of the mass of the white dwarf if our value of $`M_2`$ is correct. We employ a modified form of Kepler’s third law,

$$\mathrm{sin}^3i=\frac{K_1^3P_{\mathrm{orb}}}{2\pi GM_2}\left(\frac{q+1}{q}\right)^2$$
(2)
and solve for $`q`$ in order to determine $`M_{WD}`$. For $`i=78\pm 1^{}`$ (Hellier et al., 1987), $`M_2=0.152M_{}`$, and $`K_1=59.6\pm 2.6\mathrm{km}\mathrm{s}^1`$, we determine a white dwarf mass, $`M_{WD}=1.33\pm 0.11M_{}`$. This value is considerably higher than $`0.78\pm 0.17M_{}`$, which was determined by Hellier et al. (1987), and the even lower values of $`0.5M_{}`$ that have been determined from X-ray studies (e.g., Fujimoto & Ishida 1997). Although our $`K_1`$ amplitude agrees within the errors with the value measured by Hellier et al. (1987), our white dwarf mass does not agree with theirs due to the larger secondary mass used in the above calculation.

Smith, Cameron, & Tucknott (1993) applied the method of skew-mapping to detect the secondary stars in cataclysmic variables. Applying this to EX Hya, they found a velocity amplitude for the secondary motion of $`K_2=356\pm 4\mathrm{km}\mathrm{s}^1`$. This amplitude, along with our measured value of $`K_1`$ and calculated value of $`M_2`$ may be used to apply another determination of $`M_1`$ using $`K_1/K_2=M_2/M_1`$. This method gives a white dwarf mass of $`M_{WD}=0.91\pm 0.05M_{}`$, which is $`30\%`$ lower than the value determined above, although still larger than the earlier estimates (e.g., Hellier et al. 1987).

### 3.4 Spectral Modeling

Discerning exactly which components of the binary system are contributing UV continuum flux (and at which spin phases) is an important issue for thorough modeling of the HST spectra of EX Hya. Spin maximum phase would seem to be the most likely time in which we would have a direct view of the white dwarf, but are there other components also visible during this phase? The accretion poles (shock regions) are strong emitters in the X-ray and EUV (Belle et al., 2002). The irradiated material in the accretion curtain close to the pole, however, may emit reprocessed X-ray and EUV radiation as UV radiation, thus adding another UV emission component to the system. The accretion disk, although truncated, is likely to add to the UV continuum emission as well.

To model our HST spectra with a white dwarf photosphere, we must determine which spectrum is the best representation of the white dwarf by ascertaining which of the different components discussed above contribute to the spin maximum and spin minimum spectra. Visual inspection of the spin maximum and minimum spectra reveals their differences; the red-ward slopes of each spectrum are roughly the same while the blue end of the spin minimum spectrum has a lower flux than that of the spin maximum spectrum. There are two possible cases to consider; $`1`$ spin minimum is the true white dwarf spectrum and the accretion curtain is optically thin to the white dwarf continuum. Spin maximum therefore contains an additional continuum component that could be the accretion region located at either pole. $`2`$ The spin minimum spectrum is absorbed by the accretion curtains and the spin maximum spectrum is the true white dwarf spectrum plus an additional continuum component. These two options outline the occultation (1) versus absorption (2) models. For option 1, spin phased modulation is due to occultation of an emitting pole, while option 2 dictates that spin phased modulation is due to UV flux being absorbed by the accretion curtains.

For either model option, the spin maximum spectrum has added UV continuum components. This is evidenced by the fact that the blue slopes of the spin maximum and minimum spectra are not equal. If the spin modulation were simply due to absorption, the spin minimum spectrum would retain the same slope as that of the spin maximum spectrum. (This assumes a gray absorber, which is valid over the small, $`600`$Å, wavelength range covered by the HST spectra.) We show the difference spectrum of EX Hya, which is the spin maximum minus the spin minimum spectrum, in Figure 10. This spectrum represents the additional UV continuum and emission line components that are seen during spin maximum phases. Considering the above arguments, we rule out the spin maximum spectrum as representative of only the UV white dwarf continuum flux. For option 1 to be viable, the accretion curtains near the poles would have to be dense enough to absorb the flux emitted from the accretion region. While the outer curtains are optically thin, it is likely that the curtain close to the pole is optically thick and will absorb emitted radiation from the accretion pole. We therefore suggest the case that the spin minimum spectrum is the better representation of the pure white dwarf spectrum and that a blue blackbody and emission line component (or multiple components) is seen in addition to the white dwarf during spin maximum. Flux contributed by the truncated accretion disk must also be considered at all spin phases.

We first created white dwarf photosphere models using TLUSTY and SYNSPEC (Hubeny, 1988) for the two white dwarf masses computed in §3.3, $`0.90M_{}`$ ($`\mathrm{log}g=8.5`$) and $`1.33M_{}`$ ($`\mathrm{log}g=9.4`$). We modeled white dwarf photospheres of solar composition and considered temperatures ranging from $`T_{\mathrm{eff}}=19000\mathrm{K}35000`$ K (in steps of $`1000`$ K). Although it was decided that the spin minimum spectrum would be the best representation of the white dwarf spectrum, the model photospheres were fit to the spin maximum and minimum spectra in order to insure that the spin minimum spectrum was indeed the best choice. For each fit, the broad emission lines in the HST spectrum were masked, except for the wavelength range including the Si III$`\lambda 1298`$Å and C II$`\lambda 1335`$Å emission lines because of the white dwarf absorption components that appear between these emission features. Our best fitting $`0.90M_{}`$, $`T_{\mathrm{eff}}=23000`$ K white dwarf photosphere model is shown in Figure 11 fit to both the spin maximum and spin minimum spectra. The spin minimum fit has a reduced $`\chi ^2`$ of 1.7, and scaled to the HST continuum fluxes, gives a distance to EX Hya of 54 pc. The spin maximum fit has a reduced $`\chi ^2`$ value of 3.3 and scaled to the HST continuum fluxes gives a distance of 46 pc.

Note that many of the absorption features appearing between the Si III and C II in the spin minimum spectrum are fit well by the white dwarf photosphere model. The slope of the model spectrum also produces a better match to that of the spin minimum spectrum than the maximum spectrum. These facts, coupled with the lower reduced $`\chi ^2`$ value of the spin minimum fit, corroborate with our earlier statement that the spin minimum phases are the ideal phases during which to model the white dwarf photosphere. We therefore use only the spin minimum spectrum in our following model fits.

White dwarf photospheres with temperatures of $`T_{\mathrm{eff}}=2300035000`$ K (in steps of $`2000`$ K) were created for $`M_{WD}=1.33M_{}`$. The best fit model to the spin minimum spectrum has a temperature of $`T_{\mathrm{eff}}=25000`$ K, but while this model provides a good fit, it also places the distance of EX Hya at $`30`$ pc! This distance varies greatly from the distances of $`65100`$ pc that have been determined in previous studies (e.g., Hellier et al., 1987; Eisenbart et al., 2002). We explored the relationship of distance versus white dwarf photosphere temperature for both white dwarf masses over a range of temperatures and found that, for $`M_{WD}=1.33M_{}`$, a white dwarf temperature of $`T_{\mathrm{eff}}=35000`$ K would be required to achieve a distance of 55 pc (when scaled to the spin minimum spectrum). An $`M_{WD}=1.33M_{}`$, $`T_{\mathrm{eff}}=35000`$ K white dwarf photosphere model was considered for the spin minimum spectrum, but the slope did not match the HST continuum slope. However, before we completely rule out either white dwarf mass, we must also consider the accretion disk contribution to the UV continuum flux.

Accretion disk models were created using TLUSDISK, SYNSPEC, and DISKSYN (Hubeny, 1988; Wade & Hubeny, 1998) for white dwarf masses of $`M_{WD}=0.90M_{}`$ and $`M_{WD}=1.33M_{}`$. We used an accretion rate of $`10^{16}\mathrm{g}\mathrm{s}^1`$ (Fujimoto & Ishida, 1997) and an inclination of $`i=78^{}`$ (Hellier et al., 1987). In order to determine the radius at which the disk should be truncated, the inner disk radius was determined from our optical spectra. For a Keplerian disk, the highest velocities in an emission line originate from the inner disk. The optical line wings, therefore, should map to the inner edge of the disk. The H$`\alpha `$ $`\lambda 6563`$Å line profile in each of our optical spectra was used to determine the velocities at the blue and red wings. For a $`M_{WD}=0.90M_{}`$ white dwarf, these velocities map to an inner accretion disk radius of $`R_{\mathrm{inner}}=2.5R_{\mathrm{WD}}`$ ($`1.6\times 10^9\mathrm{cm}`$) and for a $`M_{WD}=1.33M_{}`$ white dwarf, the inner radius becomes $`R_{\mathrm{inner}}=9.0R_{\mathrm{WD}}`$ ($`2.3\times 10^9\mathrm{cm}`$). Previous modeling studies (e.g., Greeley et al., 1997) have assumed $`R_{\mathrm{inner}}5R_{\mathrm{WD}}`$ ($`3.5\times 10^9\mathrm{cm}`$) and optical spectroscopy by Hellier et al. (1987) gave $`R_{\mathrm{inner}}=6\times 10^9\mathrm{cm}`$ for $`M_{WD}0.8M_{}`$.

Two different accretion disk models must be created for the two different white dwarf masses we have been using. $`M_{WD}=0.90M_{}`$ and $`\dot{M}=10^{16}\mathrm{g}\mathrm{s}^1`$ correspond to Models $`n`$ and $`p`$ of Wade & Hubeny (1998) and $`M_{WD}=1.33M_{}`$ with the same accretion rate matches most closely with Models $`y`$ and $`z`$. Model $`n`$ ($`p`$) has UV flux contributing annuli that extend to $`6.2R_{\mathrm{WD}}`$ ($`9.40R_{\mathrm{WD}}`$), and Model $`y`$ ($`z`$) has contributing annuli that extend out to $`14.3R_{\mathrm{WD}}`$ ($`21.7R_{\mathrm{WD}}`$); the last annulus in each model has $`T10000`$ K. Although it would seem that the synthetic spectra for each disk model would be different because of the different truncation points, the truncated accretion disk spectra are in fact quite similar. This is due to the fact that the more massive white dwarf has a hotter accretion disk and while the disk is truncated at a greater radius than the less massive white dwarf, an approximately equal amount of UV flux is still emitted from the remaining annuli.

The truncated disk models were added to the white dwarf photosphere models using the following method: the photosphere model and disk model were scaled separately to a distance, $`d`$, and then summed together. The summed spectrum was then fit to the observed spectrum. The distance, $`d`$, was modified until the summed spectrum produced a good fit to the observed spectrum. The summed model fits were judged by visual inspection, rather than by a formal fitting routine. For $`M_{WD}=0.90M_{}`$, we found that a $`T_{\mathrm{eff}}=23000`$ K white dwarf plus the truncated disk produced the best fit at a distance $`d=60`$ pc. The $`M_{WD}=1.33M_{}`$ model showed a best fit for $`T_{\mathrm{eff}}=25000`$ K plus the truncated disk at a distance of $`d=33`$ pc. The $`M_{WD}=0.90M_{}`$ fit is shown in Figure 12, as we definitively rule out the higher mass model below.

The white dwarf photosphere plus truncated accretion disk model fits to the HST fluxes give two different distances: the $`0.90M_{}`$ white dwarf model gives $`d=60`$ pc, while the $`1.33M_{}`$ white dwarf gives $`d=33`$ pc. We can check these distances independently by using Bailey’s method (Bailey, 1982), which uses the $`K`$band magnitude and measured surface brightness values for K-M single stars to determine a distance;

$$\mathrm{log}d=\frac{K}{5}+1\frac{S_K}{5}+\mathrm{log}\left(\frac{R_2}{R_{}}\right)$$
(3)
where $`S_K`$ is the surface brightness. Using $`K=11.705(28)`$ from the 2MASS survey (Hoard et al., 2002), $`S_K=4.5\pm 0.5`$ for an M3 main sequence star (Dhillon et al., 1997), $`R_2=0.17R_{}`$, and the secondary $`K`$band flux contribution as $`75\pm 25\%`$ (Dhillon et al., 1997), we determine the distance to EX Hya to be $`d=54_7^{+12}`$ pc. This effectively rules out the $`1.33M_{}`$, $`T_{\mathrm{eff}}=25000`$ K white dwarf, as the distance determined by this model is 33 pc.

The greatest uncertainties in the distance calculation and the white dwarf mass determination are the secondary star mass and spectral type. An M3 main sequence secondary star (determined from IR spectra, Dhillon et al., 1997) has a mass of $`0.3M_{}`$. This is not consistent with the mass determined using the ZAMS mass-period relation (Equation 1), which gave the secondary mass as $`0.15M_{}`$. Returning to Equation 2, we find that a secondary mass of $`0.12M_{}`$ would be required to obtain $`M_{WD}=0.9M_{}`$. A more accurate secondary mass will be needed to aid future distance and white dwarf mass determinations of EX Hya.

## 4 CONCLUSION

We have made an accurate measurement of the $`K`$ amplitude of the white dwarf based on the narrow emission line radial velocity curve. The small velocities and high ionization temperatures required for production of the N V and O V emission lines indicate that they are formed close to the white dwarf, and therefore, the $`K_1=59.6\pm 2.6\mathrm{km}\mathrm{s}^1`$ value that we have measured is an accurate representation of the motion of the white dwarf within the binary system.

Synthetic spectral models for a white dwarf photosphere and truncated accretion disk were made using TLUSTY and SYNSPEC (Hubeny, 1988). Based on spectral fits using a white dwarf photosphere model, we determined that the spin minimum phases are the best representation of the white dwarf spectrum, i.e., occultation of (one of) the emitting poles is responsible for the continuum modulation seen over the white dwarf spin phase. The distance we derive based on spectral modeling, $`d60`$ pc, agrees well with two recent measurements; Eisenbart et al. (2002) determine $`d=65\pm 11`$ pc using Bailey’s method and a $`K`$ magnitude derived from their infrared spectra, while a recent parallax measurement has placed the distance to EX Hya at $`60`$ pc (K. Beuermann, private communication). This implies that our spectral modeling parameters for a solar composition white dwarf photosphere and truncated accretion disk - $`M_{WD}=0.9M_{}`$, $`T=23000`$ K, and $`R_{\mathrm{inner}}=2.5R_{\mathrm{WD}}`$ \- represent EX Hya well.

We thank the director of the Hubble Space Telescope for allocating DDT time for the EX Hya multi-wavelength project. We also thank K. Beuermann for use of his EX Hya parallax distance prior to publication. The anonymous referee is thanked for very constructive comments. Support for proposal 8807 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. SBH acknowledges partial support of this work from HST Grant GO-08807.01-A.


text:
# Errata to Local contribution to the Lefschetz fixed point formula

## Remarks

Let $`\overline{R}(x_0,y_0)=[0,x_0]\times [0,y_0]`$ denote the closure of the rectangle defined in §4.1. The new condition (d) implies the following slightly weaker condition:

* For some (and hence for any sufficiently small) $`x_0,y_0>0`$ the set

$$S=Wc_1^1t^1(\overline{R}(x_0,y_0))c_2^1t^1(\overline{R}(x_0,y_0))C$$

is closed (and hence compact),

which is equivalent to the statement that the closure of $`S`$ is contained in $`W.`$ The (erroneous) proof of Lemma 5.10, which appears in §5.11, actually shows that (d) implies (d). In fact, this weaker condition (d) suffices for the case $`j=1`$ of Theorem 4.7, although the stronger condition (d) is needed for the other cases.


text:
# Surface versus bulk Mott transition in CaxLa1-xVO3

## Abstract

The Mott insulator LaVO3 is known to become metallic if La is replaced by small concentrations of Ca. Since surface quasi-particle spectra of various perovskites are more strongly correlated than their bulk spectra, the coexistence of a metallic bulk and an insulating surface layer might be feasible. To investigate this possibility the dynamical mean field theory is used to evaluate the quasi-particle spectra in the bulk and at the surface of CaxLa1-xVO3 for various Ca concentrations.

Significant surface contributions to quasi-particle spectra of strongly correlated transition metal oxides have recently been identified in photoemission measurements on several perovskite materials: SrRuO3 fujioka , Sr2RuO4 SrRuO1 ; SrRuO2 , CaxSr1-xVO3 aiura ; maiti ; sekiyama , and CaxLa1-xVO3 maiti00 . The latter series is particularly interesting since LaVO3 is insulating and becomes metallic upon doping with Ca. By varying the photon energy Maiti et al. maiti00  were able to show that surface correlation features are more pronounced than in the bulk, in accord with the observations on the other perovskites cited above. Thus, the coherent peak at the Fermi level is less intense and the satellite peak below the V $`t_{2g}`$ band is stronger than observed in the corresponding bulk photoemission spectra. Even more remarkably, Maiti et al. identified a Ca concentration range where the bulk has become metallic while the surface still appears insulating, before, at even higher Ca doping, it also becomes metallic. Although sample preparation and experimental resolution make the data analysis difficult, these results address a very interesting issue, namely the possible coexistence of different bulk and surface phases in semi-infinite systems. In the case of magnetism, for example, it is well known that at finite temperatures the magnetic moment can be reduced in the surface layers of a ferromagnet, but that both have the same Curie temperature unless the surface is artificially decoupled from the bulk magnetism .

The identification surface-induced enhanced correlations is clearly important in order not to interprete them erroneously as bulk correlations. For instance, earlier photoemission data on CaxSr1-xVO3 aiura  appeared to be in conflict with various low-frequency measurements onada . While the latter data showed that this perovskite series is metallic for all Ca concentrations, the photoemission spectra seemed to suggest that SrVO3 is strongly correlated and that CaVO3 is on the verge to a metal-insulator transition. Recent photoemission measurements maiti ; sekiyama  covering a wide range of photon energies demonstrated instead that the bulk spectra of both materials are consistent with metallic behavior, but that correlation features are much more pronounced at the surface. While SrVO3 is a cubic perovskite, the oxygen octahedra in CaVO3 are orthorhombically distorted. This leads to a slight reduction in $`t_{2g}`$ band width and a broadening of the van Hove singularity. The density of states at $`E_F`$, however, is only weakly affected. That the bulk properties of SrVO3 and CaVO3 should indeed be very similar was recently also found by Nekrasov et al. nekrasov  in bulk electronic structure calculations combining the the local density approximation (LDA) with the many-body dynamical mean field theory (DMFT) georges ; vollhardt ; pruschke .

As shown by Liebsch lieprl  using surface DMFT calculations, the enhancement of surface correlations can at least in part be understood in terms of the reduced atomic coordination at the surface and the concomittant narrowing of the local density of states. Thus, the ratio of local Coulomb energy and band width $`U/W`$ is effectively enhanced at the surface. Accordingly, the coherent peak in surface quasi-particle spectra is less intense and the lower Hubbard band is enhanced compared to the bulk. Other mechanisms, such as additional band narrowing due to reconstruction and distortion of oxygen octahedra, or less efficient screening of electron-electron interactions at surfaces, might further increase the importance of correlations in surface photoemission spectra. The recent work on CaxSr1-xVO3 therefore established that photoemission is consistent with low-frequency bulk probes as long as surface effects are properly taken into account.

The same conclusion applies to Sr2RuO4 for which previous photoemission spectra seemed to contradict bulk de Haas-van Alphen measurements SrRuO1 ; SrRuO3 . Recent experimental and theoretical work SrRuO2 ; SrRuO4  proved, however, that this discrepancy can be resolved by taking into account the lattice reconstruction at the surface of Sr2RuO4 which leads to significant changes in the photoemission spectra.

In the present work we investigate surface correlations in the perovskite series CaxSr1-xVO3 at Ca doping levels where metallic bulk and insulating surface behavior might coexist. Since doping with Ca leads only to minor structural changes and nearly constant Hubbard $`U`$ bocquet , the modifications of the electronic properties are almost entirely caused by the degree of band filling. According to theoretical predictions by Potthoff and Nolting potthoff , a phase with metallic bulk and insulating surface properties cannot occur unless electronic coupling between surface and bulk is assumed to be weak. These results were derived for the semi-infinite half-filled simple cubic $`s`$ band, and by adopting a linearlized version of the DMFT in which low- and high-frequency excitations are essentially decoupled. It is therefore not clear to what extent they are applicable to realistic multi-band systems at arbitrary band filling. Here we calculate surface quasi-particle spectra by using a realistic layer-dependent local density of states appropriate for CaxSr1-xVO3, combined with multi-band DMFT Quantum Monte Carlo (QMC) calculations.

Fig. 1 shows the bulk and surface density of states for LaVO3 and Ca0.5La0.5VO3. The main difference is the band filling. While the $`d^2`$ occupancy of the V $`t_{2g}`$ bands makes LaVO3 insulating solovyev , Ca0.5La0.5VO3 with its $`d^{1.5}`$ occupancy is metallic. In the bulk, the insulator-metal transition occurs at about $`x0.2`$ maiti . Since our main interest here is to study the effect of band filling on this transition in the bulk and at the surface we adopt for simplicity the ideal cubic symmetry and neglect the distortion of oxygen octahedra. According to the LDA-DMFT calculations for CaVO3 nekrasov  these structural changes have little influence on the $`t_{2g}`$ quasi-particle spectra. The bulk density of states is therefore similar to that of SrVO3. It exhibits the characteristic asymmetric shape resulting from the pronounced planar nature of the $`t_{2g}`$ states which is found for various peroskite materials takegahara .

In a cubic structure, the $`t_{2g}`$ bands are degenerate. The surface lifts this degeneracy since only the intra-planar $`d_{xy}`$ states exhibit strong dispersion within the surface layer while $`d_{xz,yz}`$ bands are modified due to the reduced atomic coordination (the $`z`$-axis specifies the surface normal). The layer dependent local density of states is calculated using a tight-binding formalism for semi-infinite systems kalkstein . Details concerning the application of this method to $`t_{2g}`$ bands are given Ref. lieprb . As shown in Fig. 1, the spectral weight is shifted from the low- and high-energy regions to intermediate energies close to $`E_F`$. Thus, although the full band width coincides with the one in the bulk, the local density of $`d_{xz,yz}`$ states in the surface layer is effectively narrowed. The surface distributions shown in Fig. 1 (a) and (b) are qualitatively similar. The small difference arises from slightly different surface potentials which are employed in the electronic structure calculation for semi-infinite CaxLa1-xVO3 in order to ensure charge neutrality. In the deeper layers, the local density of $`d_{xz,yz}`$ states approaches the bulk $`t_{2g}`$ density rather quickly, the main modification consisting in an oscillatory deviation which decreases with increasing distance from the surface. In contrast to the $`d_{xz,yz}`$ densities, the local density of $`d_{xy}`$ states in the surface layer is almost identical to the isotropic bulk density. This is a consequence of the pronounced planar character of the $`t_{2g}`$ states.

The narrowing of the local density states in the surface layer has an important influence on the quasi-particle distributions which we evaluate by using the multiband QMC-DMFT. Essentially we are dealing with coupled impurity problems where layer-dependent self-energies must be calculated and iterated until self-consistency is achieved. Since the full solution of this problem for multiband systems is computationally not yet feasible, we neglect the coupling between layer baths and treat bulk and surface layers independently. In previous single-band DMFT calculations for surfaces this coupling was found to be small potthoff . Enforcing charge neutrality at the surface presumably also diminishes the importance of inter-layer coupling between the impurity baths. The key input in the DMFT calculation is therefore the orbital- and layer-dependent local density of states. As a result of the one-electron band narrowing we can expect a similar enhancement of surface correlation effects as discussed previously for SrVO3 and SrRuO3 lieprb .

As mentioned above, LaVO3 is insulating, but hole doping by replacement of La via Ca induces an insulator-metal transition at rather low Ca concentrations. This can be seen clearly in the different shapes of the imaginary-time Green’s function for bulk CaxLa1-xVO3, as shown in Fig. 2 for several values of $`x`$. Note that the value of $`G(\tau )`$ at $`\tau \beta /2`$ is representative of the weight of the quasi-particle peak at $`E_F`$, while $`G(\beta )`$ gives the occupancy per spin-band. In the case of LaVO3 ($`x=0`$) already for $`U=4.5\mathrm{}5.0`$ eV there is little weight near $`\tau \beta /2`$, indicating the reduction of intensity near $`E_F`$ as $`U`$ approaches the critical value for the Mott transition. Of course, since the DMFT-QMC calculations are carried out at finite temperatures, $`G(\beta /2)`$ cannot approach zero, i.e., a finite spectral density at $`E_F`$ remains even in the nominally insulating phase. According to the results shown at finite Ca concentrations, however, $`G(\beta /2)`$ diminishes much less rapidly, indicating the incipient metallicity for $`x>0`$. Thus, in the range of Coulomb energies where LaVO3 is insulating, i.e., for $`U5.0`$ eV, CaxLa1-xVO3 for $`x>0.2`$ is metallic. We conclude from these calculations that within our model the borderline doping level between insulating and metallic behavior is about $`x0.1`$.

According to these results we consider $`x=0.1`$ as the most promising Ca concentration at which the bulk is already close to metallicity while the surface layer might still be insulating. Fig. 3 shows bulk quasi-particle spectra for Ca0.1La0.9VO3 for various Coulomb energies in the critical region of the metal-insulator transition. The analogous surface spectra are shown in Fig. 4. All spectra exhibit the three-peak structure characteristic of strongly correlated systems: a coherent peak at $`E_F`$ and upper and lower Hubbard bands which appear as satellites in photoemission or inverse photoemission data. As in the systems studied previously lieprb , as a consequence of the effective narrowing of the $`d_{xz,yz}`$ subbands, correlations are more pronounced at the surface than in the bulk: the Hubbard peaks are more intense and the coherent peak is a narrower. Since the QMC-DMFT calculations are performed at about 1450 K ($`\beta =8`$) the coherent peak persists even in the insulating phase. Nevertheless, in view of the temperature dependence of this feature known from one-band systems georges ; bulla  the spectra in shown in Figs. 3 and 4 suggest that at the surface the likelihood for a gap to open at $`E_F`$ at low temperatures is slightly greater than in the bulk. In the bulk there is clearly more spectral weight near $`E_F`$ compared to pure LaVO3 ($`x=0`$lieprb , reflecting the beginning metallic bulk behavior at $`x=0.1`$ Ca doping concentration. The surface spectra shown in Fig. 4, on the other hand, indicate stronger tendency towards gap formation, in particular for $`U=5`$ eV.

In spite of the more pronounced surface correlations discussed above, it is not clear whether extrapolation to low temperatures would yield different Mott transitions for the bulk and the surface layer. Even though the spectral weight at $`E_F`$ is consistently smaller at the surface, both $`N_s(E_F)`$ and $`N_b(E_F)`$ could still vanish at the same critical $`U`$. This would be analogous to the behavior found in semi-infinite Heisenberg models: the magnetization near the surface is smaller than in the bulk, but the critical temperature for the transition towards paramagnetism is the same magnetism . Only if surface and bulk are strongly decoupled their critical behavior might be different. A careful examination of this question would also require to account for the coupling between the layer dependent ‘impurity’ baths which we have neglected for computational reasons. Presumably this coupling whould diminish the differences between bulk and surface spectra.

So far we have only considered the effect of band narrowing on the $`d_{xz,yz}`$ states which hybridize preferentially in planes normal to the surface. Surface reconstruction or distortion of oxygen octahedra presumably would reduce hopping also between $`d_{xy}`$ states in the first layer, thereby leading to additional enhancement of surface correlations. Moreover, less efficient electronic screening processes near the surface could increase the local Coulomb energy. Both effects would make the surface electronically different from the bulk. It is doubtful, however, whether this decoupling would be strong enough to give rise to separate bulk/surface Mott transitions.

In summary, we have calculated bulk and surface quasi-particle distributions for CaxLa1-xVO3 using the dynamical mean field theory. At integer occupancy ($`d^2`$ for $`x=0`$) this perovskite material is a Mott insulator. At small Ca concentrations, however, the bulk material becomes metallic. As a result of the reduced atomic coordination at the surface, two of the $`t_{2g}`$ subbands exhibit appreciable band narrowing, giving rise to stronger correlation features in the quasi-particle spectra. In the range of critical local Coulomb interactions for the Ca concentration $`x=0.1`$ we have found at the surface a tendency for the formation of an excitation gap at a lower value of $`U`$ than in the bulk. Whether this implies a coexistence of metallic bulk and insulating surface behavior is, however, difficult to estimate at present because of the finite temperature used in the QMC-DMFT calculations. Also, additional surface effects on the electronic structure such as reconstruction, the charge state of the surface layer, and weaker screening of the local Coulomb energy should be considered.

Acknowledgement: I like to thank A. Bringer, E. Eisenriegler, and D.D. Sarma for discussions. I also thank A.I. Lichtenstein for the QMC-DMFT code.