text:
# Effect of interface bonding on spin-dependent tunneling from the oxidized Co surface

## I Introduction

Spin-dependent tunneling (SDT) in magnetic tunnel junctions (MTJs) is a dynamically developing area of research that attracted a lot of attention due to promising applications in non-volatile random access memories and next-generation magnetic field sensors (for a recent review of SDT see Ref. TML, ). The experimental efforts have succeeded in achieving large reproducible tunneling magnetoresistance (TMR) in MTJs Moodera  but also raised fundamental questions regarding the nature of SDT. One such question is the role of the ferromagnet/insulator interfaces in controlling the spin polarization (SP) of the tunneling conductance. The SP is usually defined as $`P=(G_{}G_{})/(G_{}+G_{})`$, where $`G_\sigma =(e^2/h)_𝐤_{}T_\sigma (𝐤_{})`$ is the conductance for spin channel $`\sigma `$, $`T_\sigma `$ is the transmission function, and $`𝐤_{}`$ is the lateral wave vector.

Commonly the expected spin dependence of the tunneling current is deduced by considering the symmetry of the Bloch states in the bulk ferromagnetic electrodes and the complex band structure of the insulator. Maclaren ; Mavro  By identifying those bands in the electrodes that are allowed by symmetry to couple efficiently to the evanescent states decaying most slowly in the barrier one can predict the SP of the conductance. However, this approach has two deficiencies. First, it assumes that the barrier is sufficiently thick so that only a small focused region of the surface Brillouin zone (SBZ) contributes to the tunneling current. For realistic MTJs with a barrier thickness of about 1 nm this assumption is usually unjustified. Second, symmetry considerations alone applied to bulk materials are not always sufficient to predict the SP. It is critical to take into account the electronic structure of the ferromagnet/barrier interfaces which, as it was shown both experimentally Leclair  and theoretically, UL  controls SDT.

An important mechanism by which the interfaces affect the SP of the tunneling current is the bonding between the ferromagnetic electrodes and the insulator. TP97  This effect was put forward to explain positive and negative values of TMR depending on the applied voltage in MTJs with Ta2O5 and Ta2O5/Al2O3 barriers Sharma  and to elucidate the inversion of the SP in Co/SrTiO3/La0.67Sr0.33MnO3 MTJs. Deteresa  So far there are no theoretical studies explaining the microscopic origin of this phenomenon.

In this paper we report the results of first-principles calculations of SDT from clean and oxidized Co surfaces through vacuum into Al and demonstrate the crucial role of the bonding between Co and O atoms. This system was chosen for investigation because the MTJs based on alumina have predominantly O-terminated Co/Al2O3 interfaces. oxidized  By replacing alumina by vacuum we can ignore the complexity of the atomic structure of the amorphous alumina and focus on effects of surface oxidation. Moreover, this system can be directly studied using spin-polarized STM. SP-STM,

We show that a monolayer of oxygen on the Co surface creates a spin-filter effect due to the Co-O bonding by producing an additional tunneling barrier in the minority-spin channel. This reverses the sign of the SP from negative for the clean Co surface to positive for the oxidized Co surface, thus revealing the crucial role of interface bonding in SDT.

## II Surface transmission functions

We approach the tunneling problem in the spirit of perturbation theory. Harrison  We consider the system consisting of “left” and “right” leads separated by a relatively thick barrier and assume that two-dimensional translational periodicity in transverse directions is preserved for both electrodes including their surfaces (although it may be different for each subsystem). Every Bloch wave with a lateral quasi-wave vector $`𝐤_{}^L`$ coming from the left lead has a decay tail in the vacuum composed of the waves with transverse wave vectors $`𝐤_{}^L+𝐆_i`$ where $`𝐆_i`$ are the reciprocal lattice vectors of the SBZ of the left lead. Mazin  At sufficient distances from the surface (typically just a few lattice parameters for low-index surfaces) all waves with $`𝐆_i0`$ vanish and may be neglected in the Fourier expansion of the Bloch wave in the vacuum. On the other hand, the quasi-wave vector $`𝐤_{}^R`$ is the good quantum number for the right lead. This means that $`𝐤_{}`$ is conserved across the entire system even if there is no common in-plane periodicity and $`𝐤_{}^L`$, $`𝐤_{}^R`$ are defined with respect to different SBZs. In this case, each tunneling eigenstate is almost identical to an evanescent plane wave in the central region of the barrier.

For a given $`𝐤_{}`$ the transmission function is the sum of the transmission coefficients for all tunneling eigenstates corresponding to all incoming Bloch waves with this $`𝐤_{}`$ in the left electrode. At the same time, each transmission coefficient for a given eigenstate contains a sum over outgoing states in the right electrode with the same $`𝐤_{}`$. Let us choose a reference plane in the vacuum region at a sufficient distance from the surface of an electrode, so that the eigenstates for all $`𝐤_{}`$ are already indistinguishable from the barrier eigenstates at this plane (see Fig. 1). For each tunneling eigenstate the amplitude of the barrier eigenstate between the reference planes is the only parameter coupling the left and right electrodes. Then, the S-matrix element coupling the states in the two electrodes may be written as

$$S_{pq}^\sigma =S_{pr}^\sigma S_{rr^{}}S_{r^{}q}^\sigma $$
(1)
where $`p`$ is the incoming and $`q`$ the outgoing Bloch states in the left and right leads respectively, $`r`$ and $`r^{}`$ denote the same vacuum eigenstate at the left and right reference planes. We omitted the dependence on $`𝐤_{}`$ for all the S-matrices in Eq. (1) for brevity. The vacuum S-matrix $`S_{rr^{}}`$ simply describes the exponential decay of the wave function in the vacuum. Note that no summation is implied in (1), because the state $`r`$ is uniquely defined by $`𝐤_{}`$. The simple product of S-matrices in Eq. (1) without any multiple scattering terms is a consequence of our assumption that the barrier is sufficiently thick. Thus, we see that the transmission function $`T(𝐤_{})`$ of the MTJ is factorized:

$$T_\sigma (𝐤_{})=t_L^\sigma (𝐤_{})\mathrm{exp}\left[2\kappa (𝐤_{})d\right]t_R^\sigma (𝐤_{}).$$
(2)
Here we replaced $`|S_{rr^{}}|^2`$ by its explicit exponential form with

$$\kappa (𝐤_{})=\left(\frac{2m\varphi }{\mathrm{}^2}+k_{}^2\right)^{1/2},$$
(3)
where $`\varphi `$ is the work function, and $`d`$ is the distance between the reference planes assigned to the electrodes as shown in Fig. 1. All the information about the properties of individual surfaces is described by the *surface transmission functions* (STF) $`t_L^\sigma `$, $`t_R^\sigma `$:

$$t_L(𝐤_{})=\underset{p}{}\left|\frac{B_r}{A_p}\right|^2,t_R(𝐤_{})=\underset{q}{}\left|\frac{C_q}{B_r^{}}\right|^2,$$
(4)
where the four amplitudes characterize the behavior of tunneling eigenstates at the two surfaces (see Fig. 1). The different definitions of $`t_L`$ and $`t_R`$ are due to the fact that they obey different boundary conditions. Specifically, $`t_L`$ and $`t_R`$ are obtained by matching the Bloch wave functions with the vacuum eigenstates, decaying and growing into the vacuum respectively. The definition of $`t_L`$ implies the solution of a scattering problem for the incoming wave with amplitude $`A_p`$, and $`B_r`$ is the amplitude of this scattering eigenstate at the reference plane in the vacuum. On the other hand, $`t_R`$ describes an inverse scattering problem in which the exponentially decaying wave in the vacuum with amplitude $`B_r^{}`$ at the reference plane is scattered on the right surface; here $`C_q`$ is the amplitude of the outgoing Bloch state $`q`$ in the right electrode for this “eigenstate.” Note that physically, this state is forbidden because it grows to infinity in the vacuum, but it is still a formal solution of the Schrödinger equation with the specified boundary condition at infinity.

Each surface may be considered as forming the left or the right interface of the tunnel junction. It is straightforward to show that, just as in the case of a transmission matrix connecting the propagating states in the two electrodes, Datta  the requirement of flux conservation results in the reciprocity condition $`t_L(𝐤_{})=t_R(𝐤_{})t(𝐤_{})`$ for any (laterally periodic) surface, as long as the appropriate normalization of the wave functions is chosen. Specifically, all Bloch waves in the electrodes should be normalized for unit flux, while the vacuum eigenstates should be normalized for unit “imaginary flux” $`\kappa /m`$. The easiest way to establish this reciprocity condition is to prove it for free electrons in a semi-infinite potential well, and then to use this free-electron system as a second electrode in a junction. Since the total transmission function for a two-electrode system is reciprocal, Datta  this proves the above reciprocity condition for the STF. Note that this proof does not rely on time reversal symmetry, because time reversal also replaces $`𝐤_{}`$ by $`𝐤_{}`$. However, this symmetry gives another useful relation, $`t(𝐤_{},𝐇)=t(𝐤_{},𝐇)`$ where $`𝐇`$ is the external magnetic field.

It is easy to see from Eq. (4) with the chosen unit flux normalization that the STF is proportional to the Fermi-level value of the $`𝐤_{}`$\- and energy-resolved charge density for the given spin, which is generated by the incoming Bloch states and taken at the reference plane (any localized surface states are excluded). Indeed, the energy-resolved charge density may be written as

$$\rho _\sigma (𝐫,E)=\frac{1}{N_{}}\underset{𝐤_{}}{}\rho _\sigma (𝐫,𝐤_{},E)$$
(5)
where the $`𝐤_{}`$\- and energy-resolved charge density is

$$\rho _\sigma (𝐫,𝐤_{},E)=\underset{p}{}\left|\psi _{\sigma 𝐤_{}p}^F(𝐫)\right|^2$$
(6)
where the summation is over the incoming Bloch states with the given $`𝐤_{}`$ and $`E`$, and the superscript $`F`$ emphasizes that the scattering eigenstates $`\psi _{\sigma 𝐤_{}p}^F`$ are normalized so that the incoming Bloch waves carry unit flux normal to the interface. Setting $`A_p=1`$ in Eq. (4), we see that $`t(𝐤_{})`$ is given by Eq. (6) where $`𝐫`$ is taken at the reference plane (since we neglect all components with $`𝐆_i0`$ at the reference plane, $`\rho (𝐫,𝐤_{},E)`$ does not depend on $`𝐫_{}`$, the location within this plane).

Although we considered a vacuum tunneling barrier, the analysis can be extended to other physically important cases of insulating barriers. Indeed, the main requirement for the validity of Eqs. (1,2) is that the complex band structure be predominantly represented by a single evanescent wave for each $`𝐤_{}`$. At a minimum, this premise must hold for the “active” regions of the SBZ that contribute appreciably to the tunneling current. This means that other tunneling states must have a notably larger imaginary part of the wave vector compared to the dominant one. In the case of a vacuum barrier, this amounts to the neglect of all vacuum eigenstates with $`𝐆_i0`$. In practice this criterion is well satisfied for sufficiently thick wide-gap tunnel barriers, such as $`sp`$-bonded oxides. However, for any particular barrier, this assumption has to be carefully verified by checking the complex band structure of the barrier for the presence of additional slowly decaying states in the active regions of the SBZ. It is important to note that the area of the active region quickly shrinks as the barrier thickness is increased,Mavro  indicating that the single evanescent state criterion will often lead to an additional requirement of “sufficient thickness” of the barrier. Note that this requirement is not related to a similar one stemming from the neglect of back-scattering.

The factorization (2) allows us to study tunneling between completely different electrodes. For practical purposes, the STF can be calculated using $`T_\sigma (𝐤_{})`$ obtained for a symmetric MTJ by factoring out the vacuum decay factor for the given choice of reference planes, and then taking the square root. The resulting STFs for different surfaces may then be convolved with an appropriate vacuum decay factor, $`\mathrm{exp}(2\kappa d)`$, to obtain the transmission functions for asymmetric MTJs. In particular, $`T_\sigma (𝐤_{})`$ for the antiparallel (AP) magnetization of the two electrodes (and the TMR) can be found from the up- and down-spin STFs calculated from the transmission function for the parallel (P) magnetization.

## III Generalization of Jullière’s formula

Tunneling magnetoresistance is often discussed in terms of Jullière’s formulaJulliere

$$\mathrm{TMR}=\frac{2P_LP_R}{1P_LP_R}$$
(7)
where $`P_L`$, $`P_R`$ are the “spin polarizations” of the left and right electrodes. Eq. (7) is derived from the empirical expression $`G^\sigma \rho _L^\sigma \rho _R^\sigma `$, where $`\rho _L^\sigma `$, $`\rho _R^\sigma `$ are the “tunneling densities of states” of the electrodes, and the spin polarization is defined as $`P=(\rho ^{}\rho ^{})/(\rho ^{}+\rho ^{})`$. The popularity of this formula is due to the fact that it usually agrees reasonably well with experiment if the spin polarizations, which are directly related to the “tunneling densities of states,” are taken from Meservey-Tedrow experiments MT  with the same barrier. However, the validity of Jullière’s formula has been debated for a long time, and the reasons for its apparent agreement with experiment are unclear. The physical meaning of the “tunneling density of states” is also unclear, but it is obvious both from elementary quantum mechanics and from experiments that the tunneling properties of a magnetic heterostructure are determined not by the ferromagnet alone, but rather by the ferromagnet/barrier combination and by the structure of the interface. A number of explicit first-principles calculations for idealized MTJs without disorder (see, e.g., Ref. Maclaren97, ) confirmed this fact. However, it was suggested MathUm  that phase decoherence due to disorder which is always present in realistic MTJs may recover the factorization of the tunneling conductance in a product of “transport densities of states,” which are essentially equal to the regular densities of states at the surfaces of the electrodes if there are no resonant localized states in the barrier. Moreover, it was shown TP98  within a single-band tight binding model that in the limit of strong disorder one recovers Jullière’s formula (7) by identifying $`P_L`$, $`P_R`$ with the measurable spin polarizations of the tunneling current for the same electrode/barrier systems. Therefore, it seems that there are good reasons for the widespread use of Jullière’s formula, and it is highly desirable to elucidate these reasons.

Let us explore the connection between Eq. (2) and Jullière’s formula (7). In Eq. (2) the simple product of the “tunneling densities of states” is replaced by a convolution of STFs, which explicitly include the effects of bulk densities of states and of the surface structure. Thus we can consider Eq. (2) as a generalization of Jullière’s formula for an ideal MTJ with no disorder.

However, we may go further and identify limiting cases where Eq. (2) can be directly related to Jullière’s formula, providing formal definitions of the “tunneling densities of states” of the electrodes appropriate for these cases. First, consider the case of a disordered insulating barrier. Such a barrier may be characterized by its eigenstates, half of which are decaying from left to right, and the other half from right to left. Although these eigenstates do not have a conserved $`𝐤_{}`$ any more, it is still clear that tunneling will be dominated by Feynman paths that do not “loop back,” because each path carries a weight decaying exponentially with its length (see below). Therefore, we may still write an expression similar to Eq. (1) neglecting back-scattering, but now we should sum up over all barrier eigenstates (now defined in real space). Within this formulation $`S_{rr^{}}`$ is still diagonal because it describes the decay of a single eigenstate.

The weight of a Feynman path in the imaginary-time functional integral often used for tunneling problemsColeman  (with Euclidean action written in its reduced Maupertuis form; see, e.g., Ref. Katsnelson, ) is given by $`\mathrm{exp}\left[\kappa (l)𝑑l\right]`$ up to a prefactor, where the integral is taken along the path, $`\kappa =[2m(VE)]^{1/2}`$ and $`V(𝐫)`$ is the potential. In an ordered insulator many paths with similar weights contribute to the path integral resulting in the formation of the complex band structure. However, in a disordered insulator the tunneling current may be dominated by Feynman paths running close to a relatively small number of “easy” paths with locally maximum weights, i.e. by imaginary-time classical paths.Coleman ; Katsnelson  If there is only one such channel or one class of channels with similar properties (e.g., due to surface roughness), Eq. (2) will produce Jullière’s formula where $`\rho ^\sigma `$ is simply the Fermi-level value of the energy-resolved charge density given by Eq. (6) integrated over $`𝐤_{}`$ and taken at some reference point within the channel. (Now each term describes the scattering eigenstate corresponding to the single incoming Bloch wave with the given $`𝐤_{}`$.) Like STF, this quantity does not depend on the properties of the other side of the barrier. This conclusion agrees with the results of Ref. TP98,  showing that the tunneling current through a strongly disordered barrier is dominated by a small number of random configurations, and that Jullière’s formula is also recovered in this limit.

Now consider the case when disorder is weak close to the interfaces, but remains strong in the insulator. Obviously, the S-matrix of the disordered insulator in $`𝐤_{}`$ representation will be essentially a random matrix, and after averaging Eq. (2) thus yields Jullière’s formula with $`\rho ^\sigma _𝐤_{}t^\sigma (𝐤_{})`$. This case is the easiest from the computational point of view, because the STFs may be directly calculated for a $`𝐤_{}`$-conserving MTJ.

It is instructive to compare this result with the conclusions of Mathon and Umerski MathUm  on the applicability of Jullière’s formula obtained using the transfer Hamiltonian formalism. Our approach shares in common with Ref. MathUm,  the neglect of multiple reflections across the junction. However, the assumption of constant matrix elements (hopping integrals) for all Bloch waves made in Ref. MathUm,  completely removes all physical effects connected with orbital- and spin-dependent bonding at interfaces. This obviously contradicts the experimental findings showing that the spin polarization of the tunneling current and magnetoresistance strongly depend on the type of barrier used.TML  In our approach, the STFs for the electrodes allow us to encapsulate the effects of the interface structure and provide the proper dependence of the tunneling current on barrier type. Jullière’s formula obtained in the limiting case of full decoherence inside the insulator is expressed in terms of the spin polarization actually measured in the Meservey-Tedrow experiment (assuming that the superconductor acts as an ideal, non-biased spin detector).

Finally, for very thick $`𝐤_{}`$-conserving barriers the tunneling current may be carried predominantly by a close vicinity of some special $`𝐤_{}`$-points in the SBZ (e.g., the $`\mathrm{\Gamma }`$ point). In this case, the tunneling density of states is simply equal to the value of $`T(𝐤_{})`$ at this $`𝐤_{}`$.

It is not clear *a priori* whether any one of these three limiting cases is directly applicable to realistic MTJs, although it seems that disorder in the insulator together with the “channelization” of the tunneling current are both likely to play a major role. However, the emergence of Jullière’s formula in these different scenarios suggests that it may actually have a rather wide range of applicability. In general, the “tunneling density of states” should be identified with some appropriately averaged energy-resolved charge density taken at the Fermi level at a sufficient distance from the interface within the barrier. Unlike the bulk density of states, this function fully takes into account the relevant properties of the surface.

## IV Tunneling from clean and oxidized Co (111) surfaces through vacuum into Al

We calculated the transmission functions using the principal-layer Green’s function approach Turek-book  based on the tight-binding linear muffin-tin orbital method (TB-LMTO) in the atomic sphere approximation (ASA) and the transmission matrix formulation of Ref. Kudr, . All the atomic potentials were determined self-consistently within the supercell approach using the TB-LMTO-ASA method. The vacuum barrier was modeled using empty spheres in the positions corresponding to the continuation of the crystal lattice of the electrodes. We have also performed full-potential LMTO calculations FLMTO  which confirmed all main features of the band structure of the oxidized Co (111) surface discussed below.

We checked the validity of factorization (2) by calculating $`T_\sigma (𝐤_{})`$ for (100) and (111)-oriented fcc Co electrodes with parallel magnetizations, taking the square root, and convolving $`t_{}(𝐤_{})`$ with $`t_{}(𝐤_{})`$. Then, the result was compared with the independent calculation for the antiparallel configuration in a range of energies. The agreement was always excellent (better than 1%), except for a couple of specific energies for a (100) MTJ with four vacuum “monolayers” where narrow resonances appear in the minority channel. Wunnicke  If the vacuum barrier is extended to 8 monolayers, excellent agreement is restored.

Using the factorization (2) we investigated the SP of the conductance from ferromagnetic electrodes to a non-magnetic material, Al (111), which served as a detector of the tunneling SP in the spirit of the Meservey-Tedrow experiments. MT  As expected, the calculated STF of Al is free-electron-like, having almost perfect Gaussian shape originating from the vacuum decay factor up to the reference plane. Therefore, this surface may be considered as equally transparent for all Bloch waves, and the total transmission function for a MTJ with Al spin-detector electrode is essentially a product of the other electrode’s STF and the vacuum decay factor.

First, we discuss the properties of a Co/vacuum/Al MTJ with a clean Co (111) surface. Figs. 2a,b show the $`𝐤_{}`$-resolved transmission for the majority- and minority-spin electrons within the SBZ of Co (111). The Fermi surface of Co viewed along the  direction has holes close to the $`\overline{\mathrm{\Gamma }}`$ point with no bulk states in both spin channels, which results in zero conductance in this area. The majority-spin transmission (Fig. 2a) varies rather smoothly and is appreciable over a relatively large area of the SBZ. On the other hand, the minority-spin transmission (Fig. 2b) has a narrow crown-shaped “hot ring” around the edge of the Fermi surface hole. The analysis of layer and $`𝐤_{}`$-resolved density of states (DOS) shows that it is not associated with surface states, Wunnicke  but rather with an enhancement of bulk $`𝐤_{}`$-resolved DOS near the Fermi surface edge (compare Fig. 2b with Fig. 4c).

As seen from Figs. 2a,b, the Fermi surface hole is smaller for majority spins. Therefore, the conductance should become fully majority-spin polarized in the limit of very thick barriers. However, since the Fermi surface hole is also quite narrow for minority spins, positive SP is only achieved at very large barrier thicknesses $`d10`$ nm, while for typical values of $`d2`$ nm the SP is about $`60\%`$ and depends weakly on $`d`$.

The oxidized Co surface was modeled by an O monolayer placed on top of the Co(111) electrode. The equilibrium atomic structure of this surface was found using the pseudopotential plane-wave method CASTEP  within the generalized gradient approximation. We used both types of stacking: ABCA and ABCB, where the last symbol designates the position of the O monolayer. The O atoms are assumed to lie in symmetric positions above the second (S2) or first (S1) sub-surface Co layer, respectively (these layers are shown in Fig. 1 in the absence of the O monolayer). We allowed the O layer and two Co layers (S1 and S2) to relax in the direction normal to the surface, while the positions of atoms in deeper layers (S3,…) were kept fixed. The energies of both equilibrium configurations of oxygen monolayers were found to be very close to each other. All results of interest in the present context are quite similar for these two stackings. Below all specific data are given for the ABCB stacking. The equilibrium interlayer distances were found to be 2.14 Å between layers S3 and S2, 2.18 Å between S2 and S1, and 1.08 Å between S1 and O layers, compared to 2.07 Å between the adjacent Co layers in the bulk. The Co-O bond length is equal to 1.82 Å.

Presence of oxygen at the surface of cobalt raises the question of whether electron correlations similar to those characteristic for transition-metal oxides may be strong enough to induce significant changes in the band structure at the oxidized Co surface. However, the enhancement of correlations in oxides is due to much weaker screening of Coulomb interaction compared to the metallic state. On the other hand, cobalt atoms below the oxygen monolayer preserve the close-packed configuration of bulk cobalt except for the top three nearest neighbors being absent out of twelve. Therefore, it is reasonable to expect that screening of Coulomb interaction in the $`3d`$ shell is not much weaker compared to the bulk. For this reason, we believe that LDA electronic structure of the oxidized Co surface is correct as far as the main features are concerned.

The oxygen monolayer dramatically changes the electronic structure of the underlying Co layer making this layer almost magnetically-dead. This change can be understood from band dispersion plots shown in Fig. 3. For each spin, the free-standing oxygen monolayer has three energy bands deriving from $`2p`$ states, each doubly degenerate due to $`\sigma _z`$ reflection symmetry ($`z`$ is the axis normal to the surface). When the monolayer is deposited onto the Co surface, the degeneracy is lifted, and two sets of three bands each are formed corresponding to bonding and antibonding mixing of oxygen and cobalt orbitals. The three bonding bands (marked B in Fig. 3) lie well below the bulk Co $`3d`$ band, whereas the antibonding states are close to the Fermi energy $`E_F`$. As a result of this bonding the local DOS for the S1 layer at $`E_F`$ is strongly reduced, so that, according to the Stoner criterion, magnetism in this layer is almost completely suppressed. The magnetic moment of Co atoms in the S1 layer is only 0.17 $`\mu _B`$.

Transmission of propagating bulk states from the electrode through the barrier is very sensitive to the degree of mixing of these states with the antibonding surface states. This mixing is controlled by a selection rule which follows from the fact that all bands can be classified as “even” or “odd” according to their symmetry with respect to $`\overline{\mathrm{\Gamma }}\overline{M}`$ reflection ($`\overline{\mathrm{\Gamma }}\overline{M}`$ is the projection of a reflection plane). Although this classification is exact only along the $`\overline{\mathrm{\Gamma }}\overline{M}`$ direction, it is approximately valid throughout the entire SBZ.

According to this classification, two of the three surface bands are even, and one is odd. On the other hand, the free-electron-like band of bulk Co which forms the only Fermi surface sheet for majority-spin electrons is even, while the minority-spin states on the Fermi surface sheet closest to the $`\overline{\mathrm{\Gamma }}`$ point are odd. Even and odd bands are orthogonal and can not mix.

This selection rule results in the principal difference between the majority- and minority-spin transmission. The majority-spin bands are shown in Fig. 3a. One even and one odd antibonding surface bands (marked A) are degenerate at the $`\overline{\mathrm{\Gamma }}`$ point (at about 0.8 eV above $`E_F`$). At a short distance from the $`\overline{\mathrm{\Gamma }}`$ point both bands enter the continuum of bulk states. The odd band does not mix with the bulk states along the $`\overline{\mathrm{\Gamma }}\overline{M}`$ line and remains almost flat due to repulsion from a lower-lying band. On the other hand, the even band readily mixes with the free-electron-like majority-spin band crossing the Fermi level and completely loses its surface localization. This is evident from Figs. 4a,b which show the $`𝐤_{}`$-resolved DOS of majority-spin electrons for the bulk-like S6 and surface S1 Co layers. The DOS for the S1 layer (Fig. 4b) is appreciable in the entire area of the SBZ where bulk states are available and has no sharp features that might indicate localized surface states. This implies that the bulk majority-spin states extend to the very surface of the electrode and therefore can readily tunnel through the barrier.

The situation is very different for minority-spin states. Although the odd surface band is again almost flat and lies above $`E_F`$, the even surface band crossing the Fermi level does not mix with the odd minority-spin band. As a consequence, the $`𝐤_{}`$-resolved DOS for Co(S1) layer (Fig. 4d) is large only along the curve lying at the periphery of the SBZ where the Co-O antibonding surface band crosses the Fermi level (oxygen DOS looks very similar). As a result, the bulk minority-spin states responsible for most tunneling transmission from the clean surface only extend up to the S2 layer, encountering a band gap in the S1 and oxygen layers. Thus, an additional tunneling barrier is introduced in the minority-spin channel, and the total SP of the tunneling transmission becomes almost 100% positive, which is evident from Figs. 2c,d.

The predicted effect of interface bonding is not limited only to the Co(111) surface. We have also calculated the transmission from clean and oxidized Co(100) and Ni(111) electrodes and found that surface oxidation reverses the SP due to the bonding between Co or Ni with O. As it was shown earlier, the reversal of the SP also occurs for the Fe(100) surface. TOP

## V Conclusion

We have shown that the problem of calculating of the transmission function for a sufficiently thick insulating barrier is reduced to the solution of three separate problems, namely the penetration of the bulk wave functions into the barrier from both sides, and the behavior of the evanescent barrier eigenstates. This separation provides a natural generalization of Jullière’s formula. We identified three limiting cases when the original Jullière’s formula is recovered. The “tunneling density of states” in this formula is identified with an appropriately averaged energy-resolved charge density generated by the bulk Bloch states within the barrier and taken at the Fermi level.

Using the factorization of the transmission function into a product of surface transmission functions and a barrier decay factor we calculated the spin polarization of the tunneling current from clean and oxidized Co (111) surfaces through vacuum into Al. We showed that the bonding between Co and O atoms at the oxidized surface controls SDT by creating an additional barrier for minority-spin electrons, which results in a reversal of the spin polarization.

Experimentally, the reversal of the SP associated with surface oxidation may be detected using spin-polarized STM measurements. SP-STM  Since the ferromagnetic tip is sensitive to the SP of the total local DOS above the surface (see, e.g., Ref. Heinze, ), the TMR in the system surface/vacuum/tip should change sign when the Co surface is oxidized. In other words, for the clean Co (111) surface the tunneling current should be higher when the magnetizations of the tip and the surface are aligned parallel (the dominating minority channel is then open), but for the oxidized surface it should be higher for the antiparallel configuration.

###### Acknowledgements.

We are grateful to V. P. Antropov for useful discussions. This work is supported by NSF (DMR-0203359 and MRSEC DMR-0213808) and the Nebraska Research Initiative. MvS was supported by DARPA and ONR. Work at Sandia National Laboratories is supported by U.S. Department of Energy under contract No. DE-AC0494AL85000.


format: text

text:
# No naked singularities in homogeneous, spherically symmetric bubble spacetimes?

## I Introduction

Within five-dimensional Kaluza-Klein theory, negative energy initial data configuration at a moment of time symmetry have been presented in the past brillpfister ; brillhorowitz . Among the reasons for considering negative energy solutions is that naked singularities are associated to them. Therefore studying such solutions is an attractive problem to further test the cosmic censorship conjecture penrose . Additionally, as these bubbles can be obtained via semiclassical tunneling (following Witten’s construction witten ; aharony ), it is important to understand their dynamical behavior. To obtain a better feeling for the dynamics of the solutions presented in brillhorowitz , which describe a bubble spacetime, Corley and Jacobson corleyjacobson  analyzed the initial acceleration of the bubble (which can be calculated analytically with the information contained on a single hypersurface.) They found that negative mass bubbles start out expanding and consequently argued a naked singularity would be unlikely. This conjecture is backed by the calculations at the initial hypersurface together with the intuition gained from a time symmetric scenario. Namely, that if the bubble were to reverse its behavior it will go through another time-symmetric phase111However, we will see that one can encounter situations in which a bubble reverses its behavior without passing through a moment of time symmetry., which would also suggest the bubble should expand (assuming the solution at this time of symmetry still belongs to the family presented in brillhorowitz ).

Recently however, a numerical study was presented shinkai  which computed explicitly the solution to the future of the initial, time-symmetric, initial data. In the negative mass case, the solution found appeared to indicate a very puzzling development. Namely, that as the bubble expanded, it would encounter a naked singularity on its way.

In the present work we reexamine the dynamical behavior of the data presented in brillhorowitz  via numerical simulations, and further study a more general case where the initial data describes an initially collapsing bubble with negative mass. We restrict our simulations to the context of homogeneous spherically symmetric spacetimes. To summarize our findings, we see that negative mass bubbles will not lead to a naked singularity. By choosing different members of the initial data family, these could start out either expanding or collapsing. In the former case, the bubble continues on expanding while in the latter case, the bubble bounces back without ever collapsing. We additionally study the case of a positive mass bubble which starts out collapsing. Here, depending on the strength of the gauge field at the initial hypersurface, the bubble collapses to a black string or bounces back to expand “forever” (i.e. for as long as we let the simulation run). Furthermore, these two distinct possibilities give rise to critical behavior on the threshold of collapse or expansion.

## II Set up: Bubble data

We consider a generalization of the time symmetric family of initial data presented in brillhorowitz . We start with a spacetime endowed with the metric

$$ds^2=dt^2+U(r)dz^2+\frac{dr^2}{U(r)}+r^2d\mathrm{\Omega }^2,$$
(1)
where $`d\mathrm{\Omega }^2=d\vartheta ^2+\mathrm{sin}^2\vartheta d\phi ^2`$ is the standard metric on the unit two-sphere $`S^2`$ and $`U(r)`$ is a smooth function that has a regular root at some $`r=r_+>0`$, is everywhere positive for $`r>r_+`$ and converges to one as $`r\mathrm{}`$. The coordinate $`z`$ parametrizes the extra dimension $`S^1`$ which has the period $`P=4\pi /U^{}(r_+)`$. The resulting spacetime $`\{t,z,rr_+,\vartheta ,\phi \}`$ constitutes a regular manifold with the topology $`\times ^2\times S^2`$. The bubble is located where the circumference of the extra dimension shrinks to zero, that is, at $`r=r_+`$.

Additionally, we consider the presence of a U(1) gauge field in the following form,

$$A_adx^a=k(r_+^nr^n)dz,$$
(2)
with $`k`$ an arbitrary constant and $`n`$ an integer greater than one. It is important to point out two things. First, the gauge field here considered has a more generic dependence than those previously considered in brillhorowitz , where only the case $`n=2`$ was discussed. Second, the symmetries of the problem would a priori allow for the gauge field to have non-trivial $`t`$ and $`r`$ components. However, it follows from Maxwell’s equations that those contributions would lead to a Coulomb-like electric field which is divergent at the location of the bubble. For this reason we set those components of $`A_a`$ to zero.

In the time-symmetric case, initial data satisfying the Hamiltonian constraint obeys

$$U(r)=1\frac{m}{r}+\frac{b}{r^2}\frac{\stackrel{~}{k}^2}{r^{2n}},$$
(3)
where $`\stackrel{~}{k}kn/\sqrt{(n1)(2n1)}`$. The parameter $`m`$ is related to the ADM mass via $`M_{ADM}=m/4`$.

The fact that the bubble be located at $`r=r_+`$ requires that $`0=U(r_+)=1\overline{m}+\overline{b}\overline{k}^2`$, where $`\overline{m}m/r_+`$, $`\overline{b}b/r_+^2`$, $`\overline{k}\stackrel{~}{k}/r_+^n`$. We also require

$$0<r_+U^{}(r_+)=2\overline{m}+2(n1)\overline{k}^2$$
(4)
and avoid the conical singularity at $`r=r_+`$ by fixing the period of $`z`$ to $`P=4\pi /U^{}(r_+)`$. Repeating the analysis in corleyjacobson  one finds that the initial acceleration of the bubble area is given by

$$\ddot{A}=8\pi \left[1\overline{m}\frac{4\overline{k}^2}{3}(n1)(n2)\right].$$
(5)
For $`n=2`$ this agrees with the formula found by Corley and Jacobson, and as they discussed, for negative mass bubbles, the initial acceleration is positive, and only if chosen appropriately can a positive mass bubble have negative initial acceleration. As we will see later, in the vacuum case, this acceleration remains negative leading to a collapse of the bubble and the formation of a black string. In the non-vacuum case however, the strength of the gauge field can modify this behavior completely. For weak enough $`k`$ the bubble continues to collapse whereas when $`k`$ is large the bubble area bounces back and expands. Furthermore, by fine-tuning the parameter $`k`$ we will show that the solution approaches a static bubble solution to the five-dimensional Einstein-Maxwell equations.

We also consider non-vacuum initial data with $`n>2`$ which allow for richer scenarios. A particularly interesting case is that describing an initially negative acceleration with negative mass. If this were possible, and the bubble continued to collapse, this would likely give rise to a naked singularity. Thus, we are interested in data for which $`\ddot{A}<0`$. In view of the inequality (4) and Eq. (5) this means that

$$2+2(n1)\overline{k}^2>\overline{m}>1\frac{4\overline{k}^2}{3}(n1)(n2).$$
(6)
From this we see that for $`n=2`$ the mass can only be positive. In contrast to this, for $`n>2`$ and $`\overline{k}`$ large enough such that $`4\overline{k}^2(n1)(n2)>3`$, negative mass bubbles can be considered.

## III Numerical implementation

A detailed description of the numerical techniques used in our code will be presented elsewhere threetenors . Here we point out its salient features.

1. The equations are written in first oder symmetric hyperbolic form, with a gauge choice that is such that the characteristic speeds are given by $`0`$, $`1`$ or $`1`$. These properties facilitate achieving a stable implementation of the equations.
2. Regularity conditions are enforced at the bubble and employed to design a numerical scheme that satisfies a discrete energy estimate lsu1 ; lsu2 . This discrete energy estimate guarantees numerical stability for a related toy model problem.
3. Constraint preserving boundary conditions are enforced at the outer boundary to ensure no constraint violating modes are fed through it calabreselehnertiglio ; sarbachetal .
4. The equations are implemented using second order accurate finite difference techniques. A non-uniform radial coordinate is employed to improve the resolution near the bubble. For all cases presented here, the radial grid covered $`r[1,20]`$ with a total of $`8000`$ points. The outer boundary’s location is chosen so as to have always less than a crossing time in our evolutions, guaranteeing outer boundary conditions (though consistent with the constraints) will not affect quantities at the bubble. Second order convergence is checked in all runs. As an example, the convergence of some of the simulations representing expanding bubbles is illustrated in Fig. 2. Similar behavior is observed in all cases222As mentioned, in the cases with negative mass, our results disagree with those presented in shinkai . We believe that the reason for the disagreement is that the numerical simulations presented in shinkai  were contaminated by ill-posed modes present in the weakly hyperbolic form of the equations and the implementation of the boundary conditions used there. The understanding of these issues and how to deal with them has advanced considerably in the past few years. In particular, points 1-3 mentioned in section III are major differences between the two implementations..

## IV Results

### IV.1 Brill-Horowitz initially expanding case: $`(n=2)`$

Here we evolve the Brill-Horowitz initial data in the case of vanishing gauge field. The bubble area $`A`$ as a function of the proper time $`\tau `$ at the bubble is shown in Fig. 1 for different values of the mass parameter $`m`$. (Figure 2 illustrates the convergence of the solution throughout the evolution). As expected, the lower the mass of the initial configuration, the faster the expansion. Empirically, and for the parameter ranges used in our runs, we found that at late times the expansion rate is given by

$$\frac{\dot{A}}{A}\frac{2\overline{m}}{r_+(\tau =0)},$$
(7)
where a dot denotes the derivative with respect to proper time $`\tau `$. In particular this approximation is valid for the bubble solution exhibited by Witten witten  which describes the time evolution in the case $`\overline{m}=0`$. We monitored several curvature invariant quantities (the Kretschmann invariant in particular) for our numerical spacetimes and found no sign of divergence.

### IV.2 Brill-Horowitz initially collapsing case: $`(n=2)`$

Next, we analyze the Brill-Horowitz initial data for the case in which the bubble is initially collapsing (notice that for $`n=2`$ this implies that the ADM mass is positive). While our numerical simulations reveal that in the absence of the gauge field such a bubble continues to collapse, we also show that when the gauge field is strong enough, the bubble shrinks at a rate which decreases with time and then bounces back.

Obviously, if the collapse trend were not halted, a singularity should form at the origin. Since the ADM mass is positive, one expects this singularity to be hidden behind an event horizon, and one should obtain a black string. In fact, for the solutions which are initially collapsing and which have vanishing gauge field, we observe the formation of an apparent horizon. Furthermore, we compute the curvature invariant quantity $`IR_{AH}^4`$ at the apparent horizon (as discussed in choptuikblackstring ), where $`I=R_{abcd}R^{abcd}`$ is the Kretschmann invariant and $`R_{AH}`$ the areal radius of the horizon. For a neutral black string, this invariant is $`12`$. Figure 3 shows how this value is attained after the apparent horizon forms for representative vacuum cases (with $`m=1.1`$ and $`m=1.99`$), thus providing strong evidence for the formation of a black string. (We also find apparent horizons forming, or already present in the initial data, when a nontrivial gauge field is considered and the parameters are appropriately chosen. The corresponding bubbles will likely collapse to charged black strings and will be analyzed elsewhere threetenors .)

As mentioned, for strong enough gauge fields, the previously described dynamics is severely affected. Figure 4 shows the bubble area vs. proper time for different values of $`k`$. For large values of $`k`$ the bubble “bounces” back and expands while for small ones the bubble collapses. There is a natural transition region separating these two possibilities. Tuning the value of $`k`$ one can reveal an associated critical phenomena, the ‘critical solution’ being a member of the family of static solutions given by

$`ds^2=V(r)dt^2+{\displaystyle \frac{V(r)}{U(r)}}dr^2+{\displaystyle \frac{U(r)}{V(r)^2}}dz^2+r^2V(r)d\mathrm{\Omega }^2,`$
$`A_adx^a=\pm {\displaystyle \frac{1}{2}}\sqrt{3\left({\displaystyle \frac{r_+}{r_{}}}1\right)}{\displaystyle \frac{dz}{V(r)}},`$
where $`V(r)=1r_{}/r`$ and $`U(r)=1r_+/r`$. These solutions can be obtained by a double analytic expansion of the black hole solutions described by Eqs. (2.13) and (2.15) of Ref. horowitzmaeda\_solution  followed by the transformations $`r_{}r_{}`$, $`rrr_{}`$, $`r_+r_+r_{}`$. The parameters $`r_{}`$ and $`r_+`$ $`(>r_{})`$ are related to the period of the $`z`$ coordinate and to the ADM mass via $`P=4\pi r_+(1r_{}/r_+)^{3/2}`$ and $`M_{ADM}=r_+/4`$. Since these parameters are conserved throughout evolution, the initial data fixes the member of the static family the dynamical solution approaches to. Details of the critical behavior exhibited by this system will be presented in threetenors .

### IV.3 Collapsing negative mass case

We here restrict to cases with negative masses that start out collapsing. Interestingly enough we find that even when starting with large initial negative accelerations, which in turn make the bubble shrink in size to very tiny values, it bounces back without ever collapsing into a naked singularity. As an example, Fig. 5 shows the bubble’s area versus time for different values of $`n`$ and $`k`$. The initially collapsing bubbles decrease in size in a noticeable way but this trend is halted and the bubbles bounce back and expand. Although we have not found a simple law as that in Eq. (7), clearly the bubbles expand exponentially fast. Therefore, it seems not to be possible to “destroy” the bubble and create a naked singularity. This situation is somewhat similar to the scenarios where one tries to “destroy” an extremal Reissner-Nordström black hole by attempting to drop into it a test particle with high charge to mass ratio. There, the electrostatic repulsion prevents the particle from entering the hole wald .

## V Conclusions

We have performed numerical simulations to investigate the dynamics of spherically symmetric, homogeneous bubble spacetimes. These have revealed several interesting features.

In the vacuum case, negative mass bubbles do not form singularities, rather they continue on expanding and the asymptotic rate of expansion obeys the simple relation (7). Positive mass cases, as pointed out previously can start out either expanding or contracting. In the former case, they continue on expanding and its asymptotic rate obeys the same simple relation. For the collapsing case we have evidences that the bubble will collapse to a black string. It is interesting to note here that the ratio of a circle at infinity over the ADM Mass is given by

$$\frac{P}{M_{ADM}}=\frac{16\pi }{\overline{m}(2\overline{m})}.$$
(8)
Since $`0<\overline{m}<2`$ this value is always larger than that needed for the Gregory-Laflamme “instability”gregorylaflamme . Thus, a slight $`z`$-dependent perturbation of the scenario considered here could indeed give rise to a naked singularity. This would require infinite time to appear as shown in horowitzblackstring , but certainly would give rise to rich dynamics as illustrated in choptuikblackstring .

In the non-vacuum case, the presence of a gauge field gives rise to additional features. First, for the positive mass case, it can “halt” the collapse and subsequently induce an expansion. Furthermore, by fine-tuning the strength of the field, the solution approaches a static bubble. This critical solution separates the collapsing bubbles (which, likely, collapse to a charged black string) from the one that expand at late times. In a subsequent paper threetenors  we will show that the dynamical behavior near the static bubble exhibits critical phenomena. Note that clearly the solution does encounter another moment of time symmetry at the critical value. However, for field strengths larger than the critical one, this is not the case. The bubble bounces but the solution does not go through a time-symmetric phase.

Last, the gauge field can also make negative mass bubbles start out collapsing if its radial dependence is steeper than that originally considered in brillhorowitz . Intuitively this happens due to the contribution of the energy in the gauge field overcoming the repulsive gravitational effect of the negative mass bubble. However, as part of the field is radiated away, the collapse is halted and a subsequent expansion results. By choosing both the strength of the field or its radial dependence appropriately (with $`n>2`$), we could make the bubble start out collapsing at a very rapid rate. In these cases, the bubble shrinks to very small sizes but eventually bounces back without forming a naked singularity.

###### Acknowledgements.

We thank Robert Myers for calling this problem to our attention and for pointing out the influence a steeper radial dependence in the gauge field could have. Additionally, we thank him and Gary Horowitz for many interesting discussions. We acknowlege the LSU gravity group for comments on the manuscript. We benefited from the hospitality of the Kavli Institute for Theoretical Physics. L.L. thanks the Perimeter Institute for Theoretical Physics for hospitality. This research was supported in part by the NSF under grant numbers NSF-PHY0244335 to Louisiana State University, PHY99-07949 to the University of California at Santa Barbara and by the Horace Hearne Institute for Theoretical Physics. L.L. has been partially supported by the Alfred P. Sloan Foundation.


text:
# Estimating the distribution of dynamic invariants: illustrated with an application to human photo-plethysmographic time series

## Acknowledgments

This research was supported by a Hong Kong Polytechnic University Research Grant (NO. A-PE46). The authors wish to thank J. Bhattacharya for supplying the photo-plethysmographic time series.


text:
# Two-step melting of the vortex solid in layered superconductors with random columnar pins

## Abstract

We consider the melting of the vortex solid in highly anisotropic layered superconductors with a small concentration of random columnar pinning centers. Using large-scale numerical minimization of a free-energy functional, we find that melting of the low-temperature, nearly crystalline vortex solid (Bragg glass) into a vortex liquid occurs in two steps as the temperature increases: the Bragg glass and liquid phases are separated by an intermediate Bose glass phase. A suitably defined local melting temperature exhibits spatial variation similar to that observed in experiments.

The mixed phase of type-II superconductors with random pinning constitutes an excellent test system for studies of the effects of quenched disorder on the structure and melting of crystalline solids. In systems with weak random point pinning, the existence of a low-temperature topologically ordered Bragg glass (BrG) phase with quasi-long-range translational order is now well established brg1 ; brg2 . A variety of fascinating “glassy” behavior has been experimentally observed tifr  near the first-order melting transition of the BrG phase in both conventional and high-$`T_c`$ superconductors. It has been suggested tifr ; menon  that these observations can be understood if it is assumed that the melting of the BrG phase occurs in two steps: the BrG first transforms into a “multidomain” glassy phase which melts into the usual vortex liquid at a slightly higher temperature.

In the presence of random columnar pinning, a “strong” Bose glass (BoG) phase nelson  without quasi-long-range translational order occurs at low temperatures if the concentration of pins is larger than that of vortex lines. In the opposite limit of dilute pins, one expects radz  a “weak” BoG phase at low temperatures which would melt into an interstitial liquid (IL) as the temperature is increased. In the IL phase, some of the vortices remain pinned at the strong pinning centers, while the other, interstitial ones form a liquid. A recent numerical study bbrgl  suggests that a topologically ordered BrG phase is also possible in such systems if the pin concentration is sufficiently small. It is also found experimentally, for both point expt1  and columnar expt2  pinning, that the melting of the solid phase is “broadened”: the local transition temperature, measured by a discontinuity of the local magnetization, is different in different regions of the sample.

Here we report results of a numerical study that provides insights and explanations for some of the observations described above. From minimization of an appropriate free energy functional, we find that the vortex system in an extremely anisotropic, layered, superconductor with a random dilute array of strong columnar pins (with both pins and magnetic field perpendicular to the layers) forms a BrG phase at low temperatures. As $`T`$ is increased, this phase undergoes a first order transition into a glassy phase which we identify as a polycrystalline BoG. This phase then transforms, at a slightly higher $`T`$, into the IL phase via a second, more strongly first order transition. We also show that the local transition temperatures, obtained from the temperature-dependence of a quantity that measures the degree of localization of the vortices in a small region of the sample, exhibit substantial spatial variation correlated with the local arrangement of the pinning centers.

The model and methods we use are similar to those in our earlier work cdotv  for a periodic array of columnar pins: the main difference is that the pin array here is taken to be random. Thus, we study a layered superconductor with vanishingly small Josephson interlayer coupling (vortices on different layers are coupled via the electromagnetic interaction only). In this limit, appropriate gautam  for extremely anisotropic Bi- and Tl-based high-$`T_c`$ materials, the energy of a system of “pancake” vortices residing on the superconducting layers may be written as a sum of anisotropic two-body interactions. We use the Ramakrishnan-Yussouff (RY) free energy functional ry . Since the potential produced by a set of straight columnar pins perpendicular to the layers is the same on every layer, $`\rho (𝐫)`$, the time-averaged local areal density of vortices at point $`𝐫`$ on a layer, must be the same on all layers. The free energy $`F`$ per layer may then be written as:

$`\beta (F[\rho ]F_0)`$ $`=`$ $`{\displaystyle d^2r\left[\rho (𝐫)\{\mathrm{ln}(\rho (𝐫))\mathrm{ln}(\rho _0)\}\delta \rho (𝐫)\right]}`$ (1)
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle d^2rd^2r^{}\stackrel{~}{C}(|𝐫𝐫^{}|)\delta \rho (𝐫)\delta \rho (𝐫^{})}`$
$`+`$ $`\beta {\displaystyle d^2rV_p(𝐫)\delta \rho (𝐫)}.`$
Here, $`\delta \rho (𝐫)\rho (𝐫)\rho _0`$, $`F_0`$ is the free energy of the uniform liquid of areal density $`\rho _0`$ (= $`B/\mathrm{\Phi }_0`$ where $`B`$ is the magnetic induction and $`\mathrm{\Phi }_0`$ the flux quantum), $`\beta =1/k_BT`$, $`V_p(𝐫)`$ is the pinning potential, and $`\stackrel{~}{C}(r)_nC(n,r)`$, where $`C(n,r)`$ is the direct pair correlation function of a layered liquid of pancake vortices ($`n`$ is the layer separation and $`r`$ is the separation in the layer plane). We use the results for $`C(n,r)`$ obtained gautam  from a hypernetted chain calculation. The RY functional yields cdotv ; gautam  a correct quantitative description of the melting transition in the absence of pinning. The potential $`V_0(r)`$ at $`𝐫`$ due to a pinning center at the origin is assumed to have the form $`V_0(r)=\alpha \mathrm{\Gamma }(1r^2/r_0^2)`$ for $`rr_0`$ and $`V_0(r)=0`$ if $`r>r_0`$. Here, $`\mathrm{\Gamma }\beta d\mathrm{\Phi }_0^2/8\pi ^2\lambda ^2(T)`$, $`d`$ is the layer spacing, $`\lambda (T)`$ is the penetration depth in the layer plane, $`r_0`$ is a range parameter and $`\alpha `$ is a strength parameter. The net pinning potential $`V_p(𝐫)`$ is the sum of the potentials due to $`N_p`$ randomly placed pinning centers. We use parameter values appropriate to BSCCO i.e. $`\lambda (T=0)=1500\AA `$ and $`d=15\AA `$, and assume a two-fluid $`T`$-dependence of $`\lambda (T)`$ with $`T_c(0)=85`$K. Defining $`a_0`$ (which we will use as our unit of length) via the relation $`\pi a_0^2\rho _0=1`$, we set $`r_0=0.1a_0`$ and $`\alpha =0.05`$. For these values, each pinning center traps one vortex cdotv  in the temperature range of interest.

We discretize space by defining density variables $`\{\rho _j\}`$ at the sites of a triangular grid of size $`(Nh)^2`$ with periodic boundary conditions. The grid spacing $`h`$ is taken to be $`a/16`$ where $`a1.988a_0`$ is the equilibrium spacing cdotv  of the pure vortex lattice at melting for the value of the magnetic field ($`B`$ = 2kG) used here. The $`N_p`$ pinning centers are randomly put on computational lattice sites. Local minima of the discretized $`F[\rho ]`$ Eq.(1), written as a function of the $`\{\rho _j\}`$, are then obtained numerically using a methodology quite similar to that in Ref. cdotv . We report here primarily results for $`N`$ = 1024 (corresponding, for the chosen value of $`h`$, to including $`N_v`$ = 4096 vortices in the calculation), and relative pin concentration (the ratio of the number of pins $`N_p`$ to $`N_v`$) $`c`$ = 1/64. A larger pin concentration, $`c`$ =1/32, and samples of size $`N`$ = 512 (1024 vortices) were also studied prep . Results for the glassy phases depend somewhat, see below, on the placement of the $`N_p`$ random pins. Over twenty different random pin configurations were studied and averages taken where appropriate.

Different local minima of the discretized $`F`$, corresponding to the several possible phases of the system, are reached when the numerical minimization is performed starting from different initial states cdotv . The free energies of these different minima at a given $`T`$ determine the (mean-field) phase diagram. A great advantage of our method in identifying minima representing different phases is that at each minimum the values of the full $`\{\rho _j\}`$ set are available. We can then calculate e.g. the structure factor $`S(𝐤)=|\rho (𝐤)|^2/N_v`$. We can further characterize the structure of a minimum by analyzing additional information prep  contained in the full set $`\{\rho _j\}`$. In particular, we have calculated the local peak densities, defined as the values of the density at local peaks. The density is considered to locally peak at a mesh point $`j`$ if $`\rho _j`$ is higher than those at all other mesh points within a distance $`a/2`$ from $`j`$. At low-temperature minima with localized vortices, these local density peaks lie where the vortices are localized: their number matches the number of vortices $`N_v`$. Thus, the positions of these peaks define a “vortex lattice”. To elucidate the degree of order in this “vortex lattice”, we have found it particularly useful to carry out a Voronoi construction, thereby determining the number of nearest neighbors of each vortex and hence the defect structure of the minimum. To address questions about orientational order, we have obtained the bond-orientational vortex correlation function $`g_6(r)`$, defined as the correlation function of the field $`\psi (𝐫)=_j\mathrm{exp}[6i\theta _j(𝐫)]/n_n`$, where $`\theta _j(𝐫)`$ is the angle that the bond connecting a vortex at $`𝐫`$ to its $`j`$th neighbor makes with a fixed reference direction, and $`n_n`$ is the number of neighbors of the vortex at $`𝐫`$.

We find three different kinds of local minima of $`F`$. The simplest kind is found by quenching the system, starting from uniform initial conditions, to temperatures somewhat higher than the equilibrium melting temperature of the pure vortex lattice ($`T_m^018.4`$cdotv ). The minimum can then be slowly cooled to below $`T_m^0`$ until it becomes unstable. This phase has local vortex densities close to the uniform liquid density everywhere except near the pinning centers, each of which traps one vortex. These minima clearly correspond to the interstitial liquid (IL) state with very small ($`5`$) peak values of $`S(𝐤)`$ (see panel (c) of Fig.1).

The second kind of local minima are obtained by quenching to temperatures below $`T_m^0`$ with initial conditions corresponding to a perfectly crystalline initial state (we use the crystalline state for which the pinning energy is minimum). They can then be cooled down, or warmed up to above $`T_m^0`$. This phase is nearly crystalline: the densities at local peaks turn out to be large (5-10 times $`\rho _0`$) nearly everywhere, except at the pins where they are much higher. There are also a few small regions of lower peak density, indicating weakly localized vortices. The Voronoi plots for such minima (see panel (a) in Fig. 1) clearly illustrate the defect structure. A pair of adjacent 5- and 7-coordinated sites, shown as (red) triangles and large (green) dots respectively, corresponds to a dislocation. These dislocations form tightly bound clusters in this case. These clusters are located near pinning sites, shown by black circles in panels (a) and (b) of Fig.1. The local peak densities near the defect clusters are lower, indicating weaker localization of the vortices. The structure factor plot, shown in panel (c), exhibits six sharp Bragg peaks for these minima. The “crystalline order parameter” extracted from the peak value of $`S(𝐤)`$ is large at low temperatures ($`0.55`$ at 17K for $`N_v=1024`$) and decreases very slowly with sample size (by $``$ 6% as the sample size is doubled). The vortex bond-orientational correlation function $`g_6(r)`$ for such minima (see Fig. 2, upper inset) saturates to a large value for large $`r`$. Thus, these minima exhibit all the characteristics brg1  of a BrG phase and we conclude that they can be so identified. However, our numerical study can not rule out the occurrence of unpaired dislocations at much longer length scales. If this happens, then these minima may correspond to a “hexatic glass” hexatic  phase. In any case, it is clear that these minima represent a phase that is distinct from the polycrystalline one described below.

The third kind of minima are obtained either by slowly cooling a liquid-like minimum to below the temperature where it becomes unstable, or by quenching with a uniform initial density to a temperature well below $`T_m^0`$. The Voronoi construction results for this case (panel (b) of Fig. 1) clearly show a polycrystalline structure with the dislocations lining up along grain boundaries, which lie mainly in regions without any pinning center. As a liquid-like initial state is cooled from relatively high $`T`$, the vortices arrange themselves in triangular crystalline patches around the pins. The orientation of a crystalline patch depends on the local pin arrangement. As the temperature is lowered further, misaligned patches join each other at grain boundaries to form a new minimum of this kind. The local peak density is substantially lower near the grain boundaries. As these minima are warmed upprep , the regions near the grain boundaries begin to “melt” before the other parts of the sample. As shown in panel (c) of Fig.1, the structure factor for these minima exhibits several (typically more than six) peaks of height much lower than that of the six Bragg peaks found for the BrG minimum. The function $`g_6(r)`$ for such minima goes to zero at large $`r`$ (see Fig. 2, upper inset). We conclude, therefore, by considering all the evidence, that such minima correspond to polycrystalline BoG states. While for any given pin configuration the IL and BrG states reached upon the minimization are, within numerical uncertainty, unique, different BoG type states can be reached using different starting states and cooling rates. This is characteristic of a glassy phase. When more than one BoG minima are found at a given $`T`$, we consider the one with the lowest free energy. The polycrystalline nature of the BoG minima is consistent with the results of experiments expt3  and simulations nandini  of the Bose glass phase in the dilute-pin regime.

Fig. 2 illustrates the main result of our study. There we show the temperature dependence of $`\beta (FF_0)`$ of the BrG, IL and BoG minima for the same pin configuration at $`c`$=1/64. The BrG minimum has the lowest free energy at low temperatures. Its free energy crosses that of the BoG minimum near 17.9K, indicating a first-order BrG-BoG transition at this temperature. The BoG phase then melts into the IL phase at a slightly higher temperature, near 18.3K, as indicated by the crossing of the free energies of these two minima. Thus, the melting of the low-temperature BrG phase with increasing temperature occurs in two steps, with a small region of BoG phase separating the BrG and IL phases. The grain boundaries in the BoG minima survive thermal cycling expt3  across the BoG-IL transition temperature, whereas similar thermal cycling produces grain boundaries in the BrG minima, indicating that the BoG phase is thermodynamically stable near the BoG-IL transition.

For samples with $`c=1/64`$, the average value of the temperature interval in which the BoG is the equilibrium state is 0.42 K. This width varies from sample to sample between less than 0.1 in one single case, to 1.2K. Thus, the two-step transition is a generic feature. The entropy jump at the BrG-BoG (lower) transition ($`0.1k_B`$ per vortex) is smaller than that at the BoG-IL transition ($`0.15k_B`$ per vortex). These values do not depend significantly on the system size and their sum is slightly smaller than the value ($`0.29k_B`$) at the single melting transition in the pure system. The size of these jumps makes it very unlikely that fluctuations would change the nature of the transitions in our 3D system. The value of the upper transition temperature is between 18.2K and 18.3K for all samples at $`c`$ = 1/64. These values are close to the first-order melting temperature of the pure system cdotv . The weak dependence of the transition temperatures on $`c`$ is consistent with experiments expt2 ; expt3 . Also, a first-order freezing transition of the IL to a polycrystalline solid has been observed expt3  in experiments on BSCCO with a small concentration of columnar pins. A narrow “two-phase” region found near the BrG melting transition in Ref.bbrgl  may correspond to an intermediate BoG phase. Alternatively, the sample size ($`100`$ vortex lines) in Ref. bbrgl  may be too small for the detection of an intermediate BoG phase with large crystalline domains.

The BrG-BoG transition occurs as a result of a competition between elastic and pinning parts of the free energy. The BrG minimum has a lower elastic (free) energy than the BoG minimum, but a higher pinning energy: the vortices adjust better to the pinning potential in the BoG minimum. The softening of the lattice near melting decreases the relative importance of the elastic component, thus causing a crossing of the two free energies.

We have also calculated a space-dependent “local transition temperature” by monitoring the temperature-dependence of $`\rho _{av}^p`$, the average of the local peak density in small regions containing $``$ 100 vortices. Vortices localized at pinning centers are not included in the calculation of $`\rho _{av}^p`$. Values of $`\rho _{av}^p`$ much larger than $`\rho _0`$ indicate a solid-like local structure, while values close to $`\rho _0`$ indicate liquid-like behavior. In Fig.2, lower inset, we show the $`T`$-dependence of $`\rho _{av}^p`$ for three regions centered at points A, B and C in panel (a) of Fig.1. The local transition temperature, defined as the temperature at which $`\rho _{av}^p`$ crosses $`3\rho _0`$ cdotv , is different in the three regions. The lowest local $`T_c`$ corresponds to the BrG-BoG transition and reflects the local melting near a grain boundary of the BoG minimum (point C). The highest local $`T_c`$ is higher than the BoG-IL transition temperature, reflecting solid-like local structure near a cluster of pinning sites in the IL minimum (point A). The range of variation of the local $`T_c`$’s is comparable to that found in experiments expt1 ; expt2 .

Thus, we have shown that a layered superconductor with a small concentration of columnar pins exhibits a two-step melting transition from a low-$`T`$ BrG phase to a high-$`T`$ IL phase via an intermediate BoG phase. A suitably defined local transition temperature exhibits spatial variations correlated with the local arrangement of pinning centers. Our results are consistent with experiment, and support the suggestion tifr ; menon  of similar behavior in systems with point pinning.


text:
# Cold atoms probe the magnetic field near a wire

## 1 Introduction

The ability to control cold atom clouds in microscopic magnetic traps  and waveguides  has created the new field of miniaturized atom optics . With the use of microstructured surfaces (atom chips) it becomes possible to control cold atoms on the $`\mu `$m length scale and to anticipate the construction of integrated atom interferometers . Ultimately there is the possibility of controlling the quantum coherences within arrays of individual atoms for use in quantum information processing . For these kinds of applications it is important to avoid fluctuating or inhomogeneous perturbations, which tend to destroy the quantum coherences.

Cold atom clouds cooled below a few $`\mu `$K have recently been used to probe the magnetic field fluctuations within 100 $`\mu `$m of a current-carrying wire. Two frequency domains have been studied. Audiofrequency fluctuations of the currents that form the microtraps can excite centre-of-mass vibrations of the atoms , causing the cloud to heat up. Radiofrequency noise in the wire currents can drive spin flips, which cause atoms to be ejected from the trap . Usually this noise is due to technical imperfections of the apparatus, as elucidated by Leanhardt et al. through a comparison of magnetic and optical traps near a surface . However, there is also a fundamental component of the magnetic field noise in the near-field of the wire, which is due to thermal fluctuations of the current and has recently been measured in our laboratory .

The dc behaviour of the field is less well understood. As the atom cloud is brought close to the main current-carrying wire of the trap it breaks into fragments along the length of the wire , indicating an unexpected variation in the trapping potential. Kraft et al. have recently shown that this is due to the presence of a magnetic field component $`\mathrm{\Delta }B_z`$ parallel to the wire that varies along the length of the wire. In this paper we use cold atoms held in a microtrap close to a wire to measure $`\mathrm{\Delta }B_z`$ and to explore how it varies in the vicinity of the wire.

## 2 Preparation of cold atoms

The arrangement of wires used to form our microtrap is shown in figure 1. The main wire is a 500 $`\mu `$m diameter guide wire along the z-direction. In cross section it has a $`370\mu `$m diameter copper core, surrounded by an aluminium layer $`55\mu `$m thick with a $`10\mu `$m thick ceramic outer coating. This wire is glued by high-vacuum epoxy (Bylapox 7285) into the $`200\mu `$m-deep channel formed by a glass substrate and two glass cover slips. Below the guide wire there are four transverse wires, each $`800\mu `$m in diameter. The cover slips are coated with 60 nm of gold so that they reflect 780 nm light. In order to load cold atoms into the microtrap we first collect 87Rb atoms using a magneto-optical trap whose beams are reflected from the gold surface. This MOT collects $`1\times \mathrm{\hspace{0.17em}10}^8`$ atoms at a height of 4 mm above the surface and cools them to $`50\mu `$K. The MOT is pulled down to a height of 1.3 mm by passing a current of 3.2 A through the guide wire and adding a uniform magnetic field $`B_x^{bias}`$ of 6 G along the x-direction. This compresses the cloud into a cylindrical shape and increases the phase space density of the atoms to $`2\times 10^6`$ .

The light and the anti-Helmholtz coils of the MOT are then switched off and the atoms are optically pumped into the $`|F,m=|2,2`$ state. We collect $`2\times \mathrm{\hspace{0.17em}10}^7`$ of these atoms in the magnetic guide formed by the guide wire (8 A along z) and the transverse bias field $`B_x^{bias}`$ (10 G along x). Axial confinement is provided by the inner transverse wires (15 A each along -x ) and the outer transverse wires (15 A each along x ). The field at the centre of this trap is partly cancelled by an axial bias field $`B_z^{bias}`$ ($`6`$G along z). Next, the trap is adiabatically compressed over 0.5 s by increasing $`B_x^{bias}`$ and $`B_z^{bias}`$ to $`29`$G and $`11`$G respectively, and reducing the guide current to $`6.9`$A. This brings the trap to a height of 225 $`\mu `$m above the top of wire and reduces the net field along $`z`$ at the centre of the trap to $`1`$G. In this way, the radial and axial trap frequencies are raised to 840 Hz and 26 Hz. The elastic collision rate is now $`54`$ s-1, which is high enough for forced rf evaporative cooling to be efficient. We sweep the rf frequency logarithmically over 6 s from 13 MHz to a final frequency near 2.8 MHz, cooling the cloud down to approximately 6 $`\mu `$K. This temperature is a good choice because it gives the atoms a thermal kinetic energy comparable with the potential created by $`\mathrm{\Delta }B_z`$ that we wish to probe. The atoms are brought closer to the surface by a smooth reduction of the current flowing in the guide wire during the last $`1`$s of evaporation. The arrival of the atoms at the desired height coincides with the end of the evaporation ramp.

The cloud is viewed using a ccd camera to record the absorption of a resonant probe beam propagating along the x-direction. Two images are recorded, one with the atoms present and one without, and we compute the logarithm of their ratio. Since the cloud is optically thin, this procedure yields an image of the column density of the atoms, viewed in the $`yz`$ plane. Integrating this over $`y`$, gives the probability distribution of atoms along the z-direction. Let us call it $`\rho (z)`$. Assuming that the gas is in thermal equilibrium at temperature $`T`$, $`\rho (z)`$ is related to the confining potential $`U(z)`$ by the Boltzmann factor $`\mathrm{exp}(U(z)/kT)`$. In the absence of $`\mathrm{\Delta }B_z`$, the potential $`U(z)`$ is harmonic with an axial frequency of 26-27 Hz over the range of heights studied. Fitting a parabola to $`\mathrm{ln}\rho (z)`$ then determines the temperature $`T`$. For example, the top curve in figure 2 taken 97 $`\mu `$m above the surface of the wire gives $`T=5.8\mu `$K.

## 3 Fragmentation

When the cloud is brought closer to the surface, the distribution $`\rho (z)`$ develops additional structure, as can be seen in figure 2. Similar behaviour has been reported by other groups studying atoms near copper wires . This indicates that the atoms experience some other interaction near the surface in addition to the expected simple harmonic potential of the trap. Kraft et al. have recently shown that this is due to a magnetic field component $`\mathrm{\Delta }B_z`$ parallel to the wire . Where this field adds to (reduces) the axial bias field of the magnetic trap there is an increase (decrease) of $`\mu _B\mathrm{\Delta }B_z`$ in the potential energy of the trapped atoms. To obtain $`\mathrm{\Delta }B_z`$ for each curve in figure 2, we first derive the temperature $`T`$ from the parabola fitted to -ln$`\rho (z)`$. This parabola is chosen to pass through the centre of the $`\mathrm{\Delta }B_z`$ oscillations. For the 8 curves shown the temperature is nearly constant, being (5.8, 6.1, 6.4, 7.0, 7.8, 7.8, 7.8, 7.8) $`\mu `$K, with the variations arising from slight differences in the final stage of the evaporation trajectories. The difference between $`\mathrm{ln}\rho (z)`$ and the fitted harmonic potential is the potential of interest $`\mathrm{\Delta }B_z/kT`$, and since we know $`T`$ from the harmonic fit, we obtain $`\mathrm{\Delta }B_z`$. In this way we have measured how $`\mathrm{\Delta }B_z`$ varies along $`z`$ at the 8 different heights $`h`$ of the cloud above the wire.

The height of the cloud above the surface is set by adjusting the current in the guide wire, keeping the bias field $`B_x^{bias}`$ fixed at 29 G. For the largest and smallest heights presented in figure 2 (97 $`\mu `$m and 7.2 $`\mu `$m) the guide wire currents are 5 A and 3.7 A respectively. $`\mathrm{\Delta }B_z`$ is proportional to the current in the guide wire , so to compare $`\mathrm{\Delta }B_z`$ at different heights we scale the result at each height to the current used. Figure 3 is a map of $`\mathrm{\Delta }B_z(y,z)`$ scaled to a fixed 3.7 A in the guide wire. The surface shown interpolates between the 8 heights where measurements were made.

We see $`\mathrm{\Delta }B_z`$ undergoing 2 full periods of oscillation along the wire. The phase of the oscillation is fixed with respect to the wire and did not change over several months of experiments. The field has many more oscillations along the wire, with an average wavelength $`\lambda `$ of $`230\pm 10\mu `$m, but in this particular experiment the atoms are confined to these two periods by the trap potential. It would seem that the current, or perhaps the electron spin , follows an oscillatory or helical path along the wire with wavelength $`\lambda `$. This could be due to some fundamental physics, or perhaps to a regular arrangement of impurities or defects in the wire. In any such case, the decay of the field component $`\mathrm{\Delta }B_z`$ with distance $`y`$ should be well approximated by the modified Bessel function $`K_1(ky)\mathrm{exp}(ky)/\sqrt{ky}`$, provided the distance $`y`$ from the current to the atoms is much greater than the transverse excursions of the current. The quantity $`k`$ is $`2\pi /\lambda `$. In order to test this idea we have measured the amplitude of the $`\mathrm{\Delta }B_z`$ oscillations at each height $`h`$ above the wire and these points are plotted in figure 4. The solid line is a least-squares fit to the function $`aK_1(2\pi (h+\delta )/\lambda )`$, in which $`a,\delta `$, and $`\lambda `$ were allowed to vary freely as fitting parameters. The best fit has a $`\chi ^2`$ of 4.8 for 5 degrees of freedom and gives the results $`\lambda =217\pm 10\mu `$m and $`\delta =251\pm 12\mu `$m. The close agreement between this value of $`\lambda `$ and the value of $`\lambda =230\pm 10\mu `$m obtained directly from the images strongly suggests that our model is correct. The value of $`\delta `$ is equal to the wire radius indicating that the centre of the oscillating current coincides with the centre of the wire.

The dashed curve in figure 4 shows the usual magnetic field produced by 3.7 A in a wire, which varies inversely with the distance from the centre of the wire. Expressed as a fraction of this field, the amplitude of the $`\mathrm{\Delta }B_z`$ oscillations ranges from $`3\times 10^3`$ to $`3\times 10^4`$ over the heights studied. This ratio provides information about the transverse excursion of the effective current. For example, if we suppose that the current is helical with an effective radius of $`r`$, the ratio of field components has the value $`k^2ryK_1(ky)`$, which implies that $`r=50\mu `$m. A similarly large amplitude is found for other possible motions (e.g. planar or random walk with characteristic length scale). The Lorentz forces due to the applied magnetic fields do not generate enough transverse current to account for this. Nor do the variations in the shape and centre of the wire. Microscopic structure within the wire was revealed by polishing and etching a longitudinal section of the wire following standard metallographic methods. In the copper we found grains of $`10\mu `$m and in the aluminium we found small defects typically spaced by $`5\mu `$m. Neither of these provides a plausible explanation for the transverse current. We therefore remain unable to propose a mechanism for the observed variation of the current.

The exponential character of the decay means that $`\mathrm{\Delta }B_z`$ decreases too rapidly with distance to be described by a power law. The Tübingen group, who proposed a power law ansatz , have kindly provided us with their data  taken above a 90 $`\mu `$m diameter copper wire and we find that it is also well described by the $`K_1`$ Bessel function.

## 4 Conclusions

We have investigated the behaviour of a cloud of cold trapped atoms brought close to the surface of a current-carrying wire. As the atoms approach the wire, the cloud breaks into fragments due to the presence of an anomalous magnetic field parallel to the wire. Analysis of the density distribution of the atoms has allowed us to make a quantitative map of this anomalous field. Along the wire it oscillates with a period of 230$`\pm 10\mu `$m and with increasing radius it decays according to a modified Bessel function $`K_1`$. We point out that this behaviour is characteristic of the field produced by a periodic transverse current distribution in the wire. We have found a characteristic decay length of $`217\pm 10\mu `$m, in good agreement with the observed period and we have shown that the anomalous current is centred on the middle of the wire. Despite some effort, it remains unclear why the current should move in this way.

We are indebted to Alan Butler and John Knight for expert technical assistance. We also thank József Fortágh for helpful discussions and for providing us with data from the Tübingen experiment. This work was supported by the UK Engineering and Physical Sciences Research Council, and by the FASTNET and Cold Quantum Gases networks of the European Union.

## References


text:
# KANAZAWA-03-17 ITEP-LAT-2003-12 Entropy of monopoles from percolating cluster in quenched SU(2) QCDPresented by K. I. at Lattice’03.

## 1 INTRODUCTION

The dual superconductor picture  of the QCD vacuum is based on the existence of Abelian monopoles which appear naturally in an Abelian gauge. There are various indications that the monopoles are responsible for the color confinement (for a review, see Ref. ). One of the most important results is the observation of the monopole condensate in the confinement phase .

The monopole trajectories form clusters. Typically, each lattice configuration contains many finite-sized (ultraviolet) clusters and one large percolating (infrared) cluster . The IR cluster – which occupies the whole lattice – represents the monopole condensate . The tension of the confining string gets a dominant contribution from the IR cluster . In the deconfinement phase the IR cluster disappears , as expected.

The UV and IR monopole clusters were studied previously in Refs. . Below we investigate the action and length distribution of the infrared monopole cluster for various lattice sizes and scales $`b`$ at which the magnetic charge is defined. Using the action and length distribution we investigate the entropy of the IR clusters.

## 2 MODEL

In our simulations we use the Wilson action, $`S(U)=\frac{\beta }{2}\mathrm{Tr}U_P`$. We work in the Maximal Abelian (MA) gauge . The Abelian gauge field, $`\theta `$, is determined as a phase the diagonal component of the $`SU(2)`$ link variable, $`\theta _\mu (s)=\mathrm{arg}U_\mu ^{11}(s)`$. The Abelian field strength, $`\theta _{\mu \nu }(s)(4\pi ,4\pi )`$, is decomposed into two parts, $`\theta _{\mu \nu }(s)=\overline{\theta }_{\mu \nu }(s)+2\pi m_{\mu \nu }(s)`$, where $`\overline{\theta }_{\mu \nu }(s)[\pi ,\pi )`$ is interpreted as the electromagnetic flux through the plaquette and $`m_{\mu \nu }(s)`$ can be regarded as a number of the Dirac strings piercing the plaquette.

The elementary monopole currents are defined as follows , $`k_\mu (s)=1/2ϵ_{\mu \nu \rho \sigma }_\nu m_{\rho \sigma }(s+\widehat{\mu })`$, where $``$ is the forward lattice derivative. To study the monopole charges at various scales we use the $`n^3`$ extended monopole construction ,

$`k_\mu ^{(n)}(s)={\displaystyle \underset{i,j,l=0}{\overset{n1}{}}}k_\mu (ns+(n1)\widehat{\mu }+i\widehat{\nu }+j\widehat{\rho }+l\widehat{\sigma }).`$
The extended monopole is defined on a sublattice with the spacing $`b=na`$, where $`a`$ is the spacing of the original lattice. Both elementary and extended monopole charges are conserved and quantized.

We used the standard Monte–Carlo procedure to generate 1000-3000 configurations of the gauge field for each value of $`\beta =2.12.6`$. To fix the MA gauge we used either the usual iterative algorithm (for $`24^4`$ and $`48^4`$ lattices) or the Simulated Annealing method  with five Gribov copies (for $`32^4`$ lattice).

## 3 MONOPOLE ACTION

To get the monopole action we integrate out all degrees of freedom but the monopole ones. Following Ref.  we generate the SU(2) configurations, fix the MA gauge, then get the configurations of the IR monopole clusters as it was described in the previous section. Then we use the inverse Monte–Carlo method to get the monopole action.

The monopole action can be represented  as a sum of the $`n`$–point ($`n2`$) operators $`S_i`$:

$$S_{\mathrm{mon}}[k]=\underset{i}{}f_iS_i[k],$$

where $`f_i`$ are the coupling constants. We adopt only the two–point interactions in the monopole action ($`i.e.`$ interactions of the form $`S_ik_\mu (s)k_\mu ^{}(s^{})`$). A detailed description of the $`S_i[k]`$ interactions can be found in Ref. .

One can find that the monopole action is proportional with a good accuracy to the length of the monopole loop, $`S_{\mathrm{mon}}=fL+const`$. The dominant term in the monopole action corresponds to the most local self-interaction of the monopole currents, $`S_1[k]=_{s,\mu }k_\mu ^2(s)`$. We compare the parameter $`f_1`$ for the monopole action of the IR cluster and for the action associated with the whole monopole ensemble in Figure 1.

One can notice that the coupling constant $`f_1`$ shows scaling for large $`b=na`$ in agreement with Ref. . This coupling is independent of the lattice volume and for large blocking scales $`b`$ the type of the ensemble (the IR cluster or the whole ensemble) is not essential for determination of $`f_1`$. However, at small $`b`$ values, $`b\sqrt{\sigma }0.5`$, the type of the lattice ensemble becomes important.

## 4 LENGTH DISTRIBUTION

The distribution of the ultraviolet clusters was studied both numerically  and analytically . The distribution can be described by a power law $`D_{UV}L^\tau `$, where the power $`\tau `$ is very close to 3, Ref. . This behaviour indicates that the monopoles in UV clusters are randomly walking objects . In our simulations we are concentrated on the IR cluster because the IR cluster is important for the confinement of quarks.

The length distribution function, $`D(L)`$, is proportional to the weight with which the particular trajectory of the length $`L`$ contributes to the partition function. From the previous section it is clear that monopole action contributes to $`D(L)`$ in a form of an exponential factor, $`e^{fL}`$. The entropy of the monopole trajectory also contributes to the monopole length distribution. The entropy contribution is proportional to $`\mu ^L`$ (with $`\mu >0`$) for sufficiently large monopole lengths, $`L`$. Thus the distribution of the monopole trajectories in infinite volume must be described by the function

$`D_{inf}^{IR}(L)\mu ^Le^{fL}=e^{\gamma L},\gamma =\mathrm{ln}\mu f.`$ (1)
In this equation we neglect a power-law prefactor which is essential for the distribution of the ultraviolet clusters. In the finite volume there appears a finite–volume suppression factor and the total distribution can be described as a Gaussian :

$`D^{IR}(L)=const.\mathrm{exp}\{\alpha (b,V)L^2+\gamma (b)L\}.`$ (2)
The peak of this distribution, $`L_{max}=\gamma (b)/2\alpha (b,V)`$, is expected to be proportional to the volume of the system, $`V`$, to insure the finiteness of the IR monopole density, $`\rho _{IR}=L_{max}/V`$, in the thermodynamic limit, $`V\mathrm{}`$. Thus we conclude that $`\alpha (b,V)=A(b)/V`$, where $`A(b)`$ is a certain function. Therefore in the thermodynamic limit the parameter $`\alpha `$ vanishes and the distribution (2) is reduced to Eq.(1), as expected.

To get the monopole entropy we need to know the parameter $`\gamma `$, Eq. (1). We fit the monopole loops distributions by the function (2) and use the bootstrap method to estimate the statistical errors. We obtain that the parameter $`\gamma `$ — shown in Figure 2 — scales with $`b=na`$ and is independent on the volume of the system.

## 5 MONOPOLE ENTROPY

The knowledge of the distribution and the monopole action allows us to determine the entropy of the monopole currents. According to Eq. (1), $`\mu =e^{\gamma +f}`$. If the monopoles make a simple random walk on a $`4D`$ hypercubic lattice then we get $`\mu =7`$ since there are seven choices at each site for the monopole current to go further.

The entropy factor $`\mu `$ is shown in Figure 3. It scales with $`b`$ and is independent of the volume of the lattice. Note that in a small $`b`$ region $`\mu >7`$ because the monopole action can not be reliably described by the quadratic terms only .

At large $`b`$ the entropy factor is smaller than seven. We have fitted the entropy by a function:

$`\mu ^{\mathrm{fit}}=\mu _{\mathrm{}}+C\mu ^q,`$ (3)
where $`\mu _{\mathrm{}}`$, $`C`$ and $`q`$ are the fitting parameters. The fit gives $`\mu _{\mathrm{}}=1.6(4)`$, $`C=1.7(5)`$ and $`q=1.2(2)`$. Fixing $`q=1`$ in Eq. (3) we get $`\mu _{\mathrm{}}=1.15(25)`$ and $`C=2.2(1)`$. The corresponding best fit curves are shown in Figure 3.

The fact that the asymptotic value of the entropy is very close to the unity in the large $`b`$ limit may have a simple explanation. A monopole with a large blocking size $`b`$ behaves as a classical object and its motion is no more a simple random walk. The predominant motion of the large–$`b`$ monopole is close to a straight line.