text:
# Quantum Films Adsorbed on Graphite: Third and Fourth Helium Layers

## I Introduction

Helium films adsorbed on graphite exhibit a number of phases and have proven to be a rich source for both experimental and theoretical studies of two-dimensional phenomena. The graphite substrate is ordered on atomic length scales and offers a potential well for helium that is relatively strong for physical adsorption but is short ranged perpendicular to the substrate. As a result, a number of distinct, atomically thin layers occur, each with its own phase diagram. Near the graphite surface, the layers tend to solidify, with both commensurate and incommensurate solids occuring in the first two layers. In both layers, this solidification occurs before promotion to the next layer. The second layer exhibits a coverage region with superfluidity as well,  while the first layer apparently favors the formation of solid clusters over liquid droplets at low densities,  a lthough there is debate on this issue. A general review of physically adsorbed films such as helium on graphite can be found in the book of Bruch, Cole, and Zaremba. In this paper, we will focus on the liquid third and fourth layers.

Detailed information on the film structure at low temperatures for the third and higher layers has come from a number of experiments, including heat capacity, torsional oscillator, and third sound measurements. Ranges in which the heat capacity depends linearly on the coverage suggest that gas-liquid coexistence regions exist in the third and fourth layers. The transition from liquid droplets to a uniform liquid phase in each layer is signaled by a peak in the isothermal compressibility. Unlike the first and second layers, these higher layers do not solidify before layer promotion, since torsional oscillator measurements detect superfluidity for all coverages beginning at intermediate third layer densities. This apparently rules out earlier suggestions that the third layer may solidify, although it still may be possible for solidificatino to take place under compression of higher layers.

Perhaps the most unusual feature of the higher layers is the step-like behavior of superfluidity with increasing density. The superfluid coverage in layered films will not grow continuously because of the layering transitions. Plateaus in the superfluid density immediately after layer promotion are expected, and have been observed, when the particles are in the droplet region. These plateaus occur because the droplets in the new layer lack the connectivity to exhibit superflow across the entire surface. The interesting observation is that the plateaus actually begin before layer promotion for the third through sixth layers. At 500 mK, superfluidity even exhibits a decrease with increasing coverage near the promotion to the fourth layer. These effects have been discussed in the context of the Bose-Hubbard model. The suggestion is that the plateaus are produced by the increased localization of particles in the dense liquid.

Theoretical tools applied in the study of quantum films on realistically treated substrates include the hypernetted chain Euler-Lagrange (HNC-EL) method, density functional theory, and quantum Monte Carlo.  The HNC-EL approach has been used extensively in studies of the third and higher layers of helium films on graphite. This approach determines stable coverages of the helium layers. Because the theory requires a uniform film, calculations cannot be made for all coverages. Breakdowns in the theory are interpreted as occuring at coverages where the film is unstable to the formation of droplet patches or to layer promotion. The theory predicts at least three liquid layers will form on top of the solid first and second layers, and yields a maximum coverage value before promotion of 0.065 atom/$`\AA ^2`$ for each layer, in good agreement with, but somewhat below, the experimental value of 0.076.

In this paper we present results for the third and fourth helium layers using the path-integral Monte Carlo (PIMC) method. Our simulation is able to take into account the effect of the second layer’s corrugations and zero-point motion on the third layer. We are also able to allow for the possibility of promotion and demotion of particles between the second, third, and fourth layers, and these effects are observed. Finally, our simulation method can be applied to the entire range of possible phases in a layer, from liquid droplets to full solidification. Thus we are able to probe both the low density and high density phases of a layer.

———————————————————————–

## II Simulation method

Our simulation is performed using a path-integral Monte Carlo method that includes particle permutations and substrate effects. The general method for bulk simulation has been reviewed elsewhere,  and our modifications for simulating layered systems are given in a previous publication. Some of the present calculations have required minor changes to this method, so we will briefly outline the procedure now in order to explain the modifications.

### A Overview of procedure

The partition function $`Z`$ for a system of $`N`$ bosons at the inverse temperature $`\beta `$ can be expanded as a path-integral by inserting $`M`$ intermediate configurations:

$`Z=`$ $`{\displaystyle \frac{1}{N!}}{\displaystyle \underset{P}{}}{\displaystyle \mathrm{}d^3R_1\mathrm{}d^3R_Md^3R}`$ (2)
$`\times \rho (𝐑_1,𝐑_2;\tau )\rho (𝐑_2,𝐑_3;\tau )\mathrm{}\rho (𝐑_M,P𝐑_1;\tau ),`$
where $`\rho `$ is the density matrix, $`𝐑_i`$ is a configuration of $`N`$ particles, and $`\tau =\beta /M`$. The sum over $`P`$ is the sum of all possible permutations of particle labels. Both the configurations and the permutations are sampled by our PIMC method. The number of intermediate configurations used, $`M`$, is referred to as the number of inverse-temperature slices.

To implement this method, two important ingredients are needed. First, a starting approximation for the density matrix at $`\tau `$ is required. The simplest of these is the so-called semiclassical approximation. The drawback to this approach is that very high starting temperatures are needed to obtain an accurate approximation. Methods for improving the starting approximation can lower the starting temperature to $`\tau ^1=40K`$.

The second important feature of the method is multilevel sampling. The path-integral, Eq. (2), can be thought of as a system of $`N`$ ring polymers, each with $`M`$ beads. Particle permutations correspond to splicing the ring polymers together. In multilevel sampling, we update particle positions over a section of the polymer chain. This allows us to implement permutations by spreading them over several inverse-temperature slices. The number of slices to be updated is $`2^l`$, where $`l`$ is the level of the move. The value of $`l`$ must be chosen to balance the acceptance of particle moves and particle permutations. Increasing the value of $`l`$ increases the rate of accepting permutations but decreases the acceptance of new particle positions, while decreasing $`l`$ has the opposite effect. Typically, we take $`l=3`$ since this gives the optimum balance between the acceptance of permutations and new particle positions.

### B Third and higher layers

In our calculations for the third and higher layers, we use a simulation cell that is designed to accommodate the first layer helium solid. Periodic boundary conditions are applied in the plane of the substrate. We approximate the effects of the solid first layer on the second and higher layer atoms by placing frozen atoms at triangular lattice sites on the substrate. These frozen atoms are located 2.8$`\AA `$ above the graphite surface, the height of the first layer as indicated by neutron scattering. The second layer also solidifies before promotion to the third layer occurs. In a previous simulation, we determined that the highest second layer coverage before layer promotion was 0.2117 atom/$`\AA ^2`$. This simulation was performed with 20 second layer atoms above 30 frozen first layer atoms in a cell of dimensions $`15.075\AA \times 15.67\AA `$. The majority of our simulation results for the third and higher layers were obtained using a simulation cell of this size.

In order to investigate the effects of an inert versus active second layer, we have performed two sets of calculations. In the first, we equilibrate the second layer at 0.2117 atom/$`\AA ^2`$. The positions of these atoms are then frozen and are no longer included in the sampling. Additional atoms are then placed above this inert substrate and have their positions and permutations sampled. We refer to these additional atoms as being “active”. The bisection level used in this simulation was $`l=3`$. All errors that we report for these calculations are statistical errors arising from the Monte Carlo simulation of the active atoms above a particular frozen second layer. There will be additional systematic errors that arise from our particular choice for the frozen second layer configuration.

In the second set of calculations, both the second and third layer particles are included in the sampling. There is a potential problem with doing Monte Carlo calculations on such a system. The third layer is liquid and permutations will occur at low temperatures. This favors using $`l=3`$ for this layer. The second layer, on the other hand, is solid and furthermore increases in coverage as the overall coverage is increased. Using $`l=3`$ produces a low acceptance rate for particle moves in the compressed second layer. Sampling efficiency can be improved by using $`l=2`$ for this layer. Thus a single value of $`l`$ for the entire system is not optimum.

To give each layer the best value of $`l`$, we have partitioned atoms into “second layer” and “third layer” atoms. “Second layer” atoms have their positions initially taken from an equilibrated second layer solid. “Third layer” atoms are started from an initial configuration of atoms placed at various heights above the second layer. Promotion and demotion between the layers are allowed, but the layer label of the atoms does not change. “Third layer” atoms are sampled with $`l_{3rd}=3`$, while “second layer” atoms are sampled at $`l_{2nd}=2`$. Permutations are allowed between atoms with the same layer label, but we do not allow atoms with different layer labels to permute. This limitation is not a problem for atoms demoted to the second layer after being initially placed on the third layer, since exchanges are uncommon in the solid second layer. The promotion of a particle from the second layer to the third occasionally occurs after a particularly long simulation. At most, only one particle is promoted, so this will have little effect on the permutations in the third layer in all but the lowest coverages. This approach gives a higher acceptance rate for second layer particle moves, while allowing permutations to occur reasonably often in the third layer. We have found that this partitioning of $`l`$ lowers the energy by an amount ranging from 0.0 K to 0.1 K per atom relative to energy calculations with a single value $`l=3`$ for the entire system. The smallest energy shift occurred at the lowest coverage tested, 0.2286 atom/$`\AA ^2`$. At this coverage the two values were within error bars. The greatest energy shift occurred near the equilibrium coverage of the third layer liquid. This was the highest coverage we tested. We also compared the energy values at selected coverages to calculations performed with a single value $`l=2`$ for the entire system and found them to agree, but third layer superfluidity was suppressed, as expected.

## III RESULTS FOR THE THIRD AND FOURTH LAYERS

Before presenting our simulation results, we wish to clarify the convention that we use to report adsorbed helium coverages. Normally, the values we give for the density are for the total adsorbed helium (first two solid layers plus any additional coverages for the higher layers). Relative coverages within a layer are prefaced with a reference to that layer. For example, the coverage 0.2583 atom/$`\AA ^2`$ using our standard simulation cell corresponds to a third layer coverage of 0.0466 atom/$`\AA ^2`$ plus two solid layers with a combined coverage of 0.2117 atom/$`\AA ^2`$.

### A Results with the inert second layer

#### 1 Layer promotion and demotion

In order to study multiple layers of the helium film, it is first necessary to establish that our simulation does, in fact, produce distinct layers with increasing coverage. This is illustrated in Fig. 1, which shows the growth of the density profile perpendicular to the substrate for the third and higher layers with increasing density. The peaks associated with the third and fourth layers can be clearly observed, as can the beginning of the fifth layer peak at the highest simulated coverage. Also as can be seen in the figure, the continued growth of the third layer peak for all coverages indicates compression of this layer even after atoms are promoted to the higher layers.

Promotion to the fourth layer occurs for coverages greater than 0.2837 atom/$`\AA ^2`$ (17 active atoms). This coverage may be determined from the following observation, which is illustrated in Fig. 2. For coverages at and just below promotion, increasing the coverage increases the peak height but does not appreciably change the profile’s width. Just above this density, the peak height does not change, but an abrupt increase in the width is observed. With increasing coverage, the tail of the profile just above layer promotion evolves into the fourth layer peak. This is also illustrated in the figure. This value for layer promotion is in agreement with experiment. Heat capacity measurements show layer promotion at 0.288 atom/$`\AA ^2`$, while isothermal compressibility measurements  give a somewhat lower value 0.280 atom/$`\AA ^2`$ for the promotion density.

Layer promotion is also signaled by a change in the density dependence of the chemical potential. By differencing our calculated total energy values, we can obtain the chemical potential $`\mu `$. This is plotted in Fig. 3. The values for the energy per particle are given in Table I. As can be seen from the figure, the chemical potential increases rapidly just below layer promotion until it becomes favorable to promote an atom to the next layer. Above layer promotion, the chemical potential remains roughly constant with increasing density, as it should for the liquid-gas coexistence region. This change in the density dependence of $`\mu `$, together with the density profiles, determines the promotion density.

Layer promotion in general is related to the chemical potential. As the density of a particular layer increases above the equilibrium coverage, so does $`\mu `$. Each atom added to the system must choose whether to go in the dense layer or sit above it in the unfilled layer. This choice is governed by the chemical potentials for the two layers. Above the equilibrium coverage of a layer but below promotion, the $`\mu `$ of the dense layer is lower than $`\mu `$ for the unoccupied layer. However, the chemical potential for the dense layer increases rapidly with increasing density and at some coverage, the chemical potential for adding an atom to the dense layer exceeds the chemical potential for adding it to the unoccupied layer. The above discussion is illustrated for the third and fourth layers by Fig. 3. As can be seen in the range 0.26 to 0.28 atom/$`\AA ^2`$ each additional atom causes a rapid increase in $`\mu `$. Around 0.284 atom/$`\AA ^2`$, $`\mu _3`$ of the third layer exceeds $`\mu _4`$ of the fourth layer, and so the next atom added to the system will be preferentially promoted to the fourth layer.

The beginning of particle promotion to a new layer does not signal the end of the filling of the old layer. As can be seen in Fig. 3, adding atoms to the system at the densities 0.306 and 0.334 atom/$`\AA ^2`$, marked by the letter “D” and the dark arrowheads, causes the chemical potential to increase above the chemical potential at layer promotion. Immediately above both of these densities, $`\mu `$ drops back to values equal to $`\mu `$ at promotion. This signals particle demotion to the third layer. Figure 4 illustrates this for densities around 0.306. Notably, each additional demoted atom increases the energy required to further increase the layer density, so $`\mu `$ at promotion is only a rough estimate of the chemical potential required to demote additional atoms.

Particle demotion may be understood as a balance between many factors. Atoms are initially promoted above a layer because this is energetically favorable. The first promoted atom loses the large energy benefit for being close to the substrate, but gains kinetic energy since it is free to move about on the open surface, with its wave function no longer constrained by the other atoms. It also retains some of the potential energy advantage gained from having neighboring helium atoms. Adding more atoms to the system will favor the formation of droplets. However, once droplets have formed in the new layer, the next atom added to the system faces a different choice than the first promoted atom. This is illustrated in Fig. 5. If it goes into the less dense outer layer, it gains some attraction from the other atoms in this layer. However, it no longer gains the kinetic energy advantage that the first promoted atom had. On the other hand, if the added atom goes into the dense lower layer, it regains the benefit of being closer to the strongly attractive substrate. Furthermore, it has the energy advantage for having more helium neighbors: those in its layer, those in the layer above, and those in the layer below (not shown). In contrast, if the added atom goes into the outer layer, it has less neighbors. It should be noted that each demotion to the dense layer does significantly increase the chemical potential for this layer, so at some point it will be more favorable for atoms to be added to the outer layer again.

#### 2 Third and fourth layer phases

We expect four principal regions in the third layer before fourth layer promotion. These are a low density gas phase (which will have a negligible density at low temperatures), a droplet region, a metastable liquid region, and an equilibrium liquid phase. The droplet region consists of a liquid phase separated from a low-density gas by an interface. In the metastable region, the droplet phase is replaced by a stretched uniform phase which has negative spreading pressure. The crossover from the droplet region to the stretched liquid phase occurs at the spinodal point. Direct evidence that the layer is liquid comes from torsional oscillator measurements, which detect superfluidity up to layer completion. The isothermal compressibility has been measured for the third and higher layers and exhibits a divergence that is associated with the spinodal point. The equilibrium liquid coverage can be inferred from heat capacity measurements. Below, we present evidence of each of these phases using several different observables.

First, we can establish the existence of droplets and a uniform liquid phase at different densities with the radial distribution function, g(r), which provides a direct probe of short and long range behavior. Calculations for the third layer are shown in Fig. 6. These are plotted as functions of the magnitude of the distance vector between pairs of atoms, projected onto the plane of the substrate. These calculations of the averaged $`g(r)`$ smooth out possible anisotropies induced by the corrugations of the underlying solid helium layer. The g(r) for the three coverages shown in Fig. 6 are representative of the droplet region, the equilibrium liquid, and the liquid near layer promotion. At the lowest coverage, 0.233 atom/$`\AA ^2`$ (5 active atoms), the radial distribution function drops below unity at large distances, as would be expected for a droplet phase. The actual dimensions of the droplet for a given density depend on the size of the simulation cell. For the intermediate coverage, 0.2624 atom/$`\AA ^2`$ (12 active atoms), the first peak has changed only slightly, but the long range behavior is noticeably different, rising again past unity instead of dropping continuously. At the highest coverage, 0.284 atom/$`\AA ^2`$, the system shows evidence of increased correlation, but the long range behavior cannot be determined due to the small size of the simulation cell. This largest coverage also has a different short range behavior, showing an increased probability that the projected distance between two atoms will be less than 2.0 $`\AA `$. In part, this is a result of the thickening of the layer, as can be seen in Fig. 1. Because the atoms can be at different heights above the substrate, the projected distance between them can become smaller than would be possible in strictly two-dimensional calculations.

We can also gain some insight into the layer phases by examining contour plots of the probability distribution of atoms in the plane of the substrate. Plots near equilibrium and layer promotion are shown in Fig. 7. The high density liquid shows noticeably more localization than the equilibrium fluid. An increased correlation can also be seen in the radial distribution function at high density. This suggests that the system may be nearing a liquid-solid coexistence phase.

Having established that the third layer has gas-liquid and uniform liquid phases, we can next determine the equilibrium coverage ranges of these phases. This can be done by using the Maxwell construction. At low temperatures, the total energy and the total free energy are nearly equal, so coexistence regions may be identified by applying the Maxwell construction to the energy. The results for the low temperature (400 mK) temperature scans are shown in Fig. 8. We have verified that the energy values shown are effectively zero temperature results by recalculating some values at 500 mK. In all cases the calculations at the two temperatures agreed within error bars. The values shown in the figure have been shifted by the amount $`N_{act}e_0`$, where $`e_0=15.897\pm 0.024`$ is the minimum energy per particle, and $`N_{act}`$ is the number of active atoms in the simulation. At low temperatures, the gas phase will have zero coverage and thus zero total energy. We can thus draw a coexistence line between the beginning of the third layer, 0.2117 atom/$`\AA ^2`$, and the coverage with minimum energy per particle. This higher coverage is the equilibrium liquid density. We find this coverage to be 0.2625 atom/$`\AA ^2`$ ($`N_{act}=12`$). The best chi-squared parabolic fit around this minimum gives 0.2645(9) atom/$`\AA ^2`$. The number in parenthesis is the error in the last digit. At this density the layer is completely covered by a uniform liquid. Below this value, the system enters the gas-liquid coexistence region. The energy values in the coexistence region lie above the coexistence line, either because the liquid phase is unphysically uniform, or else because of the appreciable cost for creating a phase boundary in a finite-sized system. The third layer equilibrium coverage, 0.0528(9) atom/$`\AA ^2`$, is comparable to (and slightly higher than) the equilibrium coverage found for the second layer, 0.0480(6) atom/$`\AA ^2`$. For both layers, simulated with the same size cell, the energy minimum occurs when the system contains 12 atoms.

The equilibrium density that we determine is in good agreement with both heat capacity and torsional oscillator measurements. In the measurements of Greywall,  the low temperature heat capacity depends linearly on density from the beginning of the third layer to 0.260 atom/$`\AA ^2`$. This linear dependence is a signal of phase coexistence. The torsional oscillator measurements provide evidence of a similar region. The temperature of the oscillator’s dissipation peak, which gives a rough estimate for the superfluid transition temperature, is independent of coverage from 0.22 to 0.26 atom/$`\AA ^2`$. This is characteristic of a surface covered by liquid droplets. Increasing the coverage in the droplet region increases the size of the droplets, not their density, and so the transition temperature remains constant.

For densities above the equilibrium coverage and below layer promotion, we have not been able to rule out the onset of solid-liquid coexistence in the third layer. Our procedure for identifying such regions in the second layer was to look for the instabilities in the total energy that signal phase coexistence. This was possible because that layer solidified before layer promotion. In the present case, though, the energy cannot be used to find solid-liquid coexistence because any instabilities that might signal solidification are inextricably entangled with the liquid-vapor phase coexistence that occurs in the fourth layer at the same coverages. Specifically, layer promotion begins at the third layer coverage of 0.0720 atom/$`\AA ^2`$. The next increment in coverage that we can simulate using the cell described in Sec. II is 0.0762 atom/$`\AA ^2`$. In our previous simulation of the second layer using the same sized cell, we determined that an incommensurate solid begins to form at the second layer coverage 0.0762 atom/$`\AA ^2`$. Apparently, promotion preempts solidification. Also of note is that the density profile of the third layer at 0.0720 is much less peaked than the second layer at the same layer coverage, so any third layer solid phase must have very large zero-point motion. One possible consequence of the third layer entering solid-liquid coexistence is that it may fully solidify under compression of further adsorbed layers. As we have discussed in the previous section, the third layer density continues to increase even after fourth layer promotion. Layer promotion is not a phase transition and can occur whenever the chemical potential of the system favors it. Thus it is possible to have layer promotion in the middle of a phase transition. The third layer solid-liquid phase transition can be completed after fourth layer promotion as additional atoms added to the system get demoted to the third layer. We note finally that there is some experimental evidence suggesting that third layer solidification may occur. In the heat capacity results presented by Greywall,  a small peak at about 1.8 K can be observed for coverages beginning at 0.3100 atom/$`\AA ^2`$, between the rounded heat capacity feature associated with the fourth layer liquid and the sharp peak associated with the melting of the second layer solid.

The fourth layer also exhibits a self-bound liquid coverage. We can identify this in the same manner as before. Since promotion to the fourth layer occurs above 0.2837 atom/$`\AA ^2`$, we can consider this coverage to correspond to a zero density fourth layer gas. We apply the Maxwell construction again to the fourth layer and determine a bound liquid phase at 0.3345 atom/$`\AA ^2`$ (29 active atoms), giving a fourth layer equilibrium density of 0.0508 atom/$`\AA ^2`$. The solid line in the figure gives the maximum possible range for gas-liquid coexistence in the fourth layer. All energy values in this region are on or above the coexistence line. This coexistence region agrees with heat capacity measurements, which exhibit linear isotherms in the fourth layer up to 0.3300 atom/$`\AA ^2`$ at low temperatures.

A trend can be noticed in the equilibrium density as one progresses from the inner to outer layers. For the first layer of helium adsorbed on a flat substrate, we have calculated the equilibrium density to be 0.0450(6) atom/$`\AA ^2`$, close to the two-dimensional value. Successive increases are observed in the second and third layers, as the zero-point motion of the layers perpendicular to the plane of the substrate becomes larger. From a two-dimensional point of view, this motion has the effect of softening the hard cores of the helium atoms.

In Fig. 8, we do not plot energy values for the coverages between 0.2413 and 0.2583 atom/$`\AA ^2`$ at 400 mK. Calculations that we performed for these densities at 500 mK actually showed an energy decrease, outside of error bars. We attribute this to variations inherent in using frozen configurations for the second layer: the two calculations required two different configurations, and so there will be some systematic difference in the energy values.

We have repeated the third layer calculations at 400 mK using a larger simulation cell ($`21.105\AA \times 20.889\AA `$). This will allow us to examine finite-size effects and the effects of different configurations. We will also be able to examine the intermediate region that we excluded in the calculations using the smaller cell. The total energy values are given in Fig. 9. These energy values have been shifted by subtracting out the gas-liquid coexistence line, as described above. The energy per particle values used to construct this figure are given in Table II. The minimum energy per particle is $`16.005\pm 0.027`$ K, which occurs at 0.2676 atom/$`\AA ^2`$ (24 atoms). We fit these energy values to a polynomial in the form

$`E/V=e_0+B({\displaystyle \frac{\rho \rho _0}{\rho _0}})^2+C({\displaystyle \frac{\rho \rho _0}{\rho _0}})^3,`$ (3)
where $`\rho `$ is the coverage. Fitted parameters are given in Table III. Not all digits are significant. Notable from the table is that the equilibrium liquid coverage occurs at the same density as in previous calculation with the smaller cell, despite the fact that there is a small discrepancy in the energy values at similar coverages between the cells. In our previous calculations on the second layer, we observed that finite-size effects on the energy were negligible. Thus we believe that the difference in the minimum energy values obtained with the two cells is attributable to the particular frozen second layer configurations that we used. We also observed that a calculation using the larger cell but at 500 mK and 0.2540 atom/$`\AA ^2`$ (18 active atoms) showed an energy increase over the value obtained at 400 mK.

By differencing the (unshifted) total energy values of the larger cell, we can obtain the chemical potential $`\mu `$. Values obtained from both the fit and the actual energy values are shown in Fig. 10. The size of the chemical potential near the minimum ($`16.88K\pm 0.28`$) is in general agreement with the value range of -9 to -16 K obtained for the third layer by Clements, et al. As noted by those authors, their values are sensitive to the attractiveness of their model substrate.

We can determine the spinodal point by taking the second derivative of Eq. (3). The isothermal compressibility is given by $`\kappa _T=\rho ^2(\mu /\rho )_T^1`$, where $`\rho `$ is the coverage. This diverges when the derivative of the chemical potential is zero. Below this coverage, $`\mu /\rho `$ is negative and the speed of sound becomes imaginary. This is the spinodal point. From the fit, we determine that this occurs at 0.2518 atom/$`\AA ^2`$, a third layer coverage of 0.0386 atom/$`\AA ^2`$. This is in agreement with the HNC-EL calculations for the third layer. Experimentally, the spinodal point can be determined from the divergence of the isothermal compressibility. In the third layer, this occurs near the third layer coverage of 0.03 atom/$`\AA ^2`$ at 1.2 K.

The droplet region and the uniform coverage region for this simulation cell can be directly identified by the short wavelength behavior of the static structure function, S(k). Figure 11 shows S(k) for the coverages 0.2472 (below the spinodal point) and 0.2676 atom/$`\AA ^2`$ (near equilibrium). The upward swing of $`S(k)`$ for the lower density for small values of $`k`$ indicates the presence of a droplet. Near the equilibrium density, the droplet has entirely covered the substrate and so $`S(k0)0`$.

At the temperature of our simulation, superfluidity can be observed in the third and fourth layers. Superfluid density values, ranging from intermediate third layer to low fourth layer coverages, are shown in Fig. 12. The values shown are the ratios of the superfluid density to the third layer coverage. A slight suppression in superfluidity occurs for coverages just below and just above the promotion density to the fourth layer, although the noise of our calculations obscures this effect. Suppression of superfluidity has been observed at similar temperatures and coverage values in both the torsional oscillator measurements of Crowell and Reppy and a simulation of helium on hydrogen. Suppression just after fourth layer promotion is a consequence of gas-liquid coexistence. In the low coverage region, the fourth layer consists of droplets that lack the connectivity to exhibit superfluidity until they have percolated across the surface. The suppression of superfluidity before promotion to the fourth layer can be observed more clearly when we allow the second layer to respond to the third layer film, as described in the next section.

### B Results with active second layer

#### 1 Layer promotion and demotion

In the previous calculations we assumed that the solid second layer could be treated as inert; that is, the response of the second layer to the third was ignored. It is known experimentally, however, that the second layer solid continues to be compressed after third layer promotion. Evidence of this can be found in the heat capacity peak associated with second layer melting. This peak continues to increase in temperature for the entire range of the third layer, indicating an increase in the second layer density. Neutron scattering also detects compression of the second layer by the third, with an abrupt restructuring at intermediate third layer coverages. Second layer compression is not unexpected, since the compression of the first layer by the growth of the second is well established. To allow for the effects of second layer compression, as well as the response of this layer to the growth of the third, we have performed calculations that include both second and third layer atoms in the Monte Carlo sampling.

Density profiles illustrate the compression of the second layer in our simulation. These are shown in Fig. 13 for selected coverages. The left and right peaks in the figure are for the second and third layers, respectively. As can be seen, the height of the second layer peak grows between the coverage 0.2296 and 0.2879 atom/$`\AA ^2`$. The peak height then remains constant up to and including the coverage 0.2964, where promotion to the fourth layer is visible. Finally, the second layer peak increases again at the highest coverage examined, 0.3049 atom/$`\AA ^2`$.

By integrating the profiles up to the minimum between the two peaks (at approximately 7.5 $`\AA `$), we obtain the second layer coverages of 0.0845, 0.0885, 0.0894, 0.0890, and 0.0929 atom/$`\AA ^2`$ for the coverages shown in the figure. Thus at the lowest third layer densities, no demotion occurs. Beginning at intermediate coverages and up to layer promotion, a single atom is demoted. At low fourth layer coverages (the highest coverage examined) two atoms are demoted. As a consequence of demotion, promotion to the fourth layer is pushed to a higher overall coverage. As is illustrated in Fig. 13, we do not observe promotion to the fourth layer until the density has exceeded 0.2879 atom/$`\AA ^2`$, compared to 0.2837 atom/$`\AA ^2`$ when the second layer was frozen. This value will be further established below by examining the density dependence of the chemical potential. This value is still in agreement with the heat capacity and isothermal compressibility measurements. We note that Zimanyi, et al have proposed a slightly higher completion coverage of 0.293 atom/$`\AA ^2`$.

The layer promotion illustrated above is accompanied by a discontinuity in the chemical potential. We can obtain this quantity by differencing our calculated total energy values, given in Table IV, and plot the results in Fig. 14. Included in this figure are results obtained for the second layer from a previous calculation. As can be seen, there is a distinctive change in the density dependence of the chemical potential $`\mu `$ around the promotion density. Below the promotion density, the energy changes very rapidly with increasing coverage. Near promotion, an added atom has the choice of going to the unoccupied third layer or the dense second layer, and will choose the layer that is energetically favorable. When the promotion coverage is reached, the chemical potential for adding the atom to the second layer exceeds that for adding it to the third layer, and so the atom is added to the unoccupied layer. It can be seen from the figure that the change in $`\mu `$ associated with the promotion is quite large. In this case, the second layer is solid and relatively dense, so we expect an energy gap to be associated with promotion. As can be seen in the density profiles (see Fig. 13) for this layer, there is a range of heights above the second layer that is forbidden for the third layer atom. This is in contrast to what happens in the case of fourth layer promotion (see Fig. 2), in which there is significantly more overlap between the third and fourth layers. Each additional atom added to the system will also go the outer layer, but the rate of the energy change will decrease, since these atoms are attracted to each other and form a droplet. Thus we see that the chemical potential just after layer promotion decreases.

#### 2 Layer phases

We can determine phase boundaries for the third layer by again using the Maxwell construction. Figure 15 gives the total energy values for the three layer system, with both second and third layers active. The energy values per particle are given in Table IV. As before, a coexistence line can be drawn between the beginning of the third layer, 0.2117 atom/$`\AA ^2`$ and the equilibrium liquid coverage. The upper endpoint of the coexistence region is 0.2667 atom/$`\AA ^2`$. This is the equilibrium coverage for the third layer liquid. The intermediate energy values are in gas-liquid coexistence. Following our usual procedure, we have subtracted out the gas-liquid coexistence line from all the energy values, using the values $`30.111\pm 0.019`$ K and $`23.910\pm 0.013`$ K for the beginning (0.2117 atom/$`\AA ^2`$, 20 atoms) and end (0.2667 atom/$`\AA ^2`$, 33 atoms) of the coexistence region.

As with the results found using the frozen second layer, we can fit the energy values to a polynomial in the form of Eq. (3). This is the solid curve in Fig. 15. Energy values for the very low and high coverages of the layer were not included in the fit. The equilibrium coverage determined is $`0.2653\pm 0.0005`$ atom/$`\AA ^2`$, and the spinodal point is found to be near the third layer coverage of 0.040 atom/$`\AA ^2`$.

Suppression of superfluidity before promotion to the fourth layer is even more prominent in these calculations. Figure 16 shows the superfluid coverage (instead of the superfluid fraction) versus the coverage for the third layer. For superfluid films adsorbed on heterogeneous surfaces (such as Vycor), the superfluid coverage increases continuously with increasing coverage. This does not occur in the atomically thin fluid layer we are simulating. Superfluidity instead increases with increasing coverage until the equilibrium density is reached, after which there is a plateau and then a drop. The drop in $`\rho _s`$ is seen experimentally as well, for coverages greater than 0.280.

The reason for this suppression before layer promotion is not clear, but three related factors can contribute. First, increasing the density in the liquid layer produces excitations that should disrupt the superfluid state. This energy increase can be seen in Fig. 8 and Fig. 15. Second, increased correlation between the atoms will suppress the permutations required for superfluidity. We observe this increase in the radial distribution function, Fig. 6, and in the probability contours. The calculated probability distributions also suggests the third mechanism, namely that the layer is partially solidifying before layer promotion. As discussed above, we have not found clear evidence that this occurs.

## IV SUMMARY

This paper has presented path-integral Monte Carlo calculations for the third and fourth layers of helium on graphite using two different treatments of the underlying solid layers. In the first approach, we treated the first and second layers as inert. The second layer solid was initially equilibrated, and then atoms were frozen in some particular configuration. This simplification neglects the response of the second layer to growth of the third layer, as well effects such as particle demotion, but allows us to simulate larger systems. Using this approach, we have produced as many as three additional layers on top of the solid first and second layers. The third and fourth layers are found to possess self-bound liquid phases with layer densities of 0.0528 and 0.0508 atoms/$`\AA ^2`$, respectively. These densities are in agreement with values that may be inferred from heat capacity and torsional oscillator measurements. In the third layer, below the equilibrium density, we determine the spinodal point to be 0.0386 atoms/$`\AA ^2`$. This coverage separates the unstable region, in which droplets form on the surface, from the metastable “stretched” phase just before equilibrium. Static structure calculations confirm that droplets occur at low densities. The chemical potential was also calculated, and third layer promotion was found to coincide with a discontinuity in the density dependence of this quantity. The superfluid density has been calculated for coverages from the third layer equilibrium to low fourth layer coverages. Promotion to the fourth layer is observed at the (total) coverage 0.2837 atoms/$`\AA ^2`$. The continued growth of the third layer is also observed after fourth layer promotion. The question of partial solidification in this layer remains open for further investigation. In particular, it is possible that particle promotion preempts the formation of a third layer solid, but particle demotions at higher coverages may increase the density of this layer to the point of solidification. Experimentally, there is some suggestive evidence for third layer solidification through compression by higher layers. We find that the heat capacity feature that emerges above 0.3100 atom/$`\AA ^2`$ might be related to the melting of this solid. Solidification also provides a natural explanation for the plateaus seen in torsional oscillator measurements of the superfluid density of the higher layers.

The effects of including the second layer atoms in the Monte Carlo updating procedure have also been examined. These calculations incorporate the effects of zero-point motion in the second layer solid and allow this layer to respond to the third and higher layers. We observe particle demotion from the third to the second layer near equilibrium and again near the fourth layer promotion density. The overall coverage at which the third layer liquid phase forms increases slightly from 0.2645 to 0.2653 atoms/$`\AA ^2`$ when the second layer is active. The completion density is also slightly higher, occuring at 0.2879 atoms/$`\AA ^2`$. A very notable effect produced by the active second layer was the reduction in the third layer superfluid density at high coverages. This effect has also been seen experimentally.


format: text

text:
# Fermi Surface, Surface States, and Surface Reconstruction in Sr2RuO4

\[

## Abstract

The electronic structure of Sr2RuO4 is investigated by high angular resolution ARPES at several incident photon energies. We address the controversial issues of the Fermi surface (FS) topology and of the van Hove singularity at the M point, showing that a surface state and the replica of the primary FS due to $`\sqrt{2}\times \sqrt{2}`$ surface reconstruction are responsible for previous conflicting interpretations. The FS thus determined by ARPES is consistent with the de Haas-van Alphen results, and it provides additional information on the detailed shape of the $`\alpha `$, $`\beta `$, and $`\gamma `$ sheets.

\]

Angle resolved photoemission spectroscopy (ARPES) has proven itself to be an extremely powerful tool in studying the electronic structure of correlated electron systems. In particular, in the case of the high-temperature superconductors, it has been very successful in measuring the superconducting gap, determining the symmetry of the order parameter, and characterizing the pseudo-gap regime. On the other hand one of the fundamental issues, namely the determination of the Fermi surface (FS) topology, has been controversial, as in the case of Bi2Sr2CaCu2O8, raising doubts concerning the reliability of the ARPES results. A similar controversy has also plagued the fermiology of Sr2RuO4. In this context, the latter system is particularly interesting because it can also be investigated with de Haas-van Alphen (dHvA) experiments, contrary to the cuprates, thus providing a direct test for the reliability of ARPES.

Whereas dHvA analysis, in agreement with LDA band-structure calculations, indicates two electronlike FS ($`\beta `$, and $`\gamma `$) centered at the $`\mathrm{\Gamma }`$ point, and a hole pocket ($`\alpha `$) at the X point, early ARPES measurements suggested a different picture: one electronlike FS ($`\beta `$) at the $`\mathrm{\Gamma }`$ point and two hole pockets ($`\gamma `$, and $`\alpha `$) at the X point. The difference comes from the detection by ARPES of an intense, weakly dispersive feature at the M point just below $`E_F`$, that was interpreted as an extended van Hove singularity (evHs) pushed down below $`E_F`$ by electron-electron correlations. With the evHs below $`E_F`$, rather than above (LDA band-structure calculations place it 60 meV above $`E_F`$), the $`\gamma `$ pocket is converted from electronlike to holelike. The existence of the evHs was questioned in a later photoemission paper, where it was suggested that dHvA and ARPES results could be reconciled by assuming that the feature detected by ARPES at the M point was due to a surface state (SS). Recently, two possible explanations were proposed: first, the evHs at the M point could be only slightly above $`E_F`$ (e.g., 10 meV), so that considerable spectral weight would be detected just below $`E_F`$. Alternatively, ARPES could be probing ferromagnetic (FM) correlations reflected by the existence of two different $`\gamma `$-FS (hole and electronlike, respectively, for majority and minority spin direction), which escaped detection in dHvA experiments. Lastly, the surface reconstruction as detected by LEED, which has been proposed to be indicative of a FM surface, would also complicate the ARPES data. The resolution of this controversy is important not only for the physics of Sr2RuO4 per se, but also as a reliability test for FS’s determined by ARPES, especially on those correlated systems where photoemission is the only available probe.

In this Letter, we present a detailed investigation of the electronic structure of Sr2RuO4. By varying the incident photon energy and the temperature at which the samples were cleaved, we confirm the SS nature of the near-$`E_F`$ peak detected at the M point, and we identify an additional dispersive feature related to the ‘missing’ electronlike FS ($`\gamma `$). A complete understanding of the data can be achieved only by recognizing the presence of shadow bands (SB), due to the $`\sqrt{2}\times \sqrt{2}`$ surface reconstruction which takes place on cleaved Sr2RuO4 (as confirmed by LEED). Despite the surface complications, the FS as determined by ARPES is consistent with the dHvA results, and provides additional information on the detailed shape of the $`\alpha `$, $`\beta `$, and $`\gamma `$ sheets.

ARPES data was taken at SSRL on the normal incidence monochromator beam line equipped with a SES-200 electron analyzer in angle resolved mode. With this configuration it is possible to simultaneously measure multiple energy distribution curves (EDC’s) in an angular window of $`12^{}`$, obtaining energy-momentum information not at one single $`k`$-point but along an extended cut in $`k`$-space. The angular resolution was 0.5x0.3 \[0.3 along the cut, corresponding to a $`k`$ resolution of 1.5% of the Brillouin zone (BZ), with 28 eV photons\]. The energy resolution was 14 meV for high-symmetry cuts and photon energy dependence, and 21 meV for the FS mappings. Sr2RuO4 single crystals were oriented by Laue diffraction, and then cleaved in situ with a base pressure better than 5x10-11 torr. In order to compensate for the angular response of the analyzer, EDC’s in a single cut were normalized against those measured on polycrystalline gold. Different cuts were then normalized with respect to each other on the basis of the spectral weight above $`E_F`$ integrated in both momentum and energy.

In order to begin the discussion of our experimental results, let us first give an overview of the very rich photoemission spectra, introducing all the features we will focus on throughout the paper. Fig. 1 presents EDC’s along the high-symmetry directions for Sr2RuO4 cleaved and measured at 10 K. The $`\alpha `$, $`\beta `$, and $`\gamma `$ sheets of FS expected on the basis of LDA calculations, and dHvA experiments are indicated together with the experimentally measured cuts in the sketch depicting 1/4 of the 2D projected zone (PZ). All the features in the data are labeled following the assignments which will be given in the paper. Along $`\mathrm{\Gamma }`$-M, two peaks emerge from the background, disperse towards $`E_F`$, and cross it before the M point, defining $`\beta `$ and $`\gamma `$ electronlike pockets (Fig. 1a). Along M-X, a peak approaches and crosses $`E_F`$ before the X point defining, in this case, the $`\alpha `$ hole pocket centered at X (Fig. 1b). Similar results were obtained along $`\mathrm{\Gamma }`$-X (Fig. 1c): the $`\beta `$ pocket is clearly resolved, while $`\alpha `$ and $`\gamma `$ crossings are almost coincident. In addition, we identify a weak feature which shows a dispersion opposite to the primary peaks along $`\mathrm{\Gamma }`$-M and $`\mathrm{\Gamma }`$-X (SB, see below). Around the M point a sharp peak is observed (SS), whose weak dispersion along M-X can be followed until it crosses $`E_F`$ and loses intensity. The highest binding energy along the dispersion (which is symmetric with respect to $`\mathrm{\Gamma }`$-M) is found at the M point. Because of this behavior, the SS peak was initially associated to a sheet of FS centered at X, and holelike in character.

In the following discussion, we will concentrate on the features observed near the M point, which are relevant to the controversy concerning the character of the $`\gamma `$ sheet of FS. We will show that by working at 28 eV photon energy with sufficient momentum resolution both the $`\beta `$ and $`\gamma `$ electronlike pockets (predicted by LDA calculations) are clearly resolved in the ARPES spectra. In order to address this issue, we measured the M point region (with cuts along $`\mathrm{\Gamma }`$-M-$`\mathrm{\Gamma }`$) varying the incident photon energy between 16 and 29 eV, in steps of 1 eV. Here, we covered the location of $`\beta `$ and $`\gamma `$ pockets in both first and second zone for better illustration (i.e., four $`E_F`$ crossings will be observed). From the EDC’s shown in Fig. 2 for 16, 22, and 28 eV, we can see that the cross sections of SS, $`\beta `$, and in particular $`\gamma `$ exhibit a strong (and different) dependence on photon energy (note that by working at lower photon energies the momentum resolution increases and, in turns, the number of EDC’s becomes progressively larger in going from panel c to a). At 28 eV, $`\beta `$ and $`\gamma `$ crossings can be individually identified in the EDC’s. Owing to the high momentum resolution we can now follow the dispersion of the $`\gamma `$ peaks until the leading edge midpoints are located beyond $`E_F`$. After that, the peaks lose weight and disappear, defining the $`k_F`$ vectors for the electronlike $`\gamma `$ pockets. Right at $`k_F`$ we can resolve a double structure which then reduces to the non dispersive feature (SS) located 11 meV below $`E_F`$. The difference between the 28 eV results and those obtained at 16 or 22 eV is striking (the latter are consistent with those previously reported at similar photon energies). At low photon energies, the $`\beta `$ crossings are still clearly visible. On the other hand, we can follow only the initial dispersion of the $`\gamma `$ peaks (now broad and weak), before they merge with the SS feature, giving the impression of an evHs. At 16 eV it is impossible to identify the $`\gamma `$ crossings. At 22 eV the location of the leading edge midpoints is at best suggestive of the presence of the $`\gamma `$ crossings.

We have showed that the $`\gamma `$ electron pocket had so far escaped detection in ARPES because it is indistinguishable from the SS feature at low photon energy and/or low angular resolution. In order to have a full picture of the relevant issues to be addressed, let us proceed to the discussion of the FS mapping. Fig. 3b shows the $`E_F`$ intensity map obtained at 28 eV on a Sr2RuO4 single crystal cleaved and measured at 10 K. The actual EDC’s were taken over more than a full quadrant of the PZ with a resolution of 0.3 (1) in the horizontal (vertical) direction. The EDC’s were then integrated over an energy window of $`\pm 10`$ meV about the chemical potential. The resulting map of 73x22 points was then symmetrized with respect to the diagonal $`\mathrm{\Gamma }`$-X (to compensate for the different resolutions along horizontal and vertical directions). The $`\alpha `$, $`\beta `$, and $`\gamma `$ sheets of FS are clearly resolved, and are marked by red lines in Fig. 3b. In addition, Fig. 3b shows some unexpected features: besides the diffuse intensity around the M point due to the presence of the SS band, there are weak but yet well defined profiles (marked in yellow). They can be recognized as a replica of the primary FS, and are related to the weak SB features detected in the EDC’s along the high-symmetry lines (Fig. 1a and 1c). This result is reminiscent of the situation found in Bi2Sr2CaCu2O8 where similar shadow bands are possibly related to AF correlations, or to the presence of two formula units in the unit cell. On the other hand, in Sr2RuO4 the origin of the SB is completely different: inspection with LEED reveals superlattice reflections corresponding to a $`\sqrt{2}\times \sqrt{2}`$ surface reconstruction (see Fig. 3a), which is responsible for the folding of the primary electronic structure with respect to the ($`\pi `$,0)-(0,$`\pi `$) direction. This reconstruction, which was found on all the Sr2RuO4 samples, is absent in the cuprates. Quantitative LEED analysis of the surface structure shows a 9 rotation of the RuO6 octahedra around the surface normal. This leads to the 45 rotation of the in-plane unit cell and to the enlargement of its dimensions by $`\sqrt{2}\times \sqrt{2}`$ over that of the bulk. The reconstruction, which reveals an intrinsic instability of the cleaved surface, should be taken into account as the origin of possible artifacts in all surface sensitive measurements.

By inspecting the M point (Fig. 3b), it now becomes clear why the investigation of this $`k`$-space region with ARPES has been so controversial: in addition to the weakly dispersive SS feature (Fig. 1 and 2), there are several sheets of FS (primary and ‘folded’). At this point, the obvious question is: what is the nature of the SS feature? It has been proposed that it could be a surface state, and in order to verify this hypothesis, we investigated its sensitivity to surface degradation by cycling the temperature between 10 and 200 K. We observed that the SS peak is suppressed much faster than all other features. Furthermore, by cleaving the crystals at 180 K and immediately cooling to 10 K we suppressed the SS, affecting the intensity of the other electronic states only weakly. A more sizable effect is observed on the SB, confirming a certain degree of surface degradation. However, this degradation was not too severe, as demonstrated by the LEED pattern taken after the measurements which still clearly shows the surface reconstruction (Fig. 4d). M-region EDC’s measured at 10 K (on a sample cleaved at 180 K), and corresponding intensity plots $`I(𝐤,\omega )`$ are shown in Fig. 4a and 4b. No signature of the SS is detected, and the identification of the Fermi vectors of $`\alpha `$, $`\beta `$, and $`\gamma `$ pockets is now straightforward. Performing a complete mapping on a sample cleaved at 180 K, we obtained an extremely well defined FS (Fig. 4c). With the surface slightly degraded, we expect to see less of the relative intensity coming from SB and SS (note that the intensity scales in Figs. 3b and 4c, although not displayed, are identical). At the same time, we might expect also the primary FS to be less well defined, which is precisely opposite to what is observed. The FS shown in Fig. 4c is in very good agreement with LDA calculations and dHvA experiments. The number of electrons contained in the FS adds up to a total of 4, in accordance with the Luttinger theorem, within an accuracy of 1% (as a matter of fact, for the FS determined on samples cleaved at 10 K the accuracy in the electron counting reduces to 3% due to the additional intensity of folded bands and surface state). As a last remark, we can confirm that $`\alpha `$ and $`\beta `$ FS present the nested topology which has been suggested as the origin of the incommensurate magnetic spin fluctuations later observed in inelastic neutron scattering experiments at the incommensurate wave vectors $`𝐐(\pm 2\pi /3a,\pm 2\pi /3a,0)`$.

Our results confirm the surface state nature of the SS peak detected at the M point. The comparison of Figs. 3b and 4c suggests that a surface contribution to the total intensity is responsible also for the less well defined FS observed on samples cleaved at 10 K. At this point, one might speculate that these findings are a signature of the surface ferromagnetism (FM) recently proposed for Sr2RuO4. In this case, two different FS’s should be expected for the two spin directions, resulting in: (i) additional $`E_F`$-weight near M due to the presence of a holelike $`\gamma `$ pocket for the majority spin; (ii) overall momentum broadening of the FS contours because the $`\alpha `$, $`\beta `$, and $`\gamma `$ sheets for the two spin populations are slightly displaced from each other in the rest of the BZ. Moreover, due to the surface-related nature of this effect, it would have escaped detection in dHvA experiments. In this scenario, a slight degradation of the surface would significantly suppress the signal related to FM correlations, due to the introduced disorder. The resulting FS would be representative of the non-magnetic electronic structure of the bulk (Fig. 4c). The hypothesis of a FM surface seems plausible because the instability of a non magnetic surface against FM order is not only indicated by ab initio calculations, but it may also be related to the lattice instability evidenced by the surface reconstruction. To further test this hypothesis, we suggest spin-polarized photoemission measurements, and both linear and nonlinear magneto-optical spectroscopy experiments \[i.e., magneto-optical Kerr effect (MOKE), and magnetic second-harmonic generation (MSHG), respectively\].

In summary, our investigation confirms the SS nature of the weakly dispersive feature detected at the M point (possibly a fingerprint of a FM surface). On the basis of both ARPES and LEED, we found that a $`\sqrt{2}\times \sqrt{2}`$ surface reconstruction occurs in cleaved Sr2RuO4, resulting in the folding of the primary electronic structure. Taking these findings into account, the FS determined by ARPES is consistent with the dHvA results and provides detailed information on the shape of $`\alpha `$, $`\beta `$, and $`\gamma `$ pockets.

We gratefully acknowledge C. Bergemann, M. Braden, T. Mizokawa, Ismail, E.W. Plummer, and A.P. Mackenzie for stimulating discussions and many useful comments. Andrea Damascelli is grateful to M. Picchietto and B. Top$`\stackrel{´}{\mathrm{ı}}`$ for their valuable support. SSRL is operated by the DOE office of Basic Energy Research, Division of Chemical Sciences. The office’s division of Material Science provided support for this research. The Stanford work is also supported by NSF grant DMR9705210 and ONR grant N00014-98-1-0195.


text:
# Clusters AgeS Experiment. Hot subdwarfs and luminous white dwarf candidates in the field of the globular cluster M41footnote 11footnote 1Based on observations collected at the Las Campanas Observatory 2.5m du Pont telescope.

## 1 Introduction

The Clusters AgeS Experiment (CASE) is a long term project with the main goal of determining accurate ages and distances to globular clusters using detached eclipsing binaries (Paczyński 1997). The photometry of new eclipsing binaries in M4 has already been presented by Kaluzny et al. (1997). In this paper we investigate the presence of blue stars in the field of the M4 globular cluster: hot blue subdwarfs and young bright white dwarfs (WD).

Recently several metal-rich globular clusters have been found to possess a handful of blue subdwarfs (Rich et al. 1997). Several scenarios have been proposed to explain the formation of such objects (Moehler 2001). The study of the properties of these stars could shed some light on the mechanisms of their creation. Only one such star, Y453 (Cudworth & Rees 1990), has been reported to date in M4.

The globular cluster M4 = NGC 6121 ($`\alpha _{2000}=16^h23^m36^s,\delta _{2000}=26^{}31\mathrm{}32\mathrm{}`$) is the nearest globular cluster. Due to its proximity M4 has been the subject of searches for WDs in the past. Several WD candidates were found by Chan and Richer (1986) in a blink survey designed to select ultraviolet-excess objects. Drukier, Fahlman & Richer (1989) presented the spectra of two of those objects. One was found to be a very hot subdwarf (sdO) and the other a field DA with weak hydrogen lines. Recently a large sample of WDs has been identified with the HST (Richer et al. 1995; 1997).

Our $`UBV`$ photometry, although inferior in depth to that obtained with the HST, covers a much larger area of M4 and thus offers a possibility to uncover new young bright WD candidates at the top of the cooling sequence, located in less crowded regions. Such objects would be within the reach of ground-based spectroscopy. The study of these objects could provide further constraints on stellar evolution models of WD progenitors. Our photometry also covers the entire central part of the cluster, not investigated in the Chan & Richer blink survey.

This paper is organized as follows: in Section 2 we describe the observations. Section 3 outlines the data reduction procedure. In Section 4 we provide details of the calibration of our instrumental magnitudes to the standard system. In Section 5 we present and examine the color-magnitude diagrams (CMDs) constructed from our $`UBV`$ photometry. In Subsections 5.1 and 5.2 we describe the luminous hot subdwarfs and the WD candidates, respectively. A differential reddening map and a fiducial sequence for the $`V/BV`$ CMD are presented in Subsections 5.3 and 5.4. In Subsection 5.5 we determine the value of the reddening $`E(BV)`$ towards M4. We discuss the variable V46 in Section 6 and present two new variables in Section 7. Section 8 contains the concluding remarks.

## 2 Observations

The photometric data were obtained with the 2.5 m du Pont telescope at Las Campanas Observatory, equipped with a thinned $`2048\times 2048`$ Tektronix CCD with a scale of $`0\stackrel{}{\mathrm{.}}26`$ pixel-1. An $`8\stackrel{}{\mathrm{.}}8\times 8\stackrel{}{\mathrm{.}}8`$ field centered on the cluster was monitored during five consecutive nights, from June 1st to 5th, 1995. Observations were made through Johnson $`UBV`$ filters. Most of the frames were collected in the $`B`$ band, with exposure times varying slightly with the observing conditions, typically 120s in $`B`$ and 60s in $`V`$. The average seeing was $`1\stackrel{}{\mathrm{.}}1`$. Newly discovered variables were presented in Kaluzny et al. (1997). A more detailed description of the collected data is presented therein.

Spectra of the blue subdwarfs B1, B2 and B3 were obtained at the Baade 6.5 m telescope with the Boller & Chivens spectrograph equipped with Tek#1 CCD. The wavelength coverage extends from 3774 to 4492 Åat a resolution of 0.7 Å. Two spectra were taken for each of the stars, on May 8th, 2001 for B1 (500 s) and B2 (600 s) and on June 16th and 19th, 2001 for B3 (1200 s). The spectra were combined to increase the signal to noise ratio (S/N) which ranged from $`16`$ to $`23`$ for single exposures.

A 30-minute 6.3 Å resolution spectrum of V46, extending from 3753 to 6961 Å was obtained on the night of April 25th, 1996 with the 3.6m ESO telescope using the EFOSC1 spectrograph.

## 3 Data reduction

The photometry was carried out using Daophot/ Allstar (Stetson 1994). The positional variability of the point-spread function was modeled with a second order approximation. The extraction of profile photometry closely followed the procedure described in Mochejska & Kaluzny (1999).

The $`V`$ band photometry was obtained from an image combined from three 300 s exposures with subarcsecond seeing. The images were transformed to a common coordinate system using the ISIS image subtraction package (Alard 2000) and median combined with IRAF222IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the NSF. using an averaged sigma clipping algorithm to reject cosmic rays. To recover some of the saturated stars, the photometry from an averaged $`3\times 60`$ s image, a 20 s and a 6 s exposure was included. A similar approach was followed for the $`U`$ and $`B`$ band data. The final $`U`$ band photometry was derived from a median combined $`4\times 900`$ s image and a 480 s exposure. For the $`B`$ band we used an averaged $`2\times 800`$ s image, a median combined $`3\times 120`$ s image and a 30 s exposure. The full list of exposures used in the construction of the CMD is provided in Tab. 1.

The photometry lists derived from the images for each filter were transformed to a common instrumental system defined by the longest exposure in each filter. Saturated stars were removed from the long exposures. From the shorter exposure lists we removed stars fainter than the magnitude at which the average formal error on the short exposure was twice the error on the longest exposure. In the final step the lists were added, with the best of the available measurements (i.e., having the smallest formal error) retained for stars present on more than one list. In addition we removed 24 RR Lyrae variables from the photometry lists.

## 4 Calibration

The color terms of the transformations of instrumental magnitudes to the standard system were derived from 19 observations of Landolt (1992) fields collected on 5 nights spread between April 22/23 and June 3/4, 1995. The same instrumental setup (CCD+filters) was used on all nights. The following values were obtained by averaging color terms derived for separate nights:

$`v=V0.0160(BV)+c_1`$ (1)
$`bv=0.9440(BV)+c_2`$ (2)
$`ub=0.9480(UB)+c_3`$ (3)
$`uv=0.9450(UV)+c_4`$ (4)
The constant offsets for $`V`$ and $`BV`$ were determined relative to M4 data collected with the same telescope on the night of May 31/June 1 1998 at airmasses of 1.07-1.11. The transformation coefficients were derived from 56 stars in 6 Landolt fields observed that night at airmasses ranging from 1.07 to 1.75. Figure 1 shows the residuals of the standard stars resulting from the following transformations:

$`v`$ $`=`$ $`V+0.69030.0176(BV)`$
$`+0.1284(X1.25)`$
$`bv`$ $`=`$ $`0.3057+0.9507(BV)`$
$`+0.0731(X1.25)`$
The offsets for $`UB`$ and $`UV`$ were determined from a comparison with M4 data collected with the 1m Swope telescope at LCO at airmasses less than 1.06. We observed 39 stars in 9 Landolt fields at airmasses ranging from 1.07 to 1.36. The extinction coefficients were assumed as the average values for the winter season at LCO. The following transformations were adopted:

$`v=V+2.78370.0189(BV)+0.15(X1.25)`$ (7)
$`bv=0.3058+0.9359(BV)+0.12(X1.25)`$ (8)
$`ub=2.0252+0.9377(UB)+0.31(X1.25)`$ (9)
$`uv=2.3269+0.9313(UV)+0.44(X1.25)`$ (10)
Figure 2 shows the $`V`$, $`UB`$ and $`UV`$ transformation residuals of the standard stars.

As a check of consistency of the $`V`$ and $`BV`$ transformations we have compared the May 31/June 1 1998 calibrated M4 data (Eqs. (4) and (4)) with the June 1995 data transformed to the standard system using the color terms in Eqs. (1) and (2) and offsets relative to the transformations given by Eqs. (7) and (8). For stars with $`V<16`$ mag the offsets in $`V`$ and $`BV`$ are -0.026 and -0.020 mag, respectively. There is no correlation between the residuals and color.

We have compared our $`V/BV`$ diagram with a CMD constructed from the photometry obtained by Alcaino et al. (1997) in a $`6\mathrm{}\times 6\mathrm{}`$ field bordering our field on the west (data courtesy of Dr N. Samus). As Alcaino et al. report the existence of a reddening gradient across the face of the cluster, with the reddening increasing to the west, we selected for the comparison the westernmost $`2\stackrel{}{\mathrm{.}}2\times 8\stackrel{}{\mathrm{.}}8`$ and easternmost $`2\mathrm{}\times 6\mathrm{}`$ sections of our and Alcaino et al. fields, respectively. We fitted fiducial points to the two datasets in the magnitude range $`16.8<V<20`$. The mean offset between the two sequences is -0.007 mag. In Figure 3 we overplot the Alcaino et al. $`V/BV`$ CMD with fiducial points fitted to our data.

In Figure 4 we compare our $`UV`$ photometry with that obtained with the HST (Ibata et al. 1999; Richer et al. 1997). The offsets in the photometry are 0.095 mag in $`V`$ and 0.129 in $`U`$ for stars with $`19<V<21`$ mag and $`18<U<22`$ mag, respectively. A similar systematic offset between ground-based and HST $`V`$ magnitudes has been noted by Richer et al. (1997; Figure 2 therein). For stars above $`V=20`$ mag the difference in photometry dramatically increases with magnitude.

## 5 Color-magnitude diagrams

The $`V/BV`$, $`V/UB`$ and $`V/UV`$ color-magnitude diagrams (CMD) constructed from our data are presented in Figure 5. The CMDs exhibit a well defined main sequence that can be traced down to $`V22`$ mag on the $`V/BV`$ CMD and $`V21`$ mag on the remaining two diagrams. The subgiant, red giant and horizontal branches are well populated, but there are few stars on the asymptotic giant branch. The horizontal branch exhibits a clear RR Lyrae gap. A sequence of four luminous blue stars is present on the CMDs (plotted as open triangles). These objects are discussed in more detail in Subsection 5.1. There is also a population of faint blue stars with magnitudes and colors ($`V<21`$ mag, $`UV<0.4`$ mag) consistent with the values expected for the upper part of the white dwarf cooling sequence. A brief discussion of these objects is presented in Subsection 5.2. In Subsections 5.3 and 5.4 we present a differential reddening map and a sequence of fiducial points determined for the $`V/BV`$ CMD. In Subsection 5.5 we estimate the reddening $`E(BV)`$ towards M4.

### 5.1 Hot subdwarfs

We note a sequence of four luminous blue stars on the CMDs ($`V<20`$, $`UV<0.6`$), most probably hot subdwarfs belonging to the cluster (marked with triangles on the CMDs in Figure 5). The brightest of these, B1, was initially discovered by Cudworth & Rees (1990; Y453 in their catalog). No other stars as blue as these were found in a blink survey aimed at the detection of UV excess objects in an area between $`3\stackrel{}{\mathrm{.}}9`$ and $`23\stackrel{}{\mathrm{.}}7`$ from the cluster center (Chan & Richer 1986). Four bright blue objects with $`V<20`$ and $`UV<0.1`$ were found in that study: 1710 (reported as non-stellar), 774, 1285 and 710 (see Table 2 therein). The star 774 was investigated spectroscopically by Drukier et al. (1989) and found to be most likely a field WD. One fainter object, 831 with $`V=20.33`$ and $`UV=0.97`$ was found to be a very hot subdwarf by Drukier et al. (1989).

In Figure 6 we present the spectra of B1, B2 and B3, smoothed with a boxcar filter 5 units in length. Helium lines are indicated with continuous lines and hydrogen lines with dotted ones. B1 (Y453) has been previously studied spectroscopically by Moehler, Landsman & Napiwotzki (1998). Their analysis indicates that it could be a post-early asymptotic giant branch (PEAGB) star with $`T_{eff}=58800`$ K. Our spectrum shows broad Balmer lines, but the He I and He II lines reported in the Moehler et al. (1998) data are barely discernible, most likely due to our higher resolution and moderate S/N. B2 shows singly and doubly ionized helium lines. The lines at wavelengths corresponding to the hydrogen Balmer series are most likely due to doubly ionized helium, judging by the presence of the He II line at 4200 Å. It appears to be a very hot helium rich subdwarf, approximately sdO5:He4 in the Jeffery et al. (1997) classification scheme. The precise classification is somewhat uncertain, as it is based on only three lines in the 300 Å overlap between our and their spectra: the He II line at 4339 Å and He I lines at 4388 Å and 4471 Å. The spectrum of B3 exhibits broad Balmer lines and two lines from singly ionized helium (4026 Å and 4471 Å) and none from He II. It appears to be an sdB:He1 subdwarf.

The spectrum of B2, which appears to he helium-rich, is quite unusual. Hot subdwarfs in globular clusters are almost always found to be helium-poor (e.g. Moehler et al. 1997a). Only one hot helium-rich subdwarf in a globular cluster has been reported to date (Moehler et al. 1997b). B2 seems to be only the second such star known. Brown et al. (2001) suggest that such stars might have undergone late helium-core flash while descending the WD cooling sequence. A convection zone produced during the helium flash would mix hydrogen into the hot interior, where it would be consumed, thus greatly enhancing the envelope helium abundance. However, B2 seems to be too luminous for their scenario.

The hot subdwarf B3 and possibly B2 could belong to a vertical extended horizontal branch (EHB), not reported in previous studies of this cluster. Generally the globular cluster horizontal branch (HB) morphology correlates with metallicity - the HB becomes redder as the mean metal abundance increases. However, variations in HB morphology at a given metallicity are observed and a second, still unknown parameter has to be introduced. It has been suggested that mass loss on the red giant branch may be the second parameter (Catelan 2000). The stars which suffered strong mass loss would become EHB stars (D’Cruz et al. 1996, D’Cruz 1998). Recently the existence of high-metallicity hot blue horizontal branch stars has been reported in the relatively metal-rich globular clusters NGC 1851 (Bellazzini et al. 2001) and 47 Tuc (Ferraro et al. 2001), the very old metal-rich open cluster NGC 6791 (Kaluzny & Udalski 1992) and the Milky Way bulge (Peterson et al. 2000). Detailed study of the nature of such extreme horizontal branch stars could provide valuable clues regarding the unknown second parameter.

An alternative explanation for B4 is that it is a low-mass helium WD. Similar stars have been identified in the globular cluster NGC 6397 (Cool et al. 1998) and the spectrum of one of them was found to be consistent with a He WD (Edmonds et al. 1999). Using differential surface density data from HEQS (Köhler et al. 1997), LBQS (Hewett et al. 1995) and the 2dF QSO survey (Boyle et al. 2000) we estimate the probability of B4 being a quasar to be between 2% and 48%.

We have checked the subdwarf candidates for variability in our observations spanning 5 days and found them to be constant within the errors of the photometry.

In Table 2 we list the photometry and equatorial coordinates for the hot subdwarfs. Finding charts for these stars are are presented in Figure 14. Each chart covers $`33\mathrm{}\times 33\mathrm{}`$ on the sky.

### 5.2 White dwarf candidates

A population of faint blue stars is present on our CMDs. In Figure 7 we overplot these stars (filled circles) with the WDs observed with the HST (open circles, Richer et al. 1997). Their location is consistent with the upper part of the HST white dwarf cooling sequence. Figure 4 is a comparison of our $`UV`$ photometry and that of the HST, open circles indicate WDs that are common to the two studies.

The photometry and equatorial coordinates for the five brightest objects with $`V22`$ (marked with open circles in Figure 5) are presented in Table 3333The full table is available from the authors via anonymous ftp from cfa-ftp.harvard.edu in the directory /pub/bmochejs/M4/  and $`33\mathrm{}\times 33\mathrm{}`$ $`V`$ and $`U`$ band finding charts are presented in Figure 13. These objects are located in the outer, less crowded regions of the cluster (compared to the HST WD sample) and are bright enough for ground-based spectroscopy.

### 5.3 Differential Extinction Map

All of the principal sequences on the CMDs display a broadening due to the differential reddening across the face of the cluster, as discussed by Cudworth & Rees (1990), Fahlman et al. (1996), and Alcaino et al. (1997). The differential reddening has been investigated by Lyons et al. (1995) using Ki column density measurements towards 16 stars and by Ivans et al. (1999; 2000) on the basis of spectroscopic analysis of 22 red giant branch stars. These studies revealed the existence of a patchy reddening variation across the cluster face. In addition, Peterson, Rees & Cudworth (1995) found evidence for an unusually high value of $`R_V=A_V/E(BV)4`$.

We have constructed an $`8\times 8`$ element differential reddening map following closely the procedure presented by von Braun & Mateo (2000), based on the approach of Kaluzny & Krzeminski (1993). The errors were also derived as described therein, taking into account the correlation between the errors in color and magnitude. We adopted $`R_V=4`$ in our analysis. Figure 8 presents the derived map of differential reddening, measured relative to the section with the lowest reddening. The extinction appears to be weakest in the east and increases to the west, in accordance with results obtained in previous studies (eg. Alcaino et al. 1997). The patchy nature of the reddening is also confirmed.

We have compared our differential extinction map with the Schlegel, Finkbeiner & Davis (1998; hereafter SFD) dust map. There is some qualitative similarity between the two maps: both show lowest reddening in the northeastern corner, intermediate in the southwestern and higher in the northwestern corners. The maps disagree as to the reddening in the southeastern section: the SFD map shows the strongest reddening there, while it is rather small on our map. Another difference is that the differential reddening is larger on our map (up to 0.074) than on the SFD map (up to 0.039). We have checked a larger area around M4 on the SFD map and noted that it is located in a region of highly variable extinction. The cluster coincides with the center of a void with patchy lower extinction ($`E(BV)0.5`$) some $`1^{}\times 1^{}`$ surrounded by thick filaments of higher extinction ($`E(BV)>1`$ mag). Moreover, we notice that the $`E(BV)`$ values from the SFD map range from 0.472 to 0.511 mag for the cluster area, which is substantially higher than the reddening value of $`E(BV)=0.36`$ mag reported for M4 by Harris (1996). This discrepancy might be due to the fact that the SFD map also takes into account the extinction caused by dust located behind the cluster. This interpretation would invalidate to large extent any comparison between the two reddening maps.

### 5.4 $`V/BV`$ CMD Fiducial Sequence

To derive the $`V/BV`$ fiducial sequence we have selected a region of uniform reddening using the differential reddening map constructed in Subsection 5.3. This region is shown in Figure 8 surrounded with a thick line. The fiducial points were determined separately for the main sequence (MS; $`16.75<V<21.0`$), the subgiant branch (SGB; $`16.2<V<16.75`$) and the red giant branch (RGB; $`13.0<V<16.2`$). For the MS and SGB, the fiducial points were determined by finding the mode of the color distribution, smoothed with a boxcar filter 0.03 mag in length, in magnitude bins 0.21 and 0.135 mag wide, respectively. The magnitude for each bin was computed as the mean. In case of the RGB we adopted a different procedure, as this branch is more sparsely populated than the previous two. To avoid contamination of the RGB sample, especially with the AGB and HB stars, the following magnitude-dependent color cuts were applied: $`(29.7V)/14.5<BV<(32.6V)/14.5`$. For the RGB, the fiducial points were computed as the mean color in magnitude bins 0.55 mag wide. The magnitudes within each bin were determined as the mean, just as for the MS and SGB. The resultant fiducial sequence is presented in Figure 9, superimposed on the CMD which was used in their derivation. The fiducial sequence is listed in Table 4.

### 5.5 Foreground Extinction E(B-V)

We have obtained an independent estimate of foreground extinction in B-V towards M4 from comparison of the fiducial sequence constructed in the previous Subsection with the one given by Sandquist et al. (1996) for the M5 globular cluster. We have chosen M5 for this purpose due to its low reddening (E(B-V) = 0.03; Harris 1996), similar metallicity (\[Fe/H\] = -1.11 dex as compared with -1.19 dex for M4; Carretta & Gratton 1997) and similar age ($`\tau _{M4}\tau _{M5}=0.6`$ Gyr; Rosenberg et al. 1999).

To determine the $`BV`$ colors of the cluster turnoffs (TO) we fitted the fiducial points in the turnoff regions with a fifth order Chebyshev polynomial. The color of the turnoff was taken as the color of the bluest point on the fitted curve. We obtained the values of 0.854 and 0.470 for M4 and M5, respectively. From the isochrones of Girardi et al. (2000) we estimated the corrections in $`(BV)_{TO}`$ for the difference in metallicity (-0.08 dex) and age (0.6 Gyr) between M4 and M5 to be -0.014 and +0.006 mag, respectively. Assuming an extinction of 0.03 mag towards M5, we obtain a reddening value of $`E(BV)=0.41`$ for M4. This is consistent with the value of 0.36 listed by Harris (1996), taking into account the variability of the extinction across the face of the cluster. It is in disagreement with the total extinction estimate of 0.5 mag based on the SFD dust maps.

## 6 The variable hot subdwarf V46

V46 was detected to be variable by Kaluzny et al. (1997). The star exhibits coherent modulations with a period of 0.04358 d. It is indicated on the CMDs in Figure 5 with an open square and mean photometry is listed in Table 5. We have measured the light curve of V46 using the ISIS image subtraction package (Alard & Lupton 1998; Alard 2000). The upper panel of Figure 10 shows the raw light curve for the first night of observations. The phased light curve from all five nights is shown in the lower panel. The shape of the light curve appears to be unstable. We have searched for other coherent periodicities and have found none. Nevertheless, we cannot rule out multiple periodicity, due to the short timespan of our observations.

The spectrum of V46, shown in Figure 11, exhibits strong and broad Balmer lines, while the He I and He II lines are absent. We have compared it with spectra for stars of types B through F, taken from the spectral atlas of Pickles (1998). Our spectrum is consistent with types from late B through late A. The dereddened colors $`(BV)_0=0.13`$ and $`(UB)_0=0.32`$ suggest V46 is a late B type star. In order to be a main sequence late B star ($`M_V1.5`$), at $`V_017.1`$ it would have to be located at a distance modulus of $`18.6`$ mag or 52 kpc (assuming the larger $`E(BV)=0.5`$ from Schlegel et al. (1998)). At $`b=15.97`$ this corresponds to $`Z=14`$ kpc. This suggests that V46 is rather a hot subdwarf, type sdB1 in the Jeffery et al. (1997) classification scheme.

Although the spectrum of V46 indicates that it is an sdB, its period of brightness modulations (0.04358d=3765s) is much longer than the typical periods of 100-200s (500s in extreme cases) for pulsating subdwarfs (Koen et al. 1998).

Another possibility to consider are AM CVn type interacting binary WD systems. AM CVn itself exhibits a low amplitude flickering light curve with a period of 0.01216d=1051s (Krzeminski 1972; Provencal et al. 1995). GP Com (G61-29) is the longest period variable of this class, with P=0.03231d=2791s (Marsh 1999). However, AM CVn type variables show He I features in their spectra, accompanied by a lack of hydrogen lines. This is in disagreement with the V46 spectrum and thus renders this interpretation very unlikely.

A different interpretation, also suggested more by the light curve than by the spectrum, is that V46 is a field SX Phe star. Its period is typical for these variables. The light variations appear remarkably similar to those of SX Phe stars (i.e. Breger et al. 1995, Handler et al. 1996, Pych et al. 2001, Kaluzny & Thompson 2001). The location of V46 on the CMD (Fig. 5) would suggest that it should be located behind the cluster.

An alternative interpretation could be that V46 is a binary star with the brightness modulations due to the reflection effect or to ellipsoidal variations (P=0.08716d=2.09h). A study of the radial velocities of a sample of 70 sdBs indicates that 45% of them are post-common envelope binaries with periods of the order of a few hours to a few days (Saffer, Green & Bowers 2000). The period for V46 is consistent with this period range. Spectra with better time resolution are needed to verify this interpretation.

## 7 Other variables

The data set was reanalyzed using the ISIS image subtraction package (Alard & Lupton 1998, Alard 2000) mainly with the aim of obtaining a better light curve for V46. The light curves of all previously found variables were also extracted and two new variables were found. There is notable improvement in the quality of the photometry extracted with ISIS, as compared to traditional profile photometry. In Figure 12 we show the new ISIS light curve for the previously discovered variable V52, phased with its newly derived period of 0.7765 days (Tab. 5). The light curves of the two new variables V56 and V57 are also presented in Figure 12, their photometry is listed in Tab. 5 and finding charts are presented in Figure 13. On the CMDs (Figure 5) these variables are located on the RGB (star symbols).

## 8 Conclusions

Our investigation of the M4 $`UBV`$ color magnitude diagrams has resulted in the identification of four luminous hot subdwarfs (one previously known) and the tip of the white dwarf cooling sequence. The hot blue stars could possibly form a vertical extended horizontal branch of the cluster. The spectra for the three brightest blue stars confirm that they are hot subdwarfs. Further study would be required to elucidate the nature of the fourth object. We have selected five luminous WD candidates above $`V22`$ mag, located in the outer, less crowded regions of M4. These objects would make good targets for follow-up ground-based spectroscopy. In addition we present a fiducial sequence for the $`V/BV`$ CMD and a differential reddening map. We obtain an estimate of $`E(BV)=0.41`$ for the foreground reddening towards M4 based on a comparison with M5, a value consistent with previous determinations.

We also studied the variable star V46, discovered by Kaluzny et al. (1997) in the same data set. We present and discuss its spectrum, which indicates that it is a hot subdwarf and present a new light curve, derived with the ISIS image subtraction package (Alard 2000). We also identify two new variables, both lying on the cluster RGB.

W. Landsman, G. Preston, and K.Z. Stanek provided us with helpful comments on the manuscript. We would like to thank H. Duerbeck for help with the spectroscopic observations, G. Pojmański for $`lc`$ \- the light curve analysis utility and N. Samus for the comparison $`BV`$ photometry. We also thank the referee for a very useful report, which significantly improved this paper. BJM, JK and WP were supported by the Polish KBN grant 5P03D004.21. BJM was also supported by the Foundation for Polish Science stipend for young scientists. JK was also supported by the NSF grant AST-9819787. IT was supported by the NSF grant AST-9819786.