text:
# Phase Transitions of Soft Disks in External Periodic Potentials: A Monte Carlo Study

## I Introduction

The liquid-solid transition in two dimensional systems of particles under the influence of external modulating potentials has recently attracted a fair amount of attention from experiments , theory  and computer simulations . This is partly due to the fact that well controlled, clean experiments can be performed using colloidal particles  confined between glass plates (producing essentially a two-dimensional system) and subjected to a spatially periodic electromagnetic field generated by interfering two, or more, crossed laser beams. One of the more surprising results of these studies, where a commensurate, one dimensional, modulating potential is imposed, is the fact that there exist regions in the phase diagram over which one observes reentrant  freezing/melting behavior. As a function of the laser field intensity the system first freezes from a modulated liquid to a two dimensional triangular solid. A further increase of the intensity confines the particles strongly within the troughs of the external potential, suppressing fluctuations perpendicular to the troughs, which leads to an uncoupling of neighboring troughs and to re-melting.

Our present understanding of this curious phenomenon has come from early mean field density functional  and more recent dislocation unbinding  calculations. The mean field theories neglect fluctuations and therefore cannot explain reentrant behavior. The order of the transition is predicted to be first order for small laser field intensities, though for certain combinations of external potentials (which includes the specific geometry studied in the experiments and in this paper) the transition may become second order after going through a tricritical point. In general, though mean field theories are applicable in any dimension, the results are expected to be accurate only for higher dimensions and long ranged potentials. The validity of the predictions of such theories for the system under consideration is, therefore, in doubt.

A more recent theory  extends the dislocation unbinding mechanism for two-dimensional melting  to systems under external potentials. For a two-dimensional triangular solid subjected to an external one-dimensional modulating potential, the only dislocations involved are those which have their Burger’s vectors parallel to the troughs of the potential. The system, therefore, maps onto an anisotropic, scalar Coulomb gas (or XY model)  in contrast to a vector Coulomb gas  for the pure $`2D`$ melting problem. Once bound dislocation pairs are integrated out, the melting temperature is obtained as a function of the renormalized or “effective” elastic constants which depend on external parameters like the strength of the potential, temperature and/or density. Though explicit calculations are possible only near the two extreme limits of zero and infinite field intensities, one can argue effectively that a reentrant melting transition is expected on general grounds quite independent of the detailed nature of the interaction potential for any two-dimensional system subject to such external potentials. The actual extent of this region could, of course, vary from system to system. In addition, these authors predict that the autocorrelation function of the Fourier components of the density (the Debye-Waller correlation function) decays algebraically in the solid phase at the transition with an universal exponent which depends only on the geometry and the magnitude of the reciprocal lattice vector.

Computer simulation results in this field have so far been inconclusive. Early simulations  involving colloidal particles interacting via the Derjaguin, Landau, Verwey and Overbeek (DLVO) potential  found a large reentrant region in apparent agreement with later experiments. On closer scrutiny, though, quantitative agreement between simulation and experiments on the same system (but with slightly different parameters) appears to be poor . Subsequent simulations  have questioned the findings of the earlier computation and the calculated phase diagram does not show a significant reentrant liquid phase.

Motivated, in part, by this controversy, in Ref.  we have recently investigated the freezing/melting behavior of a two dimensional hard disk system in an external potential. The pure hard disk system is rather well studied  by now and the nature of the melting transition in the absence of external potentials reasonably well explored. Also, there exist colloidal systems with hard interactions  so that, at least in principle, actual experiments using this system are possible. Finally, a hard disk simulation is relatively cheap to implement and one can make detailed studies of large systems without straining computational resources. By these calculations we obtained a clear signature of a reentrant liquid phase showing that this phenomenon is indeed a general one as indicated in Ref. .

In the present paper we studied a system of particles interacting by a DLVO potential in a external periodic potential, motivated on one hand whether this reentrance scenario is dependent on the range of interaction, and on the other hand to compare it with experimental results .

The phase diagram has been computed by an application of finite size scaling methods similar to the methods used in our study of the hard disk systems in external potentials .

The rest of our paper is organized as follows. In Section II we specify the model and the simulation method including details of the finite size analysis used. In Section III we present our results for both zero and non-zero external potential, in particular results for the order parameter and its cumulants with a discussion on finite size effects. We also present other quantities like order parameter susceptibility, correlation functions and heat capacity, which further illustrates the nature of the phase transitions in this system. In Section IV we discuss our work in relation to the existing literature on this subject, summarize and conclude.

## II Model and method

### A The Model

#### 1 Potentials

We study a system of N soft disks in a two dimensional box of fixed volume interacting with the DLVO pair potential $`\varphi (r_{ij})`$  between particles $`i`$ and $`j`$ with distance $`r_{ij}`$,

$$\varphi (r_{ij})=\frac{(Z^{}e)^2}{4\pi ϵ_0ϵ_r}\left(\frac{\mathrm{exp}(\kappa R)}{1+\kappa R}\right)^2\frac{\mathrm{exp}(\kappa r_{ij})}{r_{ij}},$$
(1)
where $`R`$ is the radius of the particles, $`\kappa =\sqrt{\frac{e^2}{ϵ_0ϵ_rk_BT}_in_iz_i^2}`$ is the inverse debye screening length, $`Z^{}`$ is the effective surface charge, and $`ϵ_r`$ is the dielectric constant of water. We used $`ϵ_r=78`$ and, unless otherwise indicated, $`2R=1.07\mu m`$ and $`Z^{}=7800`$. Additionally we chose a temperature of $`T=293.15K`$, and the particle density such that the particle spacing of an ideal lattice is $`a_s=2.52578\mu m`$. We then obtain different values for the reduced inverse screening length $`\kappa a_s`$ by varying $`\kappa `$ as needed. In our simulation we set $`2R`$ to be the unit length. The potential in eq. (1) mainly depends on the value of $`\kappa a_s`$, so all features found for this system should be valid also for slightly different values of the other parameters mentioned above.

In addition a particle with coordinates $`(x,y)`$ is exposed to an external periodic potential of the form:

$$V(x,y)=V_0\mathrm{sin}(2\pi x/d_0)$$
(2)
The constant $`d_0`$ in Eq.(2) is chosen such that, for a density $`\rho =N/S_xS_y`$, the modulation is commensurate to a triangular lattice of disks with nearest neighbor distance $`a_s`$: $`d_0=a_s\sqrt{3}/2`$.

The main parameters which define our system are $`\kappa a_s`$ and the reduced potential strength $`V_0/k_BT=V_0^{}`$, where $`k_B`$ is the Boltzmann constant.

#### 2 box geometry

All of the data (unless otherwise indicated) presented for $`V_0^{}<0.2`$ are obtained by a simulation in a rectangular box of size $`S_xS_y`$ ($`S_x/S_y=\sqrt{3}/2`$) and periodic boundary conditions in x- and y-direction, i.e. exactly as in . We will refer to it as ’fixed box geometry’ in the rest of the paper.

For $`V_0^{}0`$ the external potential modulates the structure of the fluid and the particles form troughs oriented in the $`y`$-direction. In order to avoid unphysical results, for $`V_0^{}0.2`$ we mainly used a box with periodic boundary conditions in x-direction and movable walls in y-direction, see Fig. 1, and we will call this ’variable box geometry’. The simulation box volume is fixed as well as the side length $`S_x`$, but in $`x`$-direction the box is divided into slabs of width $`d_0`$, centered around the minima of the external periodic potential. The wall at the end of each slab can move at most $`a_s`$ upwards or downwards around its equilibrium position: $`|b|<a_s`$, such that each slab has variable length between $`S_y2a_s`$ and $`S_y+2a_s`$ in $`y`$-direction. The averaged box geometry still is $`S_x/S_y=\sqrt{3}/2`$ as for the fixed one. The constraint for neighboring walls is to have a distance less than $`a_s/2`$: $`|c|<a_s/2`$. The walls are hard, so no particle can cross them. This is indicated as thick solid line in Fig. 1. To accommodate the particles in the box as well as possible, additional ’boundary’ particles were placed in center ($`f=d_0/2`$) behind each wall at a distance of $`e=a_s/2`$. The boundary particles interact with the particles in the box by the usual DLVO potential, but do not interact with each other. The motion of the walls is done with a Monte-Carlo procedure and keeps the volume constant, so we are still in the NVT (but variable shape) ensemble.

The movable walls are chosen to give the system additional degrees of freedom to relax internal stress. We were lead to this geometry by some unphysical results when using the fixed box geometry and higher particle numbers ($`N4096`$). For a detailed discussion see end of section III B.

### B The Method

#### 1 Numerical Details

We perform NVT Monte Carlo (MC) simulations  for the system with interactions given by Eqs. (1) and (2) for various values of $`V_0^{}`$ and $`\kappa a_s`$.

Averages $`<>`$ of observables have been obtained with the canonical measure. In order to obtain thermodynamic quantities for a range of system sizes, we analyzed various quantities within subsystems and used $`<>_L`$ to denote averages in it. The subsystems are of size $`L_x\times L_y`$ where $`L_x`$ and $`L_y`$ are chosen as $`L_y=La_s`$ and $`L_x=L_y\sqrt{3}/2=Ld_0`$ consistent with the geometry of the triangular lattice. A subbox of size $`L=3`$ as shown in Fig. 2 contains in average $`N_L=L^2=9`$ particles.

Most of the simulations described below have been done for a total system size of $`N=1024`$ and $`N=4096`$ particles, additional ones with $`N=400`$. Phase transitions have been studied in most cases by starting in the ordered solid and increasing $`\kappa a_s`$ for fixed $`V_0^{}`$. Runs where $`\kappa a_s`$ is decreased were also performed for comparison.

A typical simulation run with $`10^7`$ Monte Carlo steps (MCS) per particle (including 3 $`\times 10^6`$ MCS for relaxation) took about 50 CPU hours on a PII/500 MHz PC. In addition to ordinary (local) MC moves we also used ‘trough moves’, by which particle placements in neighboring troughs are tried. Besides producing faster equilibration, including such moves ensures that at high $`V_0^{}`$ the formation of dislocations is not artificially hindered since particles can in effect bypass each other more easily — this is very unlikely with purely local MC moves.

#### 2 Observables

The potential energy of the system per particle $`\epsilon `$ is computed by:

$$\epsilon ^{}=\frac{1}{Nk_BT}\underset{i=1}{\overset{N}{}}\left[\underset{j>i}{}\varphi (r_{ij})+V(x_i,y_i)\right]$$
(3)
and the heat capacity per particle from the fluctuations of $`\epsilon ^{}`$:

$$\frac{c_V}{k_B}=N<(\epsilon ^{}<\epsilon ^{}>)^2>$$
(4)
The nature of the fluid-solid phase transition in two dimensions has been a topic of controversy throughout the last forty years . It is well known that true long range positional order is absent in the infinitely large system due to low energy long wavelength excitations so that translational correlations decay algebraically. According to the dislocation unbinding mechanism  the two dimensional solid (with quasi long ranged positional and long ranged orientational order) first melts into a “hexatic” phase with no positional order but with quasi long ranged orientational order signified by an algebraic decay of bond-orientational correlation. A second KT transition, driven by disclination unbinding, leads to melting of the hexatic into the liquid, where both the orientational and positional order is short ranged. Therefore a useful order parameter in zero external field is the orientational order parameter. For a particle $`j`$ located at $`\stackrel{}{r}_j`$ we define the local orientational order:

$$\psi _{6,j}=\frac{1}{N_b}\underset{l=1}{\overset{N_b}{}}e^{i6\theta _{lj}}$$

where $`N_b`$ is the number of nearest neighbors, and $`\theta _{lj}`$ the angle between the axis $`\stackrel{}{r}_l\stackrel{}{r}_j`$ and an arbitrary reference axis. For the total system we use:

$$\psi _6=\left|\frac{1}{N}\underset{j=1}{\overset{N}{}}\psi _{6,j}\right|$$

and as orientational correlation function:

$$g_6(r_{ij})=|<\psi _{6,i}^{}\psi _{6,j}>|$$

In an external periodic field given by Eq.(2), however, the bond orientational order parameter is nonzero even in the fluid phase . This is because for $`V_0^{}0`$ we have now a “modulated” liquid, in which local hexagons consisting of the six nearest neighbors of a particle are automatically oriented by the external field. Thus $`<\psi _6>`$ is nonzero both in the (modulated) liquid and the crystalline phase and it cannot be used to study phase transitions in this system. The order parameters corresponding to a solid phase are the Fourier components of the (non-uniform) density $`\rho (\stackrel{}{r})`$ calculated at the reciprocal lattice points $`\{\stackrel{}{G}\}`$. This (infinite) set of numbers are all zero (for $`\stackrel{}{G}0`$ ) in an uniform liquid phase and nonzero in a solid. We restrict ourselves to the star consisting of the six smallest reciprocal lattice vectors of the two-dimensional triangular lattice. In the modulated liquid phase that is relevant to our system, the Fourier components corresponding to two out of these six vectors, viz. those in the direction perpendicular to the troughs of the external potential, are nonzero . The other four components of this set consisting of the one in the direction $`\stackrel{}{G}_1`$ (as defined for the ideal crystal in Fig. 2), and those equivalent to it by symmetry, are zero in the (modulated) liquid and nonzero in the solid (if there is true long range order). We therefore use the following order parameter:

$$\psi _{G_1}=\left|\frac{1}{N}\underset{j=1}{\overset{N}{}}\mathrm{exp}(i\stackrel{}{G}_1\stackrel{}{r}_j)\right|$$

where $`\stackrel{}{r}_j`$ is the position vector of the $`j^{th}`$ particle. The corresponding susceptibility $`\chi _{G_1}`$ is:

$$k_BT\chi _{G_1}=L^2\left[\left(\psi _{G_1}\right)^2\psi _{G_1}^2\right]$$
(5)
To measure the positional correlation, we chose the Debye-Waller correlation function which we define as follows:

$$C_{\stackrel{}{G}_1}(\stackrel{}{R})=|<e^{i\stackrel{}{G}_1(\stackrel{}{u}(\stackrel{}{R})\stackrel{}{u}(0))}>|$$

where $`\stackrel{}{R}`$ points to the elementary cell of the ideal lattice, and $`\stackrel{}{u}(\stackrel{}{R})`$ is the deviation of the actual particle position from the ideal lattice: $`\stackrel{}{r}=\stackrel{}{R}+\stackrel{}{u}(\stackrel{}{R})`$. In this case we have chosen the direction of $`\stackrel{}{R}`$ to lie along the $`y`$ axis (i.e. along the troughs of the potential). In the solid we expect this quantity do decay algebraically, i.e. $`C_\stackrel{}{G}(y)1/y^{\eta _\stackrel{}{G}}`$ , where $`\eta _\stackrel{}{G}`$ depends on the elastic constants.

$`\psi _{G_1}`$ is sensitive to the phase transition where positional order is lost. Therefore, when decreasing $`V_0^{}`$ we expect the phase boundary to converge to the corresponding transition at zero external potential. But in contrast to $`V_0^{}0`$ where the crystal is oriented by the external potential, at $`V_0^{}=0`$ it is only weakly fixed by the boundary conditions and can start to rotate, so we can not apply $`\psi _{G_1}`$ there. For this purpose we use a slightly modified positional order parameter $`\stackrel{~}{\psi }_{G_1}`$ at $`V_0^{}=0`$: the phase information of $`\psi _6`$ (of course, before taking the absolute value) is used to determine the orientation of the crystal, and then a tilted coordinate system is used to compute $`\stackrel{~}{\psi }_{G_1}`$. We applied the same method when calculating $`C_{\stackrel{}{G}_1}(y)`$ at $`V_0^{}=0`$.

We have determined phase transition points by the order parameter cumulant intersection method. The fourth order cumulant $`U_L`$ of the order parameter distribution is given by :

$$U_L(V_0^{},\kappa a_s)=1\frac{<\psi _x^4>_L}{3<\psi _x^2>_L^2}$$
(6)
In order to distuingish between the cumulants of $`\psi _6`$ and $`\psi _{G_1}`$, we denote them with $`U_{L,6}`$ and $`U_{L,G}`$, respectively. In the liquid (short ranged order) $`U_L1/3`$ and in the solid (long range order) $`U_L2/3`$ for $`L\mathrm{}`$. In case of a continuous transition close to the transition point the cumulant is only a function of the ratio of the system size $`La_s`$ and the correlation length $`\xi `$: $`U_L(La_s/\xi )`$. Since $`\xi `$ diverges at the critical point the cumulants for different system sizes intersect in one point: $`U_{L_1}(0)=U_{L_2}(0)=U^{}`$. $`U^{}`$ is a non-trivial value, i.e. $`U^{}1/3`$ and $`U^{}2/3`$. Even for first order transitions these cumulants intersect  though the value $`U^{}`$ of $`U_L`$ at the intersection is not universal any more. The intersection point can, therefore, be taken as the phase boundary regardless of the order of the transition. This is useful since the order of the melting transition in $`2D`$ either in the absence  or with  external potentials is not unequivocally settled. And there is also another aspect in our special case: since the positional order correlation is predicted to decay algebraically in the solid phase (quasi long ranged order), the whole solid can be seen as consisting of a line or area of critical points with temperature- dependent critical indices $`\eta _\stackrel{}{G}(T)`$ . In that case we expect the cumulants to merge at a nontrivial value at the onset of the solid phase instead of intersecting, yielding a line of intersections. We indeed observed this behavior for hard disks in high external potentials . For the same reasons the same behavior is expected for the $`\psi _6`$-cumulants at the liquid-hexatic transition and in the hexatic phase. Also the very similar 2d-XY-spin model shows this behavior when using the magnetization as order parameter .

Note that though the order parameter $`<\psi _{G_1}>`$ decays to zero with increasing system size even in the 2-d solid (assuming quasi long ranged order there), the cumulants will stay at the non-trivial value regardless of $`L`$. So for $`L\mathrm{}`$ there should be a jump from $`1/3`$ to this nontrivial value when crossing the phase boundary from liquid to solid, which underlines the usefulness of cumulants.

In order to map the phase diagram we systematically vary the system parameters $`V_0^{}`$ and $`\kappa a_s`$ to detect order parameter cumulant intersection- or merging points which are then identified with the phase boundary.

## III Results and Discussion

### A Zero external potential

In this section we analyse the system properties for zero external field. In particular we present results for the order parameter, the cumulants, the correlation functions and the heat capacity for different values of $`\kappa a_s`$. In these studies we used the fixed box geometry.

In Fig. 3 the cumulant of the $`\psi _6`$ order parameter versus $`\kappa a_s`$ is shown for different subsystem sizes. We identify the phase transition value of $`\kappa _6a_s`$ at about 14.42 by locating the cumulant intersection point. Since the positional order is not well defined in two-dimensional systems, the positional order parameter $`\stackrel{~}{\psi }_{G_1}`$ shows a strong system-size dependency, see Fig. 4. The cumulants of $`\stackrel{~}{\psi }_{G_1}`$ intersect at a value of $`\kappa _Ga_s14.25`$, which is slightly smaller than $`\kappa _6a_s`$. This is in agreement with a KTHNY two stage melting scenario, in which the solid and the fluid phase are separated by a “hexatic” region in the phase diagram, in which the positional order is short ranged and the bond-orientational order is long ranged. Both effects are detected by the two order parameters, $`\psi _6`$ being sensitive for the bond-orientational order and $`\stackrel{~}{\psi }_{G_1}`$ on the positional order. Surprisingly, however, though in the case of the hexatic phase one expects the $`\psi _6`$-cumulants to coincide, they obviously don’t in Fig. 3 (see also ).

The Debye-Waller correlation functions $`C_{G_1}(y)`$ for different values of $`\kappa a_s`$ are shown in Fig. 5. At the transition value $`\kappa _Ga_s=14.25`$ for $`N=4096`$ we find a power law dependency of $`C_{G_1}(y)`$ from $`y`$ with an exponent $`\eta _{G_1}0.28`$ which is well within the predicted range of $`[1/4,1/3]`$. In Fig. 6 the orientational correlation function versus distance is shown. This function reveals a power law dependency of the bond-orientation correlations at $`\kappa _6a_s`$, the exponent value is about 1/4. The value of the exponent $`\eta _6`$ at the transition has been predicted by the KTHNY theory to be $`1/4`$ , which is in agreement with our results.

In Fig. 7 (left) we present the heat capacity data versus $`\kappa a_s`$ for different system sizes. It is obvious that the heat capacity does not show a singularity as would be expected in case of a first order or a conventional second-order transition. The peak maxima are not very sharp, but are roughly located close to the value of $`\kappa _6a_s`$, where the $`\psi _6`$-order parameter cumulants intersect. The peak maxima thus do not agree with the cumulant intersection point of the $`\stackrel{~}{\psi }_{G_1}`$-order parameter, which is again in agreement with the KTHNY scenario. We note that the identification of the phase transition point by the heat capacity maxima may result in misleading results on the location of the transition points in the phase diagram. In particular, for a smaller system of N = 400 particles we find that configurations wherein the crystal is rotated by a tilt angle of $`\alpha =\pm 30^{}`$ ($`\alpha `$ extracted from the phase information of $`\psi _6`$) may be present. Since this is incompatible with the box geometry it leads to a higher energy and a lower measured $`\psi _{G_1}`$. The value of $`\stackrel{~}{\psi }_{G_1}`$, on the other hand, is not appreciably altered. This is shown in the time evolution of that system in Fig. 8. We also show the configuration of a tilted crystal with $`\alpha =30^{}`$ (3.000.000 MCS) in Fig. 9 (left), and of an ’correctly’ aligned crystal ($`\alpha =0^{}`$, 7.700.000 MCS) (right).

### B External periodic potential

In this section we analyse the system properties in the presence of a periodic external potential. The studies in this section are mainly done with variable box geometry. Comparative studies with fixed box geometry show that the new method does not lead to artificial features, but rather gives improved results. For more details see the discussion of the Debye-Waller correlation function at the end of this section.

Two examples for the $`\kappa a_s`$-dependency of the $`\psi _{G_1}`$-order parameter and the cumulants for an external potential amplitude of $`V_0^{}=2`$ and $`V_0^{}=1000`$ are shown in Figs. 10 and 11. We note that the cumulant intersection, which can be clearly identified for $`V_0^{}=2`$ in Fig. 10 is developing towards an intersection “line” for $`V_0^{}=1000`$, a behavior which was found in case of the hard disk system in external potentials  as well as in related systems with a KT transition like the XY-spin-model . Another example is shown in Fig. 12. There $`\kappa a_s`$ is kept fixed at $`15.3`$ and $`V_0^{}`$ is varied. The starting point at $`V_0^{}=0.2`$ is in the mod. liquid phase, crosses slightly the solid (’laser induced freezing’), and re-enters the mod. liquid at higher $`V_0^{}`$ (’reentrance’). This is already a first sign of a ’reentrant’ phase transition scenario. For the phase diagram these cumulant intersection- or merging values were used.

In Fig. 13 the cumulant intersection values are shown as a function of $`V_0^{}`$ for fixed box geometry. We observe that $`U^{}`$ is not an universal number but, nevertheless, goes to a limiting value for large $`V_0^{}`$ .

The amount of hysteresis effects on the location of the transition point has been analyzed for the case of $`V_0^{}=2`$ by a time consuming reverse density quench simulation, in which a path in phase space was chosen in the direction opposite to the standard path. The results of this study are shown in Fig. 14. Comparing these results with the ones of Fig. 10 reveals quite close agreements showing that only small hysterese effects are present in the system.

The $`\psi _{G_1}`$ order parameter susceptibilities $`\chi _{G_1}`$ are shown in Fig. 15 versus $`\kappa a_s`$ for different system sizes. We note that close to the transition a maximum develops, the value of the maximum increasing with the system size. In Fig. 16 the susceptibilities of the largest subsystems ($`L=32`$) as functions of $`\kappa a_s`$ are compared for different values of $`V_0^{}`$. Clearly the peak position for $`V_0^{}=2`$ is shifted to larger values compared to the cases with $`V_0^{}=0.2`$ and $`V_0^{}=1000`$. This feature is another sign of a “reentrant” phase transition scenario in the phase diagram. Compared to the cumulant intersection values, $`\chi _{G_1}`$ maxima are located at slightly higher $`\kappa a_s`$. This may be due to finite size effects, which often show the feature that phase transition points in finite systems are shifted to slightly different values depending on the observable under investigation. In particular one expects (and we get) a shift towards parameter values in the disordered region (here a liquid, i.e. higher $`\kappa a_s`$) for the average order parameter and the susceptibility.

In Fig. 7 (right) the heat capacity for $`V_0^{}=2`$ and different system sizes is shown. The peak is nearly independent from system size and shifted towards the liquid phase with respect to the order parameter cumulant intersection value, i.e. the same behavior as for $`V_0^{}=0`$.

The advantage of the variable box geometry, especially for large systems, can be seen best by looking at the Debye-Waller correlation function. In Fig. 17 (left) an example is shown for fixed box geometry. The crossing from the solid phase with an algebraic decay to the mod. liquid with exponential decay is not monotonic, but at $`\kappa a_s=15.7`$, $`C_{G_1}(y)`$ drops to zero at $`y=S_y/2`$, which is not physically meaningful. At a higher value $`\kappa a_s=16`$ it rises, and then falls again at $`\kappa a_s=16.4`$, showing a exponential decay as expected. In variable box geometry this feature doesn’t show up, see Fig. 17 (right). Here we have a smooth transition from solid- to liquid-like behavior. We explain this strange behavior at $`\kappa a_s=15.7`$ in fixed geometry as follows: consider a system with $`N=10000`$ particles. Without dislocations there will be an ideal lattice with $`N_t=100`$ particles in each of the $`100`$ troughs. Assuming dislocation unbinding as melting mechanism, consider the existence of some dislocations in the system with opposite burgers vectors $`\stackrel{}{b}=\pm a_s\stackrel{}{e}_y`$. One of these dislocations increases the number of particles in the troughs by one, while the other decreases it by one. We can have for example a situation where half of the troughs has $`100`$ particles, and the other half has either $`101`$ or $`99`$ particles. If we now for simplicity assume that the distance between particles in a row is $`a=L_y/N_t`$, calculating the pair correlation function $`g(y)`$ along a trough will yield two peaks around $`y=L_y/2`$: one centered exactly at $`y=L_y/2`$ from the troughs with $`N_t=100`$ ($`a=a_s`$), and the other centered at $`y=L_y/2+a_s/2`$ due to the troughs with $`N_t=99`$ or $`N_t=101`$ ($`aa_s`$). As consequence, $`C_{G_1}(y=L_y/2)`$ will be zero. The same situation in the movable-walls geometry will not show these problems: the troughs with $`101`$ particles can expand a bit, while the ones with $`99`$ particles can contract. Now in every trough is $`a=a_s`$, and $`C_{G_1}(y=L_y/2)`$ is not necessarily zero. Also, the formation of a dislocation pair costs less energy and is closer to the true infinite system value. The discussion above is in some sense similar to the 2d-XY-spin model with a vortex in the center: in an infinite sample, the spins to the left and to the right have opposite spin directions, but periodic boundary conditions in a finite system will try to align them, so that the formation of the vortex is disturbed. Free boundary conditions won’t cause this problem.

In zero external potential with fixed geometry the formation of a dislocation pair is not so problematic, since the particles are not forced into troughs and have more degrees of freedom in movement. In the above example one end of a line of $`N_t=101`$ particles could make a slight shift in $`x`$-direction to access more space in $`y`$-direction.

However, at the first data point in the solid closest to the transition we find an algebraic decay $`C_{G_1}(y)1/y^{\eta _{G_1}}`$ with $`\eta _{G_1}`$ in the range of $`0.25\mathrm{}0.34`$ for $`V_0^{}=0\mathrm{}1000`$ and $`N=1024`$ particles. In  $`\eta _{G_1}`$ is predicted to be $`1/4`$ at the transition.

### C The Phase Diagram

For each $`\kappa a_s`$ and $`V_0^{}`$ value we computed cumulants $`U_{L,G}`$ for a range of subsystem sizes $`L`$ and located intersection- or merging points which we identify with the phase boundary. We have obtained a detailed phase diagram for $`N=1024`$ particles which is shown in Fig. 18 for fixed box geometry (left) and for variable box geometry (right). We want to emphasize that there are only slight differences and the general shape of the phase diagram is the same for both box geometries. At $`V_0^{}=0`$ also the $`\psi _6`$ cumulant intersection value is plotted for comparison. The values of $`\kappa a_s`$ at the transition initially rise and subsequently drop as $`V_0^{}`$ increases. The maximum $`\kappa _Ga_s`$ values are found for $`V_0^{}12`$. These transition points separate a high density solid from a low density modulated liquid. Thus, at a properly chosen $`\kappa a_s`$, we observe an initial freezing transition followed by a reentrant melting at a higher $`V_0^{}`$ value. Such an effect had been found earlier in experiments on colloidal systems in an external laser field .

In order to quantify residual finite size effects on the phase diagram, we have computed the transition points for different total system sizes. The resulting phase diagrams are shown in Fig. 19, again for fixed (left) and variable box geometry (right). We note that with increasing system size in fixed box geometry all transition points are slightly shifted to the solid region, whereas for the variable box this shift is towards the liquid region. One can see that the shift for the fixed box is much smaller at low $`V_0^{}`$, and for the variable box it is much smaller at medium and high $`V_0^{}`$. We also found that the cumulant intersection point smears out strongly if using the fixed box, $`N4096`$ and higher $`V_0^{}`$, probably for the same reasons as those mentioned in the discussion of the Debye-Waller correlation function. In the variable box there was no such problem. These features were the reason for us to use mainly the fixed box for $`V_0^{}<0.2`$ and the variable box for $`V_0^{}0.2`$. By the way, the same ’Debye-Waller problem’ also occurs when simulating hard disks in external periodic potentials, $`N4096`$ and $`V_0^{}`$ medium or high, and can be solved again by using the variable box geometry.

However, for all system sizes the structure of the phase diagram with a pronounced minimum at intermediate values of $`V_0^{}`$ is not affected by the shifts.

The difference in the value of $`\kappa _Ga_s`$ at the transition between the infinite and zero external potential cases, we find $`\kappa _Ga_s(V_0^{}=\mathrm{})\kappa _Ga_s(V_0^{}=0)0.82`$. This is not far away from $`0.608`$ which is a value predicted by .

We have also done simulations with slightly altered parameters, i.e using particles with diameter $`2R=3\mu m`$, effective surface charge $`Z^{}=20000`$ and $`a_s=8\mu m`$, to match the experiments in . As expected, we only observe a slightly shift of the phase diagram of $`\mathrm{\Delta }(\kappa a_s)0.35`$ towards higher values of $`\kappa a_s`$. The experimental phase diagram  qualitatively has the same shape as our results, but shows larger freezing- and reentrance regions and is shifted to higher values of $`\kappa a_s`$ at about $`\mathrm{\Delta }(\kappa a_s)4.5`$ in average. The reasons for these differences are probably due to the particle interaction. We only use pairwise interaction, which is a good approximation for low particle densities . But for higher particle densities many body interactions play a role because of macroion screening, which results in an effective pair potential that has considerable deviations  from a pure Yukawa-like potential like ours. In particular, there could be an attractive part.

### D Scaling behavior

We next try to determine the order of the phase transitions encountered in this system for two values of $`V_0^{}`$. In order to investigate this issue we studied the scaling behavior of the order parameter, susceptibility and the order parameter cumulant near the phase boundary for a small ($`2`$) and a large ($`1000`$) $`V_0^{}`$.

From finite size scaling theory (for an overview see Ref. ) we expect these quantities to scale as  :

$$<\psi _{G_1}>_LL^bf(L/\xi )$$
(7)
$$\chi _Lk_BTL^cg(L/\xi )$$
(8)
$$U_Lh(L/\xi )$$
(9)
Here $`b=\beta /\nu `$, $`c=\gamma /\nu `$ (for critical scaling), and $`f`$, $`g`$, $`h`$ are scaling functions. Defining $`\stackrel{~}{\kappa }=(\kappa a_s\kappa _Ga_s)/\kappa _Ga_s`$, we expect the correlation length $`\xi `$ to diverge as $`\xi \stackrel{~}{\kappa }^\nu `$ for an ordinary critical point, while for a KT-transition we have an essential singularity and $`\xi \mathrm{exp}(a\stackrel{~}{\kappa }^{\stackrel{~}{\nu }})`$ when approaching the transition from the liquid side.

In Fig. 20 we have plotted the left hand sides of Eqs.(7), (8) and (9) versus $`L/\xi `$ for $`V_0^{}=1000`$, where data points of the variable box geometry for $`15.2\kappa a_s16.0`$ have been considered and $`\kappa _Ga_s=15.1`$, obtained by cumulant intersection. In order not to introduce an unwarranted bias, we have separately considered ordinary critical scaling (left column) and a KT scaling form (right column) and adjusted the values of the parameters $`b`$, $`c`$ and $`\nu `$, or $`a`$, $`b`$, $`c`$ and $`\stackrel{~}{\nu }`$, till we obtained collapse of our data onto a single curve determined by a least square estimator. <sup>*</sup><sup>*</sup>*One remark concerning the KT-scaling: the errors in $`\stackrel{~}{\nu }`$ and $`a`$ are relatively big because for example increasing $`\stackrel{~}{\nu }`$ and decreasing $`a`$ by an appropriate amount resulted in a nearly equally good collapse. Good collapse of our data is observed for both scaling forms, the numerical values for $`\stackrel{~}{\nu }`$, $`2b=\eta `$ and $`c=2\eta `$ for KT scaling ($`2b0.28`$, $`c1.70`$, $`\stackrel{~}{\nu }0.37`$) are relatively close to the predicted values  ($`2b=\eta =1/4`$, $`c=1.75`$, $`\stackrel{~}{\nu }=0.5`$). The situation is similar for small $`V_0^{}=2`$, the quality of the collapse comparable to the one of $`V_0^{}=1000`$. The critical parameters were obtained in this case for values $`15.4\kappa a_s16.2`$, with $`\kappa _Ga_s=15.37`$.

We have a good internal consistency between $`\eta =0.25\mathrm{}0.33`$ extracted from the Debye-Waller correlation function, and the values obtained from data collapsing: $`\eta =0.28`$ for $`V_0^{}=1000`$, and $`\eta =0.36`$ for $`V_0^{}=2`$. Our results for the numerical values of the parameters are summarized in Table I.

A more precise classification of the phase transitions with the present data and system sizes is not easy. This topic is left for future work, in particular we plan to compute the elastic properties of the system by a method recently developed for the hard disk system  and to test the KT predictions .

## IV Summary and conclusion

In summary, we have calculated the phase diagram of a two dimensional system of soft disks, interacting via a DLVO potential, in an external sinusoidal potential. We find freezing followed by reentrant melting transitions over a significant region of the phase diagram in tune with results on hard disks , previous experiments on colloids  and with the expectations of a dislocation unbinding theory . One of the main features of our calculation is the method used to locate phase boundaries. In contrast to earlier simulations  which used either the jump of the order parameter or specific heat maxima to locate the phase transition, we used the more reliable cumulant intersection method. It must be noted that the specific heat in this system does not show a strong peak at the phase transition density so that its use may lead to confusing results. This, in our opinion, may be the reason for part of the controversy in this field. It is possible that earlier simulations which used smaller systems and no systematic finite size analysis may have overlooked this feature of $`c_V`$ which becomes apparent only in computations involving large system sizes. We have shown that finite size scaling of the order parameter cumulants as obtained from subsystem or subblock analysis, on the other hand, yields an accurate phase diagram.

What is the order of the phase transition? We know that  for the pure hard disk system in two dimensions this question is quite difficult to answer and our present understanding  is that this system shows a KTHNY transition. In our system of soft disks, for zero external potential we can rule out a strong first order transition, although smaller systems show a feature (double peak in the internal energy) which mimics such a behavior. We find several features which are consistent with the KT theory, but also one which is not. Upon turning on the external periodic potential, the difference between hexatic and liquid disappears, and an (anisotropic) KT transition  from the modulated liquid into the solid is expected. Our results show several features which suggest that this is what we have, but there are still some (not large) deviations from theory. Though we have discussed these observations in the rest of the paper, we list the important ones below for clarity:

* The behavior of the cumulants near the transition is similar to an earlier work  on the XY system which shows a KT transition.
* The specific heat is relatively featureless and does not scale with system size in a fashion expected of a true first order or conventional continuous transition.
* The decay of the correlation functions is similar to what is predicted  for an anisotropic scalar Coulomb gas.
* For two test values of $`V_0^{}`$, the scaling of the order parameter, the susceptibility and the cumulant may be reasonable described by the KT theory.

Of course, in order to resolve this issue unambiguously yet larger simulations are required. Also, we need to compute elastic properties  of this system in order to compare directly with the results of Ref. . Work along these lines is in progress.

## V Acknowledgment

We are grateful for many illuminating discussions with C. Bechinger and K. Binder. One of us (S.S.) thanks the Alexander von Humboldt Foundation for a Fellowship. Support by the SFB 513 and granting of computer time from the NIC and the HLRS is gratefully acknowledged.

|  | $`\kappa _Ga_s`$ | $`b`$ | $`c`$ | $`\nu `$ | $`b`$ | $`c`$ | $`\stackrel{~}{\nu }`$ | $`a`$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $`V_0^{}=1000`$ | 15.1 | 0.130(12) | 1.61(8) | 1.5(2) | 0.14(1) | 1.70(3) | 0.37(6) | 1.45(40) |
| $`V_0^{}=2`$ | 15.37 | 0.163(15) | 1.68(5) | 1.51(25) | 0.18(3) | 1.74(3) | 0.40(5) | 1.2(3) |
| KT theory |  |  |  |  | 0.125 | 1.75 | 0.5 | $`O(1)`$ |

Table I Parameters in the scaling plots (see Fig. (20)) for $`V_0^{}=2`$ and $`V_0^{}=1000`$. The first three parameter columns are for critical scaling, the last four for KT scaling. The last line shows the predictions of KT theory.


format: text

text:
# 1 Introduction

## 1 Introduction

Recently supertubes, as special tubular D2-branes, attracted considerable attention owing to a previously unknown mechanism of stabilisation by angular momentum as a result of crossed constant electric and magnetic fields . In particular it was shown in ref.  that the finite radius of the tube results from the product of F-string and D0-brane charges, and the BPS energy is reached by a cancellation of oppositely directed Poynting and centrifugal momenta, thus requiring a nonzero angular momentum along the axis of the supertube. It was also shown that another supersymmetry preserving configuration with the same energy exists in the limit of vanishing D2-brane radius with the D0-branes aligning like beads along a IIA superstring, with vanishing angular momentum. In these two cases the directions of the $`E`$ and $`B`$ components of the electromagnetic field are reversed. It is evident that such a seemingly simple model immediately suggests its investigation with respect to electric-magnetic  and weak-strong or S dualities . Several subsequent investigations  did not venture in this direction, and this is therefore the main objective here. Earlier investigations  studied the stability of branes in relation to their possible collapse to lower dimensional ones  or the opposite, the expansion of branes into higher dimensional ones  with further investigations in refs.. In view of apparent physical shortcomings of the case of D2 (e.g. the Poynting momentum is a scalar), it is natural to prefer the case of D3 for the investigation since here the connection between electric-magnetic and weak-strong dualities is exposed in the fullest physical context, and at least weak-strong duality of this case has been considered from time to time . In fact we shall show that weak-strong duality may be understood as the electric-magnetic duality of nonlinear electrodynamics.

Thus, in order to realise this idea, we proceed in a number of definite steps. We first demonstrate how quantised electric and magnetic charges can be connected with the numbers of F-strings and D-strings dissolved in the D3-brane. We then show that the F-string is a D3-brane with minimal electric charge. In the third step this D3-brane with minimal electric charge is mapped to the D3-brane with arbitrary electric and magnetic charges using the electric-magnetic duality rotation. Then, following this, we show that the D3-brane with arbitrary charges is a composite string in the sense of Witten which can also be described as an FD string. Finally the F-string is mapped to this composite string.

## 2 The Action

The D3-brane action is

$$I=T_3d^4\xi \sqrt{det(g_{\mu \nu }+2\pi F_{\mu \nu })},T_3=\frac{1}{(2\pi )^3g_s},$$
(1)
where $`g_{\mu \nu }`$ is the induced metric tensor and $`F_{\mu \nu }`$ is the electromagnetic field strength tensor. The dilaton field is taken to be constant and target spacetime to be flat. We consider a D3-brane of cylindrical shape $`^1𝕊^2`$ which does not allow a variation of the radius $`r`$ in time and along the axis. This choice is conditioned by our intention to investigate the brane’s dualities. The variable radius violates the electric-magnetic duality invariance and therefore we keep it constant. The brane’s worldvolume is parametrized in terms of the variables $`(t,z,\theta ,\phi )`$ and the cylindrical hypersurface is given by

$`X^0=t,X^1=z,X^2=r\mathrm{cos}\theta ,`$
$`X^3=r\mathrm{sin}\theta \mathrm{cos}\phi ,X^4=r\mathrm{sin}\theta \mathrm{sin}\phi ,others=constant.`$ (2)
The induced metric is

$$ds^2=dt^2+dz^2+r^2d\theta ^2+r^2\mathrm{sin}^2\theta d\phi ^2.$$
(3)
Hereafter we use bold face letters for the 3-dimensional vectors and symbols $`i,j,k`$ for the space indices $`(z,\theta ,\phi )`$. The electric intensity $`𝐄`$ and magnetic induction $`𝐁`$ are defined by

$$𝐄_i=2\pi F_{0i},𝐁_i=2\pi \frac{ϵ_{ijk}}{2}F^{jk}.$$
(4)
The Levi-Civita form $`ϵ_{ijk}`$ in curved space is defined by $`ϵ_{ijk}=\sqrt{g}\epsilon _{ijk}`$, where $`gdet(g_{\mu \nu })=r^4\mathrm{sin}^2\theta `$ and $`\epsilon _{ijk}`$ is the Levi-Civita form in flat space.

The action reduces to

$$I=T_3\sqrt{g}𝑑t𝑑z𝑑\theta 𝑑\phi ,$$
(5)
where the Lagrangian $``$ is

$$=\sqrt{1+𝐁^2𝐄^2(𝐄𝐁)^2},$$
(6)
the brane tension is a normalization factor and the rest is the worldvolume invariant measure. We started from the D3-brane action and obtained the Non-Linear Electrodynamic(NLE) action. The reason is that we suppressed all other degrees of freedom and kept only the gauge field.

The aspects of the NLE action we are interested in are most conveniently analyzed by introducing the electric induction $`𝐃`$,

$$𝐃=\frac{}{𝐄}=\frac{𝐄+𝐁(𝐄𝐁)}{\sqrt{1+𝐁^2𝐄^2(𝐄𝐁)^2}}.$$
(7)
Performing a Legendre transform we construct the Hamiltonian $``$ as a function of $`𝐁`$ and $`𝐃`$,

$$=\sqrt{1+𝐁^2+𝐃^2+|𝐁\times 𝐃|^2}.$$
(8)
It has $`SO(2)`$ electric-magnetic duality rotation invariance

$`𝐃`$ $``$ $`𝐃\mathrm{cos}\alpha 𝐁\mathrm{sin}\alpha ,`$
$`𝐁`$ $``$ $`𝐃\mathrm{sin}\alpha +𝐁\mathrm{cos}\alpha .`$ (9)
There is Legendre discrete duality at $`\alpha =\pi /2`$ which interchanges the electric and magnetic inductions. The $`SO(2)`$ invariance of the Hamiltonian can be extended to $`SL(2,)`$ invariance by including a scalar dilaton field $`\varphi `$ and a pseudo-scalar axion field $`\chi `$. For the time being we omit these fields for simplicity and return to this later. To complete this section we introduce the magnetic intensity $`𝐇`$

$$𝐇=\frac{}{𝐁}=\frac{𝐁𝐄(𝐄𝐁)}{\sqrt{1+𝐁^2𝐄^2(𝐄𝐁)^2}}$$
(10)
and write out the equations of motion in the static case:

$$\epsilon ^{ijk}\frac{𝐄_i}{x^j}=0,\epsilon ^{ijk}\frac{𝐇_i}{x^j}=0,$$
(11)
$$\frac{}{x^k}\sqrt{g}𝐃^k=0,\frac{}{x^k}\sqrt{g}𝐁^k=0.$$
(12)
The transformations (9) can be equally described in terms of the intensities $`𝐄`$ and $`𝐇`$. They take the solution of the NLE action to another solution. Note that the duality transformations (9) are invariances of the field equations (11,12) and Hamiltonian (8) but not of the action (5).

## 3 Brane Charges

In this section we consider two static solutions of the NLE action (5). In the first case electric and magnetic fields are parallel, the brane carries no angular momentum generated by crossed electric and magnetic fields and collapses to the bound state of the strings. In the second case electric and magnetic fields are perpendicular and the nonzero Poynting vector prevents the collapse of the brane. The flux in the D3-brane acts as a source for the strings. It comes from an integral over the whole world-volume and it is impossible to localise the strings at a particular place in the brane’s worldvolume. But they can appear as remnants when the cylindrical brane is squeezed to a line.

Using Dirac charge quantization we express the energy of the solutions as a function of corresponding quantized electric and magnetic charges $`n`$ and $`m`$. In the limit of the vanishing radius $`r`$ of the brane, $`r0`$, the term proportional to the area of the surface of the cylinder $`(r^2)`$ in the expression of the energy can be neglected in leading order. Then the rest shows the strings smeared in the brane. In this way we establish that $`n`$ and $`m`$ are F-string and D-string charges of the brane, respectively .

### Bound States of Strings

First we consider parallel uniform electric and uniform magnetic fields aligned with the $`z`$ axis, the 2-form field strength being given by

$$2\pi F=Edtdz+Br^2\mathrm{sin}\theta d\theta d\phi .$$
(13)
The Dirac electric charge quantization condition is

$$n=\underset{S^2}{}\frac{\delta I}{\delta F_{0i}}𝑑s_i=2\pi ϵ_{ijk}\frac{\delta I}{\delta 𝐄_i}𝑑x^j𝑑x^k,$$
(14)
where $`n`$ is an integer. The surface perpendicular to the electric field is the basesurface of the cylinder and is parametrized in terms of the variables $`(\theta ,\phi )`$. Then

$$T_3D\sqrt{g}𝑑\theta 𝑑\phi =\frac{n}{2\pi };D=𝐃^z.$$
(15)
The measure of the basesurface area $`S_2`$ is

$$dS_2=r^2\mathrm{sin}\theta d\theta d\phi =\sqrt{g}d\theta d\phi ,$$
(16)
so that eq.(15) takes the form

$$T_3D𝑑S_2=\frac{n}{2\pi }.$$
(17)
Eq.(17) implies that the electric induction flux through the basesurface is quantized, and as $`D`$ is constant it yields

$$D=\frac{\pi ng_s}{r^2}.$$
(18)
This result can be also derived from M-theory .

The magnetic flux quantization condition is

$$F_{\theta \phi }𝑑\theta 𝑑\phi =2\pi m$$
(19)
where $`m`$ is an integer. It yields

$$B=\frac{\pi m}{r^2}.$$
(20)
The substitution of (18) and (20) into Hamiltonian (8) and the integration over $`\theta ,\phi `$ gives for the energy $``$

$$=\sqrt{(T_3S_2)^2+(nT_s)^2+(mT_1)^2}𝑑z.$$
(21)
where $`T_s=1/2\pi `$ and $`T_1=1/2\pi g_s`$ are the tensions of the F-string and D-string, respectively. To uncover the dissolved strings we consider the energy $``$ in the limit $`S_20`$. In pure electric case it can be presented as the mass of $`n`$ fundamental strings

$$(m=0)_{r0}=nT_s𝑑z.$$
(22)
From this asymptotoic relation we conclude that there are $`n`$ fundamental strings dissolved in the worldvolume of the D3-brane, and $`n`$ is the F-string charge of the D3-brane. In the pure magnetic case the energy becomes the mass of $`m`$ D-strings

$$(n=0)_{r0}=mT_1𝑑z,$$
(23)
with the obvious interpretation that there are $`m`$ D-strings dissolved in the D3-brane, and $`m`$ is the D-string charge of the D3-brane.

The energy (21) is minimal when $`r=0`$, as it should be. The brane shrinks to a line with the energy

$$=\sqrt{(nT_s)^2+(mT_1)^2}𝑑z.$$
(24)
This is the bound state of $`(n,m)`$ strings. The D3-brane with charges $`(n,m)`$ and with no net D3-brane charge collapses to the $`(n,m)`$ strings and therefore the latter can be considered as a collapsed D3-brane. The mass formula (24) for the $`SL(2,)`$ multiplet was derived in . The basic condition was the fact that the string solutions are BPS-saturating states and preserve half of supersymmetry. Here the same relation is derived from the viewpoint of the supersymmetric D3-brane. The main point is that the string bound state is considered as a collapsed superbrane. This method allows to apply Legendre discrete duality to the $`(n,m)`$ strings. The latter interchanges the S-dual partners, i.e. F-strings and D-strings, resulting in $`(m,n)`$ strings.

### Expanded FD strings

Eq.(11,12) have nontrivial solutions with the perpendicular electric and magnetic fields. Suppose we have the electric field along the $`z`$ axis and the magnetic field along meridians, so only $`𝐄_z,𝐃_z,𝐁_\theta ,𝐇_\theta `$ components are nonzero. The eq.(12) gives a simple solution for $`𝐁_\theta `$

$$𝐁_\theta =\frac{B}{\mathrm{sin}\theta },$$
(25)
where $`B`$ is an integration constant. The $`𝐄_z`$ components can be expressed via $`𝐃_z`$ and $`𝐁_\theta `$

$$𝐄_z^2=\left(1+\frac{B^2}{r^2\mathrm{sin}^2\theta }\right)\frac{𝐃_z^2}{1+𝐃_z^2}.$$
(26)
From the eq.(11) it follows that $`𝐄_z`$ can not depend on the variable $`\theta `$. The r.h.s of the eq.(26) satisfies this restriction only when

$$𝐃_z=D\mathrm{sin}\theta ,$$
(27)
where the constant $`D`$ is defined by

$$D^2B^2=r^2.$$
(28)
Then eq.(26) becomes $`𝐄_z^2=1`$ and we rewrite it in the invariant form

$$𝐄^2=1.$$
(29)
This is the critical value of the electric intensity which forces the critical value of the magnetic intensity $`𝐇^2=1`$. Now we can make $`SO(2)`$ duality rotations and obtain more general solutions. As $`𝐄`$ and $`𝐇`$ are unit mutually perpendicular vectors, they will stay such after rotations (9). Consequently the condition (29) is right for any rotated solution and defines the true minimum of the NLE action. It has another interesting interpretation. The question whether the tubular D2-brane with the nonzero Poynting vector is supersymmetric has been clarified . The condition for the preservation of 1/4 of the supersymmetry was reduced to the eq.(29). We come to the conclusion that the classical minimum of the NLE action is a BPS-saturating state.

The charge quantization requires the compact length of the variable $`z`$. Therefore we assume that the $`z`$ direction is compactified by a circle with the circumference $`l`$. We denote by $`R`$ the compactification radius, then $`l=2\pi R`$. The quantization repeats the procedure of the previous case resulting in

$$D=\frac{4ng_s}{r^2},B=\frac{m}{R}.$$
(30)
The energy expression takes the form

$$=\sqrt{(rT_3S_2^{})^2+(4nT_sR)^2\mathrm{sin}^2\theta +\frac{(mT_1r)^2}{\mathrm{sin}^2\theta }+\frac{4n^2m^2}{\pi ^2r^2}}\mathrm{sin}\theta d\theta .$$
(31)
where $`S_2^{}=2\pi rl`$ is the azimuthal lateral area of the brane. The first term under the square root is the energy of the brane with no electromagnetic flux. The second term presents $`n`$ dissolved F-strings which is easy to see by setting $`m=0`$ and taking the limit $`r0`$.

$$(m=0)_{r0}=nT_sl.$$
(32)
In the pure magnetic case the energy is

$$(n=0)=mT_1\sqrt{\left(\frac{Rr^2}{m}\right)^2\mathrm{sin}^2\theta +r^2}𝑑\theta $$
(33)
In the limit $`r0`$ this becomes

$$(n=0)_{r0}=mT_1(\pi r).$$
(34)
This is a bunch of m D-strings with the length of a semicircle. Therefore $`m`$ is the brane’s D-string charge. The integrand in (33) represents the length measure of an elliptic arc along which the energy density is distributed and the ellipse is mapped into the azimuthal lateral face of $`^1𝕊^2`$ becoming an elliptic helix with a pitch $`r`$.

The last term in (31) is generated by the Poynting vector and owing to it the energy grows infinitely when $`r0`$. The nonzero Poynting vector $`P(Pnm)`$ stabilizes the brane at nonzero radius $`r_0`$. The substitution of (30) into (28) determines the value of the radius $`r_0`$

$$r_{0}^{}{}_{}{}^{3}=\frac{4\pi g_s}{R}|nm|.$$
(35)
The brane tension forces the collapse of the brane while the momentum works oppositely. The same value $`r_0`$ can be determined from the balance condition of these two acts. Indeed, equating the first and last terms under the square root (31)

$$rT_3S_2^{}=\frac{2|nm|}{\pi r}$$
(36)
yields the eq.(35). Therefore the D3-brane with the given charges can be regarded as an expanded FD string which is self-supported from collapse by the momentum $`𝐃\times 𝐁`$.

Legendre discrete transformation maps the above D3-brane to the its electromagnetic dual D3-brane

$$𝐃_\theta =\frac{D}{\mathrm{sin}\theta },𝐁_z=B\mathrm{sin}\theta ;|BD|=r.$$
(37)
The quantization of the electromagnetic charges gives rise to the S-dual picture, i.e. D-strings are aligned along the cylinder axis and F-strings are wound around them helically. One should note that the precise form of the helix does not have any meaning because the worldvolume reparametrization invariance allows to transmute it to another form. What is definite, is that F and D strings are wound around each other and form a 3-dimensional braid. Necessarily they carry nonzero angular momentum, otherwise the system collapses to the bound state of the strings.

## 4 Dualities

The weak-strong duality may be realized as a NLE electric-magnetic duality on the D3-brane worldvolume theory. In this section we envisage the possibility that this conception is presented in the precise form. In the case of constant flux in the flat direction of the tubular D3-brane we explicitly demonstrate the equivalence of these dualities.

Equations of motion of the NLE action (5) and Hamiltonian (8) are invariant under $`SO(2)`$ electric-magnetic duality rotations. It has been shown that this invariance becomes $`SL(2,)`$ symmetry by including dilaton and axion fields . The Lagrangian takes the form

$$^{}=e^\varphi +\chi (𝐄𝐁).$$
(38)
It defines the electric induction $`𝐃^{}`$

$$𝐃^{}=e^\varphi 𝐃+\chi 𝐁$$
(39)
which shows that the axion charge of the string is provided by electric flux in the worldvolume of the brane. To restore the symmetric dependence of the Hamiltonian on the inductions $`𝐃^{}`$ and $`𝐇`$ one has also to rescale the metric tensor.

The duality transformations are conveniently described in terms of the complex field $`\lambda =\chi +ie^\varphi `$. They mix electric and magnetic fields

$`\left(\begin{array}{c}𝐃^{}\\ 𝐁\end{array}\right)S\left(\begin{array}{c}𝐃^{}\\ 𝐁\end{array}\right)`$ (44)
and induce a Möbius transformation for the $`\lambda `$ field,

$$\lambda \frac{a\lambda +b}{c\lambda +d},$$
(45)
where $`SSL(2,)`$ is

$`S=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)`$ (48)
with $`adbc=1`$.

In the static case $`𝐃^{}`$ satisfies the same Gauss law and quantization condition (14). Its normalization change is caused by the replacement $`1/g_sexp(\varphi )`$. For instance, instead of eq.(18) one gets

$$D^{}=\frac{\pi n}{r^2}.$$
(49)
Owing to this, $`SL(2,)`$ transformations can be formulated in terms of $`𝐃`$ and $`𝐁`$ introduced in the previous section

$`\left(\begin{array}{c}\frac{1}{g_s}𝐃\\ 𝐁\end{array}\right)S\left(\begin{array}{c}\frac{1}{g_s}𝐃\\ 𝐁\end{array}\right)`$ (54)
for the quantized charges. We now consider a D3-brane with minimal electric charge, $`n=1`$. This has to shrink to the F-string and we are actually considering a fundamental string. We will map it to the $`(n,m)`$ strings for any relatively prime pair $`n`$ and $`m`$. Indeed, it is well known that the Diophante equation

$$nxmy=1$$
(55)
always has a solution for the such a pair of integers. In fact, it has countably many solutions, because if $`(x=i,y=j)`$ is a solution then for any integer $`k`$ $`(x=km+i,y=kn+j)`$ is also a solution. Then for the given pair $`(n,m)`$ we have a matrix $`SSL(2,)`$

$`S=\left(\begin{array}{cc}n& j\\ m& i\end{array}\right).`$ (58)
We perform an electric-magnetic duality rotation (54) via $`S`$, i.e.

$`\left(\begin{array}{cc}n& j\\ m& i\end{array}\right){\displaystyle \frac{\pi }{r^2}}\left(\begin{array}{c}1\\ 0\end{array}\right)={\displaystyle \frac{\pi }{r^2}}\left(\begin{array}{c}n\\ m\end{array}\right).`$ (65)
It produced parallel electric and magnetic fields with corresponding charges $`n`$ and $`m`$. We already ascertained that the D3-brane with these charges collapses to the bound state of strings. Thus electric-magnetic duality mapped the fundamental string to the $`(n,m)`$ strings which is a weak-strong duality transformation .

We note that eq.(55) has no solutions when $`n`$ and $`m`$ are not relatively prime. Consequently we can not map the fundamental string to the FD strings in this case, they are not in the same $`SL(2,)`$ multiplet. This means that the bound state of strings becomes unstable when $`n`$ and $`m`$ have a common multiplier.

## 5 Conclusions

The NLE static equations (11,12) admit a wider class of solutions with nonzero Poynting vector, but we did not consider these. Our justification is that these solutions do not allow the application of $`SL(2,)`$ tranformations. In curved subspaces SO(2) duality is maintained but other dualities, in general, are broken . One can see this from the expression (30), the charges are quantized in different units and they can not be mapped to each other. There are two possibilities to get D3-brane charges quantized in the same units, that is to consider branes of the shape $`^2𝕊^1`$ or $`^1𝕊^1𝕊^1`$. In both cases there are two orthogonal equivalent directions along which the electric and magnetic fields can be directed. It will restore $`SL(2,)`$ electric-magnetic duality and describe an $`SL(2,)`$ multiplet of D3-branes which can be mapped to each other.

Acknowledgments:
S.T. acknowledges support by the Deutsche Forschungsgemeinschaft(DFG). D.K.P. was supported by grant R05-2001-000-00106-0 from the Basic Research Program of the Korea Science and Engineering Foundation.


text:
# Monopole Condensation in full QCD using the Schrödinger Functional

## 1 INTRODUCTION

To give a possible explanation of color confinement G. ’t Hooft  and S. Mandelstam  suggested long time ago that vacuum in gauge theories behaves like a magnetic (dual) superconductor. The dual superconductivity hypothesis relies upon the very general assumption that the dual superconductivity of the ground state is realized if there is condensation of Abelian magnetic monopoles.

### 1.1 The thermal partition functional

To investigate vacuum structure of lattice gauge theories at finite temperature we introduced  a lattice thermal partition functional in presence of a static background field

$$𝒵_T\left[\stackrel{}{A}^{\text{ext}}\right]=_{U_k(L_t,\stackrel{}{x})=U_k(0,\stackrel{}{x})=U_k^{\text{ext}}(\stackrel{}{x})}𝒟Ue^{S_W},$$
(1)
$`S_W`$ is the Wilson action, the physical temperature is given by $`T=1/aL_t`$. The functional integration is extended over links on a lattice $`L_s^3\times L_t`$ with hypertorus geometry and satisfying the constraints ($`x_t`$: temporal coordinate)

$$U_k(x)|_{x_t=0}=U_k^{\mathrm{ext}}(\stackrel{}{x}),(k=1,2,3),$$
(2)
$`U_k^{\mathrm{ext}}(x)`$ being the lattice version of the external continuum gauge field $`\stackrel{}{A}^{\mathrm{ext}}(x)=\stackrel{}{A}_a^{\mathrm{ext}}(x)\lambda _a/2`$.

$$F[\stackrel{}{A}^{\mathrm{ext}}]=\frac{1}{L_t}\mathrm{ln}\frac{𝒵_T[\stackrel{}{A}^{\mathrm{ext}}]}{𝒵_T[0]}.$$
(3)
To circumvent the problem of computing a partition function which is the exponential of an extensive quantity, in practical computations we consider $`F^{}[\stackrel{}{A}^{\mathrm{ext}}]`$ the $`\beta `$-derivative of $`F[\stackrel{}{A}^{\mathrm{ext}}]`$. Indeed this quantity is simply related to the plaquette $`U_{\mu \nu }`$ ($`\mathrm{\Omega }=L_s^3\times L_t`$):

$$\begin{array}{cc}\hfill F^{}[\stackrel{}{A}^{\mathrm{ext}}]& =\frac{1}{\mathrm{\Omega }}\underset{x,\mu <\nu }{}\frac{1}{3}\text{Re}\text{Tr}U_{\mu \nu }(x)_0\hfill \\ & \frac{1}{\mathrm{\Omega }}\underset{x,\mu <\nu }{}\frac{1}{3}\text{Re}\text{Tr}U_{\mu \nu }(x)_{\stackrel{}{A}^{\mathrm{ext}}},\hfill \end{array}$$
(4)
where the subscripts on the averages indicate the value of the external field.

### 1.2 Including fermions

We are interested in gauge systems at thermal equilibrium in presence of a static (time-independent) external background field. So that in presence of dynamical fermions Eq. (1) becomes:

$$\begin{array}{cc}& 𝒵_T\left[\stackrel{}{A}^{\text{ext}}\right]=\hfill \\ & _{U_k(L_t,\stackrel{}{x})=U_k(0,\stackrel{}{x})=U_k^{\text{ext}}(\stackrel{}{x})}𝒟U𝒟\psi 𝒟\overline{\psi }e^{(S_W+S_F)},\hfill \end{array}$$
(5)
where we integrate on fermionic fields without any constraint. As usual we impose on fermionic fields periodic boundary conditions in the spatial directions and antiperiodic boundary conditions in the temporal direction.

### 1.3 Detecting monopole condensation

Abelian monopole condensation can be detected using order/disorder parameters . The disorder parameter is related to the monopole free energy and is defined by means of the thermal partition function in presence of the Abelian monopole background field:

$$\mu =e^{F_{\text{mon}}/T_{\text{phys}}}=\frac{𝒵_T\left[𝐀^{\text{mon}}\right]}{𝒵_T[0]},$$
(6)
$`F_{\text{mon}}`$ is the free energy to create a monopole. If there is condensation $`F_{\text{mon}}`$ is finite and $`\mu 0`$. Note that since our disorder parameter $`\mu `$ has been defined in terms of the thermal partition functional $`𝒵_T`$ that is gauge-invariant for time-independent gauge transformations of the external background fields gauge fixing is not needed to perform the Abelian projection in the case of Abelian background fields.

## 2 ABELIAN MONOPOLES IN FULL QCD

### 2.1 Definition

In the continuum the magnetic monopole field with the Dirac string in the direction $`\stackrel{}{n}`$ is

$$e\stackrel{}{b}(\stackrel{}{x})=\frac{n_{\mathrm{mon}}}{2}\frac{\stackrel{}{x}\times \stackrel{}{n}}{|\stackrel{}{x}|(|\stackrel{}{x}|\stackrel{}{x}\stackrel{}{n})},$$
(7)
where, according to the Dirac quantization condition, $`n_{\mathrm{mon}}`$ is an integer and $`e`$ is the electric charge magnetic charge = $`n_{\mathrm{mon}}/2e`$.

For SU(3) gauge theory the maximal Abelian group is U(1)$`\times `$U(1), therefore we may introduce two independent types of Abelian monopoles associated respectively to the $`\lambda _3`$ and $`\lambda _8`$ diagonal generators of SU(3). We shall consider here the $`\lambda _3`$ Abelian monopole field. On the lattice it is given by ($`(X_1,X_2,X_3)`$ monopole coordinates)

$`U_{1,2}^{\text{ext}}(\stackrel{}{x})`$ $`=\left[\begin{array}{ccc}e^{i\theta _{1,2}^{\text{mon}}(\stackrel{}{x})}& 0& 0\\ 0& e^{i\theta _{1,2}^{\text{mon}}(\stackrel{}{x})}& 0\\ 0& 0& 1\end{array}\right],`$
$`U_3^{\text{ext}}(\stackrel{}{x})`$ $`=\mathrm{𝟏},`$ (8)
$`\theta _{1,2}^{\text{mon}}(\stackrel{}{x})`$ $`={\displaystyle \frac{n_{\text{mon}}}{4}}{\displaystyle \frac{(x_{2,1}X_{2,1})}{|\stackrel{}{x}_{\text{mon}}|}}{\displaystyle \frac{1}{|\stackrel{}{x}_{\text{mon}}|(x_3X_3)}}.`$
### 2.2 Numerical simulations

We used the standard HMC R-algorithm for two degenerate flavors of staggered fermions with a quark mass $`am=0.075`$ (at this value of the mass and $`N_T=4`$ $``$ $`\beta _c5.35`$ ). We have collected about 2000 thermalized trajectories for each value of $`\beta `$ at $`L_s=16`$ and about 500 thermalized trajectories for each value of $`\beta `$ at $`L_s=32`$. Each trajectory consists of $`125`$ molecular dynamics steps and has total length $`1`$. The computer simulations have been performed on the APEmille crate at INFN/Bari.

### 2.3 Numerical results

In Fig. 1 we compare $`F_{\text{mon}}^{}`$ in the 2 flavors full QCD case for two different spatial volumes. Data in the weak and strong coupling regions agree within statistical errors. On the other hand the signal at the peak value gets increased going to a larger spatial volume.

Observing that $`F_{\text{mon}}=0`$ at $`\beta =0`$, we may obtain $`F_{\text{mon}}`$ from $`F_{\text{mon}}^{}`$ by a numerical integration in $`\beta `$ and from Eq. (6) we are able to estimate the disorder parameter $`\mu `$. Fig. 3 shows that $`\mathrm{ln}\mu =0`$ in the confined phase, i.e. the free energy required to create an Abelian monopole in the QCD vacuum is zero and therefore Abelian monopoles condense in the confined phase.

Fig. 4 displays $`F_{\text{mon}}^{}`$ in the 2 flavors full QCD case and $`L_s=32`$ compared with the chiral condensate. The peak in $`F_{\text{mon}}^{}`$ corresponds to the drop of the chiral condensate.

## 3 CONCLUSIONS

Using a thermal partition functional in presence of an external Abelian monopole background field we simulate 2 staggered flavors full QCD on a $`L_s^3\times 4`$ lattice. We find that Abelian monopoles condense in the spontaneously broken chiral symmetry phase of 2 flavors full QCD. This is consistent with similar results found in Ref.  For a better understanding of Abelian monopole condensation in full QCD we expect to perform new simulations at different values of the masses and of the flavor numbers as well a finite size scaling analysis.


text:
# 1 Introduction and History

## 1 Introduction and History

The TEXONO Collaboration has been built up since 1997 to initiate and pursue an experimental program in Neutrino and Astroparticle Physics . The Collaboration comprises more than 40 research scientists from major institutes/universities in Taiwan (Academia Sinica, Chung-Kuo Institute of Technology, Institute of Nuclear Energy Research, National Taiwan University, National Tsing Hua University, and Kuo-Sheng Nuclear Power Station), China (Institute of High Energy Physics, Institute of Atomic Energy, Institute of Radiation Protection, Nanjing University, Tsing Hua University) and the United States (University of Maryland), with AS, IHEP and IAE (with ) being the leading groups. It is the first research collaboration of this size and magnitude among major research institutes from Taiwan and China.

The research program  emphasizes on the the unexplored and unexploited theme of adopting detectors with high-Z nuclei, such as solid state device and scintillating crystals, for low-energy low-background experiments in Neutrino and Astroparticle Physics. The “Flagship” program  is a reactor neutrino experiment to study low energy neutrino properties and interactions at the Kuo-Sheng (KS) Neutrino Laboratory. It is the first particle physics experiment performed in Taiwan. In parallel to the reactor experiment, various R&D efforts coherent with the theme are initiated and pursued.

This article focuses on the magnetic moment data taking and analysis with a germanium detector during the Reactor Period June 2001 till April 2002.

## 2 Electromagnetic Properties of the Neutrinos

The strong and positive evidence of neutrino oscillations implies the existence of neutrino masses and mixings , the physical origin, structures and experimental consequences of which are still not thoroughly known and understood. Experimental studies on the neutrino properties and interactions which may reveal some of these fundamental questions and/or constrain certain classes of models are therefore of interests. The coupling of neutrino with the photons are consequences of non-zero neutrino masses. Two of the manifestations of the finite electromagnetic form factors for neutrino interactions  are neutrino magnetic moments and radiative decays.

The searches of neutrino magnetic moments are performed in experiments on neutrino-electron scatterings :

$$\nu _{l_\mathit{1}}+\mathrm{e}^{}\nu _{l_\mathit{2}}+\mathrm{e}^{}.$$
(1)
The experimental observable is the kinetic energy of the recoil electrons(T). The differential cross section for the magnetic scattering (MS) channel can be parametrized by the neutrino magnetic moment ($`\mu _l`$), often expressed in units of the Bohr magneton($`\mu _\mathrm{B}`$). Its dependence on neutrino energy $`\mathrm{E}_\nu `$ is given by :

$$(\frac{\mathrm{d}\sigma }{\mathrm{dT}})_{\mathrm{MS}}=\frac{\pi \alpha _{\mathrm{em}}^2\mu _{l}^{}{}_{}{}^{2}}{\mathrm{m}_\mathrm{e}^2}[\frac{1\mathrm{T}/\mathrm{E}_\nu }{\mathrm{T}}].$$
(2)
The process can be due to diagonal and transition magnetic moments, for the cases where $`l_1=l_2`$ and $`l_1l_2`$, respectively. The MS interactions involve a flip of spin and therefore do not have interference with the Standard Model (SM) cross-section. The $`\mu _l`$ term has a 1/T dependence and hence dominates at low electron recoil energy over the SM process.

The quantity $`\mu _l`$ is an effective parameter which, in the case of $`\overline{\nu _e}`$ and large mixings between the mass eigenstates, can be expressed as :

$$\mu _\mathrm{e}^2=\underset{\mathrm{k}}{}\left|\underset{\mathrm{j}}{}\mathrm{U}_{\mathrm{ej}}\mu _{\mathrm{jk}}\right|^2$$
(3)
where U is the mixing matrix and $`\mu _{\mathrm{jk}}`$ are the fundamental constants that characterize the couplings between the mass eigenstates $`\nu _\mathrm{j}`$ and $`\nu _\mathrm{k}`$ with the electromagnetic field. The $`\nu `$-$`\gamma `$ couplings probed by $`\nu `$-e scatterings is the same as that giving rise to the neutrino radiative decays :

$$\nu _\mathrm{j}\nu _\mathrm{k}+\gamma $$
(4)
between $`\nu _\mathrm{j}`$ and $`\nu _\mathrm{k}`$. The decay rate $`\mathrm{\Gamma }_{\mathrm{jk}}`$ is related to $`\mu _{\mathrm{jk}}`$ by:

$$\mathrm{\Gamma }_{\mathrm{jk}}=\frac{\mu _{\mathrm{jk}}^2}{4\pi }\frac{(\mathrm{\Delta }\mathrm{m}_{\mathrm{jk}}^2)^3}{\mathrm{m}_\mathrm{j}^3}$$
(5)
where $`\mathrm{m}_\mathrm{j}`$ is the mass of $`\nu _\mathrm{j}`$ and $`\mathrm{\Delta }\mathrm{m}_{\mathrm{jk}}^2=\mathrm{m}_\mathrm{j}^2\mathrm{m}_\mathrm{k}^2`$. The total decay rate is given by

$$\mathrm{\Gamma }_{\nu _\mathrm{e}}=\underset{\mathrm{k}}{}\left|\underset{\mathrm{j}}{}\mathrm{U}_{\mathrm{ej}}\mathrm{\Gamma }_{\mathrm{jk}}\right|^2$$
(6)
Reactor neutrinos provide a sensitive probe for “laboratory” searches of $`\mu _\mathrm{e}`$, taking advantages of the high $`\overline{\nu _e}`$ flux, low $`\mathrm{E}_\nu `$ and the better experimental control via the reactor ON/OFF comparison. Neutrino-electron scatterings were first observed in pioneering experiment  at Savannah River. However, a reanalysis of the data by Ref  with improved input parameters on the reactor neutrino spectra and $`\mathrm{sin}^2\theta _\mathrm{W}`$ gave a positive signature consistent with the interpretation of a finite magnetic moment at $`(24)\times 10^{10}\mu _\mathrm{B}`$. Other results came from the Kurtchatov  and Rovno  experiments which quoted limits of $`\mu _{\overline{\nu _\mathrm{e}}}`$ less than $`2.4\times 10^{10}\mu _\mathrm{B}`$ and $`1.9\times 10^{10}\mu _\mathrm{B}`$ at 90% confidence level (CL), respectively. However, the lack of experimental details and discussions in the published work make it difficult to assess the robustness of the results.

Theoretically, a minimally-extended Standard Model would give rise to $`\mu _\nu `$’s for Dirac neutrinos too small to be of interest . However, there are models  which can produce large $`\mu _\nu `$ by incorporating new features like right-handed currents and transition moments. Neutrino flavor conversion induced by resonant or non-resonant spin-flip transitions in the Sun via its transition magnetic moments has been considered as solution to the Solar Neutrino Problem . Stimulated by the new results from the Super-Kamiokande and SNO experiments, recent detailed work  suggested that this scenario is also in excellent agreement with the existing solar neutrino data. Alternatively, the measured solar neutrino $`\nu _{}`$-e spectra has been used to set limits of $`\mu _\nu ^{}<1.5\times 10^{10}\mu _\mathrm{B}`$ at 90% CL for the “effective” $`\nu _{}`$ magnetic moment , which is in general different from that of a pure $`\nu _l`$ state derived in laboratory experiments. In addition, there are astrophysical arguments  from nucleosynthesis, stellar cooling and SN1987a which placed limits of the range $`\mu _\mathrm{e}<10^{12}10^{13}\mu _\mathrm{B}`$. However, care should be taken in their interpretations due to the model dependence and assumptions on the neutrino properties implicit in the derivations.

These discussions show that further laboratory experiments to put the current limits on more solid grounds and to improve on the sensitivities are of interest. Beside the Kuo-Sheng (KS) experiment  reported in this article, there is another on-going experiment MUNU  at the Bugey reactor using a time projection chamber with CF4 gas.

## 3 Kuo-Sheng Reactor Neutrino Laboratory

The “Kuo-Sheng Reactor Neutrino Laboratory” is located at a distance of 28 m from the core #1 of the Kuo-Sheng Nuclear Power Station at the northern shore of Taiwan . A schematic view is depicted in Figure 1.

A multi-purpose “inner target” detector space of 100 cm$`\times `$80 cm$`\times `$75 cm is enclosed by 4$`\pi `$ passive shielding materials and cosmic-ray veto scintillator panels, the schematic layout of which is shown in Figure 2. The shieldings provide attenuation to the ambient neutron and gamma background, and are made up of, inside out, 5 cm of OFHC copper, 25 cm of boron-loaded polyethylene, 5 cm of steel and 15 cm of lead.

Different detectors can be placed in the inner space for the different scientific goals. The detectors are read out by a versatile electronics and data acquisition systems  running on a VME-PCI bus. Signals are recorded by Flash Analog-to-Digital-Convertor (FADC) modules with 20 MHz clock and 8-bit resolution. The readout allows full recording of all the relevant pulse shape and timing information for as long as several ms after the initial trigger. The reactor laboratory is connected via telephone line to the home-base laboratory at AS, where remote access and monitoring are performed regularly. Data are stored and accessed in a multi-disks array with a total of 600 Gbyte memory via IDE-bus in PCs.

It is recognized recently  that the low energy part of the reactor neutrino spectra is not well modeled or experimentally checked. Consequently, the uncertainties induced in the SM $`\overline{\nu _e}`$-e cross-sections can limit the sensitivities of magnetic moment searches at the domain where they are comparable or larger than the MS interactions. Therefore, experiments intended for measure Standard Model cross sections with reactor neutrinos should focus on higher energies (T$`>`$1.5 MeV) while $`\mu _\nu `$ searches should base on measurements with T$`<`$100 keV.

Accordingly, data were taken for Period I Reactor ON/OFF from June 2001 till April 2002 with these strategies. Two detector systems are running in parallel using the same data acquisition system but independent triggers: (a) an Ultra Low Background High Purity Germanium (ULB-HPGe), with a fiducial mass of 1.06 kg, and (b) 46 kg of CsI(Tl) crystal scintillators.

Preliminary results of the HPGe system is presented in the following sections. The performance of the CsI(Tl) modules is discussed elsewhere .

## 4 Low Background Germanium Detector

The low threshold and excellent energy resolution of the germanium detector make it optimal for $`\mu _\nu `$ searches . The set-up of the KS-Ge experiment is schematically shown in Figure 3. It is a coaxial germanium detector with an active target mass of 1.06 kg. The lithium-diffused outer electrode is 0.7 mm thick. The end-cap cryostat, also 0.7 mm thick, is made of OFHC copper. Both of these features provide total suppression to ambient $`\gamma `$-background below 60 keV, such that events below this energy are either due to internal activity or external high energy $`\gamma `$’s via Compton scattering. The HPGe was surrounded by an anti-Compton (AC) detector system made up of two components: (1) an NaI(Tl) well-detector of thickness 5 cm that fit onto the end-cap cryostat, the inner wall of which is also made of OFHC copper, and (2) a 4 cm thick CsI(Tl) detector at the bottom. Both AC detectors were read out by photo-multipliers (PMTs) with low-activity glass. The assembly was surrounded by 3.7 cm of OFHC copper inner shielding. Another 10 cm of lead provided additions shieldings on the side of the liquid nitrogen dewar and pre-amplifier electronics. The inner shieldings and detectors were covered by a plastic bag connected to the exhaust line of the dewar, serving as a purge for the radioactive radon gas.

The HPGe pre-amplifier and AC PMT signals were distributed via 10 m cables to two spectroscopy amplifiers at the same 4 $`\mu `$s shaping time but with different gain factors. A very loose amplitude threshold on the SA output provided the on-line trigger, ensuring all the events down to the electronics noise edge were recorded. The SA output, as well as the PMT signals from the AC detectors, were recorded by the FADC for a duration of 25 $`\mu `$s after the trigger. The discriminator output of the Veto PMTs were also recorded. A random trigger was provided by an external clock once per 10 s for sampling the background level and evaluating the various efficiency factors. The DAQ system remains active for 2 ms after a trigger to record possible time-correlated signatures. The activities in any part of the detector systems within 10 $`\mu `$s prior to the trigger were also recorded. The typical data taking rate for the HPGe sub-system was about 1 Hz. The accurately recorded DAQ dead time was about 10-20 ms per event.

## 5 Results

The measured spectra, after cuts of cosmic and anti-Compton vetos, during 60/46 days of reactor ON/OFF data taking are displayed in Figure 4. The background level of 1 keV-1kg-1day-1 and a detector threshold of 5 keV, comparable to those of underground Dark Matter experiment, are achieved. Additional cuts based on pulse shape and timing information are expected to further reduce the background level and the threshold.

Several lines can be identified: Ga X-rays at 10.37 keV and 73Ge at 66.7 keV from internal cosmic-induced activities, and 234Th at 63.3 keV and 92.6 keV due to residual ambient radioactivity from the $`{}_{}{}^{238}\mathrm{U}`$ series in the vicinity of the target. The ON and OFF spectra differ in two significant features. The excess of OFF over ON below the noise edge of 5 keV is due to instabilities in the trigger threshold and do not affect the analysis at higher energies. The excess at the Ga X-ray peaks originates from the long-lived isotopes (68Ge and 71Ge with half-lives of 271 and 11.4 days, respectively) activated by cosmic-rays prior to installation. It has been checked that the time evolution of the peak intensity can be fit to two exponentials consistent with the two known half-lives.

The reactor neutrino spectra was evaluated from reactor operation data using the standard prescriptions on from fission $`\overline{\nu _e}`$  together with a low energy contribution due to neutron capture on $`{}_{}{}^{238}\mathrm{U}`$ . The total flux is $`5.6\times 10^{12}\mathrm{cm}^2\mathrm{s}^1`$. The electron recoil spectra from SM and MS interactions can then be calculated. Background in the energy range of 12 to 60 keV are due to Compton scatterings of higher energy $`\gamma `$’s. Accordingly, the Reactor OFF data are fitted to a smooth function

$$\varphi _{\mathrm{OFF}}=\alpha _1\mathrm{e}^{\alpha _2\mathrm{E}}+\alpha _3+\alpha _4\mathrm{E}.$$
(7)
A $`\chi ^2`$/dof of 40/46 is obtained indicating this background description is valid. A one-parameter fit is then performed to the Reactor ON data for

$$\varphi _{\mathrm{OFF}}+\varphi _{\mathrm{SM}}+\mathrm{k}^2\varphi _{\mathrm{MS}},$$
(8)
where $`\varphi _{\mathrm{SM}}`$ and $`\varphi _{\mathrm{MS}}`$ are the expected recoil spectra due to weak interactions and MS at $`\mu _{\overline{\nu _e}}=10^{10}\mu _\mathrm{B}`$, respectively. A best-fit value of

$$\mathrm{k}^2=1.1\pm 2.5$$
(9)
at a $`\chi ^2`$/dof of 55/49 is obtained. The residual plot and the best-fit 1$`\sigma `$ region is depicted in Figure 5.

A total of 180/60 days of Reactor ON/OFF are taken for this reactor period. A limit will be derived when the data analysis is improved, complete and finalized. It is expected that the data would give world-level sensitivities in the searches of $`\overline{\nu _e}`$ magnetic moments. Depicted in Figure 6 is the summary of all the results in $`\mu _{\overline{\nu _e}}`$ searches with reactor $`\overline{\nu _e}`$ versus the achieved threshold. The dotted lines are the $`\sigma _\mu /\sigma _{\mathrm{SM}}`$ ratio at a particular (T, $`\mu _{\overline{\nu _e}}`$). The KS experiment operated at a much lower threshold of 12 keV compared to previous and current measurements. The large $`\sigma _\mu /\sigma _{\mathrm{SM}}`$ ratio at low energy implies that effects due to the uncertainties in the SM cross-sections can be neglected such that the limits derived are more robust.

Indirect bounds on the neutrino radiative decay rate are inferred using Eq. 5 and displayed in Figure 7. The sensitivity region of $`\mu _\nu =10^{10}\mu _\mathrm{B}`$ is shown for illustration purpose. Superimposed are the limits from the previous direct searches of excess $`\gamma `$’s from reactor neutrinos  and from the supernova SN1987a . Also shown is the sensitivity level of proposed simulated conversion experiments at accelerator . It can be seen that $`\nu `$-e scatterings give much more stringent bounds than the direct searches.

## 6 Summary & Outlook

With the strong evidence of physics beyond the Standard Model revealed in the neutrino sector, neutrino physics and astrophysics remains a central subject in experimental particle physics in the coming decade and beyond. There are room for ground-breaking technical innovations - as well as potentials for surprises in the scientific results.

A Taiwan, China and U.S.A. collaboration has been built up with the goal of establishing a experimental program in neutrino and astro-particle physics. It is the first generation collaborative efforts in large-scale basic research between scientists from Taiwan and China. The flagship effort is to perform the first-ever particle physics experiment in Taiwan at the Kuo-Sheng Reactor Plant. From the Period I data taking, we expect to achieve world-level sensitivities in neutrino magnetic moments and radiative lifetime studies. A wide spectrum of R&D projects are being pursued. New ideas for future directions are being explored.

## 7 Acknowledgments

The author is grateful to the scientific members, technical staff and industrial partners of TEXONO Collaboration, as well as the concerned colleagues for the invaluable contributions which “make it happen” in such a short period of time. Funding are provided by the National Science Council, Taiwan and the National Science Foundation, China, as well as from the operational funds of the collaborating institutes.