format: text

text:
# Ubiquitous finite-size scaling features in 𝐼-𝑉 characteristics of various dynamic 𝑋𝑌 models in two dimensions

## I Introduction

The phase transition between the superconducting and normal states in many two-dimensional (2D) systems is of Kosterlitz-Thouless (KT) type . The thermally excited vortices in a large enough sample, interacting via a logarithmic potential, are bound in neutral pairs below the KT transition temperature $`T_{KT}`$ , and as the temperature $`T`$ is increased across $`T_{KT}`$ from below these pairs start to unbind. The KT transition has been observed in experiments on superconducting films , 2D Josephson junction arrays , and cuprate superconductors . In these experiments, the current-voltage ($`I`$-$`V`$) characteristics have been commonly measured to detect the transition. For small enough currents it has a power law form $`VI^a`$ (or equivalently, $`EJ^a`$ with the electric field $`E`$ and the current density $`J`$), where the $`I`$-$`V`$ exponent $`a`$ is known to have a universal value 3 precisely at the transition ; For $`T<T_{KT}`$ one has $`a>3`$, whereas $`a=1`$ for $`T>T_{KT}`$ . The dynamic critical exponent $`z`$ which relates the relaxation time $`\tau `$ to the vortex correlation length $`\xi `$ via $`\tau \xi ^z`$, is connected with $`a`$ through the relation $`a=1+z`$ . Consequently, $`z`$ has the value 2 at $`T_{KT}`$.

Indeed, the values $`z=2`$ and $`a=3`$ at the KT transition have been confirmed in many numerical simulations, e.g., the lattice Coulomb gas with Monte Carlo dynamics , the Langevin-type molecular dynamics of Coulomb gas particles , and the 2D XY model both with the resistively shunted junction dynamics (RSJD) and the relaxational dynamics (RD) . That is why $`z6`$ obtained by Pierson et al. in Ref.  at the resistive transition in many 2D systems is very intriguing. We in Ref.  have re-analyzed the experimental $`I`$-$`V`$ characteristics for an ultra thin YBCO sample in Ref.  and compared with the 2D RSJD model. This led to the suggestion that a novel finite-size type scaling effect of the $`I`$-$`V`$ characteristics (which can possibly be caused not only by the actual finite size of the sample but also by a finite perpendicular penetration depth or a residual weak magnetic field ) is responsible for the large value of $`z6`$ in Ref. , but at the same time it was concluded that this exponent has nothing to do with the true dynamic critical exponent which at $`T_{KT}`$ has value two . Since similar conclusions were also reached in Ref. , we believe that there now seems to emerge some consensus , although the physical mechanism causing this finite-size scaling behavior is still unclear.

In the present paper we search for a possible reason of this finite-size scaling behavior. We use the RSJD , the RD (often referred to as time-dependent Ginzburg-Landau dynamics) , and the Monte Carlo dynamics (MCD) simulations  to study the $`I`$-$`V`$ characteristics of three different models of 2D superconductors: the Villain model , the $`XY`$ model in its original form, and the $`XY`$ model modified with $`p`$-type of potential (see Ref.  for details). These three models differ from each other by the density of the thermally created vortices. We come to the conclusion that the qualitative features of the scaling suggested in Ref.  are ubiquitous and independent of both the vortex density and the type of dynamics. We also demonstrate that the $`I`$-$`V`$ curves obtained with different types of dynamics, up to constant scale factors, coincide very well in a broad range of the external current density. This opens a possibility to use MCD simulation, which from a simulation point of view is more efficient than RSJD or RD, to examine long-time dynamic properties of the models.

The layout of the paper is as follows: In Sec. II we recapitulate the Hamiltonians for the generalized 2D XY model with the $`p`$-type of potential (the usual XY model is recovered when $`p=1`$) and the Villain model, and describe details of the dynamics used (RSJD, RD, and MCD). The results from the usual XY model with $`p=1`$ subject to different dynamics (RSJD, RD, and MCD) are presented in Sec. III, while the results from RSJD simulations applied to the different types of models, i.e., the usual XY model with $`p=1`$, the XY model with $`p=2`$ and the Villain model, are described and analyzed in Sec. IV. We summarize and make final remarks in Sec. V.

## II Models and details of dynamics

### A Model Hamiltonians

The 2D $`XY`$ model defined on a square lattice, where each lattice point $`i`$ is associated with the phase $`\theta _i`$ of the superconducting order parameter, is often used for studies of the KT transition. The phase variables in this model interact via the Hamiltonian with nearest neighbor coupling, which in the absence of frustration is given by

$`H={\displaystyle \underset{ij}{}}U(\varphi _{ij}=\theta _i\theta _j),`$ (1)
where $`ij`$ denotes sum over nearest neighbor pairs, $`\varphi `$ is the angular difference between nearest neighbors, and the interaction potential $`U(\varphi )`$ is written as

$`U(\varphi )E_J(1\mathrm{cos}\varphi ),`$ (2)
with the Josephson coupling strength $`E_J`$.

The dominant characteristic physical features close to the KT transition are associated with vortex pair fluctuations. One interesting aspect is then how the density of the thermally excited vortex fluctuation effects the critical properties. To study this we generalize the interaction potential by using a parameter $`p`$ :

$`U^p(\varphi )2E_J\left[1\mathrm{cos}^{2p^2}\left({\displaystyle \frac{\varphi }{2}}\right)\right],`$ (3)
where $`U^{p=1}(\varphi )`$ corresponds to the potential of the usual $`XY`$ model \[see Eq. (2)\]. The practical point with such generalization is that the vortex density increases with increasing $`p`$ . The variation of the parameter $`p`$ can also change the nature of the transition: for $`p`$ exceeding some maximum value ($`p>p_{\mathrm{max}}5`$) the type of the phase transition changes from KT to the first order . In the present paper we choose $`p=2`$ which is well inside the KT transition region, yet is large enough to ensure substantially more vortex fluctuation over a temperature region around the phase transition in comparison with the usual $`p=1`$ $`XY`$ model given by Eq. (2). While the $`XY`$ model with the $`p`$-type potential with $`p>1`$ has more vortices than the usual $`p=1`$ $`XY`$ model, we also study the Villain model  which has less vortex-antivortex pairs . The interaction potential $`U(\varphi )`$ in the Villain model is given by

$$e^{U(\varphi )/T}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{exp}\left[\frac{E_J}{2T}(2\pi n\varphi )^2\right].$$
(4)
To simulate the dynamic behaviors of these models we use several types of dynamics: RSJD, RD, and MCD. All these dynamics should result in the same equilibrium static behaviors if we apply them to the models with the same interaction potential. However, dynamic properties of the systems can be different. Of course, different types of dynamics have their own advantages and disadvantages. The RSJD is constructed from the elementary Josephson relations for single Josephson junction that forms the array units, plus Kirchhoff’s current conservation condition at each lattice site . Therefore, this type of dynamics has a firm physical realization. On the other hand, RSJD is quite slow which leads to the limitation in the time scale one can probe in simulations. Although the RD  is much easier to implement than RSJD, it does not converge much faster than RSJD and it does not have a similar direct physical realization as RSJD. However, a superconductor has been argued to have a RD type of dynamics rather than a RSJD . The MCD simulations  are much faster than RSJD or RD, which allows one to investigate dynamic behaviors in much longer time scale (one can also study dynamic behaviors at much lower temperatures with MCD). However, since there is no direct physical realization of the MCD in practice the applicability of this dynamics to a specific physical system must then be explicitly demonstrated. In the following discussions on the details of the different dynamics used, we focus on the original 2D $`XY`$ model with the interaction potential in Eq. (2) since the extensions to a modified 2D $`XY`$ model (3) and Villain model (4) are straightforward.

### B Dynamic models

In this section we briefly review the dynamical equations of motion for RSJD, RD, and MCD, in the presence of the fluctuating twist boundary condition (FTBC) . We perform simulations of unfrustrated square $`L\times L`$ lattices with $`L=6`$, 8, and 10 at various temperatures to measure the voltage across the lattice as a function of the external current. Although the system sizes are relatively small, which is inevitable because of the low temperatures and the small external currents used here, the FTBC has been shown to be very efficient in reducing the artifact due to small system sizes , and reliable results can be established .

In the FTBC, the twist variable $`𝚫=(\mathrm{\Delta }_x,\mathrm{\Delta }_y)`$ is introduced and the phase difference $`\varphi _{ij}`$ on the bond $`(i,j)`$ is changed into $`\theta _i\theta _j𝐫_{ij}𝚫`$, with the unit vector $`𝐫_{ij}`$ from site $`i`$ to site $`j`$, while the periodicity on $`\theta _i`$ is imposed: $`\theta _i=\theta _{i+L\widehat{𝐱}}=\theta _{i+L\widehat{𝐲}}`$. The Hamiltonian of 2D $`L\times L`$ $`XY`$ model under FTBC without external current has been introduced in Ref. , and is written as \[compare with Eq. (1)\]

$$H=E_J\underset{ij}{}\mathrm{cos}(\varphi _{ij}\theta _i\theta _j𝐫_{ij}𝚫),$$
(5)
which later in Ref.  has been extended to the system in the presence of an external current and written as

$$H=E_J\underset{ij}{}\mathrm{cos}\varphi _{ij}+\frac{\mathrm{}}{2e}L^2J\mathrm{\Delta }_x,$$
(6)
where $`J`$ the current density in the $`x`$ direction.

#### 1 RSJD and RD

We introduce first the RSJD equations of motion for phase variables and twist variables, which are generated from the local (global) current conservation for the phase (twist) variables (see Ref.  for details and discussions). The net current $`I_{ij}`$ from site $`i`$ to site $`j`$ is the sum of the supercurrent $`I_{ij}^s`$, the normal resistive current $`I_{ij}^n`$, and the thermal noise current $`I_{ij}^t`$: $`I_{ij}=I_{ij}^s+I_{ij}^n+I_{ij}^t`$. The supercurrent is given by the Josephson current-phase relation, $`I_{ij}^s=I_c\mathrm{sin}\varphi _{ij}`$, where $`I_c=2eE_J/\mathrm{}`$ is the critical current of the single junction. The normal resistive current is given by $`I_{ij}^n=V_{ij}/r`$, where $`V_{ij}`$ is the potential difference across the junction, and $`r`$ is the shunt resistance. Finally the thermal noise current $`I_{ij}^t`$ in the shunt at temperature $`T`$ satisfies $`I_{ij}^t=0`$ and $`I_{ij}^t(t)I_{kl}^t(0)=(2k_BT/r)\delta (t)(\delta _{ik}\delta _{jl}\delta _{il}\delta _{jk})`$, where $`\mathrm{}`$ is thermal average, and $`\delta (t)`$ and $`\delta _{ij}`$ are Dirac and Kronecker delta, respectively. Using the current conservation law at each site of the lattice together with the Josephson relation $`\dot{\varphi }_{ij}d\varphi _{ij}/dt=2eV_{ij}/\mathrm{}`$ one can derive the RSJD equations of motion for phase variables:

$$\dot{\theta }_i=\underset{j}{}G_{ij}\underset{k}{}^{}(\mathrm{sin}\varphi _{jk}+\eta _{jk}),$$
(7)
where the primed summation is over the four nearest neighbors of $`j`$, $`G_{ij}`$ is the lattice Green function for 2D square lattice, and $`\eta _{jk}`$ is the dimensionless thermal noise current defined by $`\eta _{jk}I_{jk}^t/I_c`$. The time, the current, the distance, the energy, and the temperature are normalized in units of $`\mathrm{}/2erI_c`$, $`I_c`$, the lattice spacing $`a`$, the Josephson coupling strength $`E_J`$, and $`E_J/k_B`$, respectively. In order to get a closed set of equations we further specify the dynamics of the twist variable $`𝚫`$ from the condition of the global current conservation that the summation of the all currents through the system in each direction should vanish :

$`\dot{\mathrm{\Delta }}_x`$ $`=`$ $`{\displaystyle \frac{1}{L^2}}{\displaystyle \underset{ij_x}{}}\mathrm{sin}\varphi _{ij}+\eta _{\mathrm{\Delta }_x}J,`$ (8)
$`\dot{\mathrm{\Delta }}_y`$ $`=`$ $`{\displaystyle \frac{1}{L^2}}{\displaystyle \underset{ij_y}{}}\mathrm{sin}\varphi _{ij}+\eta _{\mathrm{\Delta }_y},`$ (9)
where $`_{ij_x}`$ denotes the summation over all nearest neighbor links in the $`x`$ direction, and we apply the external dc current with the current density $`J`$ in the $`x`$ direction. Here, the thermal noise terms $`\eta _{\mathrm{\Delta }_x}`$ and $`\eta _{\mathrm{\Delta }_y}`$ obey the conditions $`\eta _{\mathrm{\Delta }_x}=\eta _{\mathrm{\Delta }_y}=\eta _{\mathrm{\Delta }_x}\eta _{\mathrm{\Delta }_y}=0`$, and $`\eta _{\mathrm{\Delta }_x}(t)\eta _{\mathrm{\Delta }_x}(0)=\eta _{\mathrm{\Delta }_y}(t)\eta _{\mathrm{\Delta }_y}(0)=(2T/L^2)\delta (t)`$.

In the RD, the equations of motion for the phase variables are written as

$$\dot{\theta }_i=\mathrm{\Gamma }\frac{H}{\theta _i}+\eta _i=\underset{j}{}^{}\mathrm{sin}\varphi _{ij}+\eta _i,$$
(10)
where $`\mathrm{\Gamma }`$ is a dimensionless constant (we set $`\mathrm{\Gamma }=1`$ from now one), $`H`$ is in Eq. (6), $`t`$ is in units of $`\mathrm{}/\mathrm{\Gamma }E_J`$, and the thermal noise $`\eta _i(t)`$ at site $`i`$ satisfies $`\eta _i(t)=0`$ and $`\eta _i(t)\eta _j(0)=2T\delta (t)\delta _{ij}`$. The equation of motion for the twist variables in the absence of an external current is of the form (see Ref.  for more details)

$`\dot{𝚫}={\displaystyle \frac{1}{L^2}}{\displaystyle \frac{H}{𝚫}}+\eta _𝚫,`$ (11)
which is the same as Eqs. (8) and (9) for RSJD. Accordingly, to some extent the RD may be viewed as a simplified version of the RSJD where the global current conservation is kept but the local current conservation is relaxed.

Consequently, the equations for the phase variables are different for RSJD and RD \[Eqs. (7) and (10), respectively\] while the same equations (8) and (9) apply to the twist variables for both dynamics. These coupled equations of motion are discretized in time with the time step $`\mathrm{\Delta }t=0.05`$ and $`\mathrm{\Delta }t=0.01`$ for RSJD and RD, respectively, and numerically integrated using the second order Runge-Kutta-Helfand-Greenside algorithm . The voltage drop $`V`$ across the system in the $`x`$ direction is written as $`V=L\dot{\mathrm{\Delta }}_x`$ (see Ref. ) in units of $`I_cr`$ for RSJD and in units of $`\mathrm{\Gamma }E_J/2e`$ for RD, respectively. We measure the electric field $`E=V/L_t`$ to obtain $`I`$-$`V`$ characteristics, where $`\mathrm{}_t`$ denotes the time average performed over $`O(10^6)`$ time steps for large currents for both RSJD and RD, and $`O(5\times 10^7)`$ and $`O(10^9)`$ steps for small currents for RSJD and RD, respectively.

#### 2 MCD

The technique to simulate 2D XY model with MCD is based on the Hamiltonian (6) and the standard Metropolis algorithm . The one MC step, which we identify as a time unit, is composed as follows :

1. Pick one lattice site and try to rotate the phase angle at the site by an amount randomly chosen in $`[\delta \theta ,\delta \theta ]`$ (we call $`\delta \theta `$ the trial angle range). The twist variable $`𝚫`$ is kept constant during the update of the phase variables.
2. Compute the energy difference $`\mathrm{\Delta }H`$ before and after the above try; If $`\mathrm{\Delta }H<0`$ or if $`e^{\mathrm{\Delta }H/T}`$ is greater than a random number chosen on the interval $`[0,1)`$, accept the trial move.
3. Repeat steps 1 and 2 for all the lattice sites to update the phase variables.
4. Update the fluctuating twist variables $`\mathrm{\Delta }_x`$ in the similar way that $`L\mathrm{\Delta }_x`$ is tried to rotate within the angle range $`\delta \mathrm{\Delta }`$ with $`\theta _i`$ and $`\mathrm{\Delta }_y`$ kept unchanged. (For convenience, we use $`\delta \mathrm{\Delta }=\delta \theta `$).
5. Compute the energy difference $`\mathrm{\Delta }H`$ before and after the trial step 4 for $`\mathrm{\Delta }_x`$. Accept the step 4, if $`\mathrm{\Delta }H<0`$ otherwise accept it with probability $`e^{\mathrm{\Delta }H/T}`$ like in step 2.
6. Repeat steps 4 and 5 to update $`\mathrm{\Delta }_y`$.

In the MCD simulation the trial angle $`\delta \theta =\pi /6`$ has been chosen since it is sufficiently small in order to obtain the correct $`I`$-$`V`$ characteristics while it is big enough to make MCD much faster than the other dynamic methods . The time-averaged electric field is obtained after equilibration from the averages over $`O(10^9)`$ (at large currents) to $`O(10^{10})`$ (at small currents) MC steps.

## III 2D XY model subject to different types of dynamics

In this section we use three different types of dynamics, the RSJD, the RD, and the MCD to study dynamic behavior of the 2D XY model with $`p=1`$ under the FTBC. The dynamic behavior of the system can be obtained from the complex conductivity, the flux noise spectrum, as well as the $`I`$-$`V`$ characteristics which is commonly measured in experiments. We will focus on the $`I`$-$`V`$ characteristics in the present paper. As pointed out in Sec. II the RD to some extent may be considered as a simplified version of the RSJD. Thus from this point of view it is perhaps not surprising that these two models (as we will see) contain similar features of the vortex dynamics. In Ref.  from the study of a simple dynamic model of isolated magnetic particles in a uniform field it has been shown that the actual dynamics of the model and MCD are in a good agreement when the acceptance ratio of the Metropolis step is low enough. This implies that the MCD should give the same $`I`$-$`V`$ characteristics as the RSJD after an appropriately chosen normalization of time, when the trial angle $`\delta \theta `$ is sufficiently small (it was shown in Ref.  that $`\delta \theta =\pi /6`$ is sufficiently small). We will confirm this further in the present simulations.

In Fig. 1 we compare $`I`$-$`V`$ characteristics in the form of the electric field $`E`$ versus the current density $`J`$ obtained from RSJD, RD, and MCD simulations in the temperature range $`0.70T1.50`$ (the temperature interval is 0.05 if $`0.70T1.00`$ and 0.10 otherwise) in log scales for the system size $`L=8`$. In order to make $`I`$-$`V`$ data of RSJD and RD simulations coincide, $`E`$ obtained with RD ($`E_{RD}`$) is multiplied by a temperature-independent factor represented by the horizontal line in the inset of Fig. 1. Since the measured time-averaged electric field is inversely proportional to the time scale for a given dynamics, one can from the ratio $`E_{RSJD}/E_{RD}`$ infer the correspondence between times of RD and RSJD. ¿From this comparison, we find for $`0.70T1.00`$ that one unit of time in RD approximately corresponds to 0.526 time unit in RSJD, independent of the temperature. Note that the $`I`$-$`V`$ curves corresponding to the temperatures exceeding 1.00 are almost a straight lines. Therefore the collapse of data between the different dynamics is trivial for $`T>1.00`$. Also the $`I`$-$`V`$ characteristics obtained from MCD can be made to collapse on top of the corresponding curves for the RSJD and RD, as shown in Fig. 1. However, in this case the time scale factor, describing how many RSJD time steps one MCD step corresponds to, depends on the temperature, as shown in the inset of Fig. 1, where this factor is shown to be a linear function of temperature for system size $`L=6`$, 8, and 10 up to $`T=1.00`$. This is in accordance with the model studied in Ref. , where the same linear behavior in terms of the temperature has been found. Thus for the $`I`$-$`V`$ curves we have a precise relation between RSJD and MCD: For example, at $`T=0.80`$ we get $`1`$ MCS = $`10.2\times `$RSJD time unit. This opens a practical possibility to study some aspects of dynamic behavior of the $`XY`$ model by using MCD simulation, which is usually much more efficient than RSJD and RD.

We have shown in this section that the RSJD, the RD, and the MCD applied to the usual $`XY`$ model gives basically the identical $`I`$-$`V`$ characteristics up to some constant factors. From this observation, one can conclude that the dynamic critical behaviors of the $`XY`$ model inferred from the $`I`$-$`V`$ characteristics should be identical for all these dynamic models. We in next section use the RSJD to study the Villain model and $`XY`$ models with $`p=2`$ and presume from the observation in the present section that the conclusion drawn in Sec. IV for the RSJD case should be also valid in the other dynamics (RD and MCD).

## IV Various $`XY`$ models subject to RSJD

To study the critical behavior of the system in the vicinity of the transition one can use scaling relations. Fisher, Fisher, and Huse (FFH) in Ref.  proposed that the nonlinear $`I`$-$`V`$ characteristics in a $`D`$-dimensional superconductor scales as

$$E=J\xi ^{D2z}\chi _\pm (J\xi ^{D1}/T),$$
(12)
where $`\xi `$ and $`z`$ are the correlation length and the dynamic critical exponent, respectively, and $`\chi _\pm `$ is the scaling function above ($`+`$) and below ($``$) the transition. Pierson et al. in Ref.  have applied a variant of this FFH scaling approach to the $`I`$-$`V`$ data for thin ($`D=2`$) superconductors and superfluids and suggested a phase transition with $`z6`$. In Ref.  another scaling relation has been introduced and it has been shown that a certain finite-size effect which is not included in the FFH scaling may have caused the large $`z`$ in spite of the fact that the finite size effect precludes the possibility of a real phase transition. This finite-size scaling around and below the KT transition is given by the form (see Ref.  for the details)

$$\frac{E}{JR}=h_T\mathbf{\left(}JLg_L(T)\mathbf{\right)},$$
(13)
where $`R=lim_{J0}(E/J)L^{z(T)}`$ is a finite-size induced resistance without external current, $`h_T(0)=1`$, $`h(x)x^{z(T)}`$ for large $`x`$, and $`g_L(T)`$ is a function of at most $`T`$ and $`L`$ such that a finite limit function $`g_{\mathrm{}}(T)`$ exists in the large-$`L`$ limit. For small values of the variable $`x=JLg_L(T)`$ the $`T`$-dependence of the scaling is absorbed in a function $`g_L(T)`$ for each fixed size $`L`$ giving rise to the scaling form

$$\frac{E}{JR}=h\mathbf{\left(}JLg_L(T)\mathbf{\right)}.$$
(14)
In Fig. 2 we demonstrate the existence of the finite-size scaling given by Eq. (14) for $`L=8`$ within the temperature intervals $`0.70T1.10`$. The data are obtained for MCD and scaled by the appropriate factor so as to correspond to RSJD and RD (compare Fig. 1). Because of the correspondence between the three types of dynamics (see Sec. III), this also means that the existence of the finite-size scaling given by Eq. (14) is insensitive to the choice of dynamics.

The FFH scaling given by Eq. (12) is correct only in the thermodynamic limit $`\xi /L0`$. However, from a practical point of view there is a connection between Pierson method, which is based on FFH scaling, and the finite-size scaling introduced by Eq. (14). If we assume that $`\xi `$ in Eq. (12) is proportional to $`R^\alpha `$, where $`\alpha `$ is $`T`$-independent constant, the connection between the two different scaling approach becomes: $`g_L(T)=A_LR^\alpha `$ with $`A_L`$ being a constant which may depend on $`L`$. In Fig. 3 there are presented three different functions $`g_L(T)`$ corresponding to $`L=6`$, 8, and 10. These functions are determined from the condition that curves corresponding to the different temperatures should collapse when plotted as $`E/JR`$ vs $`LJg_L(T)`$. Since it is well established that the 2D $`XY`$ model on the square lattice has the KT transition at $`T_{KT}0.892`$ (Ref. ), Fig. 3 shows that $`g_L(T)`$ over a limited region in the vicinity of the KT transition is very well represented by the $`R^\alpha `$ with $`\alpha =1/6`$ for all investigated system sizes.

Next we demonstrate how the finite-size scaling given by Eq. (14) works for the $`I`$-$`V`$ data obtained by simulations of the modified $`XY`$ model. These simulations are done with RSJD. Fig. 4(a) verifies that this scaling indeed exists for the $`XY`$ model modified with the Villain type of potential introduced by Eq. (4). The inset of Fig. 4(a) shows the scaling function $`g_L(T)`$ determined by finding the best data collapse for small values of $`LJg_L(T)`$. One can see that in the vicinity of $`T_{KT}`$, which is approximately equals to $`1.4`$ for this model, $`g_L(T)`$ can be fitted by $`A_LR^\alpha `$ with $`\alpha =1/6`$. The data collapse in Fig. 4(b) shows that $`E/JR`$ is only a function of the scaling variable $`JLg_L(T)`$ (when $`JLg_L(T)`$ is small enough) for the $`XY`$ model modified with $`p`$-type of potential (Eq. (3)), where $`p=2`$. The inset shows the function $`g_L(T)`$ together with $`A_LR^\alpha `$. Since $`T_{KT}1.15`$ for this model one can again see that the scaling function $`g_L(T)`$ in the vicinity of KT transition is proportional to $`R^\alpha `$ with the same exponent $`\alpha =1/6`$ as for the original $`XY`$ and Villain models.

The crucial difference between the Villain model, the usual $`XY`$ model and the $`p=2`$ $`XY`$ model in the present context is the vortex density. The KT-transitions for these models occur at the Coulomb gas temperatures $`T_c^{CG}=0.23`$, 0.2, and 0.1, respectively ($`T^{CG}=\frac{T}{2\pi }U^{\prime \prime }`$ (see Ref. )). Lower $`T_c^{CG}`$ means higher vortex density. Thus the finite-size scaling property given by Eq. (14) appears to be independent of vortex density.

## V Discussion and conclusions

We have simulated 2D XY model with three types of dynamics: RSJD, RD, and MCD. The main conclusion of the paper is that the qualitative features of the finite-size scaling given by Eq. (14) are independent of both the vortex density and the type of dynamics. Therefore the finite-size scaling behavior given by Eq. (14) of the finite-size induced tails of the $`I`$-$`V`$ characteristics appears to be a robust feature. From the comparisons of the current-voltage characteristics obtained for each type of dynamics we found that, up to some scale factor, $`I`$-$`V`$ curves at a given temperature are identical over a broad range of external currents. This makes it possible to use MCD simulations, which are more computer efficient than RSJD and RD simulations, to obtain $`I`$-$`V`$ curves corresponding to both RSJD and RD dynamics.

The phase transition for the 2D $`XY`$-type models are of KT type with $`z=2`$. This raises the intriguing question of the origin of the large $`z6`$ obtained by Pierson et al. . In Ref.  it was argued that the Pierson scaling in relation to the finite-size scaling given by Eq. (14) corresponds to the proportionality $`g_L(T)R^\alpha `$ where $`1/\alpha `$ is the exponent which corresponds to the “$`z`$” obtained by the Pierson scaling. The reason for the existence of this scaling like behavior is still unclear.

In the present paper we have shown that, within the class of 2D $`XY`$-type models studied, a value $`1/\alpha 6`$ is obtained independently of the type of dynamics, as well as, of the vortex density. The origin of this seemingly robust behavior calls for further investigations.

## ACKNOWLEDGMENTS

This work was supported by the Swedish Natural Research Council through Contract No. FU 04040-332.


text:
# Diffusion in a strongly correlated anisotropic overlayer

## 1 Introduction

A key of understanding chemical reaction kinetics on surfaces is the understanding of diffusion of particles on surfaces. Surface diffusion also determines many other surface properties, e.g. in adsorbed films acting as lubricants. Nevertheless, our present understanding of surface diffusion is far from being complete. In particular, situations of highly anisotropic diffusion and the role of defects is far from being understood. Surface diffusion is usually described in terms of single atom diffusion (see, e.g. ). Especially for incommensurate layers, however, this description is not adequate for mass transfer by diffusion and for the very strong influence of frozen defects in this case (see ). In order to describe this situation properly, a collective “soliton diffusion” mechanism was proposed (see ). We show below that this approach can also be adapted to describe diffusion in chain structures on an anisotropic substrate.

Two-dimensional structures which are formed on surfaces with strongly anisotropic substrate potential corrugation like (112) bcc and (110) fcc surfaces have been studied by many authors (for a review, see ). Chain structures of the type $`(1\times p)`$ are common on these substrates. At coverage less than 1/2 monolayer, these structures consist of parallel chains of adsorbed atoms. A typical example is schematically shown in Fig. 2. In experiment $`p`$ can be as large as 8. The potential corrugation of the (112) bcc surface can be represented as troughs in the $`111`$ direction. The substrate potential is also modulated along the troughs with the substrate period. In a typical $`(1\times p)`$ structure, chains are oriented normal to the trough direction. For the very existence of a chain structure with a large period, the interaction along the chain $`J_0`$ should be much larger than the interaction between chains.

Chain structures demonstrate remarkable thermal stability (see ). Their transition temperatures are high, compared with the weak interaction between chains in the structures with large period. Though chain structures are known for a long time, the first studies of diffusion in these structures have appeared only recently. As a very interesting phenomenon, these experiments on surface diffusion in pure Li and Sr overlayers on the Mo(112) face , confirm the extremely strong anisotropy of the diffusion coefficients ($`10^3`$ times and more) and reveal a strong correlation between surface diffusion parameters and the equilibrium phases of the overlayer at temperatures much higher than the melting temperatures of these phases. In these experiments the macroscopic flow of atoms in the overlayer was actually measured, not the diffusion of single atoms. A typical system is Li on Mo(112). At low coverages it forms two chain structures $`(1\times 4)`$ at the coverage $`\theta =.25`$ and $`(1\times 2)`$ at $`\theta =.5`$. At higher coverages incommensurate structures appear.

Only recently  an influence of frozen point disorder (Dy adatoms) on diffusion of Li on Mo(112) at coverages $`1/4<\theta <1/2`$, which correspond to commensurate structures (at low temperatures), has been reported for a Dy coverage $`\theta _{Dy}=0.1`$. As expected, the diffusivity is reduced, but mainly by a drastic decrease of the pre-exponential factor $`D_0`$ (up $`10^3`$ times), not by changes in the activation energies.

A strong correlation between the low-temperature phases and the diffusion coefficient at high temperature has been found experimentally in many strongly interacting and strongly anisotropic overlayer systems . The pronounced influence of defects, on the contrary, has previously only been observed for incommensurate systems .

Motivated by these experimental findings, we study in this work theoretically the diffusion in chain structures on anisotropic substrates, and focus on an explanation for the correlation between diffusion and the ordered low-temperature structures and on the influence of defects. In experiment, these phenomena have been observed for diffusion along the troughs. For this reason we focus on diffusion in the same direction, although the atom motion in the simulations described below is allowed in both directions.

Since our study only intends to clarify the basic mechanism of diffusion in this class of systems, the details of interactions are not important at this point. In fact, we are interested in mass transfer on a scale larger than the correlation length along the chain. As will be seen below, there are clear indications that the soliton diffusion mechanism is indeed still effective at temperatures far above melting. Especially for this mechanism, due to its collective nature, the details of the microscopic model should not change qualitatively our results. For this reason, we formulate a simple model which preserves the most important properties - strong correlations along the chains and weak interaction between the chains. Since in experiment close relations between the ordered structures and diffusion properties have been found, we first briefly review in Sec. 2 the equilibrium properties of commensurate and incommensurate structures on these anisotropic surfaces (see, e.g., ref. and references therein) before addressing possible mechanisms of mass transport. We show that indeed soliton diffusion plays an important role on these surfaces, and it is this mechanism which makes the connection between ordered structures and diffusion at much higher temperatures than the melting temperatures understandable. In Sec. 3 we show that the results of Monte-Carlo simulations of mass transfer, using a simple anisotropic lattice gas model, corroborate most of the predictions obtained analytically.

## 2 Mass Transport in a Chain Structure.

### 2.1 Thermodynamics of the Chain Structure.

Consider first the thermal fluctuations of a single chain of atoms aligned along the $`OY`$ axis with period $`b`$. The chain is placed into a periodic potential with period $`a`$. The chain can move from one substrate potential minimum to another by creating kinks as shown in Fig. 1. If $`E_k`$ is the kink energy, then the mean distance between kinks is $`b\mathrm{exp}(E_k/T)`$. The mean square transverse displacement of a chain of length $`L`$ is:

$$(u(L)u(0))^2=\frac{2a^2L}{b\mathrm{exp}(E_k/T)}$$
(1)
Using Eq. 1 and the corresponding relation for a continuous string, one can substitute in the thermodynamic limit the chain of atoms by a continuous string with the Hamiltonian

$$_s=\frac{\epsilon }{2}𝑑y(\frac{u}{y})^2$$
(2)
where $`\epsilon =(Tb/2a^2)\mathrm{exp}(E_k/T)`$ is the rigidity.

The structure $`(1\times p)`$ can occupy $`p`$ different sub-lattices, and has $`Z_p`$ symmetry (the group of additions modulo $`p`$). The most studied models with this symmetry are the $`Z_p`$-model, the clock model, and the discrete Gaussian model (for a review see e.g.). The thermodynamics of the clock model has been studied in detail (see e.g.). For $`p>4`$ three phases were shown to exist in the model: long range order at $`0<T<T_p`$, algebraic decay of correlations at $`T_p<T<T_m`$, and a disordered phase for $`T>T_m`$, where $`T_m`$ is the melting temperature and $`T_p`$ the depinning temperature.

Above $`T_p`$ the correlation properties of the chain structure are the same as those of a lattice of parallel strings. On scales much larger than the average distance between kinks, this lattice can be described by the Hamiltonian:

$$_{ch}=\frac{1}{2}𝑑x𝑑y(K_1(\frac{u}{x})^2+K_2(\frac{u}{y})^2)$$
(3)
where

$$K_2=\frac{\epsilon }{l}$$
(4)
is the shear modulus. We assume repulsive interactions between chains, which decay as a power of the structure period $`l=pa`$. The compression modulus $`K_1`$ can be represented as a sum of two terms:

$$K_1=\frac{J_{}}{ab}\frac{1}{p^\alpha }+\frac{\pi ^2T^2}{\epsilon a^3p^3}$$
(5)
The first term in Eq. 5 is due to direct repulsive interaction between chains, and the second is due to entropic repulsion (see e.g.). In the simplest case of the repulsive interactions between adatoms, which decay as $`r^\beta `$, the exponent $`\alpha =\beta `$. For dipole-dipole repulsion or for repulsion due to substrate elastic deformations $`\beta =3`$ (i.e. $`\alpha =3`$). In the case of an anisotropic Fermi surface, for some range of distances, the value of $`\beta =1`$ . $`J_{}p^{\alpha +2}`$ is of the order of the energy change per atom in the chain due to the displacement of the chain by one lattice constant $`a`$ along the trough.

The melting temperature $`T_m`$ is given by:

$$T_m=\frac{\pi p^2a^2}{8}\sqrt{K_1K_2}$$
(6)
Using Eqs. 4 and 5

$$T_m\frac{E_k}{\mathrm{ln}(E_k/J_1)}$$
(7)
follows, where $`J_1=const.\times J_{}`$. In other words, for highly anisotropic chain structures the melting temperature is defined mainly by the kink energy $`E_k`$.

For the depinning temperature $`T_p`$ we have the exact relation :

$$T_p=(16/p^2)T_m$$
(8)
When $`p=4`$ we have $`T_p=T_m`$ and there is only one transition for $`p<4`$. Above $`T_p`$ the thermal fluctuations ”smear” the chain lattice over many potential minima of the substrate potential. The kink free energy is zero. As a result, the chains can move due to the movement of kinks which can appear simply from thermal fluctuations at these temperatures.

Above $`T_m`$ the system is in the liquid state. This liquid is highly anisotropic: the correlation radius $`r_c`$ is large along the chains: $`r_{c||}b`$, but small perpendicular to the chains $`r_cpa`$.

### 2.2 Frozen Defects.

The influence of non-equilibrium (frozen) defects of different kinds on the order and phase transitions has been studied both theoretically and experimentally for a long time in connection with magnetic systems (for review see, e.g., ref.) and particularly with the spin glass problem (for review see, e.g.,). The disorder was classified as random bond-like and random field-like. The random bond-like disorder interacts with the same strength with different domains of the pure system. On the contrary, the random field-like disorder distinguishes between different domains. Random bond-like disorder destroys the order only at some finite concentration. In the ordered phase random bond-like disorder changes the critical behavior. The lower the dimensionality of the system, the stronger are the effects of frozen defects. It is generally accepted that in two-dimensional systems the random field-like disorder destroys the long range order at arbitrarily small values of the random field strength even at T=0 (see e.g.). However, the mean size of the well ordered domains (i.e. the correlation radius) diverges exponentially when the disorder strength decreases, e.g. in the case of the two-dimensional random field Ising model (ref.).

In reality no surface, regardless of how carefully it has been prepared, is free of defects, which are of different origin and dimensionality (point or linear). The list includes impurities segregated onto the surface, adsorbed impurities, excess atoms of the crystal proper, steps, etc. At low temperatures, the mobility of these defects is extremely small, and they can be treated as frozen. These non-equilibrium defects give rise to various random fields on the crystalline surface.

The frozen point defects on a surface, depending on their position with respect to substrate and overlayer, correspond to random bonds in magnetic models (points with special symmetry) or, in general, to both random field and random bonds simultaneously.

It is generally accepted that the random field defects destroy order in the system with algebraically decaying correlation in spatial dimensions $`D2`$. This result was first obtained by Imry and Ma .

### 2.3 Mass Transfer Mechanisms.

Frozen disorder strongly affects the large scale dynamics of an overlayer. In particular, it can lead to slowing down the relaxation processes in the layer. The problem is far from being understood, but at least two different relaxation mechanisms are conceivable: (i) via elastic deformation of the chain lattice (at $`T>T_p`$) and (ii) via defect (dislocation) formation in the chain lattice in the overlayer (plastic flow). For the case of the incommensurate overlayer the latter mechanism, the soliton diffusion mechanism, has been proposed theoretically in order to describe diffusion experiments in overlayers (see ref. and ref. Ch.10) Experiments from different experimental groups and for different systems (refs.,  and ref. Ch.10) have shown the relevance of the soliton diffusion mechanism. We will apply the approach developed in refs., to analyze the mass transfer in chain structures.

#### Ideal substrate, $`T<T_p`$.

Below $`T_p`$ the chain can shift by kink diffusion. The activation energy $`E_a`$ in this case is equal to the sum of the energy needed to create a pair of kinks $`E_a=2E_k`$ and of the energy barrier $`E_{||}`$ for atom movement along the trough. The effective diffusion coefficient can be written in the form:

$$D_{eff}\omega _0b^2\mathrm{exp}(2E_k/T)\mathrm{exp}(E_{||}/T)=D_0\mathrm{exp}(E_a/T),$$
(9)
The characteristic attempt frequency $`\omega _0`$ is of the order of the Debye frequency of the overlayer $`\omega _D`$. If only nearest neighbor interactions are considered the activation barrier for kink diffusion at $`T<T_p`$ is close to the activation barrier for single atom diffusion.

#### Ideal substrate, $`T_p<T<T_m`$.

On the ideal substrate the depinning transition should result in a drastic change of the apparent activation energy of diffusion. As we have already mentioned, the kink free energy of the chain structure disappears above $`T_p`$ due to thermal fluctuations. This implies that in the limit of an infinitely long observation time the correlation properties of the chain lattice are the same as for elastic media. However, for finite times the lattice of chains needs to overcome activation barriers. As a result in the limit of infinitely small gradients (i.e. long observation time) the apparent activation energy should change from $`2E_k+E_{||}`$ to $`E_{||}`$. We are not aware of experimental observations of this situation.

#### Substrate with disorder, $`T_p<T<T_m`$.

A real substrate always has some amount of defects. We have already mentioned that frozen defects destroy the order at large scale. These defects can effectively pin the chain structure. The simplest defect is a point trap. We discuss the case when the depth of the trap is larger than the interaction between atoms in the chain. The nature of the trap and the depth of the potential well may vary. However, the energy barrier for depinning from traps of different depths will not vary very much. The main reason for this is the mechanism of depinning. E.g. for the situation described the chain can be depinned more easily by disrupting itself rather than by removing the atom in the trap. This means that the upper bound on the depinning energy is given by the energy needed to break the chain.

Consider now in more detail the mass transfer above $`T_p`$. The gradient of the surface density of adatoms induces a pressure gradient in the depinned chain lattice, which in turn generates effective deformations at chain pinning sites. Depinning of the chain from a defect is followed by a displacement by $`l=pa`$ of all chains in a region of characteristic size $`a/\sqrt{\theta _d}`$, where $`\theta _d`$ is the defect (trap) concentration per unit cell. Such a region can be treated as a particle moving under a pressure-gradient force $`f`$ that is proportional to the pressure difference $`\mathrm{\Delta }P`$ and to its transverse size:

$$f\mathrm{\Delta }P\frac{a}{\sqrt{\theta _d}};\mathrm{\Delta }P=K_1\frac{\mathrm{\Delta }l}{l},$$
(10)
where $`\mathrm{\Delta }l`$ is the average change of the structure period due to compression, and the effective compression modulus is given by Eq. 5. The average change of period, $`\mathrm{\Delta }l`$, can be evaluated from the coverage gradient $`\theta /x`$ in a region of the order of the distance between the defects. As a result we obtain:

$$\frac{\mathrm{\Delta }l}{l}\frac{l}{\theta }\frac{\theta }{x}\frac{a}{\sqrt{\theta _d}}=\frac{l}{\sqrt{\theta _d}}\frac{\theta }{x}.$$
(11)
we have used here the fact that $`l=a/\theta `$. Finally we obtain for the force:

$$f\frac{K_1la}{\theta _d}\frac{\theta }{x}$$
(12)
The conventional equation for the diffusion flow under a constant force (see e.g. ) results in the following expression for $`D_0`$:

$$D_0\omega _0l^2\frac{K_1la}{T\theta _d}$$
(13)
#### Mass transfer in the liquid ($`T>T_m`$).

At $`T>T_m`$ there is no order at large scale. Correlations along the troughs are short ranged, the correlation radius in this direction is not larger than a few substrate periods. On the other hand, the correlations along the chains can remain strong at temperatures far above $`T_m`$ with the correlation radius exceeding many substrate periods. These latter correlations will significantly influence the diffusion above $`T_m`$ so that one expects to find a close connection between the effective diffusion coefficient above $`T_m`$ and the equilibrium structures at $`T<T_m`$. At $`T>T_m`$ the overlayer can be considered as consisting of anisotropic “particles” with a length along the chains many times larger than the length along the troughs. Let us consider a single chain in such a “particle”. In the thermodynamic limit, i.e. for infinite observation times, the kink free energy is zero. As a result one has a finite single kink density. However, we discuss the system dynamics for times of the order of the relaxation times on an atomic scale. This means that kinks need to overcome local barriers to move from one point on the chain to another. As soon as a kink appears on the chain, e.g. at the end of the chain, it propagates along the chain. Without defects the barrier it should overcome is the potential barrier along the trough $`E_{}`$. As a result of kink motion, the chain moves by one substrate period along the troughs. In the case of a substrate with frozen defects the kink stops at the frozen defect. The chain will remain at the same place until the kink overcomes the defect. As a result, in the case of strong pinning by defects, the main contribution to the diffusion time is the residence time $`\tau _{dep}`$ on the defect. One can represent the ”particle” movement by a series of jumps of the order of the average distance between defects i.e. of $`a/\sqrt{\theta _d}`$ during the time interval defined by the residence time $`\tau _{dep}`$ on the defect. As a result, the effective diffusion coefficient scales with defect concentration as

$$D_{eff}(\frac{a}{\sqrt{\theta _d}})^2\frac{1}{\tau _{dep}}\frac{1}{\theta _d}\omega _0a^2\mathrm{exp}(E_a/T).$$
(14)
This result is valid as long as the correlation length along the chain is larger than the distance between defects. Eq. 14 assumes that the residence time $`\tau _{dep}`$ is much larger than the diffusion time between defects $`\tau _{dif}(\theta _d)`$. Both conditions mean that the diffusion behavior reduces to that found on a surface without defects once the concentration of defects falls below a certain value (which depends on temperature). This directly avoids the divergence of $`D_{eff}`$ in Eq. 14, i.e. Eq. 14 considers only the case when diffusion is fully limited by defects.

## 3 Monte-Carlo Simulations.

### 3.1 The Model.

In order to study the behavior of the effective diffusion coefficient, we use a procedure which simulates in a most efficient way the experimental situation. We will not consider diffusive motion of single particles. Instead, we study the particle current as a function of applied force. There is no force applied to each particle in a real experiment. However, there is a density gradient resulting in the gradient of the chemical potential. The chemical potential gradient changes very slowly on the atomic scale. This gradient plays the role of some average force which acts on the particles.

The simplest model of diffusion in a strongly interacting overlayer is the two-dimensional lattice gas Ising model with Kawasaki dynamics. The Hamiltonian of the Ising lattice gas under the applied field $`F`$ in $`x`$ direction has the form:

$$_b=\underset{𝐫,𝐚}{}J_𝐚(𝐫)n_𝐫n_{𝐫+𝐚}\underset{𝐫}{}F\frac{x}{|𝐚|}n_𝐫$$
(15)
where summation over $`𝐫=(x,y)`$ runs over the lattice sites, $`𝐚`$ labels nearest neighbors. Here $`n_𝐫=0,1`$ and $`J_𝐚(𝐫)`$ is measured in units of $`J_0>0`$.

We have used Monte-Carlo simulations with Kawasaki dynamics (see e.g. ) to study diffusion. To model the chain structure of the experimental systems  we have introduced a strong attractive nearest neighbor interaction ($`N`$) in the $`Y`$-direction, $`J_y=8J_0`$, and repulsive nearest ($`N`$) and next nearest ($`NN`$) neighbor interactions in the $`X`$-direction, $`J_N=2J_0`$ and $`J_{NN}=J_0`$. In order to simulate scenarios that are close to the experimental situation, we have measured the particle current under the applied force $`F`$ at temperatures well above the melting temperature. First we demonstrate that the particle current is linear in $`F`$. Fig. 2 shows this dependence of the current on $`F`$ for $`F.25J_0`$ at the lowest temperature used in our simulations. The linear dependence has been checked up to $`FJ_0`$ at higher temperatures. We have also checked that the applied force has no measurable influence on the order in the system. Fig. 3 shows the correlation function along the chains and normal to the chains in a $`24\times 192`$ system for the coverage $`\theta =1/2`$ after $`2\times 10^6`$ MC steps at $`T=4.765J_0`$. The two different sets of data correspond to the system with and without applied force. The data in this figure also directly demonstrate the remaining strong correlations above $`T_m`$, in particular along the chains. Even perpendicular to the chains short range order over several lattice constants still exists at this temperature, which is roughly 20% above $`T_m`$.

### 3.2 Diffusion and Equilibrium Phases.

At equilibrium the simulated system studied here has three ordered structures at low temperatures: $`(1\times 2)`$ (near a coverage $`\theta =1/2`$), $`(1\times 3)2X`$ (near a coverage $`\theta =2/3`$), and $`(1\times 3)`$ (near a coverage $`\theta =1/3`$). The melting temperature for the $`(1\times 2)`$ structure is $`T_{m,1/2}=3.86J_0`$. Those for the $`(1\times 3)2X`$ and $`(1\times 3)`$ structures are of the same order of magnitude. As mentioned in Sec. II, there is no depinning in these cases, only melting.

We have simulated diffusive particle flow for coverages corresponding to the structures $`p(1\times 2)`$ ($`\theta =1/2`$) and $`(1\times 3)2X`$ ($`\theta =2/3`$) for the temperature range $`4.55J_0<T<12.5J_0`$ ( $`1.18T_{m,1/2}<T<3.24T_{m,1/2}`$). As can be clearly seen from Fig. 4, which shows Arrhenius plots of the effective diffusion coefficient for the coverages 0.5 and 0.66, the effective diffusion coefficients are identical at the highest temperatures. Here the correlation radii in both directions are very small so that they do not play any important role. A snapshot of the situation for $`\mathrm{\Theta }=1/2`$ in a $`25\times 200`$ system after $`10^6`$ MC steps at $`T=11.35J_0=2.94T_{m,1/2}`$ is shown in Fig. 5, demonstrating the loss of correlations at this temperature. With decreasing temperature, however, the correlation radii increase, particularly along the chains, as just shown in Fig. 3. A snapshot, again for $`\theta =0.5`$, after $`10^6`$ MC steps at $`T=4.55J_0`$ ($`T=1.18T_{m,1/2}`$), shown in Fig. 6, underlines this property. As is seen in Fig. 4 these correlations now result in different diffusion coefficients for the two coverages considered. The different slopes of the Arrhenius plots in Fig. 4 are evidence for the different activation energies in these two cases. We feel that these simulations clearly demonstrate the close connection between the ordered low-temperature structures and the effective diffusion coefficient at temperatures well above the melting temperature. This connection is established by the correlations of rather long range which remain far above the melting temperature in these strongly anisotropic systems.

### 3.3 Diffusion and Point Defects.

In order to get direct information about the main mechanism of diffusion, we have also studied diffusion in the presence of a small concentration of strong pinning centers of the random field type, i.e. particles located at pinning centers have been fixed there. Due to the strong attractive interactions along the particle chains, correlation along the $`Y`$ axis remains strong well above the transition temperature also in presence of the defects. For simplicity, the same interaction parameters as for other particles have been used for these bound particles. This simplification can not change the physics we discuss in this article. To clarify this point consider the chain of adatoms in Fig. 7. If the pinned atom interacts with neighbors with the same interactions as unpinned atoms, the kink will appear by cutting the bond $`A`$ or $`B`$ on the chain and the activation energy will be $`|J_y|`$ (note that $`J_y<0`$). Let us assume that the pinned atom interacts with neighbors with interactions $`|J_{yd}||J_y|`$. In the case $`|J_{yd}||J_y|T`$, the kink will appear by cutting the bond $`C`$ or $`D`$ on the chain. This means that the activation energy of the kink will remain the same as in the case $`J_{yd}=J_y`$. but the defect “extends” its size along the chain, i.e. the pre-exponential factor in the diffusion coefficient will decrease. In the opposite case $`|J_{yd}|<|J_y|`$ the chain can also not move through the defect atom which is infinitely strongly bound with the substrate. This is true even in the case $`J_{yd}=0`$. The activation energy in this case remains $`|J_y|`$ which is the kink energy on the chain.

The above discussion shows that activation energy will not change significantly by varying strength of the defect interaction and the use of the interaction set of non-pinned atoms is a reasonable and convenient simplification. If diffusion is now limited mainly by the residence time of the chains at the defect sites, we expect that the activation energy is also not changed significantly by varying defect concentrations, but a strong influence on the pre-exponential factor is expected. If the soliton mechanism is effective, Eq. 14 should hold as long as the defect concentration does not become too small (see previous section).

The behavior just described is indeed found in our simulations, as Fig. 8 demonstrates. In this figure, Arrhenius plots of the diffusion coefficient are shown for a coverage $`\theta =1/2`$ without point defects, and for concentrations of point defect varying between (1/64) to (1/16). First, please note that already the smallest concentration of point defects (1/64) results in an observable change of the diffusion coefficient. This is an important qualitative confirmation of the collective nature of diffusion in our model because it directly contradicts the behavior expected for a single particle diffusion mechanism. If mass transport would happen mainly via single particle diffusion, these pinning centers would not be able to change significantly the mass transport. Only for extended mass carriers can these rare pinning centers have a strong effect on diffusion.

Second, as expected from our qualitative arguments given above, the activation energy in presence of defects does not differ from the case of the ideal substrate, and does not vary at different defect concentrations. To show this more clearly, we have fitted the simulation data by straight lines in Fig. 8. These lines are also included in this figure.

Third, the pre-exponential factors strongly decrease with increasing defect concentration. The resulting dependence of the ratio of the pre-exponential factors $`D(0)/D(\theta _d)`$ of the diffusion constants as a function of the point-defect concentration is presented in Fig. 9. The pre-exponential factor depends on $`\theta _d`$ as $`D(\theta _d)=D(0)/(18.56\theta _d+1)`$. The finite limit of $`D(\theta _d)`$ when $`\theta _d0`$ is due to the finite diffusion time between defects which is comparable, at lowest defect density, with the residence time $`\tau _{dep}`$ on defects. These results fully agree with our qualitative expectations outlined above. Therefore, we think that the results of our simulations can be regarded as a confirmation of the soliton diffusion model.

## 4 Conclusion.

To conclude we have studied the mass transfer in a system with highly anisotropic interactions, which result in chain structures at low temperature. We have used both an analytical approach and numerical simulations. We propose that a mechanism of mass transfer similar to soliton diffusion is effective even at temperatures much higher than the disordering (melting) temperature. Our computer simulations corroborate this proposal: for a given coverage we find a connection between the effective diffusion coefficient at temperatures above the melting temperature and the phase that is stable at this coverage for low temperatures. An important test of this diffusion mechanism is the dependence of the effective diffusion coefficient on concentration of frozen defects. For the case of defects which are deep traps we predict that the activation energy depends weakly on the concentration of these defects, whereas the pre-exponential factor is inversely proportional to it. Although the influence of frozen disorder should be most pronounced for temperatures below the melting transition, it is still observable above this transition. This prediction has also been confirmed by Monte-Carlo simulations, which show that already a small concentration of frozen point disorder strongly influences the diffusion coefficient due to the dominant mechanism of soliton-like diffusion. This strong influence of small concentrations of defects on the diffusion coefficient would not be understandable for a predominant single particle diffusion mechanism.

The above results are not sensitive to the details of the adatom-adatom interaction and to the microscopic model describing adatom jumps from one site to another. This is both a strong and a weak point of our approach. Although our general approach shows that this mechanism can be tested experimentally by introducing small concentrations of frozen defects, significantly more detailed knowledge is needed in order to be able to describe a specific experimental system quantitatively. However, our “phenomenological” model can be used much more generally for systems which have similar correlation properties, but presumably totally different microscopic mechanisms of mass transfer. One example are vicinal surfaces. In this case we have strong correlations along the steps, relatively weak interaction between them and kinks as the elementary exitations, which are responsible for step motion. However, the microscopic mechanisms of the kink motion on the step and on the adatom chain are totally different. We plan to discuss the above similarity elsewhere.

## 5 Acknowledgments

We benefitted from discussions with A.G.Naumovets. The work is supported by the Niedersächsische Ministerium für Wissenschaft and Kultur and by the Volkswagen Stiftung. One of us (I.L.) was partly supported by the grants DE-FG03-96ER45598, NSF DMR-97-05182, THECB ARP 010366-003.


text:
# Predictions for Neutrino Structure Functions

## 1 Introduction

Deep inelastic lepton-nucleon scattering experiments have provided precision information about the quark distributions in the nucleon. However, there has been a long-standing discrepancy between the $`F_2`$ structure functions extracted from neutrino and muon experiments in the small $`x`$ range. Recently, a new analysis of differential cross sections and structure functions from CCFR $`\nu _\mu `$-Fe and $`\overline{\nu }_\mu `$-Fe data was presented; in this study, the neutrino-muon difference is resolved by extracting the $`\nu _\mu `$ structure functions in a physics model independent way.

In previous analyses of $`\nu _\mu `$ data, structure functions were extracted by applying a slow rescaling correction to correct for the charm mass suppression in the final state. In addition, the $`\mathrm{\Delta }xF_3`$ term (used as input in the extraction) was calculated from a leading order charm production model. These resulted in physics model dependent (PMD) structure functions. In the new analysis, slow rescaling corrections are not applied, and $`\mathrm{\Delta }xF_3`$ and $`F_2`$ were extracted from two parameter fits to the data.

The extracted physics model independent (PMI) values for $`F_2^\nu `$ are then compared with $`F_2^\mu `$ within the framework of NLO models for massive charm production; these are found to be in agreement, thus resolving the long-standing discrepancy between the two sets of data. However, the first measurements of $`\mathrm{\Delta }xF_3`$ are higher than current theoretical predictions. The objective of this paper is to investigate the sensitivity of $`\mathrm{\Delta }xF_3`$ upon a variety of factors including renormalization scheme and scale, quark mass effects, higher twist, isospin violation, and PDF uncertainties.

In Fig. 1 we have plotted the quantity $`\mathrm{\Delta }xF_3`$ for an isoscalar target computed to order $`\alpha _s^1`$. For comparison, we also display data from the CCFR analysis. While there is much freedom in the theoretical calculation, the difference between these calculations and the data at low $`Q`$ values warrants further investigation. We will discuss and compare the different theoretical calculations, and examine the inherent uncertainty in each with respect to different input parameters. We will also examine the experimental input, and assess uncertainties in this data.

### 1.1 Measurement of $`\mathrm{\Delta }xF_3`$

The structure functions are defined in terms of the neutrino-nucleon cross section via:

$`{\displaystyle \frac{d^2\sigma ^{\nu ,\overline{\nu }}}{dxdy}}={\displaystyle \frac{G_F^2ME_\nu }{\pi }}\left[F_2(1y)+xF_1y^2\pm xF_3y(1{\displaystyle \frac{y}{2}})\right]`$
(1)
where $`G_F`$ is the Fermi weak coupling constant, $`M`$ is the nucleon mass, $`E_\nu `$ is the incident energy, $`y=E_h/E_\nu `$ is the fractional energy transfer, and $`E_h`$ is the final state hadronic energy.

The sum of $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ differential cross sections for charged current interactions on an isoscalar target is then:

$`F(ϵ)`$ $``$ $`\left[{\displaystyle \frac{d^2\sigma ^\nu }{dxdy}}+{\displaystyle \frac{d^2\sigma ^{\overline{\nu }}}{dxdy}}\right]{\displaystyle \frac{(1ϵ)\pi }{y^2G_F^2ME_\nu }}`$ (2)
$`=`$ $`2xF_1[1+ϵR]+{\displaystyle \frac{y(1y/2)}{1+(1y)^2}}\mathrm{\Delta }xF_3`$
where $`ϵ2(1y)/[1+(1y)^2]`$ is the polarization of the virtual $`W`$ boson. In the above equation, we have used the relation:

$$2xF_1(x,Q^2)=F_2(x,Q^2)\frac{1+4M^2x^2/Q^2}{1+R(x,Q^2)}$$
(3)
where $`R=\sigma _L/\sigma _T`$ is the ratio of the cross-sections of longitudinally- to transversely-polarized $`W`$-bosons, $`Q^2`$ is the square of the four-momentum transfer to the nucleon, and $`x=Q^2/2ME_h`$ is the Bjorken scaling variable. For $`x<0.1`$, $`R`$ in neutrino scattering is expected to be somewhat larger than $`R`$ for muon scattering because of the production of massive charm quarks in the final state for the charged current neutrino production.

Using Eq. (2) and Eq. (3), $`F_2`$ and $`\mathrm{\Delta }xF_3`$ can be extracted separately. Because of the positive correlation between $`R`$ and $`\mathrm{\Delta }xF_3`$, the extracted values of $`F_2`$ are rather insensitive to the input $`R`$. If a large input $`R`$ is used, a larger value of $`xF_3`$ is extracted from the $`y`$ distribution, thus yielding the same value of $`F_2`$. In contrast, the extracted values of $`\mathrm{\Delta }xF_3`$ are sensitive to the assumed value of $`R`$, which yields a larger systematic error, shown on the data. For $`R`$ the QCD-inspired fit in  is used, but corrected for charged current neutrino scattering using a leading order slow rescaling model. This gives precisely the same type of correction as a full NLO calculation including the massive charm quark, as shown in Ref. , but leads to a somewhat higher normalization than the perturbative correction. It is arguable which prescription for $`R`$ leads to the better fit to existing data, but the difference between the two leads to an uncertainty in $`\mathrm{\Delta }xF_3`$ of the order of the systematic error shown. Clearly a further reduction in the assumed value of R (even down to zero), as suggested by models including $`\mathrm{ln}(1/x)`$ resummations, would still leave a discrepancy between the lowest $`Q^2`$ data points and the theoretical predictions.

### 1.2 Quark Parton Model Relations

Now that we have outlined the experimental method used in the extraction of the $`F`$ structure functions, it is instructive to recall the simple leading-order correspondence between the $`F`$’s and the PDF’s:000To exhibit the basic structure, the above is taken in the limit of 4 quarks, a symmetric sea, and a vanishing Cabibbo angle. Of course, the actual analysis takes into account the full structure.

$`F_2^{(\nu ,\overline{\nu })N}`$ $``$ $`x\left\{u+\overline{u}+d+\overline{d}+2s+2c\right\}`$
$`xF_3^{(\nu ,\overline{\nu })N}`$ $``$ $`x\left\{u\overline{u}+d\overline{d}\pm 2s2c\right\}`$ (4)
Therefore, the combination $`\mathrm{\Delta }xF_3`$ yields:

$$\mathrm{\Delta }xF_3=xF_3^{\nu N}xF_3^{\overline{\nu }N}4x\{sc\}.$$
(5)
Note that since this quantity involves the parity violating structure function $`F_3`$, this measurement has no analogue in the neutral current photon-exchange process. Also note that since, at leading-order, $`\mathrm{\Delta }xF_3`$ is directly sensitive to the strange and charm distributions, this observable can be used to probe the heavy quark PDF’s, and to understand heavy quark (charm) production. We discuss these possibilities further in the following subsection.

### 1.3 Implications for PDF’s

We have illustrated in Eq. (5) how $`\mathrm{\Delta }xF_3`$ is closely tied to the heavy quark PDF’s. The questions is: given the present knowledge base, should we use $`\mathrm{\Delta }xF_3`$ to determine the heavy quark PDF’S, or vice verse. To answer this question, we briefly review present measurements of heavy quark PDF’s, and assess their uncertainty.

#### 1.3.1 Tevatron $`W+Q`$ Production

The precise measurement of $`W`$ plus heavy quark $`(W+Q)`$ events provides important information on heavy quark PDF’s; additionally, such signals are a background for Higgs and squark searches.

Unfortunately, a primary uncertainty for $`W+Q`$ production comes from the heavy quark PDF’s. Given that $`\mathrm{\Delta }xF_3`$ is sensitive to these heavy quark PDF’s, we see at least two scenarios. One possibility is that new analysis of present data will resolve this situation prior to Run II, and provide precise distributions as an input to the Tevatron data analysis. If the situation remains unresolved, then new data from Run II may help to finally solve this puzzle. In the far future, a neutrino experiment from the neutrino factory at a linear collider would be an ideal tool to measure any neutrino structure function.

#### 1.3.2 DIS Di-muon production

The strange distribution is directly measured by dimuon production in neutrino-nucleon scattering.111Presently, there are a number of LO analyses, and one NLO ACOT analysis.. Results of a recent LO analysis by NOMAD  are in line with these experiments. The basic sub-process is $`\nu s\mu ^{}cX`$ with a subsequent charm decay $`c\mu ^+X^{}`$.

DIS dimuon measurements have safely established the breaking of the $`SU(3)`$ flavor symmetry,

$$\kappa \frac{𝑑xxs(x)}{𝑑xx[\overline{u}(x)+\overline{d}(x)]/2}|_{Q^2=𝒪(10\mathrm{G}\mathrm{e}\mathrm{V}^2)}\frac{1}{2}$$
(6)
in the nucleon sea. Still, there remain large uncertainties for the $`s`$-quark distribution in the kinematic regions relevant for $`\mathrm{\Delta }xF_3`$, even though the $`(x,Q)`$-range of the CCFR dimuon measurements  is comparable. CCFR recorded 5044 $`\nu _\mu `$ and 1062 $`\overline{\nu }_\mu `$ $`\mu ^{}\mu ^\pm `$ events with $`P_{\mu _1}9GeV`$, $`P_{\mu _2}5GeV`$, $`30E_\nu 600GeV`$, and $`Q^2=22.2GeV^2`$. The more recent NuTeV experiment recorded a similar sample of events, and these are presently being analysed through a MC simulation based on NLO quark- and gluon-initiated corrections at differential level. A complete NLO analysis of this data together with a global analysis may help to further constraint the strange quark distribution.

At present, PDF sets take strangeness suppresssion into account by imposing the constraint in Eq. (6) as $`s(x)=\kappa [\overline{u}(x)+\overline{d}(x)]/2`$ at the PDF input scale $`Q_01\mathrm{GeV}`$,  or by evolving $`s(x)`$ from a vanishing input at a lower scale. The residual uncertainty can be large, as can be seen from the collection of strange seas in Fig. 2.

In the following, we will take the uncertainty in the strange PDF into account by relaxing the experimental constraint in Eq. (6) as implemented in the MRST partons. This will be discussed in Sec. \[3.2\].

#### 1.3.3 Charged and Neutral Current DIS

The strange distribution can also be extracted indirectly using a combination of charged ($`F_2^\nu `$) and neutral ($`F_2^\mu `$) current structure functions; however, the systematic uncertainties involved in this procedure make an accurate determination difficult. The basic idea is to use the (leading-order) relation

$$\frac{F_2^\mu }{F_2^\nu }\frac{5}{18}\left\{1\frac{3}{5}\frac{(s+\overline{s})(c+\overline{c})+\mathrm{}}{q+\overline{q}}\right\}$$
(7)
to extract the strange distribution. Here, $`q+\overline{q}`$ represents a sum over all quark flavors. This method is complicated by a number of issues including the $`xF_3`$ component which can play a crucial role in the small-$`x`$ region—precisely the region where we observe the discrepancy. From the corresponding relation:

$$\frac{5}{18}F_2^\nu F_2^\mu \frac{1}{12}\mathrm{\Delta }xF_3$$
(8)
we see that these problems are not independent; however, this information, together with the exclusive dimuon events, may provide a more precise determination of the strange quark sea, and help to resolve our puzzle.

Prior to the DIS dimuon data, the 1992 CTEQ1 analysis  found that a combination of NC structure functions from NMC and the physics-model-dependent charged current structure functions from CCFR  seem compatible with approximate $`SU(3)`$-symmetry, i.e., $`\kappa 1`$ in Eq. (6). Recent dimuon measurements now exclude an $`SU(3)`$-symmetry $`\kappa =1`$.

We shall explore the effect of $`\kappa `$ on $`\mathrm{\Delta }xF_3`$ in Sec. \[3.2\].

## 2 Dependence of $`\mathrm{\Delta }xF_3`$ on input parameters

We now systematically investigate the sensitivity of the theoretical predictions of $`\mathrm{\Delta }xF_3`$ upon a variety of factors including: renormalization scheme and scale, quark mass effects, higher twist, isospin violation, and PDF uncertainties.

To simplify this analysis, we first examine the influence of these factors on the LO expression: $`\mathrm{\Delta }xF_3=4x\{s(\xi ,\mu )c(x,\mu )\}`$; after using this as a “toy model,” we will then return to the full NLO calculation in the next section.222Note that a possible $`s\overline{s}`$ asymmetry  would average out in $`\mathrm{\Delta }xF_3`$. For most variables, the simplified LO is sufficient to display the general behavior of the full NLO result. There are two exceptions: 1) the scheme dependence, and 2) the PDF dependence. These factors depend on the interplay of both the quark-initiated LO contributions as well as the NLO gluon-initiated contributions. For this reason, we will postpone discussion of these effects until the following section.

### 2.1 Charm Mass

We start by examining the effect of the charm mass, $`m_c`$, on $`\mathrm{\Delta }xF_3`$. In Fig. 3, we plot the leading-order expression $`\mathrm{\Delta }xF_3=4x\{s(\xi ,\mu )c(x,\mu )\}`$ vs. $`Q^2`$ for values of $`m_c`$ in the range $`=[0,1.8]`$ GeV using steps of 0.2 GeV. Here, we define $`\xi =x(1+m_c^2/Q^2)`$ which is a “slow-rescaling” type of correction which (crudely) includes mass effects by shifting the $`x`$ variable. “Slow rescaling” naturally arises at LO from the single-charm threshold condition($`W>m_c`$) and is required at NLO for consistent mass factorization .

Note, the result of this correction is most significant at low $`Q^2`$. To isolate the kinematic $`m_c`$ dependence, we have chosen the scale $`\mu ^2=Q^2+1.6^2GeV^2`$. (We will separately vary the scale in the following subsection.)

Note that this exercise is only altering the charm mass in one aspect of the calculation; to be entirely consistent it would be necessary to obtain parton distributions (particularly the charm quark) using fits with different charm masses. For charm masses in the range 1.2 to 1.8 GeV, this will be a small effect; for charm masses below $`1.2`$ GeV, (which is below the experimentally allowed range), such issues become important and the curves of Fig. 3 will be modified, i.e. lower $`m_c`$ will lead to longer evolution for charm, a larger charm distribution, and a lowering of the curves in Fig. 3.

While taking $`m_c0`$ does raise the theoretical curves in the regions where we observe the discrepancy (namely, the low $`Q`$ region), varying $`m_c`$ even within the wide range of $`[1.2,1.8]`$ GeV (lower 4 curves) does not give us sufficient flexibility to match either the shape or normalization of the data.

### 2.2 $`\mu `$ Scale

Next, we investigate the variation of $`\mathrm{\Delta }xF_3=4x\{s(\xi ,\mu )c(x,\mu )\}`$ with the renormalization/factorization scale $`\mu `$. We use three choices of the $`\mu ^2`$ scale:

* $`\mu ^2=Q^2`$ ,
* $`\mu ^2=Q^2+m_c^2`$ ,
* $`\mu ^2=P_{T_{max}}^2`$ .

Of course, $`\mathrm{\Delta }xF_3`$ with $`\mu ^2`$ scales of $`Q^2`$ and $`Q^2+m_c^2`$ differ only at lower values of $`Q^2`$; $`\mathrm{\Delta }xF_3`$ with $`\mu ^2`$ scale of $`P_{T_{max}}^2`$ is comparable to $`Q^2`$ and $`Q^2+m_c^2`$ at larger $`x`$, but lies below for smaller $`x`$. The scale choice $`Q^2+m_c^2`$ leads to an improvment over $`Q^2`$ by providing a lower bound on $`\mu `$ to keep the scale in the perturbative region.

The choice of $`P_{T_{max}}`$ is motivated, in part, by some observations by Collins. To display the relationship between $`P_{T_{max}}`$ and $`Q`$, in Fig. 5 we plot both $`P_{T_{max}}`$ and $`Q`$ vs. $`x`$ for 4 choices of $`y`$; note that the $`x`$-dependence of $`P_{T_{max}}`$ is opposite that of $`Q`$. While tying the scale choice of $`\mu `$ to $`P_{T_{max}}`$ has some interesting intuitive interpretations,333Technically, the interpretation is in terms of the characteristic $`P_T`$ of the partonic subprocess, but this is unobservable; therefore $`P_{T_{max}}`$ is used instead. for small $`x`$ this scale clearly becomes too large for the relevant physics. In any case, it cannot help us with our $`\mathrm{\Delta }xF3`$ problem as the scale choice of $`P_{T_{max}}`$ moves the theory curves away from the data.

### 2.3 Higher Twist

We now illustrate the potential effects due to higher twist contributions. We parameterize such contributions by multiplying the leading-twist terms by a correction factor that will increase at low $`Q^2`$. More specifically we use:

$$F_i^{HT}(x,Q^2)=F_i^{LT}(x,Q^2)\left(1+\frac{D_i(x)}{Q^2}\right)$$
(9)
where $`i=\{1,2,3\}`$. We vary $`D_3(x)`$ over the range  $`GeV^2`$ in steps of $`1GeV^2`$. (For the purposes of our simple illustration, it is sufficient to take $`D_3(x)`$ to be independent of $`x`$.) We will find that this range is well beyond what is allowed by experiment; however, we display this exaggerated range to make the effect of the higher twist contributions evident.

To normalize this choice with the allowable range consistent with data, we compare with the MRST higher twist fit which extracted a limit on the function $`D_2(x)`$. We see from the table of $`D_2(x)`$ that the allowed contribution from the higher twist terms is quite small. In addition, we note that the sign for $`D_2(x)`$ obtained for the relevant small-$`x`$ region tends to be negative–exactly the opposite sign that is needed to move the theory toward the data.

As we have shown constraints on $`D_2(x)`$, the obvious question is should we expect $`D_3(x)`$ to be substantially different? A calculation of the power corrections using renormalons suggests that while both $`D_1(x)`$ and $`D_3(x)`$ are of the same order of magnitude as $`D_2(x)`$ and have similar $`x`$-dependence, that $`D_1(x)`$ and $`D_3(x)`$ are even more negative than $`D_2(x)`$; again, this trend would move the theory farther from the data. Therefore, we conclude that that $`D_2(x)`$ represents a conservative limit for the $`D_3(x)`$ power corrections.

Examining Fig. 6, we note that it would take an enormous higher twist contribution to bring the normalization of the curves in the range of the data points, and still the shape of the $`Q^2`$-dependence is not well matched; hence, we conclude that this is not a compelling solution.

### 2.4 Isospin Violations

The naive parton model identity in Eq. (5) is modified if the full (non-diagonal) CKM structure, NLO QCD radiative corrections, and QCD based charm production are taken into account. This expression is also modified even in leading-order if we have a violation of exact $`pn`$ isospin-symmetry, (or charge symmetry); e.g., $`u_n(x)d_p(x)`$. In deriving Eq. (5), isospin-symmetry was necessary to guarantee that the $`u,d`$-contributions cancel out in the difference, thereby leaving only the $`s(x)`$ and $`c(x)`$ contributions.444Note we have not investigated shadowing corrections. The CCFR data on Fe is converted with an isoscalar correction, and the corresponding uncertainties are included in the data.

The validity of exact charge symmetry (CS) has recently been reexamined. Residual $`u,d`$-contributions to $`\mathrm{\Delta }xF_3`$ from charge symmetry violation (CSV) would be amplified due to enhanced valence components $`\{u_v(x),d_v(x)\}`$, and because the $`du`$ transitions are not subject to slow-rescaling corrections which strongly suppress the $`sc`$ contribution to $`\mathrm{\Delta }xF_3`$, as illustrated in Fig. 3.

We will examine possible contributions to $`\mathrm{\Delta }xF_3`$ by considering a “toy” model to parameterize the CS-violations. In isospin space, we can parameterize a general transformation as a rotation:

$`|q_n^{\mathrm{CSV}}`$ $`=`$ $`N_q{\displaystyle \underset{q^{}}{}}R_{qq^{}}(\theta )|q^{}_p`$ (10)
where $`R`$ is a rotation matrix, and $`N_q`$ is the normalization factor. For example, in this model the $`u`$-distribution in the neutron would be related to the proton distributions via the relation:

$`u_n^{\mathrm{CSV}}(x,Q^2)=`$
$`N_u^2\left[\mathrm{cos}^2(\theta )u_p(x,Q^2)+\mathrm{sin}^2(\theta )d_p(x,Q^2)\right]`$ (11)
For $`\theta =\pi /2`$, we recover the symmetric limit $`u_p(x,Q^2)=d_n(x,Q^2)`$. Note that we define the normalization $`N_u`$ such that this model preserves the sum rule:

$$𝑑xu_n^{\mathrm{CSV}}(x,Q^2)\overline{u}_n^{\mathrm{CSV}}(x,Q^2)=1.$$
(12)
Note that Eq. (10) should not be considered a serious theory,555E.g., it does strictly speaking not commute with evolution. but rather a simple one-parameter ($`\theta `$) model which is flexible enough to illustrate a range of CSV effects for $`\mathrm{\Delta }xF_3`$.

In Fig. 7 we vary $`\theta `$ over the its maximum range $`[0,\pi /2]`$ in steps of $`\pi /20`$. The exact charge symmetry limit ($`\theta _0=\pi /2`$) of Eq. (2.4) corresponds to the lowest curve in Fig. 7. Note the effect of the CSV contribution monotonically increases as $`\theta `$ deviates from the charge symmetry (CS) limit $`\theta _0=\pi /2`$. From the plot we observe that a violation of $`\theta _0\theta (3/10)(\pi /2)`$ is required to bring the theory in the neighborhood of the data.

At the relevant $`x`$ values of interest for $`\mathrm{\Delta }xF_3`$, this translates (via Eq. (2.4) and its analogue for $`d_n`$) into a $`20\%`$ symmetry violation for $`u_n`$, and a $`10\%`$ symmetry violation for $`d_n`$. Specifically,

$`{\displaystyle \frac{d_p(x)}{u_n^{\mathrm{CSV}}(x)}}`$ $``$ $`1.2`$ (13)
$`{\displaystyle \frac{u_p(x)}{d_n^{\mathrm{CSV}}(x)}}`$ $``$ $`0.9.`$ (14)
Since the electroweak couplings are flavor-independent, ($`V_q=A_q=1q`$), $`\mathrm{\Delta }xF_3`$ is in principle insensitive to a re-shuffling of the CSV contributions between $`u_n`$ and $`d_n`$. In particular, if we define the shift due to CSV in $`\mathrm{\Delta }xF_3`$ as $`\delta \mathrm{\Delta }xF_3`$, we find:

$$\delta \mathrm{\Delta }xF_3\frac{x}{4}(\delta d+\delta \overline{d}\delta u\delta \overline{u})$$
(15)
where

$$\delta u(x,Q^2)u_n^{\mathrm{CSV}}(x,Q^2)d_p(x,Q^2).$$
(16)
The expression for $`\delta \mathrm{\Delta }xF_3(x,Q^2)`$ of Eq. (15) is evaluated at $`Q^2=4\mathrm{GeV}^2`$ from Eq. (2.4) with $`\theta _0\theta =(3/10)(\pi /2)`$ and is plotted in Fig. 8.666 Again, the detailed $`x`$-shape from Eq. (2.4) should not be taken as a serious model prediction.

Although this level of isospin violation certainly improves the description of $`\mathrm{\Delta }xF_3`$, it is necessary to consider precisely what level of violation is actually allowed by other experimental data. For instance, it has previously been suggested that the discrepancy between $`F_2`$ from neutrino and muon data itself may be due to isospin violation. This type of violation required that $`\overline{u}`$ be very different to $`\overline{d}`$ in the region of interest.

The measurement of the lepton charge asymmetry in W decays from the Tevatron  places tight constraints on the up and down quark distributions in the range $`0.007<x<0.24`$, constraining them to be approximately as specified in the parton sets obtained by the global analyses. While only strictly telling us about parton distributions in the proton, this data rules out an isospin violation of this type to about $`5\%`$ as demonstrated in .

However, there are other strong constraints on isospin violation. For example, we note that while the toy model above leaves the neutron singlet combination $`q+\overline{q}`$ invariant at the $`2\%`$ level in the region $`xϵ[0.01;0.1]`$, it would lower the NC observable:

$$\left[\frac{4}{9}(u+\overline{u})+\frac{1}{9}(d+\overline{d})\right]_n|_{xϵ[0.01;0.1]}$$
(17)
by about 10%. An effect of this size would definitely be visible in NMC $`F_2^n/F_2^p`$ data which has an uncertainty of order a few percent in this kinematic region, and acts as a major constraint.

At this point, one could play clever games to evade the constraints of specific experiments. For example, a re-shuffling of CSV contributions between the individual $`\delta q,\delta \overline{q}`$ in Eq. (15) according to:

$`xu_n`$ $`=xd_p{\displaystyle \frac{2}{5}}\delta \mathrm{\Delta }xF_3`$
$`xd_n`$ $`=xu_p+{\displaystyle \frac{8}{5}}\delta \mathrm{\Delta }xF_3`$
$`x\overline{q}_n`$ $`:\text{analogous}`$
would keep Eq. (17) invariant. However, this would in turn raise the CC observable Eq. (15) by $`5\%`$; though it would help to explain the excess in $`\delta xF_3`$, it would spoil the new-found compatibility between neutral current and charged current data.

In addition, there are also fixed-target Drell-Yan experiments such as NA51 and E866 which precisely measure $`\overline{d}/\overline{u}`$ in the range $`0.04<x<0.27`$, and are also sensitive to isospin violating effects.

We therefore conclude that the many precise data sets which constrain different combinations of the PDF’s probably leave no room for CSV contributions of the magnitude necessary to fully align the theory curves with the $`\mathrm{\Delta }xF_3`$ data. However, we should add the caveat that an exhaustive investigation of the interplay of these different data sets and their influence on $`\mathrm{\Delta }xF_3`$ will only be possible within a global PDF analysis

## 3 NLO calculation of $`\mathrm{\Delta }xF_3/2`$

Now that we have used the leading-order expression for $`\mathrm{\Delta }xF_3`$ to systematically investigate dependence of this observable on various parameters, in this section we now turn to the full NLO calculation.

### 3.1 Contributions to the NLO Calculation

In Fig. 9, we have plotted the LO and NLO calculations for $`\mathrm{\Delta }xF_3`$ vs. $`Q^2`$ on an isoscalar target in 3 and 4 flavor schemes. The 3-flavor LO calculation (f=3, $`\alpha _s^0`$) involves primarily the strange quark contribution, $`s(x)`$, as the charm distribution is excluded in this case. When the higher order terms are included (f=3, $`\alpha _s^1`$), this result moves (substantially) toward the predictions of the 4-flavor scheme. We note that while the 3-flavor LO calculation (f=3, $`\alpha _s^0`$) appears consistent with the data, we cannot take this result as a precise theoretical prediction as this simplistic result is highly dependent on scheme and scale choices; a result that is verified by the large shift in going from LO to NLO.

The pair of curves in the 4-flavor scheme (using the CTEQ4HQ distributions) nicely illustrates how the charm distribution $`c(x,\mu ^2)`$ evolves as $`\mathrm{ln}(Q^2/m_c^2)`$ for increasing $`Q^2`$; note, $`c(x,\mu ^2)`$ enters with a negative sign so that the 4-flavor result is below the 3-flavor curve. For the scale choice, we take $`\mu =\sqrt{Q^2+m_c^2}`$. While the scale choice $`\mu =Q`$ is useful for instructive purposes such as demonstrating the matching of the 3- and 4-flavor calculations at $`\mu =Q=m_c`$, the choice $`\mu =\sqrt{Q^2+m_c^2}`$ is more practical as it provides a lower bound on $`\mu `$ which is important for the PDF’s and $`\alpha _s(\mu )`$. (Cf., Sec. \[3.3\], and Ref. .)

Additionally we note the stability of the 4-flavor scheme in contrast to the 3-flavor scheme. The shift of the curves when including the NLO contributions is quite minimal, particularly when compared with the 3-flavor result. This suggests that organizing the calculation to include the charm quark as a proton constituent can be advantageous even at relatively low values of the energy scale.

### 3.2 PDF Uncertainties: s(x), …

In Fig. 10, we show the variation of NLO calculation of $`\mathrm{\Delta }xF_3`$ on the strange-quark PDF. To obtain a realistic assessment of the $`s(x)`$ dependence, we have use the NLO calculation with PDF’s based on the MRST set which are re-fit with the value of $`\kappa =2s/(\overline{u}+\overline{d})`$ constrained to be $`\kappa =\{0.50,0.78,1.00\}`$. Note, by re-fitting the PDF’s with the chosen value of $`\kappa `$ we are assured to have an internally consistent set of PDF’s with appropriate matching between the quarks and gluon, and with the sum-rules satisfied.777Note, $`\kappa `$ is certainly $`Q`$-dependent, and the values for $`\kappa `$ quoted above correspond to the $`Q_0`$ of the evolution. While $`\kappa `$ compares the integral of $`s(x)`$ to the seq-quarks, there is also the possibility of an $`x`$-dependent variation. This has been studied in the fits of the strange-sea; we shall find that such subtle effects can play no role in resolving the $`\mathrm{\Delta }xF_3`$ issue.

The choice $`\kappa =0.50`$ is in line with the many experimental determinations of $`\kappa `$, cf., Table. 2; as expected, this prediction lies farthest from the data points.

The choice $`\kappa =0.78`$ is taken as an extreme upper limit given the experimental constraints; actually, in light of the results of Table. 2, this is arguably beyond present experimental bounds. This prediction is marginally consistent at the outer reach of the systematic + statistical error bars.

Finally, we take an $`SU(3)`$ symmetric set ($`\kappa =1`$) purely for illustrative purposes. It is interesting to note that even this extreme value is still below the central value of the data points at the higher $`x`$ values.

In conclusion we note that increasing the strange quark distribution does succeed in moving the theory toward the data; however, our consistent NLO analysis presented here suggests that the we have only limited freedom to increase $`s(x)`$, and that this alone is not sufficient to obtain good agreement between theory and data.

### 3.3 Scheme Choice

In our final section, we present the best theoretical predictions presently available to demonstrate the scheme dependence of $`\mathrm{\Delta }xF_3(x,Q)`$. Specifically, in Fig. 11 we show predictions for:

* NLO FFS GRV.
* NLO VFS TR calculation.
* NLO VFS ACOT calculation with CTEQ4 PDF’s.
* NLO VFS ACOT calculation with CTEQ5 PDF’s.

All these calculations use NLO matrix elements, and are matched with appropriate global PDF’s which are fitted in the proper scheme.

The first observation we make is how closely these four predictions match, especially given the wide variation displayed in previous plots such as Fig. 9. In hindsight, this result is simply a consequence of the fact that while different renormalization schemes can produce different results, this difference can only be higher order.888To be precise, different renormalization schemes can differ by i) terms of higher order in the perturbation series, and ii) terms of higher twist which do not factorize. Thus, the difference between these curves is indicative of terms of order $`\alpha _s^2`$ which have yet to be calculated.999For asymptotic results at order $`\alpha _s^2`$, see Ref. . When terms of order $`\alpha _s^N`$ are included, the span of these predictions will be systematically reduced to order $`\alpha _s^{N+1}`$.

In Fig. 11, we note the very close agreement among the VFS calculations, particularly the TR calculation and the ACOT calculation with CTEQ4 PDF’s. The ACOT calculation with the two CTEQ curves show primarily the effect of the charm distribution, as CTEQ4 uses $`m_c=1.6`$ and CTEQ5 uses $`m_c=1.3`$. The GRV calculation shows the effect of using yet a different scheme, in this case a FFN scheme, with its appropriately matched PDF. Were we to use MRST or CTEQ PDF’s, the spread of these theory curves would decrease; however, this would most likely represent an underestimate of the true theoretical uncertainties arising from both the hard cross section and PDF’s.101010The computation of PDF errors is a complex subject. For some recent approaches to this topic see: Refs..

While we consider it a triumph of QCD that different schemes truly yield comparable results (higher order terms aside), we should be cautious and note that the spread of these curves can only underestimate the true theoretical uncertainty. Note that GRV has a rather different strange distribution due to a different philosophy of obtaining this distribution rather than due to a different scheme.

## 4 Conclusions and Outlook

Comprehensive analysis of the neutrino data sets can provide incisive tests of the theoretical methods, particularly in the low $`Q^2`$ regime, and enable precise predictions that will facilitate new particle searches by constraining the PDF’s. This document serves as a progress report, and work on these topics will continue in the future.

Theoretical predictions for $`\mathrm{\Delta }xF_3`$ undershoot preliminary fixed target data at the $`1\sigma `$-level at low $`x`$ and $`Q`$. The neutrino structure function $`\mathrm{\Delta }xF_3`$ is obviously sensitive to the strange sea of the nucleon and the details of deep inelastic charm production. A closer inspection reveals, however, considerable dependence upon factors such as the charm mass, factorization scale, higher twists, contributions from longitudinal $`W^\pm `$ polarization states, nuclear shadowing, charge symmetry violation, and the PDF’s. This makes $`\mathrm{\Delta }xF_3(x,Q^2)`$ an excellent tool to probe both perturbative and non-perturbative QCD.

We have explored the variation of $`\mathrm{\Delta }xF_3(x,Q^2)`$ on the above factors and found none of these to be capable of resolving the discrepancy between the data and theory.

Although we have not eliminated the possibility that the entire set of parameters conspires to align the theory with the data, we have demonstrated this to be extremely unlikely. Of course, a definitive answer can only be obtained by a global analysis which combines the neutrino data for dimuons, $`\mathrm{\Delta }xF_3`$, $`F_2^\nu (PMI)`$, and $`F_2^{\mu ,e}`$.

As the situation stands now, this $`\mathrm{\Delta }xF_3(x,Q^2)`$ puzzle poses an important challenge to our understanding of QCD and the related nuclear processes in an important kinematic region. The resolution of this puzzle is important for future data analysis, and the solution is sure to be enlightening, and allow us to expand the applicable regime of the QCD theory.

## Acknowledgments

This work is supported by the Royal Society, the U.S. Department of Energy, the National Science Foundation, the Lightner-Sams Foundation, and the ‘Bundesministerium für Bildung, Wissenschaft, Forschung und Technologie’, Bonn.

We thank J. Bluemlein, A. Bodek, J. Conrad, J. Morfin, S. Kuhlmann, R.G. Roberts, H. Schellman, M. Shaevitz, J. Smith, and W.-K. Tung, for valuable discussions.