text:
# Time decay of the remanent magnetization in the ±𝐽 spin glass model at 𝑇=0

## I INTRODUCTION

Spin glasses are systems that, at temperatures below so-called glass transition temperature $`T_g`$, find themselves in states with frozen disorder, i. e., no long-range patternlike order typical for ordered magnets is present. Their properties are determined by competing ferromagnetic and antiferromagnetic exchange interactions that are randomly distributed in the system. Both the competition among the different interactions between the magnetic moments and their random distribution all over a given system are likely to contribute significantly to such unusual glassy behavior. From the theoretical point of view, they can be discussed in terms of the coarse-grained free energy. Namely, due to accidental degeneracy present in such systems, below $`T_g`$ their free energy landscape becomes extremely rough, with many local minima corresponding to the same macroscopically observed properties but with entirely different microscopic states in the system phase space. The minima (valleys) are separated from each other by some energy barriers, and once a system finds itself in one of them, it might take a lot of time on laboratory time scales to transit to the others. Thus, the observed properties of spin glasses may only correspond to those of one single valley in which the system happens to be, and as a result, ergodicity is practically broken . In other words, spin glasses can be seen as systems whose dynamics at low temperatures is extremally slow and whose properties measured in real experiments always correspond to situations out of equilibrium. Very good evidence for this nonstationary dynamics is the response of the system ac susceptibility to an oscillating field, i.e., its dependence both on time and frequency. Another example is the slow decay of remanent magnetizations with time . The thermoremanent magnetization (TRM) is measured by cooling the sample in a nonzero magnetic field $`H`$ from above the glass-transition temperature $`T_g`$ to a temperature $`T`$ below it and then swiching off the field. The isothermal remanent magnetization (IRM) is measured by zero-field cooling of the sample in the same way as before, then turning on the field, and subsequently turning it off. In addition, in both cases the experimental results show that for small fields applied to the sample, the remanent magnetization is strongly affected by the so-called waiting time at which the sample is kept at constant temperature before the external field is changed . According to the decay of both the remanent magnetization and the energy, the experimental results depend on the observation times after the field is switched off in all measurements. The above nonstationary dynamics has been described by a fair variety of functions. The most important ones include power law, logarithmic, stretched exponential, and others, and the question of judging which is the most universal is still far from decided . Additionally, as some experiments indicate, the remanent magnetization decays so slowly with time that some nonzero remanence is still observed over macroscopic time scales, particularly at very low temperatures . Relaxation time measurments in CeNi0.8Cu0.2 below the spin glass temperature 6 K show that the decay time increases drastically with a distinct tendency to a state with nonzero magnetization, which is higher the lower the temperature is . The theoretical background of such a property is rather unclear .

One of the simplest theoretical models of spin glasses is the Edwards-Anderson (EA) model. Ising spins are located at each site of a lattice with randomly distributed ferromagnetic and antiferromagnetic interactions between the nearest neighbors. Such a model reveals most of the crucial features typical of real spin glasses including relaxation phenomena. Both the early papers  and the newer ones  on two- (2D) and three-dimensional (3D) models with a Gaussian distribution of bonds confirm that in a wide range of temperatures, a remanent magnetization occurs that decays very slowly with time $`t`$ according to a power law $`m(t)t^\alpha `$ with $`\alpha (T)T`$. It it also known that in case the couplings among spins may take only discrete values ($`\pm J`$ models) at sufficiently low temperatures, relaxation properties of such glasses are entirely different because of the existence of energy gaps in their energy spectra . As a result, for models with bimodal distribution of interactions at temperatures well below $`T_g`$, a simple linear dependence of $`\alpha `$ on temperature $`T`$ is not satisfied, and the functional form of the remanent magnetization decay is still rather far from being established. It has been mentioned that the function $`\alpha (T)`$ might go to zero faster than linear at low temperatures . It has also been suggested that at finite small temperatures, $`m(t)`$ should coincide with $`m(t)`$ at zero temperature for some region of time, and that the relaxation causes such a system to remain trapped in one metastable state with a finite remanent magnetization.

In this paper, we revisit the remanent magnetization decay of 2D $`\pm J`$ EA spin glass because we are interested in its extremely low-temperature relaxation properties where the power law with $`\alpha (T)T`$ breaks down. We consider the limiting case of zero-temperature behavior in order to find out whether the remanence phenomena observed at such conditions could reflect the low-temperature properties of discrete systems at least qualitatively. We would also like to check out whether or not results obtained at $`T=0`$ could be treated as continous extrapolations of those at low but finite temperatures, which, according to similar research done on SK model , seems rather unlikely. Actually, in  it is shown that even the zero-temperature dynamics provides a decay of the magnetization that can be fitted by a power law with a constant exponent $`\alpha `$. We carry out simulations using the zero-temperature dynamics, which has been successfully applied to various spin lattices including the persistence probability in weakly disordered Ising model , hysteresis in the random-field Ising model on the Bethe lattice , the question of avalanches in spin systems , and others . We investigate the system relaxation process towards low-energy states by plotting both energy and magnetization versus the zero-temperature Monte Carlo steps (MCS) per spin, which are treated as ”time units” $`t`$. We discuss the influence of frustration on the nonequilibrium time properties by comparing them with corresponding properties for an unfrustrated system with pure antiferromagnetic interactions. Finally, we discuss the obtained results in terms of the system random walk on its configuration space.

## II MODEL AND SIMULATION

We use the Edward-Anderson spin glass model with a random and uniform distribution of discrete interactions $`J_{ij}=\pm 1`$ between the nearest neighbors all over a 2D square lattice with $`N`$ sites and with periodic boundary conditions. The Hamilton function of such a system in an external magnetic field $`B`$ is:

$$H=\underset{i<j}{}J_{ij}S_iS_jB\underset{i}{}S_i,$$
(1)
where $`S_i,S_j=\pm 1`$ (up/down) are Ising spins, and the sum in the first term on the right-hand side runs over the nearest-neighboring lattice sites. The samples are prepared in such a way that the fraction of antiferromagnetic bonds is $`p=0.5`$. We study how the remanent magnetization decays with time in zero magnetic field starting from the system saturation state. The magnetization per spin is given by

$$m(t)=(N_+(t)N_{}(t))/N,$$
(2)
where $`N_+`$, $`N_{}`$ denote the number of spins up and down, respectively. The system relaxation process is simulated by applying a version of the zero-temperature Metropolis algorithm given by the following steps, cf. :
(i) Consider a sample in its saturation state (all spins up).
(ii) Pick a spin at random.
(iii) Flip it only if this process does not increase the energy.
(iv) Repeat steps (ii) and (iii) $`N`$ times. So, one time unit is defined as one MC step per spin.
(v) Record magnetization and energy as functions of these time units.
(vi) Start again with step (ii) to find results for the following time step.
Thus, the above algorithm determines the system random walk on the configuration space in the direction of dropping energy. It should be mentioned that a different algorithm based on the Glauber dynamics is also used in the literature, cf. . Then for $`T=0`$ step (iii) of the above mentioned procedure is replaced by the following one: Flip it if this process decreases the energy and flip it with the probability $`1/2`$ if this process does not change the energy.

## III RESULTS

First we consider a special system of the size $`N=70\times 70`$, for which the exact ground-state energy is calculated using a branch-and-cut algorithm by De Simone et al. . Its energy and magnetization decays are compared with those of an unfrustrated one of the same size with solely antiferromagnetic interactions between the nearest neighbors. No qualitative differences in the energy relaxation can be seen, except for the fact that the ground-state energies of both systems are different (Fig. 1). On the other hand, with respect to the remanent magnetization decay, both systems are entirely different and the influence of randomness and frustration on the phenomena can be observed. While the unfrustrated system moves extremely fast to the low-energy unmagnetized states, the frustrated one relaxes more slowly and the isoenergetic states it finally explores still possess nonzero magnetization (Fig. 2). In that area, the remanent magnetization fluctuates around a constant value, which refers to simply flipping idle spins.

In order to describe the system dynamics more quantitatively, we have considered $`10`$ samples of the size $`N=50\times 50`$ with different distribution of bonds. Typically, we have performed around 2000 MC relaxation steps. In Fig. 3, we show the remanent magnetization $`m(t)`$ averaged over the samples and over $`30`$ independent runs for each. The best fit for the dependence of the magnetization on time (for $`t>0`$) is obtained by

$$m(t)=m_{\mathrm{}}+\left(\frac{t}{a_0}\right)^{a_1},a_1<0.$$
(3)
The first term $`m_{\mathrm{}}`$ corresponds to the average remanent magnetization of the system trapped in a subspace of low-lying isoenergetic states after a long time of decay, whereas the second term corresponds to the nonstationary dynamics. The parameters found are $`m_{\mathrm{}}=0.059`$, $`a_0=0.23`$, $`a_1=0.83`$. The insert of Fig. 3 shows $`\mathrm{ln}[m(t)m_{\mathrm{}}]`$ versus $`\mathrm{ln}(t)`$ with $`m_{\mathrm{}}`$ taken from the power-law fit (3). The observed dispersion of the computed points for large values of time $`t`$ is due to their small fluctuations around $`m_{\mathrm{}}`$. We have also tested some other functions with four free parameters to adjust the data. Among them, the stretched exponential law $`m\left(t\right)=m_{\mathrm{}}+a_0\mathrm{exp}\left[\left(t/a_1\right)^{a_2}\right]`$ (for $`t0`$), with $`m_{\mathrm{}}=0.060`$, $`a_0=0.94`$, $`a_1=0.43`$, and $`a_2=0.35`$, has proved reasonably close to it, however a bit worse than the power-law one.

## IV Summary

We have investigated the remanent magnetization decay of a 2D $`\pm J`$ EA spin glass model at zero temperature. All the samples were initially at their saturation states and then allowed to relax towards states with lower energy using the zero-temperature Metropolis algorithm. The observed remanent magnetization decay in the frustrated system was much slower than in the unfrustrated one, whose remanent magnetization decreased rapidly to zero. After a number of MC steps, the walk of the random system on its configuration space practically became limited to subspaces of magnetized isoenergetic states. At that region, on the average, the remanent magnetization fluctuated without further decrease. Our calculations suggest that the time decay of the the remanent magnetization can be very well represented by a simple power-law formula with three fitting parameters. From a microscopic point of view, the nonexponential relaxation phenomena in spin glasses can be discussed in terms of random diffusion on the available configuration space. It is suggested by simulations that stretched exponential relaxation behavior in glassy systems appears with the exponent $`a_2`$ going to $`1/3`$ when approaching a percolation transition in the configuration space, which is a multidimensional hypercube . The same form of the magnetization decay in a $`2D`$ ferromagnetic Ising model has been reported and some dependence of the relaxation phenomena on the system dimension has been found as well . In the so-called trap model, the system evolves among various traps with random ’trapping times’. The traps are separated from each other by energy barriers that can be crossed by thermal excitations . This kind of approach has been successfully used to study low-temperature aging of a system consisting of configurations with random energies . However, since the dynamics of ours is athermal, there is no such barrier crossing in it and a purely entropic interpretation might be adapted instead . Although for simulations at finite temperatures a mixture of ’energetic’ and ’entropic’ barriers is likely to contribute to the phenomena of slow dynamics, one of the effects vanishes with respect to the zero temperature simulations. The observed slowing down of the system relaxation could be qualitatively understood with the help of so-called entropic traps themselves. While the point in the configuration space randomly searches for available paths leading to states with lower energy, their number decreases rapidly with time. As a result, we find that the lower the energy of the states is that the system has reached the longer is the time period needed to leave them. The phenomena is also very consistent with the Markov theory . It says that every finite Markov chain contains at least one absorbing set in which the system remains forever after falling in it (a good example of an absorbing set is the closure of one of the ground states). Moreover, the theory states that the probability of passing to one of them with ’time’ going to infinity tends to one. That is, after a large number of MC steps, the available configuration space of the system becomes one of the absorbing sets, where the energy is kept constant and the magnetization fluctuates. A finite value of remanent magnetization $`m_{\mathrm{}}`$ indicates that the absorbing set contains excited states rather than ground states. This is consistent with recent results on magnetic hysteresis at zero temperature . All metastable states connected by one-spin flips without raising energy form a local ensemble of metastable states. When all these states belong to a relatively large value of $`m`$, the one-spin interconnection with the local ensemble of ground states can be prohibited.

Moreover, our simulations confirm that the observed tendency of the exponent $`\alpha (T)`$ \[in the formula $`m(t)t^\alpha `$\] to go to zero with dropping temperature seems to be in direct contradiction to the remanent magnetization decay at exactly $`T=0`$, for which $`\alpha `$ is finite. This also seems to indicate that the violations of the power law at very low temperatures might be due to the characteristic time $`\tau `$ divergence  rather than to the function $`\alpha (T)`$ going to zero faster than linear. Moreover, we think that another explanation could be suggested as well. Namely, it has been experimentally found that the remanent magnetization $`m_{\mathrm{}}`$ decreases with increasing temperature . Below a certain value of temperature $`T`$, the discrete structure of the system energy spectrum is likely to affect its properties significantly and $`m_{\mathrm{}}`$ may become relevant. That means that at very low temperatures it should be contained in the relaxation law, and a plot $`\mathrm{ln}[m(t)m_{\mathrm{}}]`$ versus $`\mathrm{ln}(t)`$ (instead of $`\mathrm{ln}[m(t)]`$ versus $`\mathrm{ln}(t)`$) could still remain linear.

###### Acknowledgements.

Support from Graduiertenkolleg ”Struktur- und Korrelationseffekte in Festkörpern” is gratefully acknowledged. We thank D. Stauffer, J. Weißbarth, and M. Nogala for helpful discussions, and the authors of Ref.  for allowance to use their computer data.


format: text

text:
# 1 Introduction

## 1 Introduction

Carbon has various meta-stable polymorph, such as diamond, graphite, fullerene, or nanotube and so on. The structural transition between these polymorph structures, often induced by high pressure, has been intensively studied for years. For example, diamond is obtained from graphite by static or shock-wave compression experiments. The diamond obtained in such ways shows several distinct phases, such as, cubic diamond(CD) and hexagonal diamond(HD), which have different stacking and networking structures of sp3 tetrahedrons. In many cases, compressed graphite directly turns into cubic diamond with the aid of catalyst, accompanied with the change in stacking from so-called AB type into ABC one. However, HD phase has also been synthesized by static compression from well-crystallized graphite or shock-compression, in which AB staking is kept after quenching. However, in the compression of polycrystalline graphite, no hexagonal phase is obtained. It is also reported that the transformation to HD occurs under quasi-hydrostatic pressure, but only CD phase is quenched. Under certain conditions, ”n-diamond” phase also forms under compression . There are several theoretical studies for graphite-diamond transition from first principles. Scandolo et al. had simulated the transformation from graphite to diamond. Their simulation has created both CD and HD structure. Tateyama et al. investigate the activation barriers and the transition paths in the transition to cubic or hexagonal diamond from graphite and discussed on the origin of diamond polymorphs, which is attributed to the sliding movement of graphite layers under compression. However, from the computational limit, the ab-initio simulations use at most tens of atoms and they may not be sufficient for describing some aspect of diamond transition, since such a crystal growth process may be subject to large scale structural fluctuation. The small size simulation is likely to ignore that the graphite-diamond transition shall essentially be a cooperative phenomena caused by local re-combination of bonding structure. In addition, the realistic mechanism in choosing HD or CD structure under compression has not been clarified enough by those ab-initio studies. In the present work we adopt a tight-binding simulation using hundreds of atoms in order to study possible effects in the transition, such as fluctuation, disorder, or catalyst. At first we simulate the cooperative process of the graphite-cubic-diamond transition. Next, the simulation for catalyst effect in the transition to high-pressure polymorphs of diamond is given. There a model for catalyst-carbon combination is included into the TB Hamiltonian for carbon bonding structure. We propose a catalyst-aided-mechanism in the transition to carbon polymorphs under pressure.

## 2 Methodology

In the simulation of pressure-induced structural transition, we adopt a constant-pressure scheme proposed by Wentzcovitch et al.. In the evaluation of total energy and force in molecular dynamics, we use tight binding Hamiltonian for hydro-carbon system proposed by Winn et al.. Throughout the present work, the time interval of one MD step is set to be 1fs.

In the simulation, we start with graphite structure at 0Pa. The initial stacking of graphite is set to be AB type. In case of the simulation for catalyst-aided-effect, the catalysts are randomly distributed between graphite sheets as intercalated atoms. The pressure is raised to the desired value at the rate of 1GPa/1MD steps and kept constant after that. The temperature is also kept constant. The simulation in this way does not cause any kind of abrupt change in the cell dynamics. In the constant pressure simulation, we set the virtual mass assigned to cell deformation heavy enough so that the inner stress shall go parallel with the increase in the external pressure.

## 3 Simulation of graphite-diamond transition

In this section, the result of graphite-diamond transition is shown.

In the present work, the graphite directly turns into cubic diamond at 150-200GPa and 5000-6000K. However, hexagonal phase is not obtained.

Figure 1 shows the growth of interlayer bonding between two neighboring graphite sheets. The compression process is executed at 200GPa, 6000K. The unit cell contains 12 graphite sheets with AB stacking and each graphite sheet has 40 carbon atoms. The development of interlayer bonding is estimated by the ratio defined as follows.

$$R_{i,i+1}=\frac{2\times \text{(Number of atoms with interlayer bonding between i and (i+1)-th layer )}}{\text{ Number of atoms in a layer}}$$
(1)
If this ratio is 0, there is no bonding structure between the corresponding two layers. On the other hand, if this ratio becomes 1, the two sheets are perfectly glued with each other by interlayer bonding.

The transition process is described as follows. When the structural fluctuation has grown enough, the diamond structure grows very rapidly. The development of diamond structure along the c-axis is slower than that goes parallel with the original graphite plane. The domain boundaries of the diamond structure come to play a role of the mold so that they combine with nearby graphite sheet. In this way, the diamond structure will grow as a collective phenomenon. Even if the more wide graphite sheets are used in the unit cell of the MD simulation, the transition occurs similarly in the same range of the temperature and the pressure.

## 4 Catalyst-aided transition

In this section, the role of catalyst is investigated by the simulation including catalyst atoms, which are randomly distributed between graphite planes as intercalated atoms. In order to simulate the catalyst effect, we make use of a model, where the only s-orbital of catalyst interacts with carbon. As the tight binding parameters between carbon and catalyst, we start from those between carbon and hydrogen. We introduce a parameter $`\lambda `$ and construct the model of catalyst-carbon coupling by multiplying $`\lambda `$ to matrix elements and repulsion terms between hydrogen and carbon in the following way.

$`t_{sp\sigma }\lambda t_{sp\sigma }`$

$`t_{sp\pi }\lambda t_{sp\pi }`$

$`t_{ss\sigma }\lambda t_{ss\sigma }`$

$`V_{\text{rep}}\lambda V_{\text{rep}}`$ .

The site energy for s orbital is fixed to that of hydrogen.

The dependence of the cohesive energy and the bond length on $`\lambda `$ in a C2H6 like molecule is given in figure 2.

The transition process from graphite to diamond proceeds at lower pressure and temperature with the aid of catalyst. In the present simulation, the transition temperature to diamond is 150-200GPa in the compression of genuine diamond. However, if the catalyst is induced, at most 4-6 catalyst atoms in 240 carbon atoms, and even if the catalyst-carbon coupling is not so strong ($`\lambda 0.6`$), the transition pressure is decreased below 100GPa. It will be attributed to the catalyst-induced buckling of the graphite basal plane, in which $`\pi `$ orbitals are activated\[Figure 3\].

## 5 Formation mechanism of hexagonal diamond

A kind of stable diamond structure called as hexagonal phase is sometimes obtained by shock compression and static compression. This phase has the structural property very similar to normal cubic diamond phase. The first principle calculation shows that the cohesive energy in hexagonal phase is slightly higher than that of cubic phase. In the TB model in the present work, the difference in cohesive energy between CD and HD is also very small. From the energetic viewpoint, it is also possible for hexagonal phase to be formed in the compression of graphite. However, the simulation in the former sections, generated only cubic phase. It will be attributed to the difference of the height of activation barrier\[Figure 4\] between the transformation paths from graphite to cubic phase and hexagonal phase, as is demonstrated by a first principles calculation .

According to the molecular dynamics simulation in the present work, under the presence of catalysts, the hexagonal phase is often formed, when we set certain types of initial geometry. The figure 5 shows one of the results(150GPa,6000K,$`\lambda =1.2`$). The transition path to diamond phase is postulated in the following way. When cubic diamond is formed, the stacking should change from AB stacking of graphite into ABC one, as graphite sheets slide collectively. On the other hand, in hexagonal phase, the AB stacking must be kept. In the geometries under the presence of impurities, the collective slide may be disturbed and ”pinned” by the coupling between catalyst and carbon\[Figure 6\]. Then the short range order of hexagonal diamond is formed around the catalyst and it finally covers up whole crystal.

Figure 7 and 8 show the simulated result for a certain initial geometry with impurities. Figure 7 shows the synthesized phases after the compression as long as 2ps (2000MD steps) with various conditions for $`\lambda `$, and the temperature and the pressure. The simulations are started from the same initial geometry. In figure 8, in order to investigate microscopic mechanism of the transition, the process of interlayer bonding with different conditions is shown, using plots similar to figure 1. Fig.8(a)-(c) show the formation of synthesized phases in Fig.7. The common initial configuration of impurities used in the simulation in Fig.7 and Fig.8 is given in Fig.8(d). In the simulated result, it seems that there is a certain tendency with respect to $`\lambda `$, temperature and pressure in the transition. In a certain range of $`\lambda `$ and the temperature, the transition to hexagonal phase occurs. However, if the temperature is set to be higher while $`\lambda `$ is fixed, the hexagonal phase tends to vanish. Instead, the transition to cubic phase occurs. When the value of $`\lambda `$ becomes large, the graphite turns into cubic phase, but the long range order is often damaged and large defects remain. The result of Fig.7 and Fig.8 is obtained by a set of test-run and it does not necessarily mean the phase diagram itself. However, if the different condition, such as the different initial velocity of each atom, is given, the simulation gives similar tendency when the similar initial configuration is used.

From Fig.8, the scenario of graphite-diamond transition induced by the competition of several physical conditions, i.e., temperature, pressure, presence of impurity, and so on, is described as follows.

* In the case of genuine graphite

In this case, the formation of interlayer bonding proceeds rapidly in the direction parallel with graphite plane, while it goes slowly along c-axis. In this case, the transition path to cubic phase having the lower activation barrier will be selected\[Figure 4\]. The simulated result shows the it costs long time before the formation of diamond structure begins. It will be interpreted to be the time lag until the compressed system will find the saddle point in the potential surface.
* Under the presence of catalyst

Since the presence of catalyst causes the combination with carbon and the buckling of graphite plane, the formation of interlayer bonding begins around impurities in earlier stage of the reaction\[Figure 8(a)\]. However, the catalyst is a kind of disorder and it somewhat disturbs the formation step of interlayer bonding and causes the fluctuation in the cooperative transition toward cubic phase, letting down the transition speed\[Figure 8(b)\]. Owing to staggering in the layer-by-layer transition steps, there arises possibility of the transition to different carbon polymorph other than cubic diamond. In addition, if the sliding movement may be hindered and ”pinned” by the weak coupling between catalyst and carbon, the graphite AB stacking is locally preserved, from which local core of hexagonal structure forms and finally covers the whole crystal\[Figure 8(b)\]. In fact, the transition to HD induced by impurities are severely dependent on randomness with respect to initial configurations of impurities and it does not always occur in the simulations like this.
* Catalyst in higher temperature

In case of too high temperature, as the structural fluctuation of graphite sheets will exceed the catalyst-induced pinning in graphite sliding, the local core of hexagonal diamond may not be formed\[Figure 8(c)\]. Thus cubic diamond will form, accompanied with strain change.
* With large catalyst-carbon interaction $`\lambda `$

When the $`\lambda `$ is set to be larger, the graphite sheets are glued firmly by impurities. Then the fluctuation in the cooperative step in the transition does not become large and the compressed graphite is rapidly stabilized into cubic phase.

In fact, in the actual conversion process, the crystal will include various structural defects other than intercalated impurities. In the neighborhood of defects, atoms can move more freely and it will be possible for local hexagonal diamond structure to be formed out of structural fluctuation. In the simulation where the points defects are included, the short-range order similar to hexagonal phase is actually obtained.

## 6 Comparison to pressure-induced transition of fcc fullerene

It is demonstrated that the fullerene fcc crystal turns into transparent glassy chips of amorphous diamond phase by shock compression. This phase is constructed from short range order of sp3 bonding and it is distinguished from normal amorphous phase of carbon. In general, amorphous carbon is interpreted to be disordered phase of graphite. By MD simulation, the pressure-induced amorphization in fullerene can also be pursued. According to the MD result, fcc fullerene collapses easily above 50GPa, 2000K and turns into amorphous phase.

Figure 9 shows the electronic structure of the finally fabricated material after the compression. These figures show the case of the following conditions, (a)65GPa,3000K,(b)125GPa,5000K, (c)125GPa,5000K with presence of 30Hydrogens. The simulation uses 240 carbon atoms , i.e., 4 C60 molecules in a unit cell. In these figures, real and dotted lines show the total DOS and the contribution from threefold carbon, respectively. The rough feature of DOS is different to that of initial fcc fullerene and it is similar to that of diamond, except that gap states contributed from dangling bonds remain between valence and conduction band. If the contribution from dangling bonds between valence and conduction bands is large, the fabricated material will lose the transparent optical property and it will not be qualified to be called amorphous diamond. However, a comparison between figure 9(a) and 9(b) suggests that the contribution from dangling bonds will be reduced after the compression at sufficiently higher temperature and pressure. (In case of figure 9(a), the ratio of fourfold atom is at most 50 %. In case of figure 9(b), the ratio of fourfold atom amounts to 80 %.) The shaded zone in Fig.9(b) corresponds to the conduction-valence gap of normal diamond structure. However, gap states are still left in it. In fact, there may be other mechanism that will reduce the contribution from dangling bonds. For example, it is well-known that the presence of hydrogen atoms reduces the number of dangling bonds in case of amorphous silicon. Figure 9(c) shows the DOS under the presence of hydrogen. In this case the ratio of fourfold carbon amounts to 60%. According to the simulations in the similar way, the number of dangling bonds decreases as the more hydrogens are doped. However, even if hydrogens are included, clear valence-conduction gap does not open because of dangling bond states. It is unlikely that a small quantity of impurities will behave as catalyst and enhance the transition to amorphous diamond, as we have seen in the previous section for the case of graphite-diamond transition. In the transition of bonding state from sp2 to sp3, it seems that higher pressure plays a more relevant role than the presence of hydrogen does. If the initial configuration with hydrogen, which is same as in Fig.9(c), is compressed by higher pressure (above 100GPa), the number of fourfold atom does not increase so much compared to Fig.9(a) of carbon only case.

In case of the compression of C60, the collapsed fullerene turns into amorphous phase. Since there are number of quasi-stable configurations in amorphous phase, the transition toward diamond-like phase is very slow and the system is apt to be trapped such transient phases by quenching. As an example, figure 10 shows the two reaction path noted as (A) and (B). In both path, the pressure is increased 0Pa to 125GPa in the interval from 0.0ps to 0.5ps in the figure, and compression is executed adiabatically in the first stage. The temperature increases from initial 300K to about 2500K at the turning point of two paths. In path (A), the simulation is executed furthermore adiabatically and the temperature finally increased to 5000K. In path(B), in which the temperature is kept to 2500K, the transition to sp3 bonding becomes slow due to lower temperature. Taking after different paths of(A) and (B), the system is trapped by different semi-stable configurations.

The presence of too many intermediate phase between fullerene and diamond will be a major difference to pressure-induced graphite-diamond transition. Graphite and diamond are located in very close configurations in the potential surface and there is no stable transient state in the transition path from genuine graphite to genuine diamond. Therefore graphite rapidly transforms into the cubic diamond without being trapped by any quasi-stable structure. Thus perfect diamond phase remains alone after the quenching and, in general, amorphous diamond phase will not be obtained in graphite-diamond transition. Concerning this, the role of catalyst in graphite diamond transition can also be stated in this way. The presence of impurities induces certain quasi-stable structures in the transition path from graphite. In the present simulation, such a quasi-stable structure corresponds to the local core of hexagonal diamond. Since such a quasi-stable structure staggers the transition steps from graphite, there arises possibility of transition to carbon polymorph other than cubic diamond.

## 7 Discussion

The present work is likely to account for some aspect of non-equilibrium process in graphite-diamond transition realized by shock compression. In shock compression, the high pressure phase, which may not have arrived equilibrium state, is quenched and frozen, since the temperature decreases very rapidly while the high pressure is still kept. In the present work, the hexagonal phase is generated as a result of non-equilibrium core formation around impurities, and such a HD phase is also stable and can be quenched. However, if the high temperature environment continues sufficiently long as is realized in the actual conversion process, and if the system will finally reach the equilibrium, the system will select cubic phase as a equilibrium phase, since cubic phase is slightly stabler than hexagonal one.

Several reports has demonstrated that the hexagonal diamond phase often grows in the compression of well-crystallized crystal. It is interpreted that the hexagonal phase is formed by the inhibited collective slide of graphite, which is caused by the difficulty of sliding movement in well-crystallized graphite. However, the present work suggests that the well-crystallized property itself is not the direct origin of hexagonal phase. It would rather be attributed to weak randomness in the crystal. The presence of disorders disturbs the collective step of the transition and gives birth to the growth of the short-range order of hexagonal phase. The well-crystallized structure of graphite helps such a short-range phase easily turn into a long-range order.

The report of the growth of hexagonal phase by static compression at low temperature will be interpreted as the catalyst-induced ”pinning effect” to collective slide in graphite sheets. In that case, the formation of local core of hexagonal phase is not disturbed by the thermal fluctuation due to low temperature. Therefore the collective slide is inhibited and the hexagonal phase is allowed to extend. However, in the present work, the hexagonal phase does not generate in the picosecond-order simulation at low temperature. Concerning with graphite-diamond transition, Yagi et al. also reported that, when the graphite is heated by YAG laser under pressure, cubic diamond is formed instead of hexagonal phase. This phenomena will be explained by the present simulation. Because of too high temperature, impurities does not suppress the sliding movement of graphite and the local core of hexagonal phase will not be formed.

The transition temperature and pressure in the present work are higher than those in the actual experiment of graphite-diamond transition ($``$20Pa, $`>`$2000K). One of the reason is as follows. The present simulation uses the perfect crystal of graphite. But the actual graphite crystal has various structural defects from which interlayer bonding will easily grow. In addition, the present simulation does not follow the shortest path from graphite to diamond. The compressed graphite drifts and fluctuates on the potential surface, and it can transform the structure only if the saddle point is find. This is different to the simulation chasing after the shortest transition path as was given in Ref.. In the present work, large fluctuation should be grown for the compressed graphite to cross the saddle point in a finite, possibly very short, simulation time. If we use much longer simulation time, it is likely that the transition occurs in lower pressure and temperature. It may be possible that the tight-binding parameters are not accurate, since they are obtained by the fitting to graphite and diamond structure, but not fitted to the intermediate phase through transition. However, the energy difference between the intermediate phase and the cubic diamond is evaluated as about 0.5-0.6 eV, which is no so different to the first principle result ($``$0.4eV)

An ab-initio calculation in the reference of  suggests that, in the compression of heavily doped GIC( LiC12), hexagonal diamond forms in a way such that carbon cages will hold lithium in their centers. In it the c-axis of the generated hexagonal diamond is parallel to that of original GIC. The periodic distortion in graphite planes by intercalated atoms in heavily doped LiC12 appears to be the key to diamond transition under pressure. However, such a transition path to hexagonal phase is essentially different to that suggested by the present work. In the present work, the c-axes of synthesized hexagonal diamond and those of the initial graphite are perpendicular, as are demonstrated in several experiments of hexagonal phase formation. In addition, the present simulation treated the case where the doping density of intercalated atom is much smaller than that of Ref.. The present work suggests that randomness with respect to the initial configuration of impurities plays an essential role.

## 8 Summary

In the present work, we simulated diamond transformation from graphite. By compression, local core structures of diamond are formed in graphite at first. The transition proceeds collectively, accompanied with fluctuation, as the local diamond domain grows its volume by turning graphite sheets into the new boundary. In the present simulation, the genuine graphite turns only into cubic diamond, but not into hexagonal one. This will be attributed to the presence of higher activation barrier from graphite to hexagonal phase. The pressure-induced structural transition from graphite is subject to the presence of disorder, such as catalyst and possibly structural defect. We used a model which mimics the catalyst-aided diamond transformation from graphite. One of the role of the catalyst is the large reduction of transition pressure. Since the combination of carbon and catalyst causes the buckling of graphite plane where $`\pi `$ electrons are activated, the interlayer bonding is easily formed there. At the same time, catalyst disturbs the steps of the rapid collective transition procedure into cubic diamond. In some cases, owing to impurities, the collective slide of graphite layers is prohibited and hexagonal diamond structure is locally formed, from which the whole cell turns into hexagonal diamond. The simulated result suggest that ”weak disorder” will be one of the origin of well-crystallized carbon polymorph such as hexagonal phase.

## Acknowledgment

The computation in this work has been done using facilities of the supercomputer center, Institute for solid state physics, University of Tokyo.

## Figure captions

* Figure 1: Interlayer bonding formation in GR-CD transition. The graphite sheets are indexed as 1,2,..12 from bottom to top in the unit cell with periodic condition. The notation, for example, ”1+2”, means the interlayer bonding formation between 1st and 2nd graphite layers and so on. One graphite sheet contains 40 carbon atoms.
* Figure 2: Parameter dependence of carbon-catalyst bonding in a C2H6-like molecule. The bond angles are little dependent on $`\lambda `$.
* Figure 3: Buckling of graphite plane (gray sphere) caused by catalyst (white sphere).
* Figure 4: Activation barrier from graphite to diamond.
* Figure 5: Hexagonal diamond structure obtained by the simulation( 150GPa, 6000K, $`\lambda =1.2`$).
* Figure 6: Formation of cubic and hexagonal diamond. (a) Growth of cubic diamond.(b) Growth of hexagonal diamond.
* Figure 7: Synthesized phase after the compression with various condition, started from a certain initial geometry. (The initial geometry is given in figure 8.) The meaning of notations in the figure is as follows. GR:Graphite, CD(W): Well-crystallized cubic diamond phase, CD(D): cubic phase with large defect, HD(W/D):Hexagonal phase, HD+CD:Coexistence of hexagonal and cubic phase. The transformed structures after MD simulations of 2ps are shown there.
* Figure 8: Formation of interlayer bonding with the presence of impurities. In each simulation, the pressure is increased from 0 to 150GPa in the first 0.15ps and kept constant after that. (a)-(c) :Plot with various conditions, each of which corresponds to synthesized phases listed in figure 7. (d):Initial geometry with 6 impurities used in the simulation of figure 7 and 8. Broad lines and crosses show the relative location of 6 graphite sheets and 6 impurities in a cell. One graphite sheet includes 40 C atoms.
* Figure 9: Density of states in amorphous structures synthesized by the compression of fullerene FCC crystal. (a): DOS in the quenched structure after the compression at 3000K, 65GPa. (b): after the compression at 5500K, 125GPa. (c): after the compression at 3500K, 55GPa, with the presence of 30 hydrogen atoms among 240 carbon atoms in the unit cell. In these figures, the real and dotted lines show the total DOS and the contribution from threefold carbon, respectively.
* Figure 10: Formation of fourfold atoms in the compression of fullerene fcc cell. This figure shows the change in the ratio of fourfold carbon thorough the compression. In path (A), the simulation is executed adiabatically and the temperature increases to about 5000K. In path (B), the simulation is executes adiabatically at first, and afterwards, the temperature is kept 2500K.


text:
# Critical exponent 𝜂ᵩ of the Lattice London Superconductor and vortex loops in the 3D XY model

\[

## Abstract

The anomalous dimension of the lattice London superconductor is determined from finite size scaling of the susceptibility to be $`\eta _\varphi =0.79(1)`$. Indirect determinations of $`\eta _\varphi `$ from properties of the vortex loops in the 3D XY model are also attempted, but it is found that the results are sensitive to details in the simulations related to vortex loop intersections. It is suggested that the same value of $`\eta _\varphi `$ can at most be obtained from vortex loop properties in the limit of low vortex density.

\]

The properties of the Meissner transition is a classical problem in statistical physics. Whereas this transition was originally believed to always be of first order, the work by Dasgupta and Halperin gave strong arguments that the transition should instead be continuous. The argument is based on a duality transformation of the lattice London Superconductor (LLS) and suggests that the transition should be 3D XY-like, but with the temperature scale inverted. A direct consequence of this relation is the expectation that the correlation length exponent $`\nu `$ should be the same in the superconductor as in the 3D XY (planar rotor) model. However, with fluctuations in both the phase angle and the gauge field, it becomes possible to define two characteristic lengths, and the behavior of the magnetic screening length $`\lambda `$ has recently been a subject of some controversy. The current evidence points to a scenario where both characteristic lengths diverge with the same exponent. This is related to the presence of an anomalous dimension $`\eta _A=1`$ for the gauge fluctuations.

What is needed beside $`\nu `$ and $`\eta _A`$ for the characterization of the universality class of the transition is a knowledge of the anomalous dimension $`\eta _\varphi `$ associated with the phase correlations. This quantity was first determined to be $`\eta _\varphi =0.20`$ from a renormalization group calculation to one loop order. Determinations from simulations have so far only been indirect through the properties of vortex loops in the 3D XY model. The main motivation for that analysis was the possible connection between the sign of $`\eta _\varphi `$ and the existence of a vortex loop blowout transition in high-$`T_c`$ superconductors in applied magnetic fields.

In this Letter we report on a direct determination of the anomalous dimension $`\eta _\varphi `$ from Monte Carlo (MC) simulations of the LLS that gives a surprisingly large negative value, $`\eta _\varphi =0.79\pm 0.01`$. We also consider in detail the approach to determine this exponent from the properties of the vortex loops in the 3D XY model. The main conclusion from that study is that the results depend strongly on the details in the simulation related to vortex loop intersections. We suggest that the theoretically expected value can at most be found in the limit of low vortex density.

The Hamiltonian of the LLS is

$$H=\underset{i\mu }{}\left\{U(\varphi _{i+\widehat{\mu }}\varphi _iA_{i\mu })+\frac{1}{2}J\lambda _0^2[𝐃\times 𝐀]_{i\mu }^2\right\},$$
(1)
where $`\varphi _i`$ is the phase of the superconducting wave function on site $`i`$, and $`A_{i\mu }`$ is the vector potential on the link starting at site $`i`$ and pointing along $`\widehat{\mu }`$. The sum is over all bonds of a 3D simple cubic lattice of unit grid spacing and $`\mu =x`$, $`y`$, $`z`$. In the first term, the kinetic energy of flowing supercurrents, $`U(\phi )`$ is the Villain function

$$e^{U(\phi )/T}=\underset{p=\mathrm{}}{\overset{\mathrm{}}{}}e^{J(\phi 2\pi p)^2/(2T)}.$$

In the second term, the magnetic field energy, $`\lambda _0`$ is the bare magnetic penetration length, and $`𝐃\times 𝐀`$ is the discrete circulation of the vector potential.

In our simulations we fix $`\lambda _0=0.3`$ to be able to compare with Ref.  and fix the gauge through $`𝐃𝐀=0`$ to facilitate a simple determination of the spin correlations. To fulfill this constraint we update the $`A_{i\mu }`$ by simultaneously adding $`\delta A`$ to four $`A_{i\mu }`$ around an elementary plaquette, a procedure which was also used in Ref. . We perform our MC simulations with the standard Metropolis algorithm, first sweeping sequentially through the phase angles $`\theta _i`$ and then sweeping three times with attempts to change the circulation of $`𝐀`$ in the three different directions. The length of the runs close to $`T_c`$ were at least $`10^7`$ sweeps through the lattice, but often much longer. The measured quantities discussed here are the second and fourth moment of the magnetization,

$$m^p=\left|\frac{1}{L^3}\underset{i}{}e^{i\varphi _i}\right|^p,$$

which are used to obtain the dimensionless fraction

$$Q=\frac{m^2^2}{m^4},$$
(2)
similar to Binder’s cumulant, which is used to determine $`T_c`$. The susceptibility at the transition is obtained from the standard relation $`\chi =L^3m^2`$.

The critical temperature for the LLS with the same parameter, $`\lambda _0=0.3`$, has already been determined with high precision. Our present analysis of $`Q`$ therefore mainly serves to confirm that the critical properties may be correctly determined from our simulations with $`𝐃𝐀=0`$. Figure 1a shows $`Q`$ versus temperature for system sizes $`L=8`$, 12, 16, and 24. We find that the curves for $`L>8`$ cross at $`T0.80`$. The $`L=8`$ data does however not cross the other curves at the same temperature, a typical sign of corrections to scaling. This data is therefore not included in the scaling collapse.

To get a more precise determination of $`T_c`$ we use all our data for $`12L32`$ in a narrow range around $`T_c`$, $`|tL^{1/\nu }|<1.1`$, and assume the scaling form $`Q(t,L)=f_Q(tL^{1/\nu })`$, where $`t=T/T_c1`$, and $`f_Q(x)`$ is a scaling function. We fix $`\nu =0.672`$ and let $`f_Q(x)`$ be a fifth order polynomial in $`x`$. The data collapse nicely, Fig. 1b, and we find $`T_c=0.800\pm 0.001`$ in good agreement with $`T_c=0.8000\pm 0.0002`$ from Ref. .

To determine the anomalous dimension $`\eta _\varphi `$ we make use of the standard finite size scaling relation $`\chi (T_c,L)L^{2\eta _\varphi }`$. The susceptibility at $`T_c=0.8`$ is shown in Fig. 2a. The points do indeed to an excellent approximation fall on a straight line and our result is $`\eta _\varphi =0.79\pm 0.01`$. Panel (b) shows a determination of the exponent $`\beta `$ from the temperature dependence of $`m^2`$. This data is for temperatures at which we expect the finite size effects to be negligable. The result is $`\beta =0.069\pm 0.003`$ in good agreement with what we expect from $`\eta _\varphi `$: $`\beta =\frac{1}{2}\nu (d2+\eta _\varphi )=0.070\pm 0.003`$.

We now turn to the alternative and indirect approach to determine $`\eta _\varphi `$ through the behavior of the dual model which is the ordinary 3D XY model. The idea is that the vanishing of the line tension of the vortex loops (see below) as $`TT_c^{}`$ is governed by $`\gamma _\varphi `$,

$$ϵ(1T/T_c)^{\gamma _\varphi },$$
(3)
which is related to $`\eta _\varphi `$ through the Fisher scaling relation $`\gamma _\varphi =\nu (2\eta _\varphi )`$. The basis for Eq. (3) is the connection between $`\varphi `$ correlations in the LLS and the vortex loops in the 3D XY model. Whereas this less direct determination of $`\eta _\varphi `$ of course is more prone to errors, it is interesting to examine if the same value of $`\eta _\varphi `$ can be obtained in that way. We will in the following examine the vanishing of $`ϵ`$ for two different cases, examine the distribution of vortex loop diameters at $`T_c`$, and finally comment on the analysis in Ref.  that gave $`\eta _\varphi 0.18`$.

The standard way to locate vortices in a XY model is from the angular difference $`\phi _{ij}=\theta _i\theta _j`$ between nearest neighbors, with $`|\phi _{ij}|<\pi `$. For each plaquette the vorticity is then obtained from

$$n=\frac{1}{2\pi }(\phi _{12}+\phi _{23}+\phi _{34}+\phi _{41}).$$
(4)
With a Villain interaction it is also possible to define vortices in a different way, which is inspired by the wire model. For each link one first calculates $`\phi _{ij}^0=\theta _i\theta _j`$ and then probabilistically sets the angular difference to $`\phi _{ij}=\phi _{ij}^02\pi p_{ij}`$, where $`p_{ij}`$ is an integer chosen with the relative weight $`e^{J(\phi _{ij}^02\pi p_{ij})^2/(2T)}`$. The vorticity is finally calculated from Eq. (4). The vortices obtained in that way are exactly equivalent to the dual vortex line model, a fact that follows from a straightforward generalization of the results in Ref.  to three dimensions. The density of these probabilistic vortices (PV) is always higher then the density of standard vortices (SV).

The simulations of the 3D XY model were performed with the Wolff cluster update method and Villain interaction on a cubic lattice with $`L=128`$. The average for each data point is based on typically $`10^5`$ measurements. The simulations were performed with both SV and PV discussed above. The tracing out of vortex loops were done by always choosing the path randomly when two (or more) vortex loops meet at the same elementary cube.

After identifying the vortex loops we measure both the perimeter $`p`$ and the “diameter” $`a`$ of each such loop. Whereas $`p`$ is directly given by the number of vortex loop segments the calculation of $`a`$ is somewhat more involved. With a loop made up by $`p`$ segments with each segment starting at $`𝐫_i`$ and pointing along the unit vector $`𝐧_i`$. the equivalent to the “magnetic moment” is

$$=\underset{i=1}{\overset{p}{}}𝐫_i\times 𝐧_i,$$

and the “diameter” becomes $`a=2\sqrt{||/\pi }`$. From these values we calculate the distributions $`D(p)`$ and $`D(a)`$.

The line tension $`ϵ`$ is obtained by fitting

$$D(p)p^\alpha \mathrm{exp}(ϵp/T),$$
(5)
with $`\alpha `$, $`ϵ`$ and a prefactor as free parameters. To get good quality fits one has to restrict the fitting to large $`p`$ only. We found that good quality fits are obtained by only using data for $`p>T/ϵ`$. Figure 3 shows the line tension $`ϵ`$ determined with both SV and PV (open and solid symbols, respectively). Note that we are here examining a very narrow temperature region to be able the probe the behavior in the critical region; most data points are within a few percent below $`T_c=3.0024(1)`$. With logarithmic scales on both axes the slope as $`TT_c`$ is expected to give the exponent $`\gamma _\varphi `$. To facilitate a comparison with the above obtained $`\eta _\varphi =0.79`$ we draw a dashed line in Fig. 3a that corresponds to $`\gamma _\varphi =\nu (2\eta _\varphi )=1.87`$. The slopes of the data seem to be in reasonable agreement with that prediction.

A more detailed analysis, however, casts doubt on this agreement. Panel (b) in Fig. 3 shows the temperature dependence of the local slope where each value is calculated from two consecutive points. The curvature in the original data is here seen as a trend towards larger values as $`T_c`$ is approached. There is, however, no evidence that the local slope saturates at $`\gamma _\varphi =1.87`$ (horizontal dashed line); the last point lies significantly above the predicted value. To check for finite size effects the relevant points in panel (a) have also been analyzed with $`L=192`$ but with no change in the result. For the PV data (solid circles) this deviation is even more significant.

We will below present more evidence that suggest that, generally speaking, the properties of vortex loops at or close to $`T_c`$ tend to not agree with theoretical predictions. We will also argue that this discrepancy is related to the existence of vortex line intersections which are normally not accounted for in theoretical considerations. The quantity we focus on is $`D(a)`$ – the distribution of vortex loop diameters – at $`T_c`$. There is a simple analytical prediction for this quantity based on the scale invariance of the system at $`T_c`$,

$$D(a)a^\lambda ,\lambda =4\text{ at }T_c.$$

Assume that there is on the average $`N`$ loops with diameter in the interval $`(a,a+da)`$ in a system with volume $`L^3`$, i.e. $`D(a)da=N/L^3`$. Scale invariance implies that the distribution should not be changed by a rescaling of the system. After rescaling with a factor $`s`$ the new system size is given by $`L^{}=L/s`$ and $`N`$ is now the number of loops with diameter in the interval $`(a/s,a/s+da/s)`$. We therefore get $`D(a/s)da=sN/(L^{})^3=s^4D(a)da`$ which in turn leads to $`D(a)a^4`$.

As shown in Fig. 4 this expectation for $`\lambda `$ is not confirmed by the simulations. The values of $`\lambda `$ for SV and PV are instead found to be $`4.13(1)`$ and $`4.18(1)`$, respectively. The deviation from the expected $`\lambda =4`$ increases with increasing vortex density and appears to be proportional to the vortex density, $`0.156`$ and $`0.226`$, respectively.

We now suggest that the reason for $`\lambda 4`$ is related to the existence of vortex loop intersections. Consider the rescaling of a loop with a scale factor $`s`$. A non-trivial effect of the rescaling is related to the smaller resolution of a rescaled loop. This may lead to new vortex intersections that could change the outcome of the tracing. As an example, consider the vortex line loop in Fig. 5. The original loop (a) will always be treated as a single one, but since the rescaled loop (b) enters the same elementary cube twice the tracing can now give the loops in panel (c) or (d), each with a probability of 50%. Since a splitting will give an increased weight to smaller loops, the effect is to make $`\lambda `$ larger. Note that a higher vortex density will make this effect more important which is in accordance with our results for $`\lambda `$.

From the above it seems that the presence of vortex loop intersections is the reason for the failure to obtain the expected $`\lambda =4`$ from $`D(a)`$. One could expect the vortex loop intersections to also affect other vortex loops properties. This is in accordance with our results from the analysis of the vortex loop line tension that the high density data (PV) is far off whereas the low density data (SV) is fairly close to the expected behavior. The correct behavior of $`ϵ`$ should then at most be expected in the limit of low vortex density.

The determination of the line tension for vortex loops in a 3D XY model was first made in Ref. . Their simulations were with cosine spin interaction ($`T_c2.2`$) and a different method for tracing out the loops; self-intersecting vortex loops are always split into two. They found $`\eta _\varphi =0.18\pm 0.07`$, in agreement with the analytically obtained $`\eta _\varphi 0.2`$. However, a comparison with our more precise MC data (now obtained with cosine interaction and splitting of intersecting loops) casts strong doubt on their analysis. Figure 6 shows our values for $`D(p)`$ (dots) together with both the data (open circles) and the fitting curves (solid lines) for $`T=2.0`$ and 2.1 from Ref. . We first note that the two different sets of MC data agree well. From the large-$`p`$ part of the data it is however clear that the fitting curves from Ref.  do not agree with our data. The same discrepancy can actually also be seen for $`T=2.0`$ in Fig. 3 of Ref. . The inset in Fig. 6 shows our values for $`ϵ`$ obtained from a fit to Eq. (5) only including data for $`p>1.2T/ϵ`$. The slope differs clearly from $`\gamma _\varphi =1.45`$ shown by the dashed line. We therefore conclude that the good agreement with the analytically obtained $`\eta _\varphi 0.2`$ was only accidental. In our analysis the slope approaches $`\gamma _\varphi 1.15`$ as $`TT_c`$. If one instead allows for vortex loop intersections $`\gamma _\varphi `$ behaves much the same as in Fig. 3 (again with $`\gamma _\varphi `$ signficantly above 1.87). This again shows the central role of the intersections for the vortex loop properties.

To conclude, our main result is a direct determination of the anomalous dimension in the LLS, giving $`\eta _\varphi =0.79\pm 0.01`$. We have also attempted determinations of the same quantity from the vanishing of the line tension in the 3D XY model, but found that it was not possible to get a well-defined value for the exponent $`\gamma _\varphi `$. From considerations involving the distribution of vortex loop diameters we suggest that it is the vortex loop intersections that are responsible for the failure to find the analytically expected values of some vortex-loop related exponents. This implies that the “correct” behavior at criticality can at most be expected in the limit of low vortex density.

The author thanks P. Minnhagen, A. Sudbø, S. Teitel, and Z. Tes̆anović for valuable discussions. This work was supported by the Swedish Natural Science Research Council through Contract No. E 5106-1643/1999.


text:
# Extended radio emission after the soft X-ray maximum of the NOAA 9077 AR solar flare on July 10, 2000

## 1 Introduction

It is widely accepted that the nonthermal radio burst emission should take place before the thermal SXR maximum, for the basic process is thought to be the following: first, high-energy particles are accelerated during a reconnection of magnetic flux and propagate along field lines, exciting the radio bursts; then, these particles pour down continuously, causing chromospheric evaporation, and increasing the SXR emission to its peak (Heyvaerts et al. 1977; Mclean & Labrum 1985).

However, an event of extended nonthermal emission after the thermal maximum has been detected by several observatories on July 10, 2000. This event seems to suggest that the extended nonthermal radiation resulted from repeated triggering of a new magnetic reconnection. There are two different ways to interpret the triggering. In the emerging flux model (Heyvaerts et al. 1977), when the magnetic flux loops emerge from below the photosphere and interact with the overlying field, continuous reconnection occurs in the current sheet between the new and old fluxes. In the SQR model (Aschwanden 1998), the shear motion of the footpoint will cause the magnetic field between the sheared small-scale loop and the overlying unsheared large-scale loop to become increasingly sheared and compressed, triggering quadrupolar X-type reconnection. The latter interpretation also regards, as a consequence of the reconnection, the various radio fine structures, including type III, and U bursts, and reverse-slope (RS)-drifting bursts. Indeed, some such radio fine structures were observed with the 1.0-2.0 GHz spectrometer of BAO in the event on July 10, 2000.

This Letter first presents the observations of the event. Then we use the SQR model to analyze the radio fine structures and the related broadband radio bursts. We find that both the shear motion and emerging motion play a part during the magnetic reconnection for this event.

## 2 Observations

The observations of PALE (Palehua station) and GOES 8 were prepared by Space Environment Center (SEC) on the internet. The other observations were made with the 1.0-2.0 GHz solar radio fast dynamic spectrometer and the 2.84 GHz radio telescope at Huairou Solar Observing Station of BAO (Fu et al. 1995).

The July 10, 2000 flare, classified as M5.7/2B, was observed at 21:05-22:27 UT in Active Region 9077 at position N19E49. Maximum was at 21:42 UT. Morphologically, the flare was a two-ribbon flare. The peak flux density at 2695 MHz was 140 sfu at 21:27 UT before the SXR maximum (by SEC). However an extended and very strong radio emission (up to 3085 sfu at 2840 MHz) was detected after 21:42 UT, and continued to about 23:00 UT, by several observatories. The data are summarized in Table 1. It indicates repeated occurrence of magnetic reconnection.

A group of radio fine structures before the broadband radio bursts were observed in the range of 1.0-2.0 GHz 21:59:30-22:00:10 UT (Fig. 1). Fig. 2 shows the GOES 8 SXR data between 21:00 and 22:30 UT. Our radio data at the several frequencies, which did not begin until 21:58 UT, are displayed in the inset, and are further blown up in the grey-scale picture below. The observed features can be summarized as follows.

(1) The fine structures contained four quasi-periods, the first one about 6 s between 21:59:32–21:59:38 UT, the second about 9 s between 21:59:38–21:59:47 UT, the third about 4 s between 21:59:47–21:59:52 UT, and the fourth about 14 s between 21:59:52–22:00:05 UT. The average period was about 6.5 s.

(2) Each period included both type U and type RS bursts. The type U bursts occurred in the lower frequency range, and the type RS bursts, in the higher frequency range. All the type U bursts had a similar turnover frequency of 1.16 GHz.

(3) In each period, the flux density of the type U bursts was 3-4 times stronger than the type RS bursts.

(4) The average frequency drift rate of the type U bursts was about 100 MHz s-1 (both positive and negative), and that of the type RS bursts was about 800 MHz s-1 (negative). The former was clearly lower.

(5) The central frequencies of the type RS bursts drifted fast towards the lower frequency. The rate of drift decreased in time in the range 100-40 MHz s-1.

(6) The type U bursts occurred after the type RS bursts.

(7) At 22:04 UT, some 4 minutes after the fine structures, a broadband radio burst appeared, which lasted until 22:49 UT.

## 3 Theoretical analysis

### 3.1 SQR model

According to Aschwanden (1998, Fig. 1 of his paper), the SQR model may contain the following process. Recognizing that shear motion along a neutral line to be an important condition for flaring, a consequence of a shear motion of the footpoint along the neutral line is that the magnetic field between the top of the sheared small-scale loop and the overlying unsheared large-scale field becomes increasingly sheared and compressed. At some point the rising sheared loop will intersect with the unsheared large-scale loop and trigger a quadrupolar X-type reconnection. During the reconnection, the connectivity of the field lines switches by exchanging the polarities of equal signs. The newly configured field lines then slip back from the X-point and relax to two dipole-like field lines. Then with the curvature of the magnetic field reduced, magnetic energy will be released (Aschwanden 1998, Aschwanden et al. 1999).

### 3.2 Radio bursts

As shown in Fig. 3, upward propagating electron beams, along open and closed field lines, can be detected in the form of Type III, J and U bursts; while downward propagating electron beams are detectable in the form of RS bursts. Therefore type U bursts occur in the lower frequency range while type RS bursts occur in the higher. The following expected observational signature can also be seen in Fig. 3: a quasi-periodic sequence of type U bursts with similar turnover frequencies produced by the upward propagating electron beams (Aschwanden 1998). The observational properties (1) and (2) in Sect. 2 are in agreement with the theoretical predictions here.

According to Zhao et al. (1997), if we assume the magnetic field of the active region to be a dipole field above the photosphere, then the relation between the height of radio source ($`H`$) and its frequency ($`f`$) can be estimated by

$$H=d\left[\left(\frac{5.6B_0}{f(MHz)}\right)^{1/3}1\right].$$
(1)
Where the depth of the dipole field $`d`$ is about $`3.5\times 10^4`$ km, and $`B_0`$ is the magnetic field strength of the photosphere. For this event $`B_0`$ was about 2500 G (Huairou Station, BAO). Using the turnover frequency of 1.16 GHz, the height of the top of the unsheared large-scale loop can be estimated to be about 4.5 $`\times `$ 104 km.

Upward-accelerated electrons encounter a decreasing magnetic field, which has a focusing effect on the pitch-angle distribution of the beam and makes it more prone to plasma emission. Also, the lower density encountered in upward direction makes it easier for plasma emission to escape. These two properties may be the main reason why upward propagating electron beams are detected much more frequently (type U bursts etc.) than downward beams (RS bursts) (Aschwanden 1998). Therefore, when the upward and downward beams are both detected in radio, the flux density of type U bursts can be stronger than the related type RS bursts. This agrees with observational property (3) in Sect. 2.

As the name implies, the large-scale loop have generally a much larger size than the small-scale loop: the typical height is 1.3 $`\times `$ 105 km for the former and 2.0 $`\times `$ 104 km for the latter (Aschwanden et al. 1992). The upward velocity of the electrons along the large-scale loop, then, is less than the downward velocity along the small-scale loop. Thus, the frequency drift rate of type U bursts is lower than that of type RS, as was stated under (4) in Sect. 2. According to Eq. 1, the spatial range ($`\mathrm{\Delta }H`$) can be obtained for each of the type U and RS bursts. Dividing $`\mathrm{\Delta }H`$ by the corresponding time interval, we can obtain the electron velocity for each burst. The average upward and downward electron velocities were estimated to be of 1.8 $`\times `$ 103 km s-1 and 1.4 $`\times `$ 104 km s-1, respectively.

From the observational data, the central frequency (cf. the observed feature (5) in Sect. 2) of each RS burst can be obtained. Substituting the central frequency into Eq. 1, the corresponding height ($`H_c`$) can be estimated for each RS burst. We found $`H_c`$ increased in time. This indicates that the small-scale loop was emerging from below. In the first period the velocity of emerging was faster, but it gradually slowed down in the three succeeding periods (cf. Feature (1) in Sect. 2). Typical emerging velocity during the second period is estimated at about 1500 km s-1. Thus, besides the shear motion of the footpoint, which has been emphasized in the SQR model, there is the emerging motion of small-scale magnetic loop taking place in the course of the reconnection.

It was also expected from the SQR model that for a large flare the reconnection will continue for a long time; for example, lasting more than 3 minutes for the 1994-Jan-06 0405 flare (Aschwanden 1998). In this event, the flare was a large flare (two-ribbon flare). Moreover the broadband radio bursts are believed to be the result of thermal bremsstrahlung from the hot and dense plasma that remains in the low corona even for some tens of minutes (Dulk 1985). Thus we can understand the reason why the broadband radio bursts appeared after the fine structures for about 4 minutes.

## 4 Discussion and Conclusions

Several conclusions can be drawn from the present study:

(1) The extended radio emission was indeed caused by a new magnetic reconnection process.

(2) There was a flux loop emerging from the below. Typical emerging velocity was about 1500 km s-1.

(3) Either the shear motion or the emerging motion or both could have triggered the quasi-periodic magnetic reconnection.

(4) The radio fine structures provide clues for the understanding of the process of energy release in the initial phase of the flare.

The original SQR model only considered the shear motion as the trigger of the magnetic reconnection, but not the emerging motion. However the emerging motion was indeed reflected in the present observations. The model may thus be improved in the future. Besides, the type U bursts occurred after the type RS bursts, and the delay between the two tended to get shorter. This property may have something to do with the emerging motion, because the consequence of the small-scale loop emerging is that the tops of the two loops eventually reach the same height.

###### Acknowledgements.

This research was supported by the National Nature Science Foundation of China under grant No. 19773016 and No. 19973008, and by Ministry of Science and Technology of China under grant No. G2000078403. We wish to thank the anonymous referee for advice on the improvement of the Letter. We are also grateful to Dr T. Kiang for his revision on the language of the Letter.


text:
# 1 Introduction

## 1 Introduction

One of the main goals of the BABAR experiment is to study CP Violation in neutral $`B`$-mesons. The neutral $`B`$ meson system is similar to the neutral Kaon system in that we have two flavor eigenstates that mix. However, the phenomenology is rather different. In the Kaon system the physics is driven by a large difference in the decay widths, such that there are two states, the $`K_L^0`$ and the $`K_S^0`$, that have substantially different lifetimes. In the $`B_d`$ system the widths of the two states are very similar and instead the physics is dominated by the mass difference, $`\mathrm{\Delta }m_{B_d}`$, which controls the oscillation frequency of the $`B^0\overline{B}^0`$ system.

In the Standard Model, CP Violation arises from complex phases in the CKM matrix. By studying the interference between decays of $`B`$ mesons that decays directly to a common final state and those that mix before they decay, we can study the phases of the CKM triangle.

Unitarity of the CKM matrix allows the construction of unitarity triangles. Applying the unitarity constraint between the first and third generation gives the least degenerate triangle. The angles of this unitarity triangle are probed in CP Violation in $`B`$ decays. Figure 1 shows the normalized CKM triangle, the studies of CP Violation in $`BJ/\psi K_S^0`$ allow us to determine $`\mathrm{sin}2\beta `$.

The decays of primary interest in this presentation are $`B^0J/\psi K_S^0`$ and the very similar $`B^0\psi (2S)K_S^0`$. These decays have been dubbed the golden modes for measuring $`\mathrm{sin}2\beta `$. These modes are experimentally fairly straight forward to construct, and have branching fractions that allow us to collect samples of events that are sufficiently large to study the time dependent asymmetry. But foremost, they are the golden modes due to the very small theoretical uncertainties as there are no penguin contributions with different weak phases.

The time dependent rate for $`BJ/\psi K_S^0`$ is given by

$$f_\pm (\mathrm{\Delta }t;\mathrm{\Gamma },\mathrm{\Delta }m_{B_d},\mathrm{sin}2\beta )=\frac{1}{4}\mathrm{\Gamma }e^{\mathrm{\Gamma }|\mathrm{\Delta }t|}[1\pm \mathrm{sin}2\beta \times \mathrm{sin}\mathrm{\Delta }m_{B_d}\mathrm{\Delta }t],$$
(1)
where $`+()`$ indicates that the $`B_{\mathrm{tag}}`$, the other $`B`$ from the $`\mathrm{{\rm Y}}(4S)`$ decay, was tagged as a $`B^0`$ $`(\overline{B}^0)`$. To experimentally fit the time distribution we need to account for two effects; finite resolution in the $`\mathrm{\Delta }t`$ determination and the possibility that the wrong tag was assigned to $`B_{\mathrm{tag}}`$. The time resolution is handled by convoluting Eq. 1 with a resolution function, $`(\mathrm{\Delta }t;\widehat{a})`$. The fraction of mistags, $`w`$, dilutes the $`\mathrm{sin}2\beta `$ measurement by $`𝒟=12w`$. Experimentally we will perform a fit to

$$_\pm (\mathrm{\Delta }t;\mathrm{\Gamma },\mathrm{\Delta }m_{B_d},𝒟\mathrm{sin}2\beta ,\widehat{a})=f_\pm (\mathrm{\Delta }t;\mathrm{\Gamma },\mathrm{\Delta }m_{B_d},𝒟\mathrm{sin}2\beta )(\mathrm{\Delta }t;\widehat{a}).$$
(2)
The time dependent CP asymmetry for the $`B^0J/\psi K_S^0`$ and $`B^0\psi (2S)K_S^0`$ decays is given by

$$a_{CP}(\mathrm{\Delta }t)=\frac{N_{B^0}(\mathrm{\Delta }t)N_{\overline{B}^0}(\mathrm{\Delta }t)}{N_{B^0}(\mathrm{\Delta }t)+N_{\overline{B}^0}(\mathrm{\Delta }t)}=𝒟\mathrm{sin}2\beta \mathrm{sin}(\mathrm{\Delta }m_{B_d}\mathrm{\Delta }t),$$
(3)
where $`N_{B^0}`$ ($`N_{\overline{B}^0}`$) is the number of events where $`B_{\mathrm{tag}}`$ was assigned a $`B^0`$ ($`\overline{B}^0`$) tag. In these expressions $`\mathrm{\Delta }t=t_{CP}t_{\mathrm{tag}}`$ refers to the time difference between the decay of $`B_{CP}`$ and $`B_{tag}`$. At the $`\mathrm{{\rm Y}}(4S)`$ where the two $`B`$ mesons are produced in a coherent state, the flavor of the $`B`$ decaying to the $`CP`$ eigenstate is determined by studying the flavor of the other $`B`$, as the two $`B`$ mesons are in a coherent $`P`$-wave state we know that the flavor of $`B_{CP}`$ is of the opposite flavor of the $`B_{\mathrm{tag}}`$ at the time of the $`B_{\mathrm{tag}}`$ decay. Hence, $`\mathrm{\Delta }t`$ is the time $`B_{CP}`$ evolved from when we knew its flavor. Further, at PEP-II the collisions are asymmetric, this allows us to measure the difference between the decay times of the two $`B`$ mesons by simply measuring the separation in $`z`$ of the two $`B`$ decays, $`\mathrm{\Delta }t\mathrm{\Delta }z/c\beta \gamma `$.

The major steps in performing this analysis are

* Reconstruct the CP eigenstates.
* Measure the vertex resolution, $`(\mathrm{\Delta }t;\widehat{a})`$.
* Determine vertex separation, which gives $`\mathrm{\Delta }t`$
* Tag the flavor of the other $`B`$, $`B_{\mathrm{tag}}`$.
* Measure the wrong tag fraction, $`w`$.
* Perform likelihood fit to the $`\mathrm{\Delta }t`$ distribution to determine $`\mathrm{sin}2\beta `$.

As far as possible we try to determine the resolution function parameters and wrong tag fractions from data.

## 2 The BABAR experiment

The BABAR experiment is located at the PEP-II storage ring at SLAC. PEP-II collides electrons and positrons with energies of about 9.0 GeV/c2 and 3.1 GeV/c2 respectively at the center of mass energy of the $`\mathrm{{\rm Y}}(4S)`$ resonance. The produced $`\mathrm{{\rm Y}}(4S)`$ mesons have a boost of about $`\beta \gamma =0.56`$. Since the first recorded collisions with the BABAR experiment on May 26, 1999, PEP-II has produced excellent luminosity that have allowed BABAR to collect the worlds largest sample of $`B`$ mesons at the $`\mathrm{{\rm Y}}(4S)`$. The analysis presented here is based on 9.0 fb-1 recorded on the resonance and 0.8 fb-1 taken about 40 MeV/c2 below the resonance. This corresponds to about $`10.5\times 10^6`$ produced $`B\overline{B}`$ pairs. (The 2000 run of BABAR ended on October 30, 2000, and the total recoded luminosity in 2000 was 23 fb-1.)

The BABAR experiment is described in detail elsewhere . Here just a few key points of particular relevance to the measurement presented will be discussed.

One of the key features of this experiment is that the produced $`B`$ mesons have a boost of about $`\beta \gamma =0.56`$, and that the time difference between the $`B`$ decays is measured by the separation in $`z`$ position of the $`B`$ decay vertices. The typical vertex separation between two $`B`$ meson decays is 250 $`\mu `$m. The BABAR experiment has a 5 layered, double sided, silicon micro strip vertex detector capable of stand-alone tracking for low momentum particles, $`p_T<120`$ MeV/c, not detected in the drift chamber. Figure 2 shows the measured track impact parameter resolution in $`z`$ as a function of momentum. At higher momenta, where multiple scattering is negligible, the resolution is about 40 $`\mu `$m. For an exclusively reconstructed $`B`$ the vertex resolution in $`z`$ is typically about 40-60 $`\mu `$m.

Another unique feature of the BABAR experiment is the charged hadron identification system, the DIRC, Detection of Internally Reflected Cherenkov light. The DIRC detects Cherenkov photons that are produced in quartz bars and reflected out to a water tank instrumented with photo multipliers, see Figure 3. The DIRC has proved to work very well, we have achieved better than 3$`\sigma `$ $`K\pi `$ separation at momenta of 3 GeV/c. Further improvements are possible with a better understanding of the alignment of the DIRC with respect to the tracking system.

Overall the BABAR experiment has performed very well, important for our ability to record data is our efficiency, typically BABAR is live 97% of the time when PEP delivers luminosity. Losses are due to background spikes that cause trips, ramp-up of voltages, and sporadic outages, e.g., of computing resources. After the data is recorded, it is typically reconstructed within 24 hours.

### 2.1 Reconstruction of event samples

The CP sample consists of events reconstructed in the following modes

| $`B^0J/\psi K_S^0`$ | ($`K_S^0\pi ^+\pi ^{}`$), |
| --- | --- |
| $`B^0J/\psi K_S^0`$ | ($`K_S^0\pi ^0\pi ^0`$), |
| $`B^0\psi (2S)K_S^0`$ | ($`K_S^0\pi ^+\pi ^{}`$), |

where the $`J/\psi `$ and the $`\psi (2S)`$ are reconstructed in both the $`e^+e^{}`$ and $`\mu ^+\mu ^{}`$ channels, the $`\psi (2S)`$ is also reconstructed in the $`J/\psi \pi ^+\pi ^{}`$ channel. In the $`e^+e^{}`$ decays, bremsstrahlung photon recovery is attempted. For more details about the event reconstruction see Ref. . We obtain a total of 168 events in these modes, as shown in Figures 4-6. The yields and purities are summarized in Table 1.

To measure the performance of the tagging, and in particular to determine the wrong tag fraction, a sample of exclusively reconstructed $`B`$ mesons is used. The modes used and the yields are summarized in Table 2. The reconstruction of these samples are detailed in Ref.  and .

### 2.2 Vertex separation

For events where one exclusively reconstructed $`B`$ meson has been found, the vertex of the other $`B`$ decay is determined by trying to combine all other tracks in the event. Candidate tracks that form a good separate vertex, e.g., a $`K_S^0`$, are combined to form a neutral candidate, which is used instead of the daughter tracks in the vertex determination. Tracks are removed if they contribute more than 6 to the $`\chi ^2`$ of the vertex fit. Events are also rejected if $`|\mathrm{\Delta }z|>3`$ mm or if $`\sigma _{\mathrm{\Delta }z}>`$ 400 $`\mu `$m.

The time resolution function is parametrized as a sum of two Gaussians,

$$(\mathrm{\Delta }t;\widehat{a})=\underset{i=1}{\overset{2}{}}\frac{f_i}{\sigma _i\sqrt{2\pi }}\mathrm{exp}((\mathrm{\Delta }t\delta _i)^2/2\sigma _i^2).$$
(4)
The resolution parameters, $`\sigma _i`$, are taken as a scale factor, $`𝒮_i`$, times the calculated resolution based on the the tracking errors. The parameters for the second, wider, Gaussian is fixed from Monte Carlo, the parameters for the first Gaussian is determined from the combined fit for mixing and the wrong tag fractions, see Section 2.3. We also allow for a wide term, $`f_w`$, with a resolution of 1.8 ps, and no bias, to handle outliers. The parameters of the resolution function are summarized in Table 5.

### 2.3 Flavor tagging

Each event with an exclusively reconstructed $`B^0`$ decay is assigned a tag as a $`B^0`$ or a $`\overline{B}^0`$ if the rest of the event satisfies the criteria for one of several tagging categories. These tagging categories are constructed such that each event will only belong to one category. The first category uses primary leptons to determine the flavor. If the event contains an identified lepton, electron or muon, with center of mass momentum greater than 1.1 GeV/c the event is tagged as a $`B^0`$ ($`\overline{B}^0`$) if the charge of the lepton is positive (negative). The second category uses charged particles identified as Kaons. If the sum of Kaons charges is positive (negative) the event is assigned a $`B^0`$ ($`\overline{B}^0`$) tag. If the lepton tags and Kaon tags disagree no tag is assigned in these categories. The last two categories, NT1 and NT2, are assigned based on the output of a neural network. The neural network combines information about Kaons, leptons, soft pions and the stiffest track in the event to form an output that distinguishes between $`B^0`$ and $`\overline{B}^0`$ tags. The output from the neural network is shown in Figure 7. The two tagging categories are defined such that NT1 corresponds to the events which has the best separation and NT2 to the events that has slightly worse separation. The events, in the middle of Figure 7, which have very little separation are not used.

The figure of merit for each tagging category is the effective tagging efficiency, $`Q_i=ϵ_i(12w_i)^2=ϵ_i𝒟_i^2`$, where $`ϵ_i`$ is the fraction of events assigned to category $`i`$ and $`w_i`$ is the fraction that had the wrong tag assigned.

To determine the wrong tag fraction we use the sample of exclusively reconstructed hadronic $`B`$ meson decays. These decays tag the flavor of the decaying $`B`$, so by performing a combined tagging and mixing analysis we can determine the wrong tag fraction for each category. For the events in each tagging category we perform a fit to

$$_\pm (\mathrm{\Delta }t;\mathrm{\Gamma },\mathrm{\Delta }m_{B_d},𝒟,\widehat{a})=\frac{1}{4}\mathrm{\Gamma }e^{\mathrm{\Gamma }|\mathrm{\Delta }t|}[1\pm 𝒟\times \mathrm{cos}\mathrm{\Delta }m_{B_d}\mathrm{\Delta }t](\mathrm{\Delta }t;\widehat{a}),$$
(5)
where $`+`$ are unmixed events and $``$ are mixed events. A log-likelihood is formed by

$$\mathrm{ln}=\underset{i}{}[\underset{\mathrm{unmixed}}{}\mathrm{ln}_+(\mathrm{\Delta }t;\mathrm{\Gamma },\mathrm{\Delta }m_{B_d},𝒟_i,\widehat{a})+\underset{\mathrm{mixed}}{}\mathrm{ln}_{}(\mathrm{\Delta }t;\mathrm{\Gamma },\mathrm{\Delta }m_{B_d},𝒟_i,\widehat{a})]$$
(6)
where $`i`$ runs over the tagging categories. Additional terms are added to the probability density functions to describe the contributions from backgrounds, details are given in .

The results of the fit are shown in Table 3. The total tagging efficiency is $`76.7\pm 0.5`$% with an effective tagging efficiency, $`Q`$, of $`27.9\pm 0.5`$% (statistical errors only).

When the tagging algorithm is applied to the sample of 168 CP events 120 events were assigned a flavor tag. Table 4 shows a break down per mode and per tagging category of the events in the CP sample. Of the 120 events tagged, 70 were $`B^0`$ and 50 were $`\overline{B}^0`$ tags.

### 2.4 Systematics

Systematic errors were considered from many different sources; input parameters to the likelihood fit, uncertainties in the time resolution function and wrong tag fractions. The $`B^0`$ lifetime was fixed to the PDG  value $`\tau _{B^0}=1.548`$ ps, $`\mathrm{\Delta }m_{B_d}`$ was also fixed to the PDG value, $`\mathrm{\Delta }m_{B_d}=0.472\mathrm{}`$ps-1. Varying the values of these parameters give the uncertainties on $`\mathrm{sin}2\beta `$ for $`\tau _{B^0}`$ and $`\mathrm{\Delta }m_{B_d}`$ of 0.002 and 0.015 respectively.

The time resolution function was determined in a high statistics sample of fully reconstructed $`B^0`$ events. We vary the parameters in the time resolution function by 1 statistical standard deviation and assign a systematic error on $`\mathrm{sin}2\beta `$ of 0.019. To study the sensitivity on the bias in $`\mathrm{\Delta }t`$, we allowed the bias of the second Gaussian to increase to 0.5 ps. This results in a change of 0.047 on $`\mathrm{sin}2\beta `$ and is assigned as a systematic uncertainty. The sensitivity to this bias is due to the different number of events that are tagged as $`B^0`$ and $`\overline{B}^0`$.

The mistag fractions are determined from exclusively reconstructed $`B^0`$ and $`\overline{B}^0`$ mesons. Sources of systematic uncertainty come from the presence of backgrounds in these samples and possible differences between the tagging performance in the CP sample and the hadronic samples. The details about the accounting of backgrounds in the hadronic samples are given in Ref. . The systematic uncertainty on $`\mathrm{sin}2\beta `$ from the measured mistag fractions is estimated to be 0.053. A rather conservative systematic error of 0.050 on $`\mathrm{sin}2\beta `$ is assigned for a possible difference between the tagging performance between the CP sample and the exclusive hadronic sample.

The CP sample is estimated to contain a background fraction of $`(5\pm 3)`$%. Backgrounds from, e.g., $`u`$, $`d`$, and $`s`$ continuum events contributes primarily at small values of $`\mathrm{\Delta }t`$ and hence do not contribute much to the determination of $`\mathrm{sin}2\beta `$. We estimate that the effective background is 3% and correct for the background by increasing the apparent asymmetry by 1.03. A fractional systematic error of 3% is assigned on the asymmetry to cover the uncertainty in the size of the background as well as any possible CP asymmetry in the background.

The systematic errors are summarized in Table 7, and a total systematic uncertainty of 0.09 on $`\mathrm{sin}2\beta `$ is obtained.

### 2.5 Results

The analysis was carried out blind to eliminate any possible experimenters bias. The blinding technique hid both the result of the likelihood fit for $`\mathrm{sin}2\beta `$, as well as the CP asymmetry in the $`\mathrm{\Delta }t`$ distribution. The error on the asymmetry was not hidden. The value of $`\mathrm{sin}2\beta `$ was hidden by adding an arbitrary offset and flipping the sign. The CP asymmetry in the $`\mathrm{\Delta }t`$ distribution was hidden by adding an offset and, on an event by event basis, multiply $`\mathrm{\Delta }t`$ with the sign of the tag.

This allowed us to carry out many systematic studies while keeping the value of $`\mathrm{sin}2\beta `$ hidden. In particular, the whole analysis procedure, including event selection, was fixed prior to unblinding the value of $`\mathrm{sin}2\beta `$.

Using Eq. 2, we perform a likelihood fit to determine one single parameter, $`\mathrm{sin}2\beta `$. The mistag fractions, $`w_i`$, and the resolution function parameters, $`\widehat{a}`$ are taken from Tables 3 and 5 respectively. We obtain

$$\mathrm{sin}2\beta =0.12\pm 0.37\pm 0.09.$$
(7)
Table 6 shows a breakdown of the fit for different modes and tagging categories. Figure 8 shows the $`\mathrm{\Delta }t`$ distributions for $`B^0`$ and $`\overline{B}^0`$ tagged events, the different yields of $`B^0`$ and $`\overline{B}^0`$ tags is apparent in this plot as well as in Figure 9, which shows the raw asymmetry.

To validate the analysis several cross checks have been made. In particular, we have used the charmonium non-CP samples such as $`B^+J/\psi K^+`$ and the self tagging $`B^0`$ mode, $`J/\psi K^0`$, $`K^0K^+\pi ^{}`$, to establish that we do not observe any time dependent asymmetry in these modes. We also use the large samples of exclusively reconstructed charged and neutral $`B`$ mesons in this check. Table 8 summarizes the fits for an apparent CP asymmetry in these modes. The result is consistent with no time dependent asymmetry in these modes.

Other important validation measurements are the lifetime and mixing measurements , also discused in detail in Rainer Bartoldus contributions to this conference.

Figure 10 shows the constraints in the $`\rho \eta `$ plane and the preliminary BABAR measurement of $`\mathrm{sin}2\beta `$. The allowed area in the plot does not include any constraint from our measurement of $`\mathrm{sin}2\beta `$.

## 3 Conclusion

BABAR has reported preliminary results based on a dataset of 9.0 fb-1 recorded at the $`\mathrm{{\rm Y}}(4S)`$ resonance from January to July in 2000. The analysis used the modes $`B^0J/\psi K_S^0`$ and $`B^0\psi (2S)K_S^0`$ and studied the time dependent asymmetry in events where the flavor of the other $`B`$ meson was determined. The preliminary result

$$\mathrm{sin}2\beta =0.12\pm 0.37\pm 0.09$$
(8)
was obtained. An analysis of the full 2000 dataset is now in progress, other modes, including $`B^0J/\psi K_L^0`$, is to be included to further improve the precision of the measurement of $`\mathrm{sin}2\beta `$. Data collected over the next few years should significantly improve the precision of this measurement. BABAR and PEP-II will start running again in February of 2001.