text:
# Magnetic-domain-controlled vortex pinning in a superconductor/ferromagnet bilayer

## Abstract

Vortex pinning in a type-II superconducting Pb film covering a Co/Pt multilayer with perpendicular magnetic anisotropy is investigated. Different stable magnetic domain patterns like band and bubble domains can be created in the Co/Pt multilayer, clearly influencing the vortex pinning in the superconducting Pb layer. Most effective pinning is observed for the bubble domain state. We demonstrate that the pinning properties of the superconductor/ferromagnet bilayer can be controlled by tuning the size, density and magnetization direction of the bubbles.

Magnetic flux penetrates type-II superconductors in the form of quantized vortices, which have the tendency to form a periodic lattice. These vortices move when a current is sent through the superconductor, thus causing dissipation and limiting the critical current density $`j_c`$ of the superconductor. To enhance $`j_c`$, vortex motion must be prevented. The latter can be achieved by introducing different pinning centers such as defects created by ion Civale et al. (1991) or neutron irradiation vanDover et al. (1989), lithographically introduced holes Fiory et al. (1978); Moshchalkov et al. (1999, 1998), or magnetic dots Otani et al. (1993); Martín et al. (1997); Morgan and Ketterson (1998); Van Bael et al. (1999).
Recently it was proposed that vortices in a superconductor/ferromagnet (SC/FM) multilayer are strongly pinned if the FM has a stripe domain structure, which typically exists in thin magnetic films with perpendicular anisotropy Bulaevskii et al. (2000). Enhanced pinning was observed in SC/FM bilayers compared to SC reference films García-Santiago et al. (2000); Zhang et al. (2001). However, in these studies the domain structure of the FM was unknown. In this Letter, we systematically investigate the influence of the domains in a ferromagnetic Co/Pt multilayer with perpendicular anisotropy on a type-II superconducting Pb film grown on top of the FM. Stable domain structures can be produced in a controlled way by magnetizing the sample before measuring the superconducting properties of the Pb film. We show that the strongest vortex pinning is obtained when localized domains (bubble domains) are present in the FM, and that by changing size and density of the bubbles one can control and optimize the vortex pinning.
The Co/Pt multilayer is deposited on a Si substrate with amorphous SiO2 top layer in an MBE apparatus by e-beam evaporation. The multilayer has a \[Co(0.4 nm)/Pt(1.0 nm)\]10 structure on a 2.8 nm Pt base layer. A 10 nm Ge film, a Pb film with thickness $`d_{Pb}=50`$ nm and a 30 nm Ge capping layer are subsequently evaporated on the Co/Pt multilayer (see Fig. 1) at a substrate temperature of 77 K.

The Pb film has a critical temperature of $`T_c=7.23`$ K. The penetration depth $`\lambda (0)=42`$ nm and the coherence length $`\xi (0)=41`$ nm are estimated from measurements of the temperature dependence of the upper critical field. The Ge film between Pb and Co/Pt is insulating at low temperatures, so that the proximity effects between Pb and Co/Pt are suppressed.
The Co/Pt multilayer has perpendicular magnetic anisotropy Zeper et al. (1989), which is confirmed by hysteresis loop measurements using the magneto-optical Kerr effect (MOKE). Fig. 2 shows the magnetization $`M_{fm}`$ of the Co/Pt multilayer, normalized to the saturation magnetization $`M_s`$, as a function of the magnetic field $`H`$ applied perpendicular to the surface.

The loop has an almost rectangular shape with $`H_n=0.6`$ kOe, $`H_c=0.93`$ kOe, and $`H_s=1.45`$ kOe, where $`H_n`$, $`H_c`$ and $`H_s`$ are the nucleation, coercive and saturation field, respectively.
The microscopic domain structure of the Co/Pt multilayer has been investigated by magnetic force microscopy (MFM) at room temperature and $`H=0`$. Figs. 3(a)-(c) show MFM images of the sample in different stable remanent magnetic states.

In Fig. 3(a), the sample has been demagnetized by oscillating $`H`$ around zero with decreasing amplitude perpendicular to the surface. The dark and bright regions in the image are domains with local magnetic moments $`m`$ either pointing up ($`m>0`$) or down ($`m<0`$), perpendicular to the sample surface. The observed band domain structure is typical for Co/Pt multilayers Valera et al. (1995). The width of the band domains is about $`0.20.4`$ $`\mu `$m.
The images shown in Figs. 3(b) and (c) are obtained after applying a negative field of $`10`$ kOe, sweeping $`H`$ to a positive value between $`H_n`$ and $`H_s`$ as indicated in Fig. 2, and then removing $`H`$. When $`H`$ is increased to a value slightly higher than $`H_n`$, small localized domains (bubble domains) with magnetic moments $`m>0`$ nucleate in the film. These bubbles are stable after removing $`H`$, as can be seen in Fig. 3(b). When $`H`$ is further increased, the bubbles grow and aggregate. At $`H=H_c`$, there is an equal number of positive and negative $`m`$. If $`H`$ is further increased, the number of negative $`m`$ decreases and the domains with $`m<0`$ become localized (see Fig. 3(c)). At $`H=H_s`$, the film is saturated with all $`m>0`$. We introduce the parameter $`s`$, which is defined by the percentage of magnetic moments that are pointing up ($`m>0`$), to describe the different remanent magnetic states obtained after this magnetization procedure. For the images shown in Figs. 3(b) and (c), $`s`$ is given by the dark area divided by the total area, which is in good agreement with $`s`$ calculated from Fig. 2.
The bubble domains in Fig. 3(b) typically have a size of about $`0.35`$ $`\mu `$m, whereas the average size of the bubbles in Fig. 3(c) is about $`0.25`$ $`\mu `$m. Another parameter changing during the reversal process is the number of bubbles per unit area $`\delta `$. For $`H`$ only slightly higher than $`H_n`$, $`\delta `$ is very small, but once the bubbles have nucleated in the FM, the reversal process is dominated by domain wall movement and $`\delta `$ does not change significantly during this process. From the MFM images, typical values of $`\delta =12`$ bubbles/$`\mu `$m2 are derived.
After MOKE and MFM measurements, the superconducting Pb film is deposited on the Co/Pt multilayer and the vortex pinning properties are studied by using a Quantum Design SQUID magnetometer. First the FM is brought into a specific magnetic state by applying the same magnetization procedure as for the MFM measurements at 300 K. The magnetization of the FM $`M_{fm}`$ is determined from $`M(T)`$ data as $`M_{fm}=M(T>T_c)`$, from which the parameter $`s`$ of each domain structure is calculated. Repeating the $`M(T)`$ measurements in several applied fields ($`|H|<120`$ Oe $`H_n`$, see Fig. 2) did not change the value of $`M_{fm}`$, implying that the domain structures are not altered by these small fields. The pinning strength related to each domain structure is determined at $`6.9`$ K by measuring $`M`$ as a function of $`H`$ applied perpendicular to the surface from $`120`$ Oe to +120 Oe. In order to compare the magnetization curves measured in the different magnetic states, $`M_{fm}`$ is subtracted from all curves shown in Fig. 3 (d)-(f), i.e. only the irreversible magnetization of the SC $`M_{sc}=MM_{fm}`$ is plotted. At $`6.9K`$, $`\xi (6.9\mathrm{K})=0.19`$ $`\mu `$m and the effective penetration depth $`\mathrm{\Lambda }(6.9\mathrm{K})=\lambda ^2(6.9\mathrm{K})/d_{Pb}=0.76`$ $`\mu `$m.
The upper $`M_{sc}(H)`$ curve shown in Fig. 3(d) is obtained after demagnetization (corresponding to the MFM image shown in Fig. 3(a)) and the lower one in the $`s=0`$ state, i.e. with the FM fully magnetized in the negative direction. For both magnetic states, a symmetric curve is found. The amplitude of $`M_{sc}(H)`$ (which is a measure for the magnitude of $`j_c(H)`$ Moshchalkov et al. (1998)) is lower for all $`H`$ when the FM is in the demagnetized state compared to the $`s=0`$ state. This unambiguously shows that the domain structure influences the SC, but does not enhance the vortex pinning in this sample. In the $`s=0`$ state, the stray field of the magnetic film has its largest amplitude at the sample boundary and is much smaller above the center of the FM. This means that the interior of the SC is only weakly influenced by the stray field. Contrary to that, the alternating stray field above the band domains suppresses the order parameter in the SC and may even lead to the creation of vortex-antivortex pairs connecting the domains with opposite magnetization and resulting in a reduction of $`j_c`$ Bulaevskii et al. (2000). This is also reflected by a decrease of $`T_c`$ of the SC to 7.04 K when the band domains are present in the FM.
When bubble domains are present in the FM, the $`M_{sc}(H)`$ curve becomes strongly asymmetric with respect to $`H`$. The direction of the asymmetry depends on the mutual orientation of $`H`$ and the magnetic moments $`m`$ of the bubbles. Fig. 3(e) shows two $`M_{sc}(H)`$ curves with $`m`$ of the bubbles pointing in the positive direction as in Fig. 3(b). Both curves have a larger amplitude of $`M_{sc}`$ for positive H. The upper curve in Fig. 3(f) is obtained after magnetizing the FM in $`H=20`$ kOe $`+H_c0`$, i.e. $`s=0.5`$. $`M_{sc}(H)`$ is almost completely suppressed for this domain state. The other two curves of Fig. 3(f) are obtained with $`m`$ of the bubbles in negative direction as in Fig. 3(c). Also for this configuration the $`M_{sc}(H)`$ curves are asymmetric, but reversed with respect to the field polarity compared to Fig. 3(e).
The asymmetry of $`M_{sc}(H)`$ can be explained by a magnetic interaction between vortices and bubble domains. The magnetic interaction can be caused by (i) the interaction between the magnetic field of the vortices $`B`$ and the magnetic moments $`m`$ in the FM (interaction energy $`E=mB`$, resulting in attractive interaction for $`m`$ and $`B`$ having the same polarity, and in repulsive interaction for $`m`$ and $`B`$ having opposite polarity), and (ii) the interaction between the vortices and the supercurrents generated in the SC by the stray field of the domain structure. During the measurement of $`M_{sc}(H)`$ shown in Fig. 3(e), the FM has had bubbles with $`m>0`$. Hence, evaluating the interaction term $`E=mB`$, the bubbles act as pinning sites for vortices when $`H>0`$, resulting in the enhanced amplitude of $`M_{sc}`$ for positive $`H`$. Consequently, the bubble domains with $`m<0`$ in Fig. 3(f) are pinning centers for vortices when $`H<0`$, giving a larger amplitude of $`M_{sc}`$ for negative $`H`$. This behavior is consistent with the vortex pinning properties of out-of-plane magnetized dots Morgan and Ketterson (1998). Assuming that one bubble can pin one vortex, we can estimate a ”matching” field $`H_{match}`$ at which all bubbles pin a vortex: $`H_{match}=\varphi _0\delta 2040`$ Oe, using $`\delta =12/\mu `$m2 and the flux quantum $`\varphi _0`$. In this field range, $`M_{sc}`$ indeed has a strongly enhanced amplitude for $`H>0`$ in the $`s=0.1`$ state and for $`H<0`$ in the $`s=0.85`$ state, indicating strong pinning of vortices.
At small fields $`|H|<3`$ Oe, $`j_c`$ is not increased by any of the investigated domain structures compared to the ”reference” $`s=0`$ state in which domains are absent. This can be attributed to the stray field above the domains, which results in an overall reduction of the superconducting order parameter $`\psi `$. Figs. 3(e) and (f) also show that more effective vortex pinning is achieved for bubbles with smaller size and lower density. The amplitude of $`M_{sc}`$ can be increased up to a factor of 3 when comparing the $`s=0`$ state with the lower curve of Fig. 3(f), and a factor of 2 compared to the middle curve of Fig. 3(f). This indicates that a reduced $`\psi `$ associated with bubbles with larger size and higher density counteracts efficient vortex pinning, even leading to the extremely weak amplitude of $`M_{sc}(H)`$ for the $`s=0.5`$ state shown in the upper curve of Fig. 3(f).
In conclusion, we have shown that the domain structure in a magnetic film can significantly influence the vortex pinning in a type-II superconductor. By tuning the density of bubble domains in the magnetic film, one can control the vortex pinning. This could be used to design logical devices based on the control of superconductivity by magnetic domain structures Bulaevskii et al. (2000).
The authors acknowledge L. Van Look, K. Temst and G. Güntherodt for help with sample preparation, and J. Swerts for MOKE measurements. This work was supported by the Belgian IUAP and the Flemish GOA programs, by the ESF ”VORTEX” program and by the Fund for Scientific Research (F.W.O.) - Flanders. MJVB is a Postdoctoral Research Fellow of the F.W.O. - Flanders.


format: text

text:
# Lagrangian evolution of global strings

\[

## Abstract

We establish a method to trace the Lagrangian evolution of extended objects consisting of a multicomponent scalar field in terms of a numerical calculation of field equations in three dimensional Eulerian meshes. We apply our method to the cosmological evolution of global strings and evaluate the energy density, peculiar velocity, Lorentz factor, formation rate of loops, and emission rate of Nambu-Goldstone (NG) bosons. We confirm the scaling behavior with a number of long strings per horizon volume smaller than the case of local strings by a factor of $``$ 10. The strategy and the method established here are applicable to a variety of fields in physics.

\]

Extended objects consisting of scalar fields such as topological defects or solitons play important roles in many fields of physical science ranging from condensed matter physics to cosmology. Much work has been done on the cosmological formation of topological defects by thermal  or nonthermal  phase transitions, their subsequent evolution , and cosmological implications of Q-balls . On the other hand, it has been pointed out that cosmological defect formation may be realizable in a laboratory . The phase transition of liquid helium can also be described in terms of Ginzburg-Landau theory as in the case of cosmological phase transitions. The order parameter, which is an expectation value of a scalar field in cosmological phase transition, is just replaced by a Bose condensate wave function. Following Zurek’s suggestion, the formation of topological defects has been studied in 3He  and 4He . In fact, vortex formation is observed in 3He. Thus, it is very important to investigate the formation and subsequent evolution of extended objects, such as topological defects, produced both in cosmology and laboratories.

Among various types of topological defects predicted in high energy physics, strings hold a unique position, which, unlike magnetic monopoles and domain walls, do not overclose the universe because they settle down to the scaling solution . The key mechanism to achieve it is intercommutation of infinite strings to lose their energy by producing closed loops, which decay through radiating gravitational waves for local strings***The possibility has also been pointed out that infinite strings lose their energy by directly emitting massive particles for local strings . or Nambu-Goldstone (NG) bosons for global strings.

While many grand unified theories may contain local strings, there are no particle physics motivations for their existence. Nonetheless they were extensively studied because it was expected that they might be the origin of density fluctuations which seeded cosmic large-scale structures. In particular, time evolution of local cosmic strings have been investigated by a number of authors  who used the Nambu-Goto action  that applies for infinitely thin strings. These numerical analyses confirmed the scaling behavior and the scaling parameter, $`\xi `$, which is defined as $`\xi =\rho _st^2/\mu `$ with $`\rho _s`$ energy density of strings and $`\mu `$ string tension per unit length, converged to $`\xi 10`$ in the radiation dominated universe . Recent observations of anisotropies of cosmic microwave background radiation, however, disfavor cosmic-string scenario of structure formation , and motivations to study local strings have somewhat diminished.

Global strings, on the contrary, are well motivated in particle physics, whose existence is predicted by the Peccei-Quinn solution to the strong CP problem . They are produced when the Peccei-Quinn U(1) symmetry is broken . They radiate axions as associated NG bosons  which are one of the best candidates of cold dark matter. Furthermore, they are also realized in condensed matter physics. Vortices of 3He and 4He studied in a laboratory are analogous to global strings .

In contrast to local strings, energy density of global strings is dominated by the gradient energy of the NG scalar field rather than the potential energy of the core, because there are no gauge fields to cancel the former. As a result they have a long-range interaction and their dynamics cannot be analyzed with the Nambu-Goto action. It should be described with the Kalb-Ramond action instead , which incorporates NG bosons and their coupling with the string core in addition to the Nambu-Goto action. In this action, however, the line density of the string is not well defined which has a logarithmic dependence on the radius. It is also very difficult to analyze it numerically. Hence the dynamics of strings governed by the Kalb-Ramond action has not been clarified in cosmological context yet even though evolution of axionic strings has very important implications in cosmology and particle physics. Instead, the results of numerical simulations of local strings have been used to discuss cosmology of global strings as well, which may not be justified because evolution of global defects can be significantly different from that of local counterparts as seen in the case of local  and global  monopoles.

In order to surmount such a situation, we have investigated evolution of global strings by solving equation of motion of the complex scalar field numerically in a three dimensional lattice in the expanding universe instead of dealing with the Kalb-Ramond action in collaboration with Kawasaki . In these fully Eulerian simulations, it is difficult to identify strings and follow their dynamics. In the previous publications, whether a lattice contains a string core was judged by comparing potential energy density at each vertex with a static cylindrically symmetric solution of a global string. Counting number of lattices thus identified as containing a string, we estimated the scaling parameter $`\xi `$ and found that $`\xi 1`$, which is significantly smaller than the case of local strings. This result has been criticized by Martins and Shellard  who claim that we have obtained incorrectly smaller value of the scaling parameter because the energy loss of long strings due to direct emission of NG bosons is much larger in simulation with small dynamic range than in the actual cosmological context, and that both types of long strings should behave quantitatively in the same way on cosmological scales.

In order to characterize cosmological evolution of global strings completely, we need to establish the method to follow cosmological evolution of global strings from a Lagrangian point of view. The energy density of long strings, $`\rho _{\mathrm{}}`$ satisfies the following equation:

$$\frac{d\rho _{\mathrm{}}}{dt}=2H(1+v^2)\rho _{\mathrm{}}\mathrm{\Gamma }_{\mathrm{loop}}\rho _{\mathrm{}}\mathrm{\Gamma }_{\mathrm{NG}}\rho _{\mathrm{}},$$
(1)
where the second and the third terms on the right-hand-side represent energy loss due to loop formation and direct emission of NG bosons or axions, respectively, and $`v^2`$ denotes average square velocity of string segments. In the scaling regime, the string network is characterized by a scale $`L`$, which is defined as $`L\sqrt{\mu /\rho _{\mathrm{}}}`$ and grows with the horizon scale $`Lt`$. Here $`\mu `$ is the effective line energy density of a global string. We can introduce a loop production coefficient $`c`$ and an emission coefficient $`\kappa `$ by

$`\mathrm{\Gamma }_{\mathrm{loop}}\rho _{\mathrm{}}=c{\displaystyle \frac{\rho _{\mathrm{}}}{L}},\mathrm{\Gamma }_{\mathrm{NG}}\rho _{\mathrm{}}=\kappa {\displaystyle \frac{\rho _{\mathrm{}}}{L}}.`$ (2)
Then, the scaling parameter $`\xi `$ is characterized by these coefficients as

$$\xi =\left(\frac{1v^2}{c+\kappa }\right)^2.$$
(3)
Indeed if $`\kappa `$ incorrectly turned out to be much larger than $`c`$, we would find a smaller value of $`\xi `$ than it should really be. Hence it is essential to evaluate contribution of each term in (1) with the help of numerical analysis. For this purpose we must know how each string segment moves and intercommutes with each other to trace formation of loops, which is impossible in conventional Eulerian simulations. Notice that this problem is ubiquitous because in many fields of physics we must solve equation of motion in Eulerian meshes, for example, in order to trace time evolution of solitons like Q-balls .

In the present Rapid Communication, we establish a method to follow Lagrangian evolution of extended objects in the context of Eulerian calculation of scalar field equations on the lattice. As a specific example, we concentrate on global strings here, but the strategy and the method established here directly apply to other situations in physics with some proper modifications.

We express the complex scalar field, which constitute global strings, in terms of two real scalar fields $`\varphi _a(x)`$ ($`a=1,2`$) and define their Lagrangian density as

$$[\varphi _a]=\frac{1}{2}g_{\mu \nu }(^\mu \varphi _a)(^\nu \varphi _a)V[\varphi _a,T],$$
(4)
where $`g_{\mu \nu }`$ is the flat Robertson-Walker metric. The finite-temperature effective potential $`V[\varphi _a,T]`$ is taken as

$`V[\varphi _a,T]`$ $`=`$ $`{\displaystyle \frac{\lambda }{4}}(\varphi ^2\sigma ^2)^2+{\displaystyle \frac{\lambda }{6}}T^2\varphi ^2,\varphi ^2\varphi _1^2+\varphi _2^2`$ (5)
which exhibits a second-order phase transition with the critical temperature $`T_c=\sqrt{3}\sigma `$ which produces global strings by breaking a global U(1) symmetry.

Equations of motion for the scalar fields are given by

$$\ddot{\varphi _a}(x)+3H\dot{\varphi _a}(x)\frac{1}{R(t)^2}^2\varphi _a(x)=\frac{V}{\varphi _a},$$
(6)
where a dot denotes time derivative and $`R(t)`$ is the cosmic scale factor. In the radiation dominated universe, the Hubble parameter $`H=\dot{R}(t)/R(t)`$ and cosmic time $`t`$ are given by

$`H^2={\displaystyle \frac{8\pi }{3M_{\mathrm{pl}}^2}}{\displaystyle \frac{\pi ^2}{30}}g_{}T^4,t={\displaystyle \frac{1}{2H}},`$ (7)
respectively, where $`M_{\mathrm{pl}}=1.2\times 10^{19}`$GeV is the Planck mass and $`g_{}`$ is the total number of degrees of freedom for the relativistic particles. For the ease of numerical calculations we take $`\sigma =0.1(45/16\pi ^3g_{})^{1/2}M_{\mathrm{pl}}`$ and $`\lambda =0.08`$ but the results are insensitive to these choices. We start numerical simulation at the temperature $`T_i=2T_c`$ corresponding to $`t_i=t_c/4`$ and adopt as an initial condition the thermal equilibrium state with a mass squared equal to the inverse curvature of the potential at that time. We have simulated the system from ten different thermal initial conditions under the periodic boundary condition. The number of lattices is $`128^3`$ with the lattice spacing $`\delta x=\sqrt{3tt_i}/8`$ and the time step $`\delta t=0.01t_i`$. Thus, the box size is nearly equal to horizon volume $`(H^1)^3`$ and the lattice spacing to a core size of a string $`\delta =1/(\sqrt{2\lambda }\sigma )`$ at the final time $`t_f=200t_i`$.

The above setup of calculation is the same as our previous fully Eulerian simulations . In order to reproduce Lagrangian evolution of string segments in the Eulerian meshes, we develop a different method to identify string cores from scalar field configurations, because the previous method  was inadequate for configurations with large curvature such as small loops and it was impossible to find more correct position of string core in a box beyond the lattice spacing. Here we use a new two-step algorithm to identify the locus of a string more accurately. The first step is to find plaquettes which a string penetrates. This is done in terms of the Vachaspati-Vilenkin algorithm  by monitoring phase rotation around each square. Then using the value of $`\varphi _a`$ at each vertex of a plaquette penetrated by a string, we linearly interpolate $`\varphi _a(x)`$ to calculate the position where both $`\varphi _1(x)`$ and $`\varphi _2(x)`$ vanish. We can thus find where a string penetrate in each plaquette and find more accurate trajectory of a string by connecting these points.

The next and the more important task to obtain Lagrangian evolution of string network, that has been missing in the previous Eulerian analyses, is to find velocity of each string segment. Since motion tangential to a string is a gauge mode, we should evaluate velocity normal to it. Suppose that a string core exists at a point $`𝒙_0`$ at time $`t_0`$, namely, $`\varphi _1(𝒙_0,t_0)=\varphi _2(𝒙_0,t_0)=0`$. In order to estimate where this string segment moves $`\mathrm{\Delta }t`$ later, we expand scalar fields $`\varphi _a(𝒙,t_0+\mathrm{\Delta }t)`$ around $`\varphi _a(𝒙_0,t_0)`$ up to the first order,

$`\varphi _a(𝒙,t_0+\mathrm{\Delta }t)`$ $``$ $`\varphi _a(𝒙_0,t_0)+\varphi _a(𝒙_0,t_0)(𝒙𝒙_0)`$ (9)
$`+\dot{\varphi }_a(𝒙_0,t_0)\mathrm{\Delta }t(a=1,2).`$
Again the loci of string core at $`t=t_0+\mathrm{\Delta }t`$ are given by the line where both scalar fields $`\varphi _a(𝒙,t_0+\mathrm{\Delta }t)`$ vanish. That is, under the first-order approximation (9) it lies on the intersection of two planes given by

$$\varphi _a(𝒙_0,t_0)(𝒙𝒙_0)+\dot{\varphi }_a(𝒙_0,t_0)\mathrm{\Delta }t=0,$$
(10)
with $`a=1,2`$. Now suppose that the line normal to the string segment at $`(𝒙,t_0)`$ intersect with the trajectory of string at $`t=t_0+\mathrm{\Delta }t`$ at $`𝒙=𝒙_l(𝒙_0,t_0,\mathrm{\Delta }t)`$. Then, making use of the fact that $`𝒙_l`$ satisfies Eq. (10), we can easily calculate the velocity of this string segment, and its magnitude is given by

$`𝒗(𝒙_0,t_0)`$ $`=`$ $`\underset{\mathrm{\Delta }t0}{lim}{\displaystyle \frac{𝒙_l(𝒙_0,t_0,\mathrm{\Delta }t)𝒙_0}{\mathrm{\Delta }t}}`$ (11)
$`=`$ $`{\displaystyle \frac{\dot{\varphi }_1\varphi _2\dot{\varphi }_2\varphi _1}{|\varphi _1\times \varphi _2|}}.`$ (12)
One should note that the formula obtained here applies to all soliton-like objects consisting of a field configuration with some proper modifications.

Now that we have developed a method to identify the location and velocity of strings at each time, the Lagrangian evolution of the string network can be traced just as in the case of numerical simulation of the evolution of local strings based on the Nambu-Goto action. In fact, our simulation, being based on scalar field equations, contains even more information, that is, we can find fate of intersecting string segments, namely, whether they reconnect or simply pass through each other, without assigning the probability of reconnection by hand unlike in the case of simulations based on the Nambu-Goto action. We can therefore calculate formation rate and spectrum of string loops without ambiguities.

Thus we can calculate all the terms of eq. (1) directly from simulation data except for the last term, which can be evaluated by use of this equation itself. As pointed out in , we should use $`\mu =\gamma \mu _s`$ as the energy per length to calculate the energy density of long strings, $`\rho _{\mathrm{}}`$. Here $`\gamma `$ is the average Lorentz factor and $`\mu _s`$ is the line density of a static string given by $`\mu _s2\pi \sigma ^2\mathrm{ln}(t/(\delta \xi ^{1/2}))`$.

Figure 1 depicts time evolution of $`\xi `$ with three different identification methods of strings. Filled squares represent results of our new identification scheme, while blank circles and blank squares correspond to our previous method  and Vachaspati-Vilenkin algorithm, which was adopted in , respectively. Both methods estimate total string length by simply counting the number of boxes penetrated by a string. As is easily seen, our previous method slightly overestimates $`\xi `$ and Vachaspati-Vilenkin approach overestimates $`\xi `$ by a factor of $`1.4`$. We confirm that $`\xi `$ becomes constant after some relaxation period and find that the more correct value of the scaling parameter in the radiation dominated era is given by $`\xi 0.80`$.

Figures 2 depict time evolution of average velocity, average velocity squared, and average Lorentz factor. The first two quantities are also found to relax to constant values after some relaxation time, namely, $`v0.65`$ and $`v^20.50v^2`$. The average Lorentz factor, however, has a large scatter in time, although long-time average is fairly constant with $`\overline{\gamma }1.8`$. Since string segments moving with a speed close to light velocity have extremely large Lorentz factor and push up the average dramatically, fluctuation in the number of such string segments results in such large scatter. This also explains why the average Lorentz factor $`\gamma =1.62.0`$ is larger than $`1/\sqrt{1v^2}`$ and $`1/\sqrt{1v^2}`$. Thus the energy per unit length of a string is enhanced by a factor $`\overline{\gamma }1.8`$.

Calculating the length of all infinite strings and that of loops at each time step and comparing with those at the preceding time step, the loop production parameter is found to be $`c=0.430.53`$. Then taking account of the relation (3), the emission parameter $`\kappa `$ is calculated as $`\kappa =0.030.13`$. We therefore conclude that the direct emission of NG bosons is a subdominant channel of energy loss of long strings even in the case cosmic horizon scale is not much larger than the string width due to the limitation of the dynamic range of numerical simulations. Thus, even if $`\kappa `$ might become significantly smaller in the cosmological situation because it is proportional to $`\mathrm{ln}(t/\delta )`$ as pointed out in , the scaling parameter $`\xi `$ would not increase to coincide with that of local strings. The result of our simulation with $`\mathrm{ln}(t/\delta )5`$ implies that $`\kappa 0.01`$ in the cosmological situations. Inserting this new value of $`\kappa `$ into the relation (3) together with $`c`$, $`\kappa `$, and $`v^2`$ yields $`\xi =0.91.4`$ in the cosmological situations, which is still significantly smaller than the case of local strings. Hence, it is confirmed that although global strings relax to the scaling solution just as local strings, their number density is much smaller than the local ones. This is mainly because global strings intercommute more often and $`v^2`$ is larger due to long-range attractive forces between strings. Finally, we apply our results into the constraint of the breaking scale $`f_a`$ of the Peccei-Quinn U(1) symmetry. Then, in accordance with Ref. , $`f_a`$ is constrained as $`f_a(0.201.6)\times 10^{12}`$ GeV for the normalized Hubble parameter $`h=0.7`$ with $`\xi =0.91.4`$ and $`\gamma 1.62.0`$.

In summary, we have developed a new method to trace Lagrangian evolution of topological defects solving scalar field equation on three dimensional Eulerian lattices and given all the quantities characterizing cosmological evolution of the global string network, namely, energy density, peculiar velocity, Lorentz factor, intercommutation rate, and emission rate of NG bosons. All these quantities can be obtained in Eulerian simulations only after establishing the method given in this Rapid Communication, which enables us to extract Lagrangian information from simulations done in Eulerian meshes. Our method is directly applicable to situations in a variety of fields in physics.

M.Y. and J.Y. are partially supported by the Japanese Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science, and Technology, Nos. 12-08555(MY) and 13640285(JY).


text:
# Switching current noise and relaxation of ferroelectric domains

## I Introduction

Polarization reversal and properties of the accompanying switching current are relevant for various applications of ferroelectric materials. In particular, it is a key process that takes place in electronic devices based on ferroelectrics thin films, where a reproducible switching process is required. Apart from the extensive experimental and theoretical work , a complete understanding of the domain-wall motion and predictability of the switching processes remains elusive in real ferroelectrics, especially when several physical parameters, determining either structure or dynamics, are varied. Domain switching and relaxation (back-switching) in ferroelectrics can be studied by hysteresis loop, switching current, and switching current bursts (electrical Barkhausen noise).

The development of domain structure and motion of domain walls during switching can be regarded in the scope of nonlinear dynamic processes occurring in driven disordered systems, which are studied extensively in the past decade. Particularly, the measurements on ferromagnets (a short review is given in ) show that Barkhausen noise exhibits interesting fractal properties, that can be characterized by a full set of scaling exponents of the sizes, durations and energies of the bursts and of the power spectrum. An example with a complete set of measured scaling exponents and corresponding fractal dimensions with derivation of the associated scaling relations can be found in Ref. . A possible grouping of the scaling exponents into different universality classes of scaling behavior offers a new and simple link between the measured Barkhausen noise and the domain structure. For example, in a recent study  such a relation was established between the universal scaling exponents of acoustic emission noise and the packing properties near the first-order structural transformation.

The electrical Barkhausen noise has been observed experimentally in $`\mathrm{BaTiO}_3`$  and $`\mathrm{Gd}_2(\mathrm{MoO}_4)_3`$ and $`\mathrm{LiTaO}_3`$  single crystals and in $`9.4/65/35(\mathrm{Pb},\mathrm{La})(\mathrm{Zr},\mathrm{Ti})\mathrm{O}_3`$ (PLZT) ceramics  and $`\mathrm{PbMg}_{1/3}\mathrm{Nb}_{2/3}\mathrm{O}_3`$ (PMN) relaxor ferroelectric . Physically, the polarization reversal and domain wall motions in ferroelectrics is accompanied with the local change of sign of bound charges, which influences values of the electric field in the sample. The depolarization field of bound charges can be partially screened by an adequate addition of charges on the surface and by other methods (see  and references therein). Hence at a given point $`x`$ in the sample at time $`t`$ un-compensated depolarization field is given by

$$\delta E_{dep}(x,t)=E_{dep}(x,t)E_{sreen}(x,t),$$
(1)
where $`E_{sreen}(x,t)`$ is a total field due to screening (surface and bulk) charges at that site. Incomplete screening of the depolarization field tends to change the electric field inside the sample, thus influencing the course of switching process itself. For instance, there is a finite probability that inverse jumps of polarization (negative Barkhausen pulses)  occur, that slow down the switching. Recently a series of sophisticated measurements were done  that reveal the effects of the depolarization field on a created plane domain wall in $`\mathrm{Gd}_2(\mathrm{MoO}_4)_3`$ and on a multi-domain structure in $`\mathrm{LiTaO}_3`$. The obtained distributions of lengths, sizes and rest time of Barkhausen jumps obey a power-law behavior with an exponential cut off, both for switching and back-switching processes. Using a suitable experimental set-up the power-spectrum in PMN relaxors was measured, which was found to decay with frequency $`\omega `$ as $`S(\omega )\omega ^2`$ .

Recently studied uniaxial relaxor ferroelectrics $`\mathrm{Sr}_{0.61}\mathrm{Ba}_{0.39}\mathrm{Nb}_2\mathrm{O}_6`$ (SBN61) and $`\mathrm{Pb}_{0.61}\mathrm{Ba}_{0.39}\mathrm{Nb}_2\mathrm{O}_6`$ are regarded as systems with the Ising symmetry of the order parameter, having the tetragonal structure in the low-temperature phase with the spontaneous polarization along c-axis . In solid solutions $`\mathrm{SBN61}:\mathrm{X}`$ doping with $`\mathrm{X}=\mathrm{Ce}^{3+},\mathrm{Cr}^{3+},\mathrm{Co}^{2+}`$ reduces the transition temperature and induces local electric random fields, which give rise to local polar clusters . Recently a comprehensive series of measurements in SBN:Ce and similar systems  gave a strong support to the presence of random electric fields in this class of relaxor ferroelectrics. In general, more common sources of uncorrelated random fields in ferroelectrics are neutral defects that can locally break the polar symmetry, as for instance in TGS (see detailed discussion in Ref. , page 149). \[Note that, in contrast to ferroelectrics, in ferromagnets local random fields are not allowed by the symmetry, but are theoretically constructed as being induced by coarse graining due to other types of disorder (see also ).\]

In this work we simulate temporal behavior of polarization driven by the electric field in three-dimensional systems of Ising spins subject to quenched local random fields and positive nearst-neighbor interaction, which is compatible with the occurrence of a ferroelectric long-range order. We consider two types of driving. First, by stepwise increase (decrease) of the external field along a hysteresis loop (see below) we study bursts in the switching current $`j(t)=dP(t)/dt`$ (electrical Barkhausen noise). We discuss effects of incomplete screening of depolarization fields on the coarse of switching and hysteresis loop properties. When the screening is complete, we study how the statistical properties of the noise depend on the strength of pinning and on the frequency with which the driving field is changed. Second, by imposing a high opposite value of the external field in the homogeneous state we let the polarization to reach its saturation at that field, and then remove the field. We determine reversal and spontaneous back-switching of the polarization with time at different strength of disorder. This type of switching is more familiar in the experiments on ferroelectrics .

Apart from allowing negative Barkhausen pulses as described below, the present study is complementary compared to earlier simulations of Barkhausen pulses in 3-dimensional Ising systems in the presence of disorder  in the following sense: In the earlier work  we considered only strong disorder that is compatible with nucleation of small domains under driving field, as for instance in relaxor ceramics, and analyzed the signal from only central part of the hysteresis loop, as it is usually done on experiments. In the present work the study is extended to entire range of disorder including the weak pinning region where extended domain walls may occur due to lack of pinning centers . In this case a potential depinning of domain walls at a critical driving field contributes significantly to the scaling behavior of Barkhausen pulses, as recent theoretical results on Bethe lattice suggest . Therefore, we include integration over entire hysteresis loop in order to incorporate the contribution of critical depinning to the universal scaling properties of the induced electrical noise in the case of weak pinning . Furthermore, here we employ the simulation algorithm at finite driving rates, that mimics closely the experimental situation with stepwise increase (decrease) of the electric field. We examine in detail the effects of driving rates in both strong and weak pinning case by varying height of the steps in the external field. The earlier simulations in Refs.  were done in the limit of zero driving rate, which appears physically inaccessible to real experiments. Whereas an entirely different aspect concentrating on the dynamic scaling during the growth of individual avalanches was considered in . It should be stressed that in this paper we implement strictly quenched disorder, as described in Section II. Note that some speed algorithms, as for instance used in , allow certain dynamic variations of pinning while the avalanche grows, which can be relevant to systems with annealed rather than quenched disorder. Theoretically, having unfixed pinning centers during the avalanche growth may alter the long-range correlations in the avalanche-triggering fields and in distances between initial points of avalanches, which were shown to be closely associated with the appearance of the fractal noise in two-dimensional systems .

The paper is organized as follows. In Sec. II we introduce the model and discuss all relevant parameters that influence the hysteresis loop properties. In Sec. III we study in detail scaling behavior of the switching current noise under varying disorder and driving conditions. Sec. IV contains the results of switching in a high field and spontaneous back switching of polarization when the pinning is varied. The paper contains a short summary of the results and discussion in Sec. V.

## II Model and relevant parameters

We consider a model with Ising-type spins situated on three-dimensional simple cubic lattice with positive interactions between neighbor spins and local random field at each lattice site :

$$=\underset{x,y}{}J_{xy}S(x)S(y)\underset{x}{}h_xS(x)E\underset{x}{}S(x).$$
(2)
We assume that the random exchange interactions $`J_{x,y}`$ are distributed around a positive value $`J`$ with a narrow distribution. The random fields $`h_x`$ are given by Gaussian distribution with zero mean and the width $`f`$ (measured in units of $`J^2`$). As mentioned above, the source of local random fields can be neutral defects that break polar symmetry  or charge disorder due to doping , which cause occurrence of polar regions in relaxor ferroelectrics. In addition, a random part of the interaction between these polar regions contributes to the observed glassy behavior of the dielectric response . Since the random fields locally break the symmetry of the order parameter, we believe that the dominant disorder effects on the dynamics that we consider are due to random fields. Therefore, here we neglect the randomness in the spin-spin interaction.

The spin dynamics that we are considering is field-assisted and takes part deeply in the ordered phase. Hence, we neglect possible temperature effects, both for the reasons of keeping the number of parameters finite and assuming that they play a secondary role in these processes. The dynamics consists of spin alignment along a locally dominant field. Thus, when the local field at site $`x`$ and time $`t`$, $`h_{loc}(x,t)`$ exceeds zero, the spin $`S(x,t)`$ at that site flips and thus aligns along the local field. The local field $`h_{loc}(x,t)`$ consists of the interaction and pinning part $`h_{ip}(x,t)=_yJ_{xy}S(y,t)+h(x)`$ and the driving field $`E(x,t)`$. The driving field itself is given by $`E(x,t)=E_{ext}(t)\delta E_{dep}(x,t)`$, where $`E_{ext}(t)`$ is the external field and $`\delta E_{dep}(x,t)`$ a non-compensated part of the depolarization field at that site . We employ a slow stepwise increase of the external field $`E_{ext}(t_i)=E_{sat}+t_i\mathrm{\Delta }E`$ at discrete time intervals $`\{t_i\}`$, while assuming that a non-compensated part of the depolarization field at site $`x`$ may assist a probabilistic (probability parameter $`b`$) back-flip of an already aligned spin at that site. It should be stressed that the parameter $`b`$ represents the probability that at a site $`x`$ in time $`t`$ the depolarization field opposes the effect of external filed. (The situation where the depolarization field is parallel to the external field is less interesting within the present numeric implementation of the dynamics where driving rate is not fixed, as described below.) Having in mind that the depolarization field at a given site depends on the actual domain structure and its kinetics, the uniform distribution of the stochastic variable $`b`$ can be regarded only as a first step towards modeling of this complex phenomena.

Implementation of the dynamics is as follows: Quenched random fields of Gaussian distribution (in double precision) are generated and stored. They are kept constant until entire hysteresis loop is completed by the driving field. The external field is increased by a small amount $`\mathrm{\Delta }E`$ and the set of local fields $`\{h^{loc}(x,t)\}`$ computed at each lattice site and stored. The system is updated in parallel, i.e., with respect to stored set of local fields, and the spin flipped when $`h^{ip}(x,t)+E_{ext}(x,t)`$ exceeds zero by a small amount $`10^{10}`$. When it is assumed that the depolarization field is not entirely compensated, the spin flip is executed with a reduced probability $`1b`$. This completes one time step of avalanche evolution. Then a new set of local fields at affected sites (neighbors of just flipped spins) is computed and stored. The conditions for spin flips at those sites examined and spins flipped according to the above rule. The process continues until no more spins are found which satisfy the conditions for flip. This completes one Barkhausen avalanche. Then the external field is increased again by the amount $`\mathrm{\Delta }E`$. When the whole loop is completed, we set another distribution of quenched random fields and repeat the driving. The results are averaged over the entire ensemble (up to 200) of different random-field realizations.

In this model the following relevant parameters appear: Width of the random-field distribution $`f`$, that determines strength of pinning and/or domain size; jump $`\mathrm{\Delta }E`$ of the external field, which defines the driving rate; probability of back-flips $`b`$, which is related to the unscreened part of the depolarization field. In practice, usually an additional parameter is due to the section of the hysteresis loop at which the signal is monitored, since the properties of the signal vary along the loop. As mentioned above, here we integrate the statistics over the entire hysteresis loop in order to capture the contribution of the depinning of large domain walls to the scaling behavior of noise by weak pinning.

The hysteresis loop dependence on pinning by random fields and on driving rates is shown in Fig. 1. By varying the strength of pinning $`f`$ the form of the hysteresis loop changes (Fig. 1a) from the rectangular loop at low disordered 3-d sample to the slim loop, when the pinning is strong (large $`f`$). In the case of low disorder, i.e., at weak pinning of domain walls, large domains may form due to lack of pinning centers. In real ferroelectrics occurrence of large domains are usually accompanied by extended plane domain walls . On the other side, when the pinning is strong, such as in the case of pinning at grain boundaries of small grain ceramics , many small domains can be nucleated at spatially distributed centers with weak pinning, and their walls can propagate along few active sites with weakest pinning, as the driving field increases. This makes effectively low-dimensional walls moving in a three-dimensional sample. In addition, a massive nucleation of domains in the central part of hysteresis imposes a spatial restriction on the growth of later domains. Hence we expect a different scaling behavior of sizes and durations of induced switched domains in the case of weak and strong disorder. Our simulations confirm this phenomenological picture. In ferroelectrics a continuous variations of the size of domains can be achieved by increasing temperature or by applying stress, as for instance described in 8/65/35 PLZT ceramics . It was suggested in Ref.  that the appropriate illumination of samples can notably reduce the strength of pinning in doped relaxor $`\mathrm{SBN61}:\mathrm{Ce}^{3+}`$ ferroelectrics.

We would like to stress that in the algorithm keeping the external field constant during the avalanche evolution and increasing it only after the avalanche has stopped does not allow to study rest periods between avalanches. However, this type of field update leads to a better resolution of the noise, while it changes the time dependence of the driving field (cf. Fig. 1c). The driving rate can be better defined as the average derivative $`r<dE(t)/dt>`$ along the loop. In the simulations we find that when back flips occur with a finite probability $`b>0`$ they effectively reduce increase of the driving field and prolong the duration of a jump, as well as inducing inverse jumps, as shown in Fig. 1 c,d.

Assuming that the back flips occur very rarely ($`b0`$), we obtain the switching current noise which is shown in Fig. 2. The train of bursts of the switching current exhibits a long-range correlations—Fourier spectrum $`j(\omega )`$ of the signal is shown in Fig. 2 (top). It is related to the power-spectrum via $`S(\omega )=|j(\omega )|^2`$, that can be directly measured (see ). Hence, the exponent of the power spectrum in the present model is given by $`2\times \varphi =1.48\pm 0.02`$ for the case of moderate pinning $`f=3.6`$. The power-law behavior of the power spectrum suggests that the switching current noise exhibits fractal character when certain conditions are met, in particular, when the depolarization fields are successfully screened.

## III Fractal properties of the switching-current noise

The fractal noise of the switching current have a range of universal scaling properties, that can be determined more clearly if the field inside the sample is fully controlled. Therefore, in this section we put $`b=0`$ in order to study in detail these scaling properties and to determine how they depend on two relevant parameters—the domain size and the driving rate.

### A Varying strength of pinning

First we fix the driving rate to a small value $`\mathrm{\Delta }E/J=1\times 10^3`$ and consider how the domain size or strength of pinning influence the properties of the switching current. In Fig. 3 we show the cumulative distributions of size (area) $`s`$ of switched domains and duration $`t`$ of switching following a single field update, for several values of pinning strength $`f`$. As the Fig. 3 shows, both switched areas and durations of their switching show a wide range of values that are described by power-law distributions with finite cut-offs ($`X`$ stands for $`Xs,t`$):

$$P(X,L)=X^{(\tau _X1)}𝒫(XL^{D_X}).$$
(3)
The scaling exponents $`\tau _s1`$ and $`\tau _t1`$ determine the respective slopes of the cumulative distributions of switched area and duration of switching (see Fig. 3). The fractal dimension $`D_s`$ of size of switched domains and the dynamic exponent $`D_tz`$ in Eq. (3) determine how the corresponding cut-offs scale with the characteristic length (here the system size $`L`$), $`X_0L^{D_X}`$.

Both the slopes $`\tau _X`$ and the cut-offs, related to $`D_X`$, of the distributions depend on the strength of pinning. First, for weak pinning (small values of $`f`$, bottom curve in Fig. 3) we notice that extended domains, bounded by presumably plane domain walls, may occur. By increasing the field these domain walls move through a weakly random medium and eventually depin at a critical field $`H_d(f)`$. Depinning of large domain walls leads to a peak in the histogram at large values of $`s`$, or a plateau in the cumulative distribution as shown in Fig. 3. Before the depinning is reached a series of small avalanches occurs that are power-law distributed. The plateau disappears at a critical value of pinning strength $`f_c(r,L)2.4`$, where a pure power-law distribution with virtually infinite cut-off (controlled only by the system size) occurs. Here the measured slopes give $`\tau _s=1.74`$ and $`\tau _t=2.14`$ for the fixed driving rate $`\mathrm{\Delta }E/J=1\times 10^3`$ (see Fig. 3).

At strong pinning (large $`f`$) many small domains are being nucleated at isolated points in the bulk and grow slowly along a few active sites per time step. The maximum sizes of domains are finite and decrease with stronger pinning (cf. top two curves in Fig. 3). Here the measured slopes give $`\tau _s=1.56`$ and $`\tau _t=1.94`$ for size and duration, respectively, suggesting another universality class of the scaling behavior, compared to the case of weak pinning. The transition to the strong pinning regime for $`f>f_c`$ is marked by disappearance of the critical depinning of large domain walls at all fields and domination of the nucleation processes. Numerical data in Fig. 3 suggest that the critical value $`f_c`$ is close to 2.4 for the used system size and driving rate.

Qualitatively similar behavior with two universality classes but with different exponents, was found in the case of two-dimensional system where an extended one-dimensional domain wall was initially prepared . For comparison, in Table 1 we summarize exponents both for 3-dimensional and for 2-dimensional case. (Notice that the listed exponents for three-dimensional system are obtained for a small finite driving rate $`\mathrm{\Delta }E/J=4\times 10^4`$, whereas the simulations at zero driving rate were employed in two dimensions.) The results for 2-dimensional case are expected to be relevant for domains on the surfaces of the bulk samples (see Ref. ) and for thin films in which the width of the domain wall exceed the film thickness. An interesting physical realization of varying strength of pinning was suggested recently  in ferromagnetic $`\mathrm{C}_o/C_oO`$ film by varying temperature through the transition temperature at which the substrate orders anti-ferromagnetically. As already mentioned, in ferroelectrics the local random fields are physical . Therefore varying the strength of pinning due to random fields in ferroelectric thin films and bulk materials can be achieved in a more direct way, offering diverse possibilities to detect the change in the scaling behavior of the switching current noise and critical properties on the hysteresis.

### B Varying driving rate

In the following we examine how the scaling properties of the switching current noise change when the driving rate is increased. Very low frequencies of the driving field are usually not accessible on experiments, both for physical and technical reasons. A finite driving rate, which is defined as $`r<dE(t)/dt>`$, can be approximated with $`\mathrm{\Delta }E/J`$ in the case of stepwise increase of the external field in the absence of back flips $`b=0`$. For finite rate $`r>0`$ the switching may start at different parts of the sample and the measured signal is a superposition of these separate events. In contrast to the mathematical limit of infinitely slow driving, by finite field jump $`\mathrm{\Delta }E`$ the domain wall may get enough energy to overcome several pinning centers at the same time. Thus either different small domains are initiated simultaneously, or different parts of a large domain wall advance in parallel. In both cases the statistical properties of the signal are changed. Merging of many individual jumps makes the occurrence of large events more probable.

For a range of values of driving rates we still find a power-law behavior with reduced values of the scaling exponents and enlarged cut-offs. In Fig. 4 we show how the distributions for sizes and durations of pulses vary when the driving rate $`\mathrm{\Delta }E/J`$ is increased, while the strength of pinning was kept constant in the strong pinning regime. In Fig. 5 the same type of analysis was done for the case of weak pinning. In the strong pinning regime a dramatic increase of the cut-offs in the distribution of switched area occurs due to merging of many small signals. Whereas, in the weak pinning the relative enlargement of the cut-off is smaller, that can be understood as an effect of parallel advancement of different sections of the same extended wall, rather than merging of different domains. Both in strong and weak pinning case the increased incidence of large events leads to decrease of the scaling exponents. We find logarithmic dependence of the exponents for the size and duration of events on $`\mathrm{\Delta }E/J`$ both in weak and strong pinning regimes. The results are shown in Fig. 6 and the fit lines are given by the general form

$$\tau _X=A_XB_X\mathrm{ln}(\mathrm{\Delta }E/J).$$
(4)
For the strong pinning both area and duration exponents decrease equally, i.e., $`B_sB_t0.08`$, with the constants $`A_s1.0`$ and $`A_t1.4`$. In the weak pinning, however, the duration exponent decreases with nearly double rate compared to the size exponent. We find $`B_t2B_s0.11`$, (see Fig. 6) and $`A_t1.37`$, $`A_s1.30`$ within numerical error bars.

It is important to note that the measuring conditions depend on the experimental set up and they can vary from one experiment to the other. The effective driving rate depends on the steps in the external field and the degree of compensation of depolarization fields inside the sample. In real experiments the driving rate is always finite, resulting in the scaling exponents that are always smaller (and the fractal dimensions larger) compared to those computed in the theoretical limit of zero driving rate $`r0`$. Therefore, in order to compare measured and predicted scaling exponents it is advisable to make several measurements at different driving rates and find out if the results, for instance for the area exponent $`\tau _s`$, belong to a line in the range 1.66—1.40, which is compatible with occurrence of small domains, or in the range 1.78—1.60, suggesting that the sample has extended plane domain walls. When the driving rate is too large (cf. top line on both panels in Fig. 5), the scaling behavior of the distributions is lost. The cumulative probability distribution shows curvature both due to the lack of small events and the excess of large events.

## IV  Polarization reversal and back-switching

In the previous section we have shown how different kinds of domain walls in the sample determine the scaling properties of switching current noise, reflecting the nature of domain-wall motion when the electric field is slowly ramped. Here we consider another type of driving. First, in the initially polarized state (all spins down) we impose a large opposite field $`E_s`$ and keep it until polarization reverses and reaches saturation in that field. Then we switch off the field to zero and monitor how the polarization decays. Both processes depend on the type of domains in the sample, and thus, on the strength of pinning. The results for polarization as a function of time are shown in Fig. 7 for various strengths of pinning $`f`$.

For weak pinning spontaneous back switching of the polarization is exponentially fast, reaching the value $`P_r`$ of remanent polarization. Whereas, the relaxation becomes increasingly slow as the strength of pinning $`f`$ increases. The limiting values of $`P_r`$ are lower, which correspond to increasingly slimmer loops for strong pinning (cf. form of the hysteresis loop in Fig. 1a). In general, polarization decay in Fig. 7 can be approximated with a stretch exponential law $`P(t)=A+B\mathrm{exp}[(x/\tau )^\sigma ]`$, with $`\sigma =\sigma (f)`$ decreasing with the pinning strength. We have $`\sigma `$1, 0.92 and 0.80 for $`f=`$ 2.4, 3.4, and 4.4, respectively. Qualitatively similar back-switching was recently observed in doped $`\mathrm{SBN}:\mathrm{Ce}^{3+}`$ relaxor ferroelectric with two different stretching exponents corresponding respectively to illuminated and non-illuminated sample .

The polarization reversal in the strong constant field shows slow increase in the presence of small domains, but with another type of time dependence (see Fig. 7). The positive part of the polarization curve is fitted by the power-law increase $`P(t)=1C(ta)^\nu `$, where $`C0.8`$, $`a`$ is the respective shift (see Fig. 7) and the exponent $`\nu `$ 2, 1.75, and 1.68, for $`f=`$ 2.4, 3.4, and 4.4, respectively. Again, the characteristic exponent $`\nu `$ decreases when the pinning becomes stronger, thus contributing to the slow dynamics but within a power-law, rather than a stretch-exponential dependence.

## V Conclusions and discussion

We have shown that, in addition to the measurements of the hysteresis loop properties and smooth component of the switching current, the domain structure in a ferroelectric sample can be indirectly determined by the non-invasive dynamic measurements of the switching current noise and back-switching processes. Fractal character of the dynamic response indicates that a nontrivial domain structure can occur under slow field driving. The measured properties in these processes can be expressed in terms of universal exponents, which turn to be related in a unique way to the underlying structure and motion of domain walls. We have shown how the universal features of the switching current noise are quantified through two distinct sets of scaling exponents. The universal scaling behavior is associated with a dominant scenario of the domain wall dynamics: either the motion and depinning of large domain walls, in the case of weak pinning, or bulk nucleation and growth of many small 3-dimensional domains, in the case of strong pinning. The pinning in our model is provided by the local electric random fields.

Similarly, the polarization reversal in a strong applied field, and back-switching of the domains when the field is set to zero, are quick when large domains occurs, whereas they are slow when many small domains are nucleated. In the latter case the polarization relaxes according to a stretch-exponential law (the stretching exponent $`\sigma <1`$) whereas it reverses with a fractal power of time (the characteristic exponent $`\nu <2`$) when a large opposite field is switched on. In the case of weak pinning, on the contrary, the relaxation is exponential ($`\sigma 1`$) and the reversal algebraic in time with $`\nu =2`$. Dependence of the stretching exponent $`\sigma (f)`$ on the pinning strength is in qualitative agreement with measurements of back-switching in doped $`\mathrm{SBN}:\mathrm{Ce}^{3+}`$ relaxor ferroelectrics , again suggesting that the local electric fields contribute significantly to slow the dynamics in this class of relaxors . (Note that the characteristic time appears unrealistically small in our lattice model.) We are not aware of measurements of temporal behavior of the polarization reversal and of the switching current noise in these systems.

The observed scaling behavior of the switching current depends on the rate at which the driving field changes. We have demonstrated that in both universality classes the scaling exponents decrease logarithmically with increased driving rate $`\mathrm{\Delta }E/J`$. This is in sharp contrast to the cases where the dominant type of disorder is due to random vacancies in spin-systems  and other avalanche-type dynamic systems driven at a finite rate , where the scaling exponents decrease and fractal dimensions increase linearly with the driving rate. Therefore, occurrence of the logarithmic dependences of the scaling exponents on the driving frequency may serve as an additional test of the presence of the real random fields in ferroelectrics.

Throughout this work we used the simulation algorithm at finite driving rates $`\mathrm{\Delta }E>0`$, that mimics the realistic situation in experiments. The obtained set of scaling exponents, together with their dependences on the driving rate, belongs to one (of the two) universality classes, which implies a definite domain structure of the system, and hence the expected character of the switching processes. Our results suggest that a crossover from one to the other universal behavior both in the switching current noise and in the polarization reversal and relaxation can be achieved by varying the strength of pinning due to random electric fields. As mentioned above, in practice there are various ways to accomplished the crossover by reducing the pinning in ferroelectrics . Here we have not discussed the properties of the dynamic phase transition separating these two scaling behaviors at a critical value $`f_c`$, nor the critical exponents describing the depinning transition along the critical field $`H_d(f)`$ for $`f<f_c`$. This questions require an additional study and proper finite size scaling analysis, for instance along the lines as done for the two-dimensional case in Ref.. Finally, by allowing a frequent appearance of the unscreened depolarization fields a finite probability of negative pulses in the switching current builds up ($`b0`$). The presence of negative pulses slows the switching process, interferes with the effects of finite driving rates, and alters the scaling properties of noise. A more complete analysis of the switching processes for non-vanishing probability $`b0`$ is left for a future study.

Finally, it should be stressed that, in spite of its simplicity, the proposed model seems reasonable for dynamic properties of three-dimensional uniaxial relaxor ferroelectrics. Whereas, in the single-crystalline normal ferroelectrics a more realistic modeling of the effects of interactions and depolarization fields on domain structure is required. Nevertheless, we believe that, when restricted only to the universal scaling behavior of the noisy component of the switching current, the presented results for the case of large domains may not depend significantly on the details of interactions. Further research on the subject is necessary.

###### Acknowledgements.

Work is supported by the Ministry of Education, Science and Sports of the Republic of Slovenia.


text:
# Effective description of the LOFF phase of QCD

## 1 Introduction

In this talk We will discuss the recently proposed crystalline color superconducting phase of QCD. This state is the QCD analogus to a state studied in QED by Larkin and Ovchinnikov and Fulde and Ferrell. Therefore it has been named after them LOFF state.
It is known that at very high density and at low temperature quark matter exhibits color superconductivity. Color superconductivity is a phenomenon analogus to the BCS superconductivity of QED. According to the BCS theory, whenever there is an arbitrary small attractive interaction between electrons, the Fermi sea becomes unstable with respect to the formation of bound pairs of electrons. Indeed the creation of a Cooper pair costs no free energy, because each electron is on its own Fermi surface. On the other hand there is an energy gain due to the attraction of electrons. We call $`\mathrm{\Delta }`$ the binding energy of each pair with respect to the Fermi level. It is possible to show that the most energetically favored arrangement of electrons is made of two electrons with opposite momenta and spins. That is the condensate has spin zero and momentum zero. This is a normal superconductor. Now let us see what is a LOFF superconductor.

## 2 LOFF phase in QED

Let us consider a system made of a ferromagnetic alloy containing paramagnetic impurities. In first approximation the action of the impurities on the electrons may be viewed as a constant self-consistent exchange field. Due to this field the Fermi surfaces of up and down electrons (with respect to the direction of polarization of the medium) split, see figure (1). The half-separation of the Fermi surfaces will be

$$I=NSa,$$
(1)
where $`N`$ is the concentration of impurities, $`S`$ is their spin and $`a`$ is the integral of the exchange potential. To determine the energetically favored state we have to compare the BCS free energy with the free energy of the unpaired state. One has:

$$F_{\mathrm{𝐵𝐶𝑆}}F_{free}I^2\mathrm{\Delta }^2/2.$$
(2)
Therefore for $`I>\mathrm{\Delta }/\sqrt{2}`$ the BCS state is no more energetically favorite. For $`I=\mathrm{\Delta }/\sqrt{2}`$ one expects that a first order phase transition takes place from the BCS to the normal state. Or at least this is what one should naively expect.

Instead it has been shown that for

$$\frac{\mathrm{\Delta }}{\sqrt{2}}<I<.75\mathrm{\Delta },$$
(3)
it is energetically favorite a superconducting phase characterized by a condensate of spin zero but with non zero total momentum. We call $`\mathrm{\Delta }_{LOFF}\mathrm{\Delta }_{LOFF}(I)`$ the electron binding energy when the distance between the Fermi surfaces is $`2I`$. Increasing the half-distance between the Fermi spheres, from $`\frac{\mathrm{\Delta }}{\sqrt{2}}`$ to $`.75\mathrm{\Delta }`$, the LOFF gap decreases to zero. Therefore, for $`I.75\mathrm{\Delta }`$, a second order phase transition from the LOFF phase to the normal phase takes place. In the LOFF phase pairs are made of fermions which remain close to their own Fermi surfaces as shown in figure (1). Therefore the creation of a pair costs no free energy. The pair has total momentum $`2\stackrel{}{q}`$, given by the formula

$$\stackrel{}{K}_u+\stackrel{}{K}_d=2\stackrel{}{q},$$
(4)
where $`\stackrel{}{K}_u`$ and $`\stackrel{}{K}_d`$ are the momenta of up and down electrons. The magnitude of $`q`$ is determined by minimizing the free energy of the LOFF state. The direction of $`\stackrel{}{q}`$ is chosen spontaneously by the condensate. We observe that not all electrons can condense in LOFF pairs. Only the electrons of a restricted region of the phase space can pair because of equation (4). Therefore $`\mathrm{\Delta }_{LOFF}\mathrm{\Delta }_{\mathrm{𝐵𝐶𝑆}}`$ because of the reduced phase space available for coupling. Moreover $`\mathrm{\Delta }_{LOFF}`$ turns out to be a periodic function of the coordinates and we have a crystalline structure described by $`\mathrm{\Delta }_{LOFF}(\stackrel{}{r})`$.

## 3 LOFF phase in QCD

Let us consider quark matter at asymptotic density and at low temperature. In the previous section we have learnt that the LOFF phase may arise when fermions have different chemical potentials, and the potential difference lies inside a certain window. Crystalline color superconductivity is also expected to occur in case of different quark masses. In what follows we shall consider matter made of massless quarks and electrons, in weak equilibrium and electrically neutral. If we impose weak equilibrium we have $`\mu _u=\overline{\mu }\frac{2}{3}\mu _e`$ and $`\mu _{d,s}=\overline{\mu }+\frac{1}{3}\mu _e`$ where $`\mu _{u,d,s,e}`$ are quarks and electrons chemical potential. Requiring electrical neutrality we get:

$$\frac{2}{3}N_u\frac{1}{3}N_d\frac{1}{3}N_sN_e=0,$$
(5)
where $`N_{u,d,s,e}`$ are the concentrations of quarks and electrons. When $`m_s=0`$ we get $`N_u=N_d=N_s`$, $`N_e=0`$ and we expect matter to be in the color-flavor locked phase (CFL) . For $`m_s\overline{\mu }`$ we have $`\mu _u=\overline{\mu }\delta \mu ,\mu _d=\overline{\mu }+\delta \mu `$, and there are just two active flavors.

## 4 Velocity dependent effective lagrangian

In this Section we give a brief review of the effective lagrangian approach for the LOFF phase, based on velocity-dependent fields. We perform a Fourier transformation of the fermion field

$$\psi (x)=\frac{d^4p}{(2\pi )^4}e^{ipx}\psi (p),$$
(6)
and we decompose the fermion momentum as

$$p_i^\mu =\mu _iv_i^\mu +\mathrm{}_i^\mu ,$$
(7)
where $`i=1,2`$ is a flavor index, $`v_i^\mu =(0,\stackrel{}{v}_i)`$ and $`\stackrel{}{v}_i`$ the Fermi velocity (for massless fermions $`|\stackrel{}{v}|=1`$); finally $`\mathrm{}_i^\mu `$ is a residual momentum. By the decomposition (7) only the positive energy component $`\psi _+`$ of the fermion field survives in the lagrangian in the $`\mu \mathrm{}`$ limit while the negative energy component $`\psi _{}`$ can be integrated out. These effective fields are velocity dependent and are related to the original fields by

$$\psi (x)=\underset{\stackrel{}{v}}{}e^{i\mu vx}\left[\psi _+(x)+\psi _{}(x)\right],$$
(8)
where

$$\psi _\pm (x)=\frac{1\pm \stackrel{}{\alpha }\stackrel{}{v}}{2}\frac{d\stackrel{}{\mathrm{}}}{(2\pi )^3}_{\mathrm{}}^+\mathrm{}\frac{d\mathrm{}_0}{2\pi }e^{i\mathrm{}x}\psi _\stackrel{}{v}(\mathrm{}).$$
(9)
Here $`{\displaystyle \underset{\stackrel{}{v}}{}}`$ means an average over the Fermi velocities and

$$\psi _\pm (x)\psi _{\pm ,\stackrel{}{v}}(x)$$
(10)
are velocity-dependent fields.
Now we write the condensates. We have a scalar condensate (sc)

$$<ϵ_{ij}ϵ_{\alpha \beta 3}\psi ^{i\alpha }(\stackrel{}{x})C\psi ^{j\beta }(\stackrel{}{x})>\mathrm{\Delta }_Ae^{2i\stackrel{}{q}\stackrel{}{x}}$$
(11)
and a spin 1 vector condensate (vc)

$$<\sigma _{ij}^1ϵ_{\alpha \beta 3}\psi ^{i\alpha }(\stackrel{}{x})C\sigma ^{03}\psi ^{j\beta }(\stackrel{}{x})>\mathrm{\Delta }_Be^{2i\stackrel{}{q}\stackrel{}{x}}.$$
(12)
The lagrangian for the (sc) is

$`_\mathrm{\Delta }^{(s)}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }_A}{2}}{\displaystyle \underset{\stackrel{}{v}_1,\stackrel{}{v}_2}{}}\mathrm{exp}\{i\stackrel{}{x}\stackrel{}{\alpha }(\stackrel{}{v}_1,\stackrel{}{v}_2,\stackrel{}{q})\}ϵ_{ij}ϵ^{\alpha \beta 3}\psi _{\stackrel{}{v}_i\alpha }(x)C\psi _{\stackrel{}{v}_j\beta }(x)`$ (14)
$`(LR)+\mathrm{h}.\mathrm{c}.,`$
where $`\stackrel{}{\alpha }(\stackrel{}{v}_1,\stackrel{}{v}_2,\stackrel{}{q})=2\stackrel{}{q}\mu _1\stackrel{}{v}_1\mu _2\stackrel{}{v}_2`$. We take $`\stackrel{}{q}=(0,0,q)`$, therefore the components of $`\stackrel{}{\alpha }`$ satisfy the equations

$`\alpha _x`$ $`=`$ $`\mu _1\mathrm{sin}\alpha _1\mathrm{cos}\varphi _1\mu _2\mathrm{sin}\alpha _2\mathrm{cos}\varphi _2,`$ (15)
$`\alpha _y`$ $`=`$ $`\mu _1\mathrm{sin}\alpha _1\mathrm{sin}\varphi _1\mu _2\mathrm{sin}\alpha _2\mathrm{sin}\varphi _2,`$ (16)
$`\alpha _z`$ $`=`$ $`2q\mu _1\mathrm{cos}\alpha _1\mu _2\mathrm{cos}\alpha _2.`$ (17)
In the limit $`\mu _1,\mu _2\mathrm{}`$ we get

$`\varphi _2`$ $`=`$ $`\varphi _1+\pi `$ (18)
$`\alpha _2`$ $``$ $`\theta _q=\mathrm{arccos}{\displaystyle \frac{\delta \mu }{q}}+𝒪\left({\displaystyle \frac{\delta \mu }{\mu }}\right),`$ (19)
$`\alpha _1`$ $`=`$ $`\alpha _2+\pi +𝒪\left({\displaystyle \frac{\delta \mu }{\mu }}\right).`$ (20)
Therefore the velocities are almost opposite $`\stackrel{}{v}_1\stackrel{}{v}_2`$, and we have to deal with a no more symmetric sum over velocities, because the angle $`\alpha _2`$ is fixed.

## 5 Phonon quark interaction

The condensates (11) and (12) explicitly break rotations and translations. We have an induced lattice structure given by parallel planes perpendicular to $`\stackrel{}{n}=\stackrel{}{q}/|\stackrel{}{q}|`$:

$$\stackrel{}{n}\stackrel{}{x}=\frac{\pi k}{q}(k=0,\pm 1,\pm 2,\mathrm{}).$$
(21)
Lattice planes are allowed to fluctuate in two ways. We describe these fluctuations by means of two fields $`\varphi `$ and $`R`$:

$$\stackrel{}{n}\stackrel{}{x}\stackrel{}{n}\stackrel{}{x}+\frac{\varphi }{2qf}\frac{\mathrm{\Phi }}{2q},$$
(22)
$$\stackrel{}{n}\stackrel{}{R}.$$
(23)
But $`\varphi `$ and $`R`$ are not independent, indeed it is possible to show that

$$\stackrel{}{R}=\frac{\stackrel{}{}\mathrm{\Phi }}{|\stackrel{}{}\mathrm{\Phi }|}.$$
(24)
Therefore there is just one independent degree of freedom, i.e. the phonon.

## 6 Effective lagrangian

By a bosonization procedure it is possible to get the effective action for the phonon and to calculate the polarization tensor. At the lowest order in the momentum $`\stackrel{}{p}`$ of the phonon we have

$$\mathrm{\Pi }(0)=\mathrm{\Pi }(0)_{s.e.}+\mathrm{\Pi }(0)_{tad}=\mathrm{\hspace{0.17em}0},$$
(25)
$$\mathrm{\Pi }(p)=\frac{\mu ^2}{4\pi ^2f^2}\left[p_0^2v_{}^2(p_x^2+p_y^2)v_{}^2p_z^2\right],$$
(26)
where

$`v_{}^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{sin}^2\theta _q+𝒪\left({\displaystyle \frac{\mathrm{\Delta }^{(v)}}{q}}\right)^2,`$ (27)
$`v_{}^2`$ $`=`$ $`\mathrm{cos}^2\theta _q,`$ (28)
are the velocities perpendicular and parallel to $`\stackrel{}{q}`$. From eq.(25) we see that the phonon is massless. From eq.(26) we have the dispersion law for the phonon:

$$E(\stackrel{}{p})=\sqrt{v_{}^2(p_x^2+p_y^2)+v_{}^2p_z^2}.$$
(29)
In conclusion we can say that the dispersion law for the phonon is anisotropic in two ways. Indeed the velocity of propagation of the phonon in the plane perpendicular to $`\stackrel{}{q}`$, i.e. $`v_{}`$, is not equal to $`v_{}`$, that is the velocity of propagation of the phonon along $`\stackrel{}{q}`$. Moreover $`p_z`$ is a quasi-momentum.


text:
# A study of accuracy and precision in oligonucleotide arrays: extracting more signal at large concentrations

### Introduction

High-density oligonucleotide arrays manufactured by Affymetrix are among the most sensitive and reliable microarray technology \[Chee et al., 1996, Lipshutz et al., 1999\] available. Based on a photolithographic oligonucleotide deposition process, labeled and amplified mRNA transcripts are probed by 22-40 (depending on chip models) short DNA sequence each 25 bases long. The probes are preferentially picked near the 3’ end of the mRNA sequence, because of the limited efficiencies of reverse transcription enzymes. In addition, the probes come in two varieties: half are perfect matches (PM) identical to templates found in databases, and the other half are single mismatches (MM), carrying a single base substitution in the middle (13th) base of the sequence. MM probes were introduced to serve as controls for non-specific hybridization, and most analysis methods postulate that the actual signal (the target’s mRNA concentration) to be proportional to the difference of match versus mismatch (PM-MM).

The purpose of this work is twofold. First, we present a detailed calibration study of GeneChips. Specifically, we apply latest analysis methods (MAS 5.0 algorithm and others) to a large yeast calibration dataset, in which a number of transcripts are hybridized at known concentrations. We investigate the concentration dependence of both the accuracy and precision of differential expression scores. By accurate, we mean that the reported numerical ratios values are close to the known expression ratios, and have therefore little bias. On the other hand, a measurement is precise if it has a low noise level, also referred to as small variance. Our results show that the ability of conventional analysis techniques to detect small changes strongly deteriorates toward high transcript concentrations. While the variance is smallest for high concentrations, it appears that the question of the bias in this regime has been neglected. In fact, the bias is strong enough that real changes of 2 (even 4) can often not be detected. This sounds at first counter-intuitive, which we believe is rooted in the following widespread interpretation of hybridization data. Namely, when examining the data from two replicated conditions (Figure 1a) most would focus on the low intensity region, and observe how noisy this regime appears to be in comparison to the high-intensity tail. However, this view is misleading, as it does not consider the question of the bias. Turning to a comparison of two different conditions (Figure 1b), we notice that the noise envelope is essentially unchanged, and that real changes appear as points lying distinctively outside the noise cloud. Looking at multiple such comparisons, we would then conclude that the high intensity data is almost always very tightly scattered about the diagonal, and that there are rarely genes in that region that show fold changes greater than, say, 1.5 or 2. The interpretation that no differential regulation occurs in highly expressed transcripts seems unlikely. In fact, we show evidence that real changes are often compressed for large concentrations. This saturation effect can actually be observed in Affymetrix’s own data111Figure 7 at
http://www.affymetrix.com/products/algorithms\_tech.html although the issue is not commented there (a qualitative report has been given in \[Chudin et al., 2001\]). The physical origin for the compression effect invokes non-linear probe affinities and chemical saturation. This is a separate issue from optical saturation (cf. Results). Chemical saturation occurs below the detector threshold and is attributed to the fact that some probes will exhaust their binding capacities at relatively low concentrations, simply because their binding affinities are high. Binding affinities are in fact very sensitive to the sequence composition, resulting in measured brightnesses that usually vary by several decades within a given probeset \[Naef et al., 2002\].

Our second goal is to present an analysis method that reduces the bias at high concentrations. Our approach uses all PM and MM probes equally, in contrast to the standard view in which the PMs are thought to carry the signal, while the MMs serve as non-specific controls. In fact, it has become clear that the MM probes also track the signal, usually with lower (although often with higher) affinities than the PM’s \[Naef et al., 2002\]. In that sense, the MM’s should be viewed as a set of on average lower affinity probes. It is then reasonable to expect that some MM probes will more accurate at high intensities, since they will be less affected by saturation than the the PM’s (cf. Figure 3).

### Methods

The existing methods for the analysis of the raw data fall into two main classes. The first methods are similar to Affymetrix’s Microarray Suite software, providing absolute intensities on a chip by chip basis, or differential expression ratios from two experiments \[Affymetrix, 2001, Naef et al., 2001, Naef et al., 2002\]. The second class are called “model-based” approaches \[Li & Wong, 2001\], and attempt to fit the probe affinities from a large number of experiments.

The method described below belongs to the second class and is specifically designed for improved accuracy in the compressive high-intensity regime. It is based on ideas borrowed from the theory of locally linear embeddings \[Roweis & Saul, 2000\].

#### Notation

We construct the following matrix

$`A_i^j`$ $`=`$ $`\{\begin{array}{cc}\mathrm{𝑃𝑀}_i^j& 1jN_p\\ \mathrm{𝑀𝑀}_i^{(jN_p)}& N_p<j2N_p\end{array}`$
or in expanded notation

$`A_i^j`$ $`=`$ $`\left[\begin{array}{cccccc}\mathrm{𝑃𝑀}_1^1& \mathrm{}& \mathrm{𝑃𝑀}_1^{N_p}& \mathrm{𝑀𝑀}_1^1& \mathrm{}& \mathrm{𝑀𝑀}_1^{N_p}\\ \mathrm{}& \mathrm{𝑃𝑀}_i^j& \mathrm{}& \mathrm{}& \mathrm{𝑀𝑀}_i^j& \mathrm{}\\ \mathrm{𝑃𝑀}_{N_e}^1& \mathrm{}& \mathrm{𝑃𝑀}_{N_p}^{N_e}& \mathrm{𝑀𝑀}_{N_e}^1& \mathrm{}& \mathrm{𝑀𝑀}_{N_p}^{N_e}\end{array}\right]`$
which contains the raw, background subtracted and normalized data. By background, we mean fluorescence background, which we identify by fitting a Gaussian distribution to the subset of all $`(PM,MM)`$ pairs satisfying the criterion that $`|PMMM|<ϵ`$, with $`ϵ=50`$. This provides us with the mean and SD in the background fluorescence. 222for details, see
http://xxx.lanl.gov/abs/physics/0102010. For a fair comparison of compression effects in various methods, we used the global normalization factors from the MAS 5.0 software in all cases, however, the technique remains applicable with other normalization schemes.

$`N_p`$ is the number of probe pairs and $`N_e`$ is the number of experiments. We introduce a set of weights $`w_i`$ such that

$$\underset{i=1}{\overset{N_e}{}}w_i=1$$

and define the column means (or center of mass)

$$m^j=\underset{i=1}{\overset{N_e}{}}w_i\mathrm{log}(A_i^j)$$

Note, we are computing the mean of the logs of the components of $`A_i^j`$.

#### Local principal component analysis

Local embeddings are adequate in situations where compression is important because non-linearities (resulting from chemical saturation) affect the one-dimensional manifold $`\{\mathrm{𝑃𝑀}^j(c),\mathrm{𝑀𝑀}^j(c)\}`$ (the concentration $`c`$ is the one-dimensional ’curve’ parameter) by giving it a non-zero local curvature. The results section, in particular Figure 3 contains ample evidence that these non-linearities are significant (cf. also Figure 2 in \[Chudin et al., 2001\]). Our method is a multidimensional generalization of the schematic depicted in Figure 2, which shows the typical situation of two probes in which one of the probes (PM2) saturates at concentrations lower than the other. If both probes were perfectly linear, the curve would be a straight line with slope 1. In the multidimensional case, the directions of largest variation (analogous to D1 or D2 in Figure 2) are computed from the principal components of the matrix

$$\sqrt{w_i}\left(\mathrm{log}(A_i^j)m^j\right)=\underset{k=1}{\overset{N_p}{}}U_{ik}D_kV_k^j$$

which can easily be done via singular value decomposition (SVD). In order to reconstruct the concentrations, one needs first to consider the unspecified sign of the vector $`V_1^j`$ (when returned by the SVD routine), which has to be chosen such that the total amount of signal comes out positive. This is easily achieved by adjusting the sign of $`V_1^j`$ such that $`_{i,j}\mathrm{log}(A_i^j)V_1^j>0`$ . The logarithm of the concentration $`s_i`$, is then computed by projecting the original matrix onto the first principal component $`V_1^j`$, corrected by a factor $`v_{max}`$. This factor accounts for the fact that the vector $`V_1^j`$ is $`L_2`$-normalized ($`_j(V_1^j)^2=1`$) by definition in the SVD. The signal then reads:

$$s_i=v_{\mathrm{𝑚𝑎𝑥}}\underset{j=1}{\overset{N_p}{}}\mathrm{log}(A_i^j)V_1^j$$

where $`v_{\mathrm{𝑚𝑎𝑥}}=\mathrm{max}_j\left|V_1^j\right|`$. In addition, the above procedure automatically yields a signal-to-noise ($`S/N`$) measure for the entire probeset

$$\frac{S}{N}=\frac{D_1}{\sqrt{_{j=2}^{N_p}D_j^2}},$$

where $`\{D_k\}`$ are the singular values. Large $`S/N`$ values imply that the probeset measurements in the $`N_e`$ experiments had a well-defined direction of variation, and can for instance be used as a filter for identifying genes that exhibit significant changes across the experiments tested. In the results shown in Figure 4, we used the following weights

$$w_i=\frac{1}{W}\{\begin{array}{cc}1,& \widehat{s}_i>\mathrm{\Gamma },\hfill \\ 1/\left(1+b(ii_s)^2\right),& \widehat{s}_i\mathrm{\Gamma }\hfill \end{array}$$

where $`W=_iw_i`$, $`\widehat{s}_i`$ is the signal obtained with uniform weights, $`i_s=20`$ (out of 28 experiments), $`b=2`$, and $`\mathrm{\Gamma }=\widehat{S}_{i_s+1}`$ with $`\widehat{S}_i`$ being the ascendantly sorted $`\{\widehat{s}_i\}`$. In other words, lower concentration points are suppressed according to their rank (computed with uniform weights) using a slowly decaying Cauchy weight function. There are of course other weight functions that could serve the same purpose.

Note that the fitting procedure used in the Li-Wong \[Li & Wong, 2001\] method is identical to an SVD decomposition, however, with different input data than was use here. The three main difference between our method and the Li-Wong technique are: (i) in the analysis here, we used log transformed PM and MM intensities, rather than the bare PM-MM values; (ii) we introduced optional weights, which can account for non-linearities of the probe response in the high concentration regime; (iii) we subtract the column mean before we compute the principal components, which is crucial for capturing the local directions of variation. Indeed, as can be seen in Figure 2, the principal component would be dominated by the mean itself without subtraction.

### Results

#### Data sets

The yeast Latin square (LS) experiment is a calibration data set produced by Affymetrix, that uses a non-commercial yeast chip. 14 groups of 8 different genes, all with different probe sets, are spiked onto 14 chips at concentrations, in pM, corresponding to all cyclic permutations of the sequence (0, 0.25, 0.5, 1, 2, $`\mathrm{}`$, 1024). Hence, each group is probed at 13 different finite concentrations, logarithmically spaced over a range of magnitudes from 1 to 4096 (in Figures 3 and 4, we refer to these concentrations as ($`1=0.25\mathrm{pM}`$, $`2=0.5\mathrm{pM}`$ $`\mathrm{}`$ $`13=1024\mathrm{pM}`$), and each group is completely absent in one array. Besides the spiked-in target cRNA’s, a human RNA background was added to mimic cross-hybridization effects that would occur in a real experiment. In addition, each experiment was hybridized twice leading to 2 groups of 14 arrays called R1 and R2.

The reason why this dataset is attractive as compared to the similar human and e. coli datat available at www.netaffx.com are several. First, the e. coli data exhibits severe optical saturation, which interfers with the chemical saturation issue we are trying to address here. The yeast dataset, on the other hand, has virtually no optically saturated cells, as can be inferred the SD in the pixel intensities reported in the raw data (.CEL) files. In total, fewer than $`0.1\%`$ of the probes have $`SD=0`$ ($`SD=0`$ characterizes optical saturated cells). Further, optical saturation is no longer an issue in GeneChips with the current scanner settings. More important, the present dataset permits far better statistics, as the number of spiked-in genes is 112 as compared to 14 for the human chip. In the latter dataset, there is only one transcript per concentration group as compared to 8 in the yeast case. In fact, we verified that the compression effects discussed below are virtually identical in the human case (not shown).

#### Summary of two-array methods

The figures in the results section show the log-ratios as function of concentration in the form of boxplots. In these plots, the central rectangle indicates the median, 1st and 3rd quartiles ($`Q1`$ and $`Q3`$). The “whiskers” extend on either side up to the last point that is not considered an outlier. Outliers are explicitly drawn and are defined as laying outside the interval $`[Q11.5IQ,Q3+1.5IQ]`$, where $`IQ=Q3Q1`$ is the interquartile distance. For each method, we show three plots, the top two measure the false negative rate for ratios of 2 and 4 fold respectively, and the last one shows the false positive rate. For the top two plots, all combinations (within R1 and R2 separately) of arrays leading to ratios of 2 and 4 were considered, and plotted as function of their baseline intensity (the lesser of the two concentrations). For the third, each gene was compared between the groups R1 and R2, at identical concentrations. Of the $`814=112`$ transcripts, 8 were left out of the analysis because they did not lead to a signal that was tracking the concentrations at all (presumably due to bad probes or transcript quality).

In Figure 3, we summarize the results obtained by the Microarray Analysis Suite 5.0 (MAS 5.0) software and the “2 chips” method discussed in \[Naef et al., 2002\]. The later method computes for each gene probed in two arrays a ratio score $`R`$ such that

$$\mathrm{log}(R)=\underset{j}{\overset{robust}{}}\mathrm{log}(R^j)$$

is a robust geometric mean (a least trimmed squares estimator) of the probe ratios $`R^j`$. Figure 3 shows the cases where

$$R^j=\frac{\mathrm{𝑃𝑀}_1^j\mathrm{𝑀𝑀}_1^j}{\mathrm{𝑃𝑀}_2^j\mathrm{𝑀𝑀}_2^j}$$

and

$$R^j=\frac{\mathrm{𝑀𝑀}_1^j}{\mathrm{𝑀𝑀}_2^j}$$

In both cases, only probes with numerator and denominator above background are retained. The first case ($`\mathrm{𝑃𝑀}\mathrm{𝑀𝑀}`$) is in essence similar to the MAS 5.0 program, differences are in the choice of the Tuckey bi-weight mean as the robust estimator, and in the treatment of negative differences. For our purpose here, we like to think of the Affymetrix method as two-array, $`\mathrm{𝑃𝑀}\mathrm{𝑀𝑀}`$ based method. In all the results presented below, the arrays were scaled according to the MAS 5.0 default settings.

The main features of Figure 3 are: there is an optimal range of baseline concentrations ($`116`$ pM) in which the ratio values from both $`\mathrm{𝑃𝑀}\mathrm{𝑀𝑀}`$ methods (the two first columns) are fairly accurate, for both ratios of 2 and 4. For both lower and higher concentrations, there is a noticeable compression effect, which is most dramatic at the high end. At the highest baseline concentration (512 pM for the ratios of 2 and 256 pM for ratios of 4), changes of 2 are basically not detected and real changes of 4 are compressed on average to values around 1.25. The analysis of the false positive rate (last row) shows that both methods yield very tight reproducibility: the log2 ratio distributions are well centered around 0 and the interquartile distances are roughly intensity independent and smaller than $`0.2`$, meaning that 50% of the measurements fall in the ratio interval $`[0.93,1.07]`$. To be fair, we should point out that as a ($`\mathrm{𝑃𝑀}\mathrm{𝑀𝑀}`$) method, the MAS 5.0 algorithm is on average a bit cleaner, having slightly fewer outliers. However, we like to emphasize that the qualitative behavior in the two ($`\mathrm{𝑃𝑀}\mathrm{𝑀𝑀}`$) methods is unchanged, especially as far as the high-intensity compression is concerned. Further, similar behavior is also found using the ($`\mathrm{𝑃𝑀}\mathrm{𝑀𝑀}`$) Li-Wong method (data not shown). The above observations are consistent with what was reported in \[Chudin et al., 2001\], confirming that these effects are independent of the chip series.

The third column in Figure 3 illustrates our contention that the MM are in essence a set of lower affinity probes. We notice that using only the MM measurements in the two-array method changes the picture qualitatively. Whereas the low concentration regime is far worse than in the ($`\mathrm{𝑃𝑀}\mathrm{𝑀𝑀}`$) methods, the behavior toward the high end has changed and the drop off occurs now at higher concentrations: approximately 256 pM for the ratios of 2 and 128 pM for ratios of 4. On the other hand, even in the optimal range, the magnitude of the medians are always a bit lower than the real ratios, and the false positive rate also suffers. To summarize, this result suggests that if one is interested in accuracy at high concentrations, then the MM-only methods offers the best two-array alternative. We have tried other variations: PM only, or the double size set consisting in the merged PM and MM’s, both being worse at high concentrations than the MM only method.

#### Multi-array methods

The data analyzed using our new method is shown in Figure 4. It is clear that both are capable of reducing the high intensity compression, as compared to existing methods. The second column explicitly shows the higher accuracy of the local method. It should be noted, however, that the precision is significantly lower than with MAS 5.0, which is the trade-off to pay for higher signal detection. As compared to the “two chip” MM method, which was previously the least compressive in this regime, the medians are systematically more accurate. Also, the method does not perform well at low-concentrations which is expected since it was not designed for that range.

#### Significance scores

Although ratio score may suffer severe compression, there remains the possibility that they would be attributed a significant increase or decrease call. Figure 5 displays the relation between the MAS 5.0 log-ratios and their associated p-values. MAS 5.0 change ‘p-values’ $`p_M`$ are symmetric about 0.5 and designed such that the ratio score is called increased when $`p_M<\gamma _1`$ and decreased (D) when $`p_M>1\gamma _1`$, with a default $`\gamma _1=0.0025`$. This definition is not well suited for plotting purposes, we therefore work with $`p_{\mathrm{𝑀𝐴𝑆}}=p_M`$ when $`p_M<0.5`$, and $`p_{\mathrm{𝑀𝐴𝑆}}=1p_M`$ otherwise. This way, both I and D genes have $`p_{\mathrm{𝑀𝐴𝑆}}<\gamma _1`$, the direction being given by the sign of the log-ratio. The results show that there remarkably few false positive calls: only 4 out of 624 for concentrations $`c8\mathrm{pM}`$, and 6 of 728 when $`c16\mathrm{pM}`$. Fold changes of 4 are also well detected despite the compression at high intensities: there are 21 false negatives (and 3 false positives having ratios with the wrong sign) out of 1248 for $`c8\mathrm{pM}`$, and 84 of 1040 false negatives for $`c16\mathrm{pM}`$. The situation deteriorates for fold changes of 2, with 124 false negatives (and 3 false positives) out of 1248 for $`c8\mathrm{pM}`$, and 425 of 1248 false negatives for $`c16\mathrm{pM}`$.

### Discussion

We have shown that high-concentration bias is a serious issue in GeneChips, which is probably related to chemical saturation in the adsorption process of the target to the probes. Exploiting the broad range of affinities of different probes (PM and MM included) offers an approach toward improvements. However, the gain in accuracy comes with an expected decrease in precision, since effectively, the weight of a measurement is transferred to a smaller set of probes. Hence, the reduction in noise levels resulting from averaging over probes is diminished.

Our method should serve as a useful complement to those who use microarrays primarily as a gene discovery tool, and are interested in maximal signal detection. In fact, we often hear that severe constraints like pharmaceutical treatments or gene knockouts appear to have no detectable transcriptional effects. While there is the possibility that transcription regulatory networks can compensate for such changes, or that some effects would be mostly post-transcriptional, real transcriptional changes may also be masked by compressive effects like those discussed above.

### Conclusion

We have summarized the performance of existing methods for analyzing Affymetrix GeneChip data, using the yeast calibration dataset from Affymetrix. The results show unambiguously the compressive tendency of GeneChip measurements in the high-intensity range, namely that fold changes as large as 4 in expression levels can be reduced to fold changes barely larger than 1 (Figure 3). Interestingly, we showed that among the standard techniques, the one using only the MM signals offers the highest accuracy at high concentrations. Additionally, we have how it is possible to achieve higher accuracy at high concentrations by exploiting the probeset’s wide affinity range. One should realize, however, that saturation problems of the sort encountered present a hard challenge in signal processing, and it is therefore expected that higher accuracy is obtained at the expense of reduced precision.

Our observations raise the following design issue for oligonucleotide arrays. Since it will likely be difficult to manufacture oligonucleotide probes with linear responses over 4 or more decades in concentration, an option would be to optimize the design of probesets such that each of its probe would be optimally linear in smaller ranges (say 2 decades at most) centered around graded concentrations. In this way, the weights of a measurement could be transferred to an appropriate subset of probes known to be optimal in a given range. Hence, one would use a different set of probes for high or low concentration values to increase the overall dynamic range of a probeset.

## Acknowledgments

The authors are thankful to Affymetrix for having provided the useful calibration data. F. N. is a Bristol-Meyers Squibb Fellow and acknowledges support from the Swiss National Science Foundation. This work was also supported by NIH/NINDS Grant NS39662 (N. D. S) and the Meyer Foundation (M. M.)

### Figures