text:
# Critical State Behaviour in a Low Dimensional Metal Induced by Strong Magnetic Fields

## I introduction

Of all the quasi-two-dimensional (Q2D) organic charge-transfer salts that exist , those of the composition $`\alpha `$-(BEDT-TTF)$`{}_{2}{}^{}M`$Hg(SCN)4, where BEDT-TTF stands for bis(ethylenedithio)tetrathiafulvalene and where $`M=`$ K, Tl, Rb or NH4, have been, perhaps, the most difficult to understand . While, the $`M=`$ NH4 salt is a superconductor with a transition temperature $`T_\mathrm{c}`$ 1 K , the $`M=`$ K, Tl and Rb salts all undergo a transition into a more resistive state with a reconstructed Fermi surface below $`T_\mathrm{p}`$ 8 (in the $`M=`$ K and Tl salts) or 10 K (in the $`M=`$ Rb salt) . The underlying physical reason for the behavioural differences between the isostructural $`M=`$ NH4 and $`M=`$ K, Tl and Rb salts remains a contemporary issue. The $`M=`$ K salt has, nevertheless, been shown to become superconducting under uniaxial stress applied perpendicular to the conducting layers , and it has further been suggested that the $`M=`$ K and Rb salts could exhibit filamentary superconductivity even at ambient pressure .

Unquestionably, the main physical traits of the $`M=`$ K, Tl and Rb salts at ambient pressure, at moderately low temperatures $`TT_\mathrm{p}`$ and at magnetic fields $`B<B_\mathrm{k}`$, where $`B_\mathrm{k}`$ ($``$ 23 T in the $`M=`$ K salt) is known as the kink transition field, are more typical of a density wave (DW) ground state . Yet, no direct evidence for a superlattice structure has been found. Most of the more recent experimental and theoretical surveys point towards a charge-density wave (CDW) rather than a spin-density wave (SDW) being the more likely candidate . The electrical transport is strongly magnetoresistive, with the increase in resistance below $`T_\mathrm{p}`$ becoming increasingly pronounced in a magnetic field . On the approach to $`B_\mathrm{k}`$, however, the magnetoresistance changes, becoming much lower above this field . The nature of the high magnetic field phase ($`B>B_\mathrm{k}`$), at low temperatures, has since been the subject of speculation . The only fully established fact is that the transition field $`B_\mathrm{k}`$ is strongly first order, with hysteresis appearing both in the magnetoresistance and the magnetisation over a few Tesla interval in field.

From the theoretical perspective, $`B_\mathrm{k}`$ lies remarkably close to the critical field $`B_\mathrm{c}`$ at which one would expect a CDW, with a mere transition temperature of 8 K, to approach the Pauli paramagnetic limit . A first order transition could be expected, with the material transforming at higher fields either into a normal metal , a spatially modulated CDW  or a mixed CDW-SDW hybrid . The latter, in particular, represents the CDW analogue of the Fulde-Ferrel phase that is anticipated to occur in strongly type II superconducters, with suppressed orbital effects, over the equivalent region of the phase diagram . Were either of these latter two phases actually to be realised at high magnetic fields, their physical properties have not yet been the subject of a thorough investigation. Experimentally, data has been published that could be considered conducive to the existence of a different thermodynamic phase above $`B_\mathrm{k}`$ at low temperatures , albeit that the magnetic data, at that time, was discussed largely in the context of SDW’s with no metion of CDW’s. Yet other data, concerning the unusual beheviour of the interlayer magnetotransport and induced currents in the magnetisation in pulsed magnetic fields, has been notionally connected with the quantum Hall effect (QHE) . While a definitive experiment has not been performed that can firmly establish either of these latter two scenarios at high magnetic fields ($`B>B_\mathrm{k}`$), one cannot dispute the fact that this high magnetic field phase is somewhat exotic.

With a view to understanding more about the high magnetic field phase (at fields above $`B_\mathrm{k}`$), in this paper we present the results of extensive measurements of the interlayer magnetotransport and the magnetic torque on several $`\alpha `$-(BEDT-TTF)2KHg(SCN)4 crystals over a broad range of temperatures, 20 mK $`<T<`$ 10 K, and in continuous magnetic fields of up to 33 T. The observed behaviour of the magnetotransport and magnetic torque is neither typical of that expected for a DW system, nor does it support recent claims of the QHE in this material. For example; as the sample is cooled in a magnetic field, the resistivity undegoes an abrupt drop at $``$ 3 K, at all fields $`BB_\mathrm{k}`$, with the drop being particularly pronounced at integer Landau level filling factors, falling exponentially by roughly two orders of magnitude. The occurrence of a change in slope in the magnetic torque at precisely the same temperature, on field-cooling the sample, lends further support to the existence of a different low temperature thermodynamic phase below $`T_\mathrm{c}`$ 3 K . At lower temperatures still, (i.e. $`T`$ 2 K), the field-dependence of the magnetic torque develops a pronounced hysteresis at all fields above $`B_\mathrm{k}`$ that increases approximately exponentially with decreasing $`T`$, but with the hysteresis being particularly strong at half integer Landau level filling factors. Only at much lower temperatures ($`T`$ 27 mK) does the hysteresis become more pronounced at integer filling factors . While the hysteresis throughout the entire high magnetic field phase is difficult to explain in terms of conventional DW ground states or the QHE, it has many of the features of a critical state model like those used to explain magnetic hysteresis in type II superconductors. In the absence of any other model specific to the high magnetic field phase of $`\alpha `$-(BEDT-TTF)2KHg(SCN)4, the irreversible hysteretis, together with the sharp drop in the resistivity at low temperatures, is remarkably similar to the behaviour of superconductor in strong magnetic fields. Such behaviour cannot be explained consistently in terms of the QHE or any conventional DW ground state.

## II experimental

The samples chosen for this study, referred to hereon as 1, 2 and 3, were grown using conventional electrochemical means , while static magnetic fields extending to $``$ 33 T were provided by the the National High Magnetic Field Laboratory, Tallahassee. Temperatures between $``$ 20 mK and $``$ 1.6 K were provided by a top-loading dilution refrigerator, with higher temperatures obtained using a 3He refrigerator. Magnetotransport measurements were made using standard four wire techniques with low frequency ($``$ 10 Hz) currents of between 1 and 10$`\mu `$A applied perpendicular to the conducting layers. For the magnetic torque measurements, the samples were mounted on the moving plate of a phosphor bronze capacitance cantilever torque magnetometer, which was itself attached to a rigid but rotatable platform in such a way that the axes of torque and rotation were parallel to each other, but perpendicular to the applied magnetic field $`B`$. The angle between $`B`$ and the normal to the capacitance plates was approximately the same as the angle $`\theta `$ between $`B`$ and the normal to the conducting planes of the sample. The capacitance of $``$ 1.3 pF was measured by means of a capacitance bridge energized at 5 KHz and 30 V (rms); the largest change in C observed throughout the experiment at the highest magnetic field was $``$ 0.04 pF, corresponding to a net angular displacement of $``$ 0.1. Torque interaction effects at $`\theta =`$ 7 were therefore not a significant factor. To eliminate possible artefacts due to time constants of the instrumentation setup or the data acquisition system, the capacitance was measured with the time constant of the lock-in amplifier set to a low value of 10 ms. Furthermore, both the output of the lock-in amplifier and the shunt voltage that determines the current flowing in the Bitter magnet coils, were measured using a digitizer with resolution of $``$ 50 $`\mu `$s, while the field was swept slowly at $``$ 8 mTs-1.

## III temperature-dependence of the magnetoresistance

Examples of the interlayer magnetoresistance measured in $`\alpha `$-(BEDT-TTF)2KHg(SCN)4 sample 1 at selected temperatures are shown in Fig. 1. The magnetoresistance displays the usual behaviour observed for this material, reaching a maximum at $``$ 10 T before falling again in a linear fashion as the kink transition field is approached . In accordance with these earlier studies, the oscillations within the low magnetic field DW phase exhibit a pronounced second harmonic. Above the kink transition, at temperatures higher than $``$ 3 K, the oscillations grow rapidly in amplitude with increasing field. The behaviour of the interlayer Shubnikov-de Haas (SdH) waveform over this range of temperatures is well understood in these materials . At integer Landau level filling factors, when the chemical potential $`\mu `$ is situated in a Landau gap, the SdH maxima increase with decreasing temperature in a insulating-like fashion. Conversely, at half-integer filling factors, the resistivity of the minima behave in a metallic fashion. This is entirely consistent with the theoretically expected behaviour of a quasi-two-dimensional metal in a magnetic field .

What is not predicted by the simple theory is the abrupt inversion of the waveform at low temperatures, resulting in the high temperature ($`T`$ 3 K) SdH maxima becoming minima at low temperatures ($`T`$ 3 K) with a resistivity several times lower than that of the sample at zero field. Figure 1 does not constitute the first observation of this effect; previous studies had reported this effect in both continuous and pulsed magnetic fields . Because no phase inversion is observed in the dHvA effect , which is entirely a thermodynamic function of state, this implies that the Landau level structure remains largely unchanged over the same temperature range. It was the occurrence of this phase inversion only in the magnetotransport that was then attributed to the effect of a chiral Fermi liquid , following a suggestion that the QHE may be taking place at high magnetic fields in this material . Because the interlayer resistance of the bulk of the sample becomes insulator-like at integer filling factors according to magnetotransport theory , higher conductivity chiral Fermi liquid states were proposed to take over the majority of the interlayer conductance. For this conjecture to explain the most recent data (i.e. Fig. 1), the interlayer conductivity of the edge states, which occupy $``$ 1 part in 104 of the sample cross section, would have to be both strongly temperature-dependent and attain a conductivity at least 106 times higher than that of the field-averaged (or background) conductivity of the bulk.

Since this conjecture was made, the interlayer conductance of chiral surface states has been measured directly in semiconductor superlattices that at the same time quite clearly exhibit the QHE . Not only is the conductance of these states found to be temperature-independent, as predicted by chiral Fermi liquid theory , but their net conductance is also observed not to be particularly high, resulting in only a weak suppression of the SdH maxima. This behaviour is therefore quite different from that observed in $`\alpha `$-(BEDT-TTF)2KHg(SCN)4. When coupled with the fact that quantized Hall plateaux proportional to $`\mathrm{}/ie^2`$ have not been observed  (in the $`M=`$ Tl salt), the chiral Fermi liquid model can no longer be considered valid.

A more distinct behaviour emerges when the resistivity, both at integer and half-integer filling factors, is plotted versus temperature in Fig. 2. For two different samples, the resistivity varies weakly with temperature for $`T`$ 3 K, but then undergoes an abrupt drop for $`T`$ 3 K. This drop is particularly pronounced at integer filling factors (i.e. when $`\mu `$ is situated in the Landau gap of the Q2D pockets), but is also clearly discernable at half-integer filling factors (i.e. when $`\mu `$ is situated in the middle of a Landau level). This behaviour is somewhat similar to that observed by Kartsovnik et al. , except that in the present case, at integer filling factors in sample 1, the resistivity is observed to fall by two orders of magnitude between 3 K and 490 mK, reaching 1.6 $`\mathrm{\Omega }`$ at the lowest temperature. While 1.6 $`\mathrm{\Omega }`$ does not sound that low, it is $``$ 2000 times lower than the room temperature resistance of 3.6 k$`\mathrm{\Omega }`$ for this sample, and this is even in the presence of a $``$ 30 T magnetic field. In spite of the fact that the samples were cooled in a dilution refrigerator on a subsequent experimental run to temperatures as low as 20 mK, zero resistivity was never observed. The abrupt change in resistivity versus temperature $`\rho (T)`$ has previously been interpreted as evidence for a phase transition , and the present study appears to support this claim. An abrupt drop in the resistance, as observed here, cannot be explained by chiral Fermi liquid theory , nor is the occurrence of a transition at half-integer filling factors (as well as integer filling factors) compatible with the QHE. Thus, it would appear to be the case that the inversion of the SdH waveform is merely an artefact of the drop in resistivity being more pronounced at integer filling factors.

From a phenomenalogical point of view, the low temperature interlayer resistivity in Fig. 2 obeys an approximate $`\rho \mathrm{exp}(T/T_0)`$ law, with $`T_0`$ 520 mK at integer filling factors for both samples 1 and 2. Such an exponential law is clearly unlike that of a normal metal, and the very observation of a sudden drop in resistivity is not expected for a DW ground state. Transitions into DW ground states are invariably followed by an immediate increase in resistivity on decreasing the temperature through the transition , owing to a net loss of carriers. A drop in resistivity, of the form observed here, more commonly precedes a superconducting state. One could argue that, while a zero resistivity state is not observed, zero resistance states are always more difficult to come by in strong magnetic fields owing to vortex dynamics . It is also interesting to note that the exponential dependence of the resistivity at low temperatures is strikingly similar to that observed in the filamentary superconductor La2CuO4-y .

## IV temperature-dependence of the magnetic torque

The temperature-dependence of the magnetotransport alone cannot be considered as proof for a transition into a new thermodynamic phase. This usually requires measurement of a thermodynamic function of state, such as specific heat or magnetisation. Specific heat measurements have thus far confirmed the existence of a second order transition into a low temperature phase within the low magnetic field regime ($`B<B_\mathrm{k}`$) . Within this same phase, magnetic moment measurements, made using a superconducting quantum interference device (SQUID) in fields of 5 T, have shown that the susceptibility drops most notably when the field is oriented within the conducting planes . Such a behaviour could be expected either for a CDW or a SDW phase following a net loss in Pauli paramagnetism with the opening of gaps on the Fermi surface . Only a very weak change, however, was observed when the magnetic field was oriented perpendicular to the conducting layers. Sasaki et al. notionally connected this anisotropy with antiferromagnetism associated with a SDW ground state , although there has since not been any compelling evidence supporting this hypothesis .

An alternative, yet rather trivial, explanation could be that when the field is oriented perpendicular to the conducting layers, the loss in Pauli paramagnetism is balanced by a loss in Landau diamagnetism following the partial opening of gaps on the Q2D Fermi surface pockets due to a CDW. Gaps on the Q2D pocket are, after all, predicted by all models of the reconstructed Fermi surface . Perhaps the strongest evidence for the latter explanation is that the magnetic susceptibility is completely isotropic for fields oriented within the planes, as demonstrated by the magnetic torque measurements of Christ et al. . This detail would be difficult to explain in terms of antiferromagnetism, but is quite easy to explain in terms of Landau diamagnetism, which in a Q2D conductor, manifests itself only perpendicular to the conducting layers.

The essential advantage of magnetic torque measurements is that they are sensitive only to the magnetic anisotropy, with the net torque being given by the product

$$\tau =𝐌\times 𝐁.$$
(1)
Thus, if we consider the susceptibility to be resolved into its components perpendicular to the conducting layers $`\chi _{}M_{}/B`$ and parallel to the conducting layers $`\chi _{}M_{}/B`$, the net magnetic torque is then given by

$$\tau =\frac{1}{2}\left[\chi _{}\chi _{}\right]B^2\mathrm{sin}2\theta ,$$
(2)
where $`\theta `$ is the angle between the magnetic field and the normal to the conducting layers (i.e. the $`𝐛`$ axis of the crystal). With the first term of Equation (2) being dominant and negative, according to Sasaki et al. , the change in anisotropy on entering the low temperature DW phase exerts a net negative torque on the sample, that attempts to align the $`𝐛`$ axis of the sample more closely with $`B`$. This is exactly the signature observed by Christ et al. . For a constant anisotropy $`\chi _{}\chi _{}`$, the negative torque should increase in a manner that is proportional to $`B^2`$. The fact that this is not the case, as shown originally by Christ et al.  and again in Fig. 3, could be considered as further supportive evidence for our explanation in terms of Landau diamagnetism, whereby the Landau diamagnetic contribution from the Q2D pocket returns gradually at higher magnetic fields due to magnetic breakdown across the CDW gaps.

Even though the anisotropy becomes reduced at high magnetic fields, a change in magnetic torque on field-cooling the sample is still visible in Fig 4. This is perhaps assisted by the increasing sensitivity ($`B^2`$) of torque magnetometry with field. While the reason for the continued anisotropy at high magnetic fields is not obvious, the abrupt change in slope in fields as high as $``$ 30.1 T clearly indicates the continuing presence of a thermodynamic phase boundary. The field of 30.1 T was chosen because, at this field, $`\mu `$ is exactly that for a half-integer filling factor, thereby eliminating any possible contribution to the torque from the dHvA signal . If our interpretation of the magnetic measurements of Christ et al. and Sasaki et al. is correct, this could implies that the loss of Landau diamagnetism on entering the low temperature, high magnetic field phase continues to be the dominant factor. Were part of the sample to become superconducting at high magnetic fields, as might be suggested from the behaviour of the resistivity, it is unclear how this would affect the results. Usually, the pure Abrikosov diamagnetism of the vortex lattice of a Q2D superconductor is largest when the field is applied perpendicular to the conducting layers . The torque resulting from the net gain in Abrikosov diamagnetism should therefore have the opposite sign of that coming from the loss in Landau diamagnetism. For our results to be consistent with the establishment of a high magnetic field-induced superconducting state, the Landau diamagnetism would have to be larger in magnitude than the Abrikosov diamagnetism at these fields. Whether this is the case depends on parameters which are presently unknown.

## V hysteresis in the magnetic torque

One of the most intriquing apsects of this material at high magnetic fields is its magnetic hysteresis. An example of the magnetic torque of sample 3 measured both on rising and falling magnetic fields at 27 mK is shown in Fig. 5. Apart from the fact that the temperature is lower in the current work, the hysteresis is very similar to that observed by Christ et al. , although this effect was largely disregarded as originating from a “complex magnetic groundstate,” with the original interpretation of the anisotropic susceptibility work of Sasaki et al. in terms of antiferromagnetism, essentially forming the foundation of their argument.

Consideration of the fact that only the anisotropy of the magnetic suscptibility gives rise to magnetic torque is the key to understanding the behaviour of the magnetic torque in this material. Hence, because the dHvA effect is the oscillatory component of the Landau diamagnetism, thereby only involving orbital effects within the conducting planes, it manifests itself as large oscillations of the magnetic torque that are rather straightforward to interpret. At low magnetic fields ($`B<B_\mathrm{k}`$), the oscillations exhibit the usual double-peaked structure that has been interpreted both as spin-splitting  and as the frequency-doubling effect that accompanies a pinned CDW or SDW phase . This is well known to be an effect observed at small angles ($`\theta `$ 20) in this material (as well as the $`M=`$ Tl and Rb salts), which occurs only within the low magnetic field DW phase. On increasing the field at these low temperatures (i.e. $`T=`$ 27 mK), the kink transition is observed to be extremely abrupt, with a small but reproducible spike feature signalling the transformation into the high magnetic field regime at $`B`$ 24.2 T in sample 3. Above this field, the dHvA oscillations develop the characteristic triangular form that has been observed within the high magnetic field phase .

The sudden change in the waveform that occurs at $``$ 24.2 T on the rising magnetic field and at $``$ 21.4 T on the falling magnetic field is indicative of a first order phase transition between two distinctly different regimes. Because the magnetic torque is entirely reversible within the low magnetic field DW phase, provided the field does not exceed $``$ 24.2 T, both the dHvA oscillations and the monotonic background torque remain unchanged when the direction of sweep of the magnetic field is reversed. Similarly, when the field sweep direction is reversed within the high magnetic field phase, the dHvA oscillations continue to have the same triangular form, provided that the field is not swept below $``$ 21.4 T. This type of behaviour indicates that the low and high magnetic field phases are immiscible, with the formation of domains being energetically unfavourable. Because it costs energy to mix the two phases, the low magnetic field DW phase is “supercooled” on increasing the magnetic field until the free energy difference between the two phases is no longer sustainable. To help understand this behaviour in $`\alpha `$-(BEDT-TTF)2KHg(SCN)4, it is instructive to consider a simplified model where the zero temperature free energies of the low and high magnetic field regimes are approximated as

$$F_0g_{1\mathrm{D}}\left[\frac{\mathrm{\Delta }_0^2}{2}h^2\right]$$
(3)
and

$$F_xg_{1\mathrm{D}}\frac{\mathrm{\Delta }_x^2}{2}$$
(4)
respectively, with $`g_{1\mathrm{D}}`$ being the density of Q1D states and $`h=g\mu _\mathrm{B}B/2`$ , and we assume, to first order, that the high magnetic field phase is not affected by magnetic field. Here, the subscript $`0`$ denotes the CDW phase that is stable at low magnetic fields, while $`x`$ denotes the spatially modulated or CDW-SDW hybrid phase . If we then assume the free energy of the domain term $`F_{\mathrm{mix}}ϵy(1y)`$ to be approximately parabolic (i.e. the lowest order even funtion of $`y`$), with $`y`$ being the fraction of the material in the low magnetic field phase, then no domain structure can ever be stable provided $`ϵ0`$. Since the total free energy is

$$F_{\mathrm{tot}}yF_0+(1y)F_x+ϵy(1y),$$
(5)
it then follows that, on increasing the magnetic field within the low magnetic field regime, within which $`y=1`$, the material cannot “snap” into the high magnetic field phase until $`F_{\mathrm{tot}}/B0`$. Conversely, on decreasing the magnetic field from within the high magnetic field regime, within which $`y=0`$, the material cannot snap into the low magnetic field phase until $`F_{\mathrm{tot}}/B0`$. The kink transitions for rising and falling magnetic field are therefore

$$h_\mathrm{k}=\sqrt{\frac{\mathrm{\Delta }_0^2\mathrm{\Delta }_x^2}{2}\pm \frac{ϵ}{g_{1\mathrm{D}}}},$$
(6)
with the $`\pm `$ sign corresponding to the direction of sweep of the magnetic field. The abruptness of the transitions should be accompanied by the release of latent heat, although this cannot be detected in the present isothermal experiment. If we assume the strength of the coupling to be similar in the two ordered regimes, then the ratio of order parameters $`\mathrm{\Delta }_x/\mathrm{\Delta }_0`$ is equal to the ratio of transition temperatures $`T_\mathrm{c}/T_\mathrm{p}3/8`$, and we obtain $`\mathrm{\Delta }_0`$ 4.1 meV and $`\mathrm{\Delta }_x`$ 1.5 meV on inserting $`B_\mathrm{c}B_\mathrm{k}`$ 23 T. Meanwhile, for $`1/g_{1\mathrm{D}}`$ 30 meV, we obtain $`ϵ`$ 24 $`\mu `$eV.

The hysteresis that starts at magnetic fields above $``$ 24.5 has quite a different form from that associated with the kink transition. Since (1) the dHvA oscillations have the same shape and size between rising and falling magnetic fields and (2) the hysteresis takes place continuously over a wide interval in field (from $``$ 24.5 T up to the highest available field of $``$ 33 T), this new type of hysteresis appears not to be a consequence of a second prolonged first order phase transition as a function of magnetic field, as was suggested by McKenzie . For a first order phase transition to take place over such an extended interval of field, either (1) the free energies of the two competing phases would have to be very similar but with slightly different dependences on $`B`$ or $`\mu `$, or (2) the mixing term $`ϵ`$ would have to be very large (and negative). Were this the case, the hysteresis would have to result from two competing CDW or CDW-SDW hybrid phases with different anisotropic magnetic susceptibilities, different nesting vectors and therefore different reconstructed Fermi surface topologies. Any difference in Fermi surface topology between rising and falling magnetic field should ultimately manifest itself on the dHvA oscillations. This is to be expected since the peak-to-peak amplitude of the dHvA waveform is directly proportional to the number of Q2D states, while the gradient on the falling side of the oscillation (i.e. at integer filling) is proportional to the density of background states . As it turns out, however, there is no detectable change in the peak-to-peak height nor in the falling slope of the dHvA oscillations between rising and falling fields in Fig. 5, at least for $`B`$ 26 T. In fact, from the fraction of the waveform $`\gamma =g_{2\mathrm{D}}/(g_{2\mathrm{D}}+g_{1\mathrm{D}})`$ over which the dHvA magnetisation increases with field , it is rather straight forward to infer that the density of background states is 68 $`\pm `$ 5 % that of the Q2D pocket, both on rising and falling fields, therefore having the same density of states, within experimental uncertaintly, as the Q1D Fermi surface sheets prior to nesting . This could imply that the Q1D states are not nested (or “ungapped”). Indeed, it has been a common conclusion of all quantum oscillation and angle-dependent magnetoresistance oscillation measurements, that the Fermi surface appears to be unreconstructed within the high magnetic field phase . Alternatively, the high magnetic field phase could be a superconductor, whereby (in the absence of a $`𝐐`$-vector) the gap automatically remains pinned to the oscillating chemical potential, enabling the Q1D Fermi surface sheet to continue affecting the waveform of the dHvA oscillations. One could also envisage a scenario whereby the $`𝐐`$ vectors of a CDW or CDW/SDW phase adjust themselves continuously in order to maintain the gap symmetrical with respect to the oscillating chemical potential; this situation would therefore contrast from that in the low magnetic field phase of $`\alpha `$-(BEDT-TTF)2KHg(SCN)4  and that in NbSe3 . In either case, the size and shape of the quantum oscillations indicates that there is no discernable difference in the density of states between rising and falling magnetic fields, therefore making it somewhat unlikely that the electronic ground states could be different between rising and falling magnetic fields. The fact that no hysteresis is observed in the magnetotransport further implies that it cannot be the electronic structure which is itself hysteretic.

The only other possibility is that the hysteresis originates from some form of dynamic magnetism involving the ferromagnetic alignment of spins or circulating currents. Whenever magnetic hysteresis occurs, whether it originates from ferromagnetism, metamagnetism, vortex pinning in a type II superconductor or even induced currents in a quantum Hall system, the tendency is always to retain magnetic flux within the sample. Using the field-cooled magnetisation as a point of reference (or the magnetisation averaged between up and down sweeps), the magnetisation is always slightly more diamagnetic on the rising field and slightly more paramagnetic on the falling field. Only the anisotropy of the magnetic susceptibility gives rise to magnetic torque: the sign of the hysteresis observed here thus implies that it involves magnetic moments that are predominantly oriented perperdicular to the conducting planes, or equivalently, currents flowing within the conducting planes. Clearly, $`\alpha `$-(BEDT-TTF)2KHg(SCN)4 contains no magnetic ions with partially filled $`d`$\- or $`f`$-electron shells, so we can eliminate this rather trivial source of ferromagnetism or metamagnetism. The hysteresis in Fig. 5 is therefore more easy to explain in terms of induced currents. The fact that the hysteresis is observed only in the magnetic torque and not in the magnetotransport could be considered consistent with this hypothesis. These cannot be quantum Hall currents, however, since the hysteresis is observed at all fields $`B`$ 24.5 T, both at integer and half-integer filling factors.

## VI irreversible processes

The behaviour that is probably most supportive of the notion of induced persistent currents, is the irreversibility of the magnetic torque on stopping or reversing the sweep direction of the magnetic field. In most metals, including those with DW phases, the magnetisation is usually reversible, as is the case, for example, within the low magnetic field DW phase (i.e. see Fig. 5). On sweeping the field back and forth within the high magnetic field phase, on the other hand, even over a very small interval in field, the magnetic torque of $`\alpha `$-(BEDT-TTF)2KHg(SCN)4 is observed to be entirely irreversible, giving rise to a hysteresis loop of the form shown in Fig. 6(a). In fact, the nature of the hysteresis has several of the characteristic features of a critical state model, such as the Bean model  that is used to describe magnetic hysteresis caused by vortex pinning in type II superconductors. This becomes particularly evident in Fig. 6(b) when the monotonic background and dHvA contributions are subtracted; the background is obtained by averaging data taken on full up and down sweeps of the magnetic field. In Fig. 6(b) we can see that, on reversing the sweep direction of the magnetic field, the magnetic hysteresis increases gradually until a critical value is reached. In type II superconductors, this would be proportional to the critical current density $`j_\mathrm{c}`$. In Fig. 6(b), the induced magnetisation corresponding to the critical state is that obtained by performing full up and down sweeps \[also shown in Fig. 6(b)\]. By “critical state” we imply that the sample has a tendancy to trap flux exactly like a type II superconductor, with the local currents determining the local magnetic field gradient. While the unit of trapped flux in type II superconductors is the vortex which contains a net flux equivalent to that of the flux quantum, it is unclear what the equivalent unit of trapped flux is in the high magnetic field phase of $`\alpha `$-(BEDT-TTF)2KHg(SCN)4.

It nevertheless follows from the critical state model as normally applied to type II superconductors, that as soon as the direction of sweep of the magnetic field is reversed, the initial slope of the magnetisation with respect to magnetic field is given by

$$\chi =\frac{M}{B}=\frac{f[1\eta ]}{\mu _0},$$
(7)
where $`1/\mu _0`$ corresponds to perfect diamagnetism, $`\eta `$ is the demagnetisation factor and $`f`$ is the volume fraction of the sample in which pinning occurs. This would then be, perhaps, as close to performing a Meissner-type experiment as one could get in extreme magnetic fields. If we consider that the anisotropic hysteresis originates from currents flowing only within the conducting planes and that $`\eta `$ 0.5 (as for a cylinder), upon taking the initial slope of the magnetic torque in Fig. 6(b), we obtain $`f`$ 1 % for sample 3 (of volume $``$ 0.8 mm3). Thus, were the hysteresis in Fig. 6 really caused by vortex pinning, this would imply either that only 1 % of the vortices are pinned, or that only 1 % of the sample is superconducting. This conclusion, that only part if the sample is pinned or superconducting, could be considered consistent with the failure of the sample to transform into a bulk zero resistance state. The value of $`f`$ 1 % would, of course, be underestimated were any of the supercurrent to flow perpendicular to the conducting planes via Josephson tunneling. The value of 1 % should therefore be considered as a lower limit. That it is necessary to sweep the magnetic field by $``$ 0.4 T before the trapped flux in the sample is completely reversed, could either imply that the local field difference in parts of the sample reaches strengths of order 0.4 T, or that shape effects, which give rise to demagnetisation, are important.

Another property of the critical state model as applied to type II superconductors, is the that they often give rise to relaxation phenomena, through thermally assisted flux flow or quantum flux creep . One way to test whether such effects occur in $`\alpha `$-(BEDT-TTF)2KHg(SCN)4 is to monitor the dependence of the width of the hysteresis loop on sweep rate. Non-linear current-voltage characteristics have already been detected in pulsed magnetic fields experiments, although it is quite likely, should $`\alpha `$-(BEDT-TTF)2KHg(SCN)4 (or the $`M=`$ Tl salt) be a superconductor, that these experiments were performed close to the flux flow regime . On sweeping the field at rates between 8 mTs-1 and 50 mTs-1 in Fig. 6, the hysteresis increases by no more than 10 %, implying that the relaxation rate is rather low. A more notable degree of relaxation of the magnetic torque is observed by stopping the sweep of the magnetic field abruptly and then observing changes over several minutes. At integer filling factors \[Fig. 7(a)\], the magnetic torque decays much more rapidly than at half-integer filling factors \[Fig. 7(b)\] and is considerably more noisy. On fitting the Anderson-Kim flux-creep model  (making the substitution of $`\tau `$ for $`M`$)

$$\tau _{\mathrm{in}}\tau _{\mathrm{crit}}\left[1\frac{\alpha _1}{2\alpha _\mathrm{c}}\mathrm{ln}t\right],$$
(8)
which is thought to be applicable at low temperatures in most type II supeconductors, the better fit is obtained at half-integer filling factors, with the characteristic parameter $`\alpha _1/\alpha _\mathrm{c}`$ for logarithmic decay being of order 1/620. This is certainly within an order of magnitude of the theoretical prediction for most type II superconductors .

## VII field and temperature-dependence of the hysteresis

Perhaps, a further analogy can be made with type II superconductors on consideration of the temperature dependence of the width of the hysteresis loop. In many type II superconductors, this is found to be strongly dependent on temperature, caused, for example, by the dependence of the vortex pinning potential $`U_\mathrm{c}(T)`$ on $`T`$ . Such a behaviour is also observable for the critical state of $`\alpha `$-(BEDT-TTF)2KHg(SCN)4, as shown in Fig. 8(a), albeit the hysteresis contains an oscillatory component. While no hysteresis could be detected at $`T`$ 4 K, it is already quite pronounced by 1.16 K, particularly at half-integer filling factors. Half-integer filling factors, realised whenever $`F/B+1/2`$ assumes an integer value, are depicted in Fig. 8 by solid vertical lines. Note that the oscillatory component of the width of the magnetic hysteresis loop oscillates in quadrature with the dHvA oscillations \[extracted in Fig. 8(b)\], but in phase with the density of states, with the hysteresis and the density of states at $`\mu `$ both exhibiting maxima at half-integer filling factors. Clearly, as discussed in the preceding section, it cannot be the dHvA oscillations themselves that are hysteretic because (1) the dHvA and oscillatory component of the hysteresis effects are at quadrature with respect to each other and (2) the actual dHvA amplitude is only relatively weakly dependent on temperature over the same range 27 mK $`<T<`$ 1.15 K, shown in Fig. 8(b). This together with the absence of any hysteresis in the SdH waveform, again, is consistent with our hypothesis in terms of induced currents.

As $`T0`$, the hysteresis becomes double-peaked, with the second sharper and more strongly temperature-dependent peak emerging at integer filling factors. Sharp peaks in the hysteresis were observed at integer filling factors in pulsed magnetic fields , and it was this “coincidence” that led to the data being interpreted in terms of the QHE. When plotted as a function of temperature in Fig. 9, the hysteresis is strongly dependent on temperature both at integer and half-integer filling factors, with the form of the $`T`$-dependence not being too different from that typically observed in high temperature superconductors .

## VIII phase diagram

It is instructive to construct a tentative phase diagram for $`\alpha `$-(BEDT-TTF)2KHg(SCN)4, from the data accumulated, at least for fields applied within $`\theta `$ 10 of the $`𝐛`$ axis of the crystals. Solid symbols have been chosen for the thermodynamic data, in Fig. 10, which is often the most reliable indicator for a thermodynamic phase transition. Solid squares denote the change in slope of the field-cooled magnetic torque extracted from Fig. 4. At high magnetic fields, this agrees well with the resistive transition at low temperatures (indicated by open circles), as, for example, can be seen in both samples 1 and 2 in Fig. 2. Note that this transition is not sample-dependent; only the extent to which the resistivity decreases at low temperatures is sample-dependent, varying somewhat between samples 1 and 2.

At low magnetic fields, the field-cooled magnetic torque data agrees rather well with the transition into the low magnetic field DW phase obtained from specific heat measurements , represented here by solid triangles. In the vicinity of the kink transition, however, the second order phase boundary, on cooling, becomes more difficult to extract from resistivity data, owing to the competition between different types of resistivity behaviour from the low magnetic field DW and exotic high magnetic field phases. For this reason, these points, depicted here as open squares, have been taken from the non-hysteretic kink observed in the magnetoresistance at higher temperatures between 2.5 K and 5 K. The hysteretic kink transition at lower temperatures, on the other hand, has been extracted from the magnetic torque measurements made in the dilution refrigerator, denoted here by vertical crosses on the rising magnetic field and diagonal crosses on the falling magnetic field.

The emerging phase diagram can be described as follows: at temperatures above $``$ 8 K, $`\alpha `$-(BEDT-TTF)2KHg(SCN)4 behaves like any other organic metal, transforming into a DW phase at low temperatures. Given the lack of evidence for static antiferromagnetically configured spins and the close similarity of our experimentally determined phase diagram to that predicted both by McKenzie  and Zanchi et al., this is more likely to be the CDW0 phase. Certainly, the first order transition into a new low temperature, high magnetic phase is mostly consistent with the theoretical phase diagram proposed by McKenzie . What is particularly interesting in the current work is that while the high magnetic field phase (CDWx) is expected, largely on theoretical grounds, to be a modulated CDW phase or mixed CDW-SDW hybrid, its transport and magnetic properties are quite unlike those observed in all known DW systems . Rather, its physical properties more closely mimic those of a type II superconducting phase in strong magnetic fields.

## IX conclusion

We have shown that the magnetic behaviour of $`\alpha `$-(BEDT-TTF)2KHg(SCN)4 at high magnetic fields has many of the characteristic features of a critical state model, strongly resembling a type II superconductor in a magnetic field . This critical state-like behaviour, together with the pronounced drop in resistivity at low temperatures, is clearly unlike that expected for a DW systems , nor can it be connected with the QHE.

While certain physical aspects of the behaviour of $`\alpha `$-(BEDT-TTF)2KHg(SCN)4 at high magnetic fields are similar to those of a type II superconductor, there are, however, many reasons for questioning any interpretation in terms of a field-induced superconducting state. While the presence of a CDW at low magnetic fields indicates that electron-phonon interactions are important in this material, a CDW ground state is invariably always more stable than an electron-phonon mediated superconductor in a strong magnetic field. Both CDW’s and singlet-paired superconductors are Pauli limited at a critical field $`B_\mathrm{c}`$. However, superconductors have the additional disadvantage of being suppressed by orbital effects at much lower magnetic fields than that at which the Pauli critical field is expected to occur . Thus, if the ground state of $`\alpha `$-(BEDT-TTF)2KHg(SCN)4 results from the competition between superconductivity and CDW’s, as one might gather from the fact that the $`M=`$ NH4 salt is superconducting and that $`\alpha `$-(BEDT-TTF)2KHg(SCN)4 becomes superconducting under uniaxial stress, the fact that the CDW has already “won” at $`B=`$ 0, makes the emergence of superconductivity in an applied field seem somewhat remote. Furthermore, while several models have been proposed for “reentrant” superconductivity in very strong magnetic fields , a common prerequisite of these models is that the material is already a stable superconductor at $`B=`$ 0.

From the CDW perspective, on the other hand, neither the modulated CDW phase nor the CDW-SDW hybrid phase have yet been the subject of detailed theoretical or experimental studies . We cannot exclude the possibility that the CDWx phase itself exists in some form of critical state, which then only mimics the inductive behaviour of a superconductor in strong magnetic fields, and has yet to be identified. With such an explanation, however, there would be some difficulty in explaining the pronounced drop in the resistivity at low temperatures in this material. One could also imagine that, within the CDWx regime, the Frölich CDW mode is no longer efficiently coupled to impurities  or to the lattice, leading to dissipationless currents as $`T`$ 0.

Clearly, there exists a challenge in this material, and similar materials within the same series (i.e. the $`M=`$ Tl and Rb salts), to find alternative ways of probing the possible presence of non-dissipative currents (or supercurrents) above $`B`$ 24.2 T. The present data indicates only that the high magnetic field phase of $`\alpha `$-(BEDT-TTF)2KHg(SCN)4 exhibits properties that are qualitatively similar to those of type II superconductors in strong magnetic fields. On the other hand, without the existence of a type I superconducting phase, it would be very difficult to experimentally prove the existence of type II superconductivity at such strong magnetic fields. What is clear from these experiments is that the high magnetic field phase of $`\alpha `$-(BEDT-TTF)2KHg(SCN)4 is very different from any CDW (or SDW) observed in any other material, nor can any of the observations in static magnetic fields be considered consistent with the quantum Hall effect.

## X acknowledgements

The work is supported by the Department of Energy, the National Science Foundation (NSF) and the State of Florida. One of us (JSB), acknowledges the provision of an NSF grant (DMR-99-71474), and NH would like to thank Albert Migliori for useful discussions. One of us (LB) would also like to acknowledge the provision of a FSU visiting scientist scolarship. We would futher like to acknowledge Eric Palm and Timothy Muphy for maintaining the dilution refrigerator during the lowest temperature experiments.


format: text

text:
# Slabs of stabilized jellium: Quantum-size and self-compression effects

## I Introduction

Thin films or slabs are systems made out of a few layers of atoms: they are finite in one direction and infinite in the other two perpendicular directions. Various theoretical models are available to calculate the electronic structure of slabs. One of the important features coming out of these calculations is the so-called quantum-size effect (QSE), i.e., the influence of the finite size on various physical properties of the slab. These effects, which can be experimentally recognized, decrease as the size of the slab increases. In fact, surface energies and work functions of the semi-infinite system are often derived from thin-slab calculations, which are simply extrapolated to this limit.

The simplest model to predict the electronic structure of simple $`sp`$-bonded metals is the jellium model, where the ions are replaced by a positive neutralizing background. Within this model, the QSE of thin films was examined by Schulte. He found an oscillatory behavior of the work function, as a function of the thickness of the slab. The same oscillatory behavior is found for the surface energy, defined as the energy required, per unit area of new surface formed, to split the solid in two along a plane.

The jellium model has been referred to as giving insight into the realistic QSE appearing in real systems. Notwithstanding important differences, an oscillatory pattern also appears in atomistic first-principles slab calculations of both the work function and the surface energy. Although the presence of the lattice may obscure the periodicity and the amplitude of the QSE, extrema were found at positions which agree with the jellium results. However, in the case of first-principles calculations difficulties arise (due to the cumbersome numerics) when one is to extract well-converged surface properties from thin films made of typically 2 to 15 layers. Hence, the clean jellium QSE, with no uncertainties in the extrapolated results, remains as a guide for more realistic investigations.

In this paper, we consider slabs in the framework of a simple modification of the jellium model which yields energy stability against changes in the background density. This so-called stabilized jellium  or structureless pseudopotential model yields realistic results, especially in the case of metals with high valence-electron density. For instance, the stabilized jellium model predicts positive surface energies that increase rapidly at high electron densities, as shown by experiment, while the jellium model predicts surface energies that are strongly negative at these densities. The stabilized jellium model, first introduced by Perdew, Tran, and Smith and similar to the ideal-metal concept developed by Shore and Rose, has been applied to the study of surfaces and clusters. In a way, the stabilized jellium is in between the jellium model and more sophisticated atomistic approaches: although it is still a continuous model (one may choose slabs with arbitrary thickness) with an analytical expression for the bulk energy, its physical predictions are in reasonable agreement with experiment. Besides including electrostatic corrections to the jellium model, the stabilized jellium model contains an averaged pseudopotential correction.

We calculate the self-consistent energetics (surface energy and work function) of slabs of stabilized jellium, with use of the local-density approximation (LDA) of density-functional theory (DFT). We take two metals: Al ($`r_s=2.07`$, $`Z=3`$) and Li ($`r_s=3.24`$, $`Z=1`$), investigate the QSE, and compare our self-consistent slab calculations with those obtained for a semi-infinite stabilized jellium. We also test an extrapolation rule, which has already been used to describe non-local surface energies of the bounded electron gas.

Although the stabilized jellium model can be tailored to give face-dependent results, it cannot describe the inhomogeneous relaxation predicted by first-principles calculations where the distances between atomic planes of the same family are optimized. However, an interesting effect displayed by the stabilized jellium model, which cannot be accounted for by the jellium model, is the so-called self-compression  (or self-expansion, in the case of charged systems ) of clusters. This effect, which can be classically viewed as the compression of a finite system due to the surface tension, is most prominent for small systems and almost negligible in the case of large clusters.

We investigate here the self-compression of thin films. We fix the size of the system along the direction perpendicular to the surface, and search for the background density which minimizes the total energy per valence electron of the slab. The equilibrium density is found to increase as the thickness of the slab decreases, and to converge to the bulk electron density in the infinite-thickness limit. Furthermore, the equilibrium electron-density parameter $`r_s^{}`$ is found to oscillate with the slab thickness, as a manifestation of the QSE, but the general trend is found to be well described within the liquid-drop model (LDM) based only on the knowledge of the bulk energy per unit volume and the surface energy. We discuss the relationship between this self-compression effect and the relaxation of metal slabs predicted by atomistic first-principles calculations.

In section II we present briefly the stabilized jellium model for slabs. In section III we discuss the results we obtained within this model. The main conclusions are drawn in section IV, where further comments on the relationship between the stabilized jellium and more elaborated models are made. Equations are written in atomic units throughout, i.e., $`e^2=\mathrm{}=m_e=1`$.

## II Slabs of stabilized jellium

The stabilized jellium model takes into account the lattice ions, but keeps the essential simplicity of the jellium model. The total energy is obtained as a functional of the electron density $`n(𝐫)`$, in the following way:

$`E_{SJ}[n,n_+]=E_J[n,n_+]+\left(e_M+\overline{w}_R\right){\displaystyle d^3rn_+(𝐫)}`$ (1)
$`+\delta v_{WS}{\displaystyle d^3r\frac{n_+(𝐫)}{\overline{n}}\left[n(𝐫)n_+(𝐫)\right]},`$ (2)
where

$$n_+=\overline{n}\mathrm{\Theta }(𝐫)$$
(3)
represents a positive neutralizing background density, $`\mathrm{\Theta }(𝐫)`$ being a function which equals 1 inside a given surface and 0 outside, and

$$\overline{n}=\frac{3}{4\pi r_{s}^{}{}_{}{}^{3}}$$
(4)
is the average valence-electron density. $`E_J`$ is the regular-jellium total energy, $`e_M`$ is the Madelung energy arising from the Coulomb interaction between a uniform negative background inside the spherical Wigner-Seitz cell and a point ion at its center,

$$e_M=\frac{9Z^{2/3}}{10r_s},$$
(5)
$`\overline{w}_R`$ is the average value of the repulsive non-Coulomb part of the Ashcroft empty-core pseudopotential,

$$\overline{w}_R=2\pi \overline{n}r_c^2,$$
(6)
and $`\delta v_{WS}`$ represents the difference between the local pseudopotential and the jellium potential, averaged over the Wigner-Seitz cell,

$$\delta v_{WS}=\frac{3r_c^2}{2r_s^3}\frac{3Z^{2/3}}{10r_s}.$$
(7)
The core radius $`r_c`$ of the Ashcroft empty-core pseudopotential is chosen to stabilize the metal for given values of the electron-density parameter $`r_s`$ and the chemical valence $`Z`$.

The two terms added to the regular-jellium energy $`E_J`$ are a volume term and a surface term. They simply account for the subtraction of the spurious self-interaction of the positive jellium background and the inclusion of a constant structureless potential inside the metal. This procedure may be understood as a first-order perturbation to a jellium system, but with the perturbation treated in an averaged manner.

The density functional of Eq. (1) represents the total energy of an arbitrary inhomogeneous system. In the case of an infinite uniform system, the equilibrium density is obtained from the bulk stability condition

$$\frac{dϵ_{SJ}^{bulk}}{dr_s}=0,$$
(8)
where

$$ϵ_{SJ}^{bulk}=e_J^{bulk}+e_M+\overline{w}_R$$
(9)
represents the average bulk energy per valence electron, $`e_J`$ being the regular-jellium contribution. Within this model any individual metal minimizes the energy at a given equilibrium density, while the jellium energy presents a single minimum at $`r_s4.2`$ close to the electron-density parameter of sodium.

We consider slabs of stabilized jellium. Slabs are translationally invariant in the plane of the surface, which is assumed to be perpendicular to the $`z`$ axis. Hence, the single-particle wave functions can be separated into a plane wave along the surface and a component $`\varphi (z)`$ describing motion normal to the surface with energy $`ϵ`$. This component is obtained by solving self-consistently the Kohn-Sham equation

$$\left[\frac{1}{2}\frac{d^2}{dz^2}+V_H(z)+V_{xc}(z)+V_{ps}(z)\right]\varphi (z)=ϵ\varphi (z),$$
(10)
where $`V_H(z)`$ represents the Hartree electrostatic potential, $`V_{xc}(z)`$ is the exchange-correlation potential, and $`V_{ps}(z)`$ accounts for the pseudopotential,

$$V_{ps}(z)=\delta v_{WS}\mathrm{\Theta }(z).$$
(11)
$`V_{xc}(z)`$ is obtained in the LDA, using the electron-gas correlation energy of Ceperley and Alder, as parameterized by Perdew and Wang. Essentially the same results are obtained from the parameterizations of Vosko, Wilk, and Nusair and of Perdew and Zunger. We have not chosen to use extensions such as the generalized gradient approximation (GGA), since the LDA has been shown to give surprisingly good results in describing the properties of jellium planar surfaces.

Outside the positive background the electron-density profile $`n(z)`$ decays rapidly from its bulk value $`\overline{n}`$. The electronic system can therefore be taken to be finite in the $`z`$ direction by assuming that $`n(z)`$ actually vanishes at a given distance $`z_0`$ from the surface. Hence, we introduce infinite potential walls at a distance $`z_0`$ from each surface, and follow Ref. to expand the wave functions $`\varphi (z)`$ in a Fourier sine series. The distance $`z_0`$ (typically 2 or 3 Fermi wavelengths) and the number of sines kept in the expansion of the wave functions $`\varphi (z)`$ have been chosen to be sufficiently large for our calculations to be insensitive to the precise values employed. These calculations have been compared with others which we have carried out for a semi-infinite electron system by using the Monnier-Perdew code for the numerical integration of Eq. (10).

For a given thickness $`L`$ of the slab, we obtain the surface energy from the difference between the total energy of Eq. (1) and the corresponding result for a homogeneous electron gas of density $`n_+`$, i.e.,

$$\sigma (L)=\frac{1}{2A}\left[E_{SJ}(L)\overline{n}LAϵ_{SJ}^{bulk}\right],$$
(12)
where A is the normalization area. The work function is obtained as the difference between the computed values for the vacuum and Fermi levels of our electron system.

## III Results and discussion

First of all, we compare jellium and stabilized-jellium electron densities $`n(z)`$ and effective potentials,

$$V_{eff}(z)=V_H(z)+V_{xc}(z)+V_{ps}(z).$$
(13)
Jellium and stabilized-jellium valence-electron densities and effective potentials for an Al slab of $`L=2\lambda _F`$ \[$`\lambda _F=\left(32\pi ^2/9\right)^{1/3}r_s`$ is the Fermi wavelength\] are shown in Fig. 1, together with the positive background density $`n_+`$. We note that the stabilized-jellium electron density is steeper at the two surfaces, so that the electronic spill-out is slightly smaller within this model. This is due to the fact that electrons feel a deeper effective potential. Both jellium and stabilized-jellium electron densities exhibit quantum oscillations inside the metal, the so-called Friedel oscillations, and an exponential decay outside.

Figs. 2 and 3 show our calculated stabilized-jellium surface energies for slabs of Al and Li, respectively, as obtained from Eq. (12) versus the thickness $`L`$ of the slab. Both curves show damped oscillations with minima occurring at the slab width $`Ln\lambda _F/2`$ ($`n=1,2,\mathrm{}`$). The same QSE, which reflects the quantization of the electronic motion along one direction, is known to occur within the jellium model.

Both the average bulk energy per valence electron $`ϵ_{SJ}^{bulk}`$ and the surface energy of the semi-infinite stabilized jellium

$$\sigma =\underset{L\mathrm{}}{lim}\sigma (L)$$
(14)
may be obtained from a linear fit of the following equation:

$$\frac{E_{SJ}(L)}{A}=2\sigma +\overline{n}Lϵ_{SJ}^{bulk},$$
(15)
where $`E_{SJ}(L)`$ represents the total energy of Eq. (1). Following this procedure, we reproduce the bulk energy of Eq. (9) and predict surface energies of $`925\mathrm{erg}/\mathrm{cm}^2`$ and $`311\mathrm{erg}/\mathrm{cm}^2`$ for Al and Li, respectively. These surface energies, represented in Figs. 2 and 3 by horizontal solid lines, agree with those reported in Ref. for semi-infinite media.

An alternative procedure to extrapolate the surface energy $`\sigma `$ of the semi-infinite medium from our calculated thin-film surface energies $`\sigma (L)`$ is to use the relation

$$\sigma =\frac{\sigma (L_n\lambda _F/4)+\sigma (L_n)+\sigma (L_n+\lambda _F/4)}{3},$$
(16)
where $`L_n`$ represents the threshold width for which the $`n`$th subband for the $`z`$ motion is first occupied. Analytical insight for this procedure is encountered within the infinite-barrier model (IBM), where the effective potential $`V_{eff}(z)`$ is replaced by an infinite square well and the one-particle wave functions $`\varphi (z)`$ are simply sines. Based on this procedure, the numerical error introduced in $`\sigma `$ by our slab calculations is found to be within $`0.1\%`$. The advantage of this algorithm is that we simply need three points to obtain the asymptotic limit, while the linear fitting may yield erroneous results if one only takes a few thin films.

Slabs with $`L<0.5\lambda _F`$ are interesting in their own, since they can be constructed in the laboratory, e.g., by joining two different semiconductors. Nevertheless, we do not give results for these ultra-thin slabs, since they fall within the two-dimensional limit where the three-dimensional LDA and GGA formulae for exchange and correlation are known to fail.

For comparison, first-principles thin-film calculations of the surface energy of the most dense faces of Al and Li \[(111) for fcc Al and (0001) for hcp Li\] are represented in Figs. 2 and 3 by solid circles and triangles, with the slab width of a $`\nu `$-layer unrelaxed crystalline film taken to be $`\nu `$ times the interplanar distance. For Al there is reasonable agreement between our stabilized-jellium results and atomistic first-principles calculations, the amplitude of the stabilized-jellium oscillations being comparable to that exhibited by first-principles calculations. For Li, however, there is a serious discrepancy between stabilized-jellium and first-principles calculations. Since lithium has been found to behave to some extent like a covalent solid rather than a free-electron gas, it is not expected to be well described by a jellium-like model.

A face-dependent approach extension of the stabilized-jellium model consists in obtaining the self-consistent electron density by adding to the constant potential $`\delta v_{WS}`$ a structure-dependent corrugation factor. This procedure yields an increased surface energy \[horizontal dashed-dotted lines of Figs. 2 and 3\], which in the case of Al is found to be close to the experimental result.

Figs. 4 and 5 exhibit our calculated stabilized-jellium work functions for slabs of Al and Li, respectively, as a function of the thickness $`L`$ of the slab, together with first-principles thin-film calculations. As in the case of the surface energy, a procedure similar to that of Eq. (16) yields a work function \[represented by horizontal solid lines\] that agrees within less than $`0.1\%`$ with the result we also obtain after solving Eq. (10) for the semi-infinite medium, a precision that is difficult to achieve by a fitting procedure. For $`L0.5\lambda _F`$, the QSE yields oscillations with relative amplitudes of $`20\%`$ and $`10\%`$ for Al and Li, respectively. For Al both the amplitude and the oscillation pattern are comparable to those exhibited by atomistic calculations. In the case of a 3-layer film of Al(111), the slab width is $`L4(\lambda _F/2)`$. Hence, the stabilized-jellium model predicts a minimum for this film, which is in reasonable agreement with the deep minimum exhibited by atomistic calculations with $`\nu =3`$. In the case of Li(0001), the stabilized-jellium model predicts a minimum for a one-layer film \[$`L1(\lambda _F/2)`$\], also in agreement with the minimum exhibited by first-principles calculations with $`\nu =1`$. Finally, we note that adding a structure-dependent corrugation factor to the slabilized-jellium $`\delta v_{WS}`$ constant potential yields a smaller value of the work function \[horizontal dashed-dotted lines of Figs. 4 and 5\], which in the case of Al is in reasonable agreement with the experiment. For Li, both the stabilized-jellium model and first-principles calculations predict work functions that are well above the experimental result.

For given values of the equilibrium-density parameter $`r_s`$ and the valence $`Z`$, all these calculations have been carried out with the core radius $`r_c`$ \[characteristic of each metal\] that is obtained from the bulk stability condition expressed by Eq. (8). However, while at the equilibrium density $`\overline{n}`$ of Eq. (4) the infinite homogeneous system is stable, at this density a finite system is not stable against changes of the background density, i.e.,

$$\frac{d(E/N)}{dr_s}0,$$
(17)
where $`N`$ represents the particle number. Instead, there is a modified equilibrium-density parameter $`r_s^{}`$, which stabilizes the finite system. This modified parameter depends on the size $`L`$ of our system and is expected to approach $`r_s`$ as $`L\mathrm{}`$.

Fig. 6 shows the result of our full self-consistent Kohn-Sham calculations of the deviation $`r_s^{}r_s`$, as a function of the thickness $`L`$ of the slab. These calculations indicate that there is a self-compression effect, which is more pronounced when the two surfaces are separated by a multiple of $`\lambda _F/2`$.

The self-compression effect exhibited in Fig. 6 may be approximately predicted with use of the LDM, a simple model to evaluate the total energy of a finite system. In this model, the energy is the sum of a volume term (the bulk energy per unit volume, $`\overline{n}ϵ_{bulk}`$, times the volume) and a surface term (the surface energy, $`\sigma `$, times the transversal area):

$$E_{LDM}=\overline{n}ϵ_{bulk}V+\sigma A.$$
(18)
For fixed $`r_c`$, and evaluated at the bulk equilibrium-density parameter $`r_s`$,

$$\frac{d(E_{LDM}/N)}{dr_s}=\frac{A}{N}\frac{d\sigma }{dr_s}+\sigma \frac{d(A/N)}{dr_s}>0.$$
(19)
The first term is positive, as can be found from the data in Table I of Ref.. For a fixed slab width $`L`$, the second term is also positive, and the surface term self-compresses, therefore, stabilized-jellium slabs. The deviation of the electron density parameter $`r_s^{}`$ obtained from the LDM stability condition

$$\frac{d(E_{LDM}/N)}{dr_s}=0$$
(20)
with respect to the bulk equilibrium density parameter $`r_s`$ is also plotted in Fig. 6, showing that the LDM provides a nice average of our self-consistent Kohn-Sham calculations, as previously demonstrated in the case of clusters.

In Ref. , thin films of Be with 1-3 layers were examined and a jellium version of a crystalline calculation was considered. The electron density parameter $`r_s^{}`$ needed to define each slab was derived from the optimized (relaxed) structural parameters. The results reported in Ref. are in agreement with the compression effect we report here, with $`r_s^{}`$ increasing with the number of layers and approaching the bulk equilibrium-density parameter $`r_s`$ as $`L\mathrm{}`$. These results show deviations of the electron density parameter, $`r_s^{}r_s`$ of $`3.2\%`$, $`1.9\%`$ and $`0.9\%`$ for thin films with 1, 2 and 3 layers, respectively. This is in agreement with our stabilized-jellium calculations, which in the case of thin films with $`2`$ layers of Li and Al predict (see Fig. 6) differences between $`r_s^{}`$ and $`r_s`$ of $`1.6\%`$ and $`1.8\%`$, respectively. The self-compression of structural parameters in ultra-thin crystalline films has also been discussed in terms of the so-called coordination model which, however, seems to fail in some cases (see, e.g., Ref.).

Finally, we note that if for each value of $`L`$ the corresponding equilibrium-density parameter $`r_s^{}`$ is taken, instead of the bulk parameter $`r_s`$, modified surface energies and work functions are obtained which are quite similar to those displayed in Figs. 2-5. This is in contrast with the discussion of Ref..

## IV Conclusions

We have modeled thin films of two simple metals, aluminum and lithium, using the stabilized-jellium model, and have studied the convergence of some physical quantities (work function and surface energy) to the semi-infinite planar-surface results. We have found the same oscillatory behavior which is typical of the QSE in jellium. Although this behavior also shows up in atomistic first-principles thin-film calculations, the clean QSE of continuous background models is obscured in the more realistic calculations. A trend consisting in surface energy minima coinciding with work function maxima was reported for first-principles crystalline calculations.  However, within the stabilized-jellium model we have found minima and maxima of both quantities at the same positions \[as also reported in Ref.  from first-principles for Al(111)\]. On the other hand, we have found that both the absolute and the relative amplitude of stabilized-jellium QSE oscillations are larger for aluminum than for lithium, in agreement with first-principles evaluations. The disagreement between our stabilized-jellium results for lithium and the more realistic atomistic all-electron calculations cannot be attributed to some property of the pseudopotential, and simply shows that this metal does not display free-electron behaviour.

Stabilized-jellium slabs of aluminum and lithium have been found not to be stable at the bulk equilibrium density, the size-dependent equilibrium density being larger. This self-compression effect, which was already known to exist for clusters, has been found to become more important as the slab width decreases. Both LDM and full self-consistent DFT calculations have shown a larger self-compression for aluminum than for lithium, which is a consequence of the larger surface energy of the former material. The self-compression of thin simple-metal films is a general rule that is also exhibited by atomistic first-principles calculations, where the unitary cell of thin films is found to be slightly smaller than that of the bulk solid.

The stabilized jellium model is computationally as simple as the jellium model; however, for the two high-density metals we have considered, it is much more realistic. In particular, we have found it to be more realistic for aluminum than for lithium. The stabilized-jellium model is adequate to obtain general qualitative conclusions and an understanding of trends of simple metals but, obviously, is unable to provide precise quantitative conclusions on particular metals. These can only be extracted from the now standard first-principles, but computationally more demanding, calculations.

## V Acknowledgments

The authors gratefully acknowledge J. P. Perdew for fruitful discussions. This project has been supported by the University of the Basque Country, the Basque Hezkuntza, Unibertsitate eta Ikerketa Saila, the Spanish Ministerio de Educación y Cultura, and the Portuguese PRAXIS XXI Program (Project PRAXIS/2/2.1/FIS/473/94).


text:
# UTCCP-P-84UTHEP-421April 2000 Large-Flavor QCD on the Lattice11footnote 1Presented by K. Kanaya at TMU-Yale Symposium on Dynamics of Gauge Fields — an External Activity of APCTP, Dec. 13–15, 1999, Hachioji, Japan

## Abstract

We study the nature of the QCD vacuum at general number of flavors $`N_F`$ by numerical simulations on the lattice. Combining the results with those of the perturbation theory, we propose the following picture: 1) For $`N_F17`$, there exists only one IR fixed point at vanishing gauge coupling, i.e., the theory in the continuum limit is trivial. 2) For $`16N_F7`$, there is a non-trivial fixed point. Therefore, the theory is non-trivial with anomalous dimensions, however, without quark confinement. 3) For $`N_F6`$, theories satisfy both quark confinement and spontaneous chiral symmetry breaking in the continuum limit.

## 1. Introduction

The beta function of QCD is universal up to two-loop in the perturbation theory: $`\stackrel{~}{\beta }(g)=b_0g^3b_1g^5,`$ with $`b_0=\frac{1}{16\pi ^2}\left(11\frac{2}{3}N_F\right)`$ and $`b_1=\frac{1}{(16\pi ^2)^2}\left(102\frac{38}{3}N_F\right).`$ When the number of flavors ($`N_F`$) exceeds $`16\frac{1}{2}`$, $`b_0`$ becomes negative, so that the asymptotic freedom of QCD is lost. In this case, the origin $`g=0`$ becomes an IR fixed point, i.e. the low energy limit of QCD is a free theory. Quarks cannot be confined and the chiral symmetry cannot be broken spontaneously. $`N_F>16\frac{1}{2}`$ is a sufficient, but not a necessary condition for these. A natural question is, whether confinement and spontaneous breaking of the chiral symmetry are satisfied at all $`N_F`$ below $`16\frac{1}{2}`$.

It is well-known that the second coefficient $`b_1`$ changes its sign already at $`N_F8.05`$. Therefore, a non-trivial IR fixed point may appear at finite $`g`$. From the two-loop beta function, this happens for $`N_F=9`$–16. At least for $`N_F16`$, the IR fixed point shows up in a perturbative region. Therefore, it is plausible that the full beta function has a non-trivial IR fixed point for $`N^{}N_F16`$ with some $`N^{}16`$. When such an IR fixed point exists, the coupling constant cannot become arbitrarily large in the IR region — this will imply that quarks are not confined. In particular, for $`N_F16`$, QCD is perturbative also in the IR limit, i.e., we cannot expect confinement.

In this paper, we study the $`N_F`$-dependence of the QCD vacuum. A non-perturbative investigation is required. Lattice QCD is the only systematic method that enables us to study non-perturbative properties of QCD. We performed a series of lattice simulations for a wide range of $`N_F`$, from 2 to 300, to study the phase structure of QCD. When the phase diagram becomes clear, we are able to see the nature of the QCD vacuum in the continuum limit, and eventually answer the question about the condition on $`N_F`$ for confinement and spontaneous breakdown of the chiral symmetry.

In 1982, Banks and Zaks published a pioneering paper on the $`N_F`$-dependence of the QCD vacuum . Based on an early result obtained on the lattice , they assumed that, in the strong coupling limit, the theory is confining and the beta function is negative, for any value of $`N_F`$. Using the perturbative results for the beta function discussed above, they conjectured Fig. 1.(a) as the simplest $`N_F`$ dependence of the beta function, and studied the phase structure of QCD based on this beta function. The assumption of negative beta function in the strong coupling region leads to an additional non-trivial UV fixed point for $`N_FN^{}`$. Due to this UV fixed point, their conjecture for the phase structure is complicated. \[For different approaches see Refs. .\]

In the arguments of Banks and Zaks, the assumption of confinement and negative beta function in the strong coupling limit plays an essential role. Here, it should be noted that the lattice study on which this assumption was based, is a study of a pure gauge theory . With dynamical quarks, the vacuum structure can be much more complicated. Actually, in the realistic cases, there exist no proofs of confinement for general $`N_F`$ even in the strong coupling limit.222We can formally rewrite the theory in terms of mesons and baryons. We may argue that, when the resulting effective action for hadrons is a well-defined action of weakly interacting particles, it is plausible that quarks are confined — the hadrons are good variables to describe the dynamical contents of the theory. In the strong coupling limit, the effective action can be computed either in the large $`N_c`$ limit , using a meanfield approximation ($`1/d`$ expansion) , or using a heavy quark mass expansion . In these cases, the effective action for mesons and baryons seems to lead to a well-defined world with spontaneously broken chiral symmetry, suggesting quark confinement. For the realistic cases of $`N_c=3`$ with light quarks, however, the effective action is quite complicated; diagrams with quark loops, which are suppressed at large $`N_c`$, large $`d`$ or large $`m_q`$, become important. Therefore, in the realistic cases, it is difficult to analytically deduce the vacuum structure as a function of $`N_F`$ even in the strong coupling limit.  A non-perturbative investigation on the lattice is required.

In a previous study, we performed simulations of QCD in the strong coupling limit at various $`N_F`$ . We found that, when $`N_F7`$, light quarks are not confined and chiral symmetry is restored, even in the strong coupling limit. Here, we extend this study to weaker couplings and to a wider range of $`N_F`$ up to 300. Combining the lattice results and perturbative arguments, we conjecture Fig. 1.(b) for the non-perturbative beta function of QCD. This leads to the following $`N_F`$-dependence of the QCD vacuum.

* When $`N_F`$ is smaller than a critical value, the beta function is negative for all values of $`g`$. Quarks are confined and the chiral symmetry is spontaneously broken at zero temperature. We find that the critical number of flavors is 6. (Corresponding critical value of $`N_F`$ is 2 for $`N_c=2`$, to be compared with the two-loop value 5 .)
* On the other hand, when $`N_F`$ is equal or larger than 17, we conjecture that the beta function is positive for all $`g`$, in contrast to the conjecture by Banks and Zaks shown in Fig. 1.(a). The theory is trivial in this case.
* When $`N_F`$ is between 7 and 16, the beta function changes sign from negative to positive with increasing $`g`$. Therefore, the theory has a non-trivial IR fixed point, i.e., the theory is non-trivial with anomalous dimensions. We conjecture that quarks are not confined in this case.

Our results are in part reported in Refs.  and .

This report is organized as follows: In Sec. 2., we describe our lattice model and simulation parameters. Results for the phase structure in the strong coupling limit are given in Sec. 3.. The phase structure at finite coupling is discussed in Sec. 4.. Finally, in Sec. 5., we summarize our conclusions and conjectures.

## 2. Model and simulation parameters

We consider a 4-dimensional Euclidian lattice with the lattice spacing $`a`$ and finite lattice volume $`(N_sa)^3\times (N_ta)`$, where $`N_s`$ and $`N_t`$ are the number of lattice sites in spatial and temporal directions. The action consists of the gauge part and the quark part, $`S=S_{gauge}+S_{quark}`$. For the gauge part, we use the standard one plaquette action.

$$S_{gauge}=\frac{\beta }{6}\underset{x,\mu \nu }{}\mathrm{Tr}\left(U_{x,\mu }U_{x+\widehat{\mu },\nu }U_{x+\widehat{\mu }+\widehat{\nu },\mu }^{}U_{x+\widehat{\nu },\nu }^{}\right),$$
(1)
where $`U_{x,\mu }`$ is the gauge connection between the site $`x`$ and the neighboring site in the positive $`\mu `$ direction, and $`\beta =6/g^2`$.

For the quark part, there are several alternatives. Two conventional choices are the staggered (Kogut-Susskind) fermion action  and the Wilson fermion action . Staggered quark has been preferably used in previous studies because a part of the chiral symmetry is preserved in this formulation. On the other hand, the flavor structure is quite complicated: The action is local only when $`N_F`$ is a multiple of 4. Off the continuum limit, there exists a lattice artifact mixing these 4 degenerate flavors. For the more realistic cases $`N_F=2`$ and 3, we modify by hand the power of the fermionic determinant in the numerical path-integration. This necessarily makes the action non-local, that sometimes poses conceptually and technically difficult problems.

Another conventional choice for the quark action is the Wilson fermion action, which we adopt in the followings.

$$S_{quark}=\underset{f,x}{}\left[\overline{\psi }_x^f\psi _x^fK\left\{\overline{\psi }_x^f(1\gamma _\mu )U_{x,\mu }\psi _{x+\widehat{\mu }}^f+\overline{\psi }_{x+\widehat{\mu }}^f(1+\gamma _\mu )U_{x,\mu }^{}\psi _x^f\right\}\right],$$
(2)
The index $`f`$ ($`=1,\mathrm{},N_F`$) is for the flavors. We assume that all $`N_F`$ quarks are degenerate with the bare mass $`m_0`$, which is related to the hopping parameter $`K`$ by $`m_0=(1/2a)(1/K1/K_c)`$, where $`K_c=1/8`$ is the point where the bare mass vanishes. In this formulation, the chiral symmetry is explicitly broken at finite $`a`$. On the other hand, the flavor symmetry as well as the C, P and T symmetries are manifestly satisfied without mixing etc. also on the lattice.

In this study, we use the Wilson fermion formalism because this is the only formalism that preserves manifest flavor symmetry. This feature will be important in a study of $`N_F`$-dependence in QCD. Concerning the chiral symmetry, a perturbative study of Ward identities shows that the effects of the $`O(a)`$ chiral breaking terms can be absorbed by appropriate renormalizations, including an additive renormalization of quark mass . As a result, the location of the chiral limit $`K_c`$ shifts from the tree value 1/8 as a function of $`\beta `$.

In order to take the additive mass renormalization into account, we define the current quark mass $`m_q`$ in terms of an axial vector Ward identity .

$$2m_q\mathrm{\hspace{0.17em}0}|P|\pi (\stackrel{}{p}=0)=m_\pi \mathrm{\hspace{0.17em}0}|A_4|\pi (\stackrel{}{p}=0)$$
(3)
where $`P`$ is the pseudoscalar density, $`A_4`$ the fourth component of the local axial vector current, and $`m_\pi `$ the pion (screening) mass. A multiplicative renormalization factor for the axial current, which is not important in this study, is absorbed into the definition of the quark mass. We then define $`K_c`$ by the condition $`m_q=0`$. Numerical simulations for small $`N_F`$ shows that $`K_c(\beta )`$ is a smooth curve connecting 1/8 at $`\beta =\mathrm{}`$ and 1/4 at $`\beta =0`$.

In QCD, even when the confinement is realized at low temperatures, we expect deconfinement at high temperatures . Therefore, we have to study the temperature dependence of the results to distinguish between the finite temperature deconfinement transition/crossover due to high temperatures and a bulk transition at zero temperature due to the effects of many flavors. On a lattice with $`N_t`$ sites in the Euclidian time direction, the temperature is given by $`T=1/N_ta`$. When the beta function is negative, as in the case of small $`N_F`$, the lattice spacing $`a`$ is a decreasing function of $`\beta =6/g^2`$. In this case, $`T`$ increases with increasing $`\beta `$ for a fixed finite $`N_t`$. On the other hand, when the beta function is positive, $`T`$ decreases with increasing $`\beta `$. Therefore, in order to study the temperature dependence for general $`N_F`$, we have to compare the results at various $`N_t`$.

Numerical study shows that the value of $`m_q`$ defined through the axial vector Ward identity does not depend on whether the system is in the high or the low temperature phase when $`\beta `$ is large; $`\beta \begin{array}{c}>\hfill \\ \hfill \end{array}5.5`$ for $`N_F=2`$ . When we define the pion decay constant $`f_\pi `$ by $`\mathrm{\hspace{0.17em}0}|A_4|\pi (\stackrel{}{p}=0)=m_\pi f_\pi ,`$ either $`m_\pi =0`$ or $`f_\pi =0`$ has to be satisfied in the chiral limit. The chiral symmetry is restored when $`f_\pi =0`$ in the chiral limit, while $`m_\pi =0`$ if the chiral symmetry is spontaneously broken .

We perform simulations on lattices $`8^2\times 10\times N_t`$ ($`N_t=4`$, 6, and 8), $`16^2\times 24\times N_t`$ ($`N_t=16`$) and $`18^2\times 24\times N_t`$ ($`N_t=18`$). We vary $`N_F`$ from 2 to 300. For each $`N_F`$, we study the phase structure in the coupling parameter space $`(\beta ,K)`$. It should be noted that, in QCD with dynamical quarks, there are no order parameters for quark confinement. We discuss about confinement by comparing the screening pion mass and the lowest Matsubara frequency, and, simultaneously, consulting the values of plaquette and the Polyakov loop. In the followings, we call the pion screening mass simply the pion mass, and similarly for the quark mass. Further details about our simulation are described in Refs.

## 3. Strong coupling limit $`\beta =0`$

In a previous paper , we have shown that, when $`N_F7`$, light quarks are deconfined and chiral symmetry is restored at zero temperature even in the strong coupling limit. For $`N_F6`$, we have only one confining phase from the heavy quark limit $`K=0`$ up to the chiral limit $`K_c=0.25`$. On the other hand, for $`N_F7`$, we find a strong first order transition at $`K=K_d`$:

* When quarks are heavy ($`K<K_d`$), both plaquette and the Polyakov loop are small, and $`m_\pi `$ satisfies the PCAC relation $`m_\pi ^2m_q`$. Therefore, quarks are confined and the chiral symmetry is spontaneously broken in this phase. We find that $`m_q`$ in the confining phase is non-zero at the transition point $`K_d`$. While an extrapolation in $`1/K`$ shows that $`m_\pi ^2`$ decreases towards zero at $`K=0.25`$, the chiral limit does not belong to this phase ($`K_d<K_c`$).
* On the other hand, when quarks are light ($`K>K_d`$), plaquette and the Polyakov loop are large. In this phase, $`m_\pi `$ remains large in the chiral limit and is almost equal to twice the lowest Matsubara frequency $`\pi /N_t`$. This implies that the pion state is a free two-quark state and, therefore, quarks are not confined in this phase. The pion mass is nearly equal to the scalar meson mass, and the rho meson mass to the axial vector meson mass. The chiral symmetry is also manifest within corrections due to finite lattice spacing.

See  for details.

We extend the study to larger $`N_F`$. We show the results of plaquette at $`N_t=4`$ for $`N_F=7`$–300 in Fig. 2.(a). Clear first order transition can be seen at $`K`$ below 0.25 ($`1/K>4`$). We then study the $`N_t`$ dependence of the results, as shown in Fig. 2.(b) for the case of $`N_F=240`$. We find that, although the transition point shows a slight shift to smaller $`1/K`$ when we decrease $`T=1/N_ta`$ from $`N_t=4`$ to 8, it stays at the same point for $`N_t8`$. The chiral limit, $`1/K=4`$, cannot be achieved in the confining phase even at $`T=0`$. Similar result for $`N_F=7`$ is given in Ref. . Therefore, we conclude that the transition is a bulk transition.

From these results, we obtain the phase diagram at $`\beta =0`$ shown in Fig. 3.. The first-order transition line separating the confining phase and the deconfining phase reaches $`1/K=4`$ (chiral limit for small $`N_F`$) between $`N_F=6`$ and 7. Note that the critical value of the inverse critical hopping parameter is greatly increased from $`1/K=4.08`$ to 7.75 when $`N_F`$ increases from 7 to 240. Recall that $`1/K=8`$ is the point for massless free quarks at $`\beta =\mathrm{}`$.

## 4. Phase structure at finite $`\beta `$

Let us now study the case $`\beta >0`$. We perform simulations at $`\beta =0.0`$–6.0 varying $`N_F`$ from 2 to 300. On a lattice with a fixed finite $`N_t`$, we have the finite temperature deconfining transition when the quark mass is sufficiently heavy. In particular, at $`K=0`$, we have the first order deconfinement transition of the pure gauge theory at finite $`\beta `$ irrespective to the value of $`N_F`$.

For $`N_F6`$, the finite temperature transition/crossover line crosses the $`K_c`$ line at finite $`\beta `$ . A schematic diagram of the phase structure for this case is shown in Fig. 4.. (In the followings, we call the transition/crossover line simply as transition line.) We have a monotonous RG flow towards the strong coupling limit also on the $`K_c`$ line. Therefore, the whole finite temperature transition line moves towards a larger $`\beta `$ as $`N_t`$ is increased. In the limit $`N_t=\mathrm{}`$ (i.e., at $`T=0`$), only the confining phase is left.

For $`N_F7`$, we have seen that there exist a bulk deconfining transition in the strong coupling limit. We have to clarify the relation between the finite temperature transition at small $`K`$ and the bulk transition at $`\beta =0`$. From our simulations, we obtain the phase diagrams sown in Fig. 5. for $`N_F=12`$ and 18, and in Fig. 6. for $`N_F=240`$. When $`N_t=4`$ or 8, the phase boundary line between confining and deconfining phase bends down at finite $`\beta `$ due to the finite temperature phase transition of the confining phase. The darker shaded lines are our conjecture for the phase boundary between the confining phase and the deconfining phase at zero temperature ($`N_t=\mathrm{}`$).

In Figs. 5. and 6., the dashed lines in the deconfining phase represents the chiral limit, $`m_q=0`$. This line also corresponds to the minimum point of $`m_\pi ^2`$. We have checked that, for the case $`N_F=240`$, the quark propagator in the Landau gauge actually shows the chiral symmetry, $`\gamma _5G(z)\gamma _5=G(z)`$, on the $`m_q=0`$ line . For $`N_F\begin{array}{c}>\hfill \\ \hfill \end{array}240`$, we find that the $`m_q=0`$ line in the deconfining phase, which starts at $`1/K=8`$ in the weak coupling limit $`\beta =\mathrm{}`$, reaches the strong coupling limit, as shown in Fig. 6. for $`N_F=240`$. For a smaller $`N_F\begin{array}{c}<\hfill \\ \hfill \end{array}100`$, because the bulk transition line $`K_d`$ in the strong coupling region shifts toward larger $`K`$ with decreasing $`N_F`$, the $`m_q=0`$ line in the deconfining phase hits the $`K_d`$ line at finite $`\beta `$. For example, in the case $`N_F=18`$, it hits at $`\beta =4.0`$ – 4.5, as shown in Fig. 5.(b).

### 4.1. RG-flow for $`N_F17`$

The $`m_q=0`$ point at $`\beta =\mathrm{}`$ is a trivial IR fixed point for $`N_F17`$. The phase diagram shown in Fig. 6. suggests that there are no other fixed points on the $`m_q=0`$ line at finite $`\beta `$.

Figure 7.(a) shows the results of $`m_\pi ^2`$ and $`2m_q`$ in the deconfining phase for $`N_F=240`$ at $`N_t=4`$. We find that the shape of $`m_\pi ^2`$ and $`2m_q`$ as a function of $`1/K`$ only slightly changes for $`1/K<8`$ when the value of $`\beta `$ decreases from $`\mathrm{}`$ down 0. We obtain similar results also for $`N_t=8`$. This suggests that quarks are almost free down to $`\beta =0`$ in the deconfining phase. The results for $`N_F=300`$ are essentially the same, except for very small shifts of the transition point.

We make a Monte Carlo renormalization group (MCRG) study along the $`m_q=0`$ at $`N_F=240`$. Performing a block transformation with scale factor 2, we estimate the quantity $`\mathrm{\Delta }\beta =\beta (2a)\beta (a)`$: We generate configurations on an $`8^4`$ lattice on the $`m_q=0`$ points at $`\beta =0`$ and 6.0 and make twice blockings. We also generate configurations on a $`4^4`$ lattice and make once a blocking. Then we calculate $`\mathrm{\Delta }\beta `$ by matching the value of the plaquette at each step. We obtain $`\mathrm{\Delta }\beta 6.5`$ at $`\beta =0`$ and $`10.5`$ at $`\beta =6.0`$.333 In order to get a more precise value of $`\mathrm{\Delta }\beta `$, one has to make many steps of block transformations, together with a fine tuning of parameters of the block transformation, and do matching using several types of Wilson loops. We reserve elaboration of this point for future works. For $`N_F=240`$, because the velocity of the RG flow is large, we will be able to obtain the correct sign and approximate value of $`\mathrm{\Delta }\beta `$ by a simple matching.  The value obtained from the two-loop perturbation theory is $`\mathrm{\Delta }\beta 8.8`$ at $`\beta =6.0`$. The signs are the same and the magnitudes are comparable. This suggests that the direction of the RG flow on the $`m_q=0`$ line at $`\beta =0`$ and 6.0 is the same as that at $`\beta =\mathrm{}`$. This further suggests that there are no fixed points at finite $`\beta `$. These imply that the theory is trivial for $`N_F=240`$.

The area of the deconfining phase decreases with decreasing $`N_F`$ in the strong coupling region. However, the shape and the values of $`m_\pi ^2`$ and $`m_q`$ are quite similar when we vary $`N_F`$ from 300 down to 17. Fig. 7.(b) shows $`m_\pi ^2`$ and $`2m_q`$ for $`N_F=18`$ – 300 at $`\beta =6.0`$ and 0. At $`\beta =6.0`$, the shapes of $`m_\pi ^2`$ are almost identical to each other, except for a small shift toward smaller $`1/K`$ as $`N_F`$ is decreased.

These facts suggest that, for $`N_F17`$, the nature of physical quantities in the deconfining phase is almost identical to that observed for $`N_F=240`$ and 300, i.e., quarks are almost free. Therefore, we conjecture that the direction of the RG flow at $`\beta \begin{array}{c}<\hfill \\ \hfill \end{array}6`$ is identical to the case $`N_F=240`$. Combining this with the perturbative result that $`\beta =\mathrm{}`$ is an IR fixed point for $`N_F17`$, we conjecture that the RG flow along the massless quark line in the deconfining phase is uniformly directing to smaller $`\beta `$, as in the case of $`N_F=240`$. When this is the case, the theory has only one IR fixed point at $`\beta =\mathrm{}`$, i.e. the theory is trivial for $`N_F17`$.

### 4.2. $`16N_F7`$

As discussed in Sec. 3., quark confinement is lost for $`N_F7`$ at $`\beta =0`$. We have performed simulations for $`N_F=12`$ and 7. The phase diagram for $`N_F=12`$ is shown in Fig. 5.(a). At $`\beta \begin{array}{c}<\hfill \\ \hfill \end{array}6.0`$ we studied, the gross feature of the phase diagrams for $`N_F=7`$ and 12 is quite similar to the case $`N_F=18`$ shown in Fig. 5.(b). Physical quantities at $`\beta \begin{array}{c}<\hfill \\ \hfill \end{array}6.0`$ also show behavior similar to that shown in Fig. 7. for $`N_F17`$. Therefore, we consider it probable that the direction of the RG flow in the deconfining phase at small $`\beta `$ is towards a larger $`\beta `$ as in the case of $`N_F17`$. However, the direction of the RG flow at $`\beta =\mathrm{}`$ is opposite to that for $`N_F17`$ because the theory is asymptotically free for $`N_F16`$. This means that we have an IR fixed point somewhere at a finite value of $`\beta `$. The continuum limit is governed by this IR fixed point.

On the $`m_q=0`$ line around $`\beta =4.5`$, we find that $`m_\pi `$ is roughly twice the lowest Matsubara frequency. This implies that the quarks are not confined and almost free. Therefore, the anomalous dimensions are small at the IR fixed point, suggesting that the IR fixed point locates at finite $`\beta `$. Unfortunately, as shown in Fig. 5., the massless line hits the boundary in the coupling parameter space. Without much space in the strong coupling region where simulation is easy, it is hard to study numerically the location of the IR fixed point. A detailed study of RG flow in a wider coupling parameter space is required. We reserve these studies for future works.

## 5. Conclusions

Summarizing our numerical results and conjectures, we propose Fig. 8. for the phase structure of QCD at $`T=0`$: When $`N_F6`$, quarks are confined for any values of the current quark mass and $`\beta `$. The chiral limit $`K_c(\beta )`$, where the current quark mass $`m_q`$ vanishes, belongs to the confining phase. When $`N_F7`$, there is no chiral limit in the confining phase. There is a line of a first order phase transition from the confining phase to a deconfining phase at a finite current quark mass for all values of $`\beta `$. The chiral limit exists only in the deconfining phase, and, therefore, can not be taken from the confining phase.

In order to see the orientation of the RG flow on the chiral limit line in the deconfining phase, we performed a MCRG study at $`N_F=240`$, because, in this case, we have sufficiently long massless quark line in the strong coupling region where numerical simulation is easy. We found that the RG flow along the massless quark line at small $`\beta `$ is consistent with the result of the perturbation theory. This suggests that the direction of the RG flow is uniform for all values of $`\beta `$, as shown in Fig. 8.. Therefore, there is only one IR fixed point at $`\beta =\mathrm{}`$, so that the theory in the continuum limit is trivial. Accordingly, we find that, in the deconfining phase, the $`K`$-dependences of pion screening mass and current quark mass at small $`\beta `$ are almost identical to those of free Wilson quarks.

Also for smaller values of $`N_F7`$, the general features of the phase diagram and the $`K`$-dependence of physical quantities are quite similar to the case of large $`N_F240`$. Therefore, we conjecture that the direction of the RG flow is the same as in the case of $`N_F=240`$, at least for $`\beta \begin{array}{c}<\hfill \\ \hfill \end{array}6.0`$ we have studied.

Together with the results of the perturbation theory, we conjecture the following picture.

* For $`N_F17`$, because the RG flow in the weak coupling limit has the same direction due to the lack of asymptotic freedom, we conjecture a uniform RG flow towards the trivial IR fixed point, as in the case of $`N_F=240`$. Thus the theory is trivial for $`N_F17`$.
* When $`7N_F16`$, because the theory is asymptotically free, the direction of the RG flow is opposite at large $`\beta `$. Therefore, we expect a non-trivial IR fixed point at finite $`\beta `$. In this case, the theory in the continuum limit is a non-trivial theory with anomalous dimensions. Because the coupling constant in the IR limit is finite, quarks are not confined in this case.
* When $`N_F6`$, quarks are confined and the chiral symmetry is spontaneously broken for any values of the current quark mass and $`\beta `$. The chiral limit belongs to the confining phase.

## ACKNOWLEDGEMENTS

Numerical simulations are performed with HITAC S820/80 at KEK, and Fujitsu VPP500/30 and QCDPAX at the University of Tsukuba. This work is in part supported by the Grants-in-Aid of Ministry of Education, Science and Culture (Nos. 09304029 and 10640248).