text:
# Brownian Dynamics Studies on DNA Gel Electrophoresis. II. ‘Defect’ Dynamics in the Elongation-Contraction Motion

## I Introduction

Gel electrophoresis is a major technique which separates polyelectrolytes according to their size. Among them gel electrophoresis of DNA is one of the most interesting problems of statistical physics of polymers. Computer simulations Deutsch (1988); Deutsch and Madden (1989) and experiments Masubuchi et al. (1993); Oana et al. (1994) reported a peculiar behavior of individual DNA that they exhibit contracted and extended forms alternatively under a constant field. Recently we have developed a Brownian dynamics (BD) method which simulates electrophoresis of DNA in a 3 dimensional (3D) space by a chain of electrolyte beads. Its details have been reported in our preceding paper (hereafter refereed to as I), where we have demonstrated ‘quasi-periodic’ time evolution of the velocity of the center of mass, $`v_\mathrm{G}(t)`$, and of the radius of the longer principal axis, $`R_l(t)`$, of the chain. The present paper is devoted to clarify the role of ‘defects’ introduced by de Gennes de Gennes (1971) on the elongation-contraction motion.

The concept of ‘defect’ proposed by de Gennes de Gennes (1971) is useful to investigate dynamics of individual polymers confined to a space with entanglement such as gel electrophoresis of DNA. A defect represents a relatively loose portion, or a slack segment, of a polymer and is assumed to behave as a particle of ideal gas. He showed that its diffusion in the tube made of other polymers leads to the famous growth law $`\varphi (t)t^{1/4}`$ in the short time regime, where $`\varphi (t)`$ is the mean square displacement of a segment in equilibrium.

Electrophoresis of a chain in an ideal regular network of other fixed chains has been extensively studied by Monte Carlo (MC) simulations. de Cruz et al. (1986); Duke (1990); Batoulis et al. (1989) Among them, the repton model proposed by Rubinstein is important. Rubinstein (1987) It simulates dynamics of a chain only by stochastic movements of reptons each of which can be regarded as embodiment of a defect. The model is applied to simulation of the constant field gel electrophoresis (CFGE) by Duke. Duke (1990) It was shown that the mobility is inversely proportional to the number of reptons in a very weak field region. The relation was consistent with what the reptation theory predicted under the vanishing field limit. Lerman and Frisch (1982); J.Lumpkin and Zimm (1982) Under finite fields, however, the results derived from the repton model qualitatively differ from those derived from the biased reptation theory Lumpkin et al. (1985): the mobility plateau is reached at lower fields in the former model than those predicted by the latter theory. Duke also observed snapshots of tubes which look like hooks, and argued that once a chain is hooked around a gel fiber, both its arms tend to stretch and rarely retract. This is because defects prefer to move in direction of the field and to create new tube segments which in turn prefer to align to the field. The observed fact that the chain rarely exhibits a retraction is considered to be an artifact of the MC method which uses local hopping rules. The incorporation of a hopping rule which allows non-local movement of defects was done by Duke and Viovy in their advanced version of the repton model. Duke and Viovy (1992a, b) Simulations of chains with over a thousand segments revealed that a chain exhibits ramified configurations where ‘hernias’ are created and annihilated under a constant field. Duke and Viovy (1992a)

Recently, we have investigated gel electrophoresis and diffusive dynamics of a chain in a space with immobile obstacles by making use of an extended version of the bond fluctuation model (BFM). Azuma and Takayama (1999) The conventional BFM (CBFM) Carmesin and Kremer (1988) is known as an efficient model to simulate general problems involved in polymers. However, it is faced with a problem which is similar to that of the repton model mentioned above: the chain motion extremely slows down when it is applied to CFGE under relatively large fields. Batoulis et al. (1989) To overcome this problem, we have empirically incorporated semi-local dynamics of slack parts of a chain. The latter correspond to ‘defects’ due to de Gennes and were named s-monomers. They are explicitly defined in each configuration appearing in the CBFM. We have called such a BFM with introduction of semi-local movements of s-monomers the extended BFM (EBFM). Both the DV model and the EBFM incorporate nonlocal movements of slack portions of DNA. There are, however, the following difference between the two methods. In the DV model any defect is not created or annihilated except at the extremity, whereas in the EBFM, s-monomers are continually created and annihilated anywhere in the chain, since it involves the CBFM trials of conventional local movements of all monomers.

Using the EBFM it has been confirmed that the trapping problem of the CBFM under large fields is overcome and the field-dependence of mobility is reproduced qualitatively as observed by the experiment. Hervet and Bean (1987) Furthermore, quasi-periodic behavior, or a damped oscillation, is observed in the correlation function of $`R_l(t)`$ in 3D Azuma and Takayama (2000) as it has been observed in the experiment Masubuchi et al. (1993); Oana et al. (1994) as well as by the BD calculations in 2D (without the excluded volume effect) Matsumoto and Doi (1994) and in 3D (with the effect). Azuma and Takayama (2001)

The defect is thus a quite important concept to understand DNA dynamics in gel. The purpose of the present paper is to clarify the microscopic bases of the defect dynamics, thereby presenting a microscopic support of the EBFM. For this purpose we examine the corresponding object in a chain, i.e., slack beads (abbreviated as s-beads) by means of the BD method. The work which is similar to this approach has been done by Masubuchi et al. Masubuchi et al. (1997) in their 2D BD study on the biased sinusoidal field gel electrophoresis (BSFGE). Looking at time evolution of snapshots of kinks which correspond to slack segments, they have found that a portion with less kinks develops near the segment which is hindered by an obstacle, and that the portion gradually shifts as the chain slides off the obstacle. They have argued that the growth of kinks plays a key role on the decrease of the mobility observed in chains with particular sizes under BSFGE.

In the present paper we examine static and dynamic properties of defects much more in details by means of our BD method for a chain in a 3D space. As described in our accompanying paper, hereafter referred to as I, Azuma and Takayama (2001) the method incorporates semi-microscopic ingredients of dynamics of a polymer, i.e., the excluded volume effect, the non-linear elastic force which keeps the distance between the two adjacent beads (segments) and the random force from a solvent. In each configuration of the chain simulated by the BD method s-beads, or ‘defects’, are explicitly specified in a quite similar way as the s-monomers have been defined in the EBFM. Then we investigate various statistics associated with the s-beads, or more explicitly, portions between the neighboring s-beads which we call extended parts of the chain.

In particular, we have concentrated on the normalized histogram of the extended parts of length $`n`$, $`r_n(t)`$, and that of their displacement along the chain direction, $`r_{\mathrm{\Delta }u}(t)`$. In equilibrium under a vanishing field it is found that both $`r_n(t)`$ and $`r_{\mathrm{\Delta }u}(t)`$ obey Poisson distributions as expected. In the elongation-contraction motion under CFGE, $`r_n(t)`$ exhibits a peculiar time evolution reflecting the elongation and contraction of the chain. The histogram $`r_{\mathrm{\Delta }u}(t)`$, on the other hand, does not fluctuate correspondingly to $`r_n(t)`$ but exhibits random hopping nature rather similar to the one observed in equilibrium. These results indicate that defects move quite stochastically even in the deterministic regime where the chain elongation-contraction motion looks apparently deterministic. These features of $`r_n(t)`$ and $`r_{\mathrm{\Delta }u}(t)`$ are confirmed to be in qualitative agreement with those observed for s-monomers in the EBFM.

The organization of the present paper is as follows. In the next section we briefly explain our numerical method, and discuss equilibrium dynamics of s-beads in Sec. III. Characteristic features of $`r_n(t)`$ and $`r_{\mathrm{\Delta }u}(t)`$ in the field range where a chain exhibits elongation-contraction motions are examined in Sec. IV. The corresponding results obtained by the EBFM are discussed in Sec. V, and the final section is devoted to the conclusion.

## II Numerical Method

The BD method we adopt in the present work is to solve the following Langevin equation of motion for a chain of spherical beads in a continuous 3D space,

$$\zeta \dot{𝒙}_i=q𝑬_\mathrm{b}+\underset{\alpha }{}\underset{j}{}𝑺_{ij}^\alpha +\underset{\alpha }{}\left[𝑭_i^\alpha 𝑭_{i1}^\alpha \right]+𝒇_i,$$
(1)
where $`𝒙_i`$ is the position of the $`i`$-th bead, and $`\zeta `$, $`q`$, and $`𝑬_\mathrm{b}`$ are the viscosity, the charge of a bead and the (bare) external field, respectively. The second and the third terms in r.h.s. are the constraining forces required for the beads to satisfy the conditions,

$`𝒙_{ij}=`$ $`𝒙_j𝒙_i`$ $`>1`$ (2)
$`𝒍_i=`$ $`𝒙_{i+1}𝒙_i`$ $`<\sqrt{2}.`$ (3)
putting the bead diameter unity. The last term in eqn. 1 stands for the random force from a solvent. Explicitly we adopt the following form for it

$$𝒇_i(t)=\underset{n}{}𝒇_i^n\delta (tn\mathrm{\Delta }t),$$
(4)
where $`\mathrm{\Delta }t`$ is regarded as the average period of random forces. The distribution of $`𝒇_i^n`$ is characterized by the averages

$`𝒇_i^n`$ $`=`$ $`\mathrm{𝟎}`$ (5)
$`𝒇_i^m𝒇_j^n`$ $`=`$ $`2\zeta k_\mathrm{B}T\mathrm{\Delta }t\delta _{ij}\delta _{mn}.`$ (6)
Here $`k_B`$ and $`T`$ are the Boltzmann constant and the absolute temperature, respectively.

The details of solution of the above set of equations have been explained in our accompanying paper I. As a representative of ‘defect’, or a loose portion of a chain in the BD calculation we define a slack bead (s-bead). For any sequences of three beads, if the distance between the first and the third is not greater than $`\sqrt{2}`$, the second one is regarded as an s-bead. This definition of ‘defect’ is the same as that of s-monomers in a 3D version of the EBFM we investigated before Azuma and Takayama (2000) (see in Sec. V). We call a sequence of beads between neighboring s-beads an extended part of the chain.

## III ‘Defects’ in equilibrium

In equilibrium under a vanishing field, ‘defects’ are considered to be distributed randomly along the polymer coordinate. Correspondingly, in our BD analysis, s-beads are randomly distributed throughout the chain. Then the ‘extended’ parts consisting of $`n`$ beads is expected to obey a distribution

$$\mathrm{\Psi }(n)=\frac{(Mn)(1\rho )^n}{_{i=1}^M(Mi)(1\rho )^i},$$
(7)
where $`M`$ is the total number of beads of the chain and $`\rho `$ the density of s-beads. For $`Mn1`$, $`\mathrm{\Psi }(n)`$ reduces to a Poisson distribution:

$$\mathrm{\Psi }(n)\rho (1\rho )^n=\rho e^{n\mathrm{log}(\frac{1}{1\rho })}.$$
(8)
Next let us consider $`\mathrm{\Delta }u`$, displacement of the tangential component, of extended parts between neighboring s-beads. In our BD calculation, it is evaluated as

$$\mathrm{\Delta }u\left|\underset{i=i_m}{\overset{i_{m+1}}{}}\mathrm{\Delta }s_i\right|,\mathrm{with}\mathrm{\Delta }s_i=\mathrm{\Delta }𝒙_i\frac{𝒙_{i1}𝒙_{i+1}}{𝒙_{i1}𝒙_{i+1}}$$
(9)
where $`\mathrm{\Delta }𝒙_i`$ is the displacement of the $`i`$-th bead in one MD step solving the Langevin equation 1 and index $`i_m`$ denotes a figure which is numbered on the s-beads. From eqn. 1 we may assume that the distribution of $`\mathrm{\Delta }𝒙_i`$ is essentially given by that of the random force, $`𝒇_i^n`$, defined by eqn. 6. Then the distribution of $`\mathrm{\Delta }s_i`$ is written as

$$\mathrm{\Xi }(\mathrm{\Delta }s_i)=\sqrt{\frac{\zeta }{2\pi k_\mathrm{B}T}}\mathrm{exp}(\zeta \mathrm{\Delta }s_i^2/2k_\mathrm{B}T).$$
(10)
Using this $`\mathrm{\Xi }(\mathrm{\Delta }s_i)`$ and $`\mathrm{\Psi }(n)`$ of eqn. 8, we can evaluate the distribution of $`\mathrm{\Delta }u`$ as

$`\mathrm{\Phi }(\mathrm{\Delta }u)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑n\mathrm{\Psi }(n){\displaystyle \underset{i=1}{\overset{n}{}}}\left[{\displaystyle _{\mathrm{}}^{\mathrm{}}}\mathrm{\Xi }(\mathrm{\Delta }s_i)𝑑\mathrm{\Delta }s_i\right]\delta ({\displaystyle \underset{i=1}{\overset{n}{}}}\mathrm{\Delta }s_i\mathrm{\Delta }u)`$ (11)
$`=`$ $`\rho \sqrt{{\displaystyle \frac{2\pi ^2\zeta }{k_\mathrm{B}T\mathrm{log}\frac{1}{1\rho }}}}\mathrm{exp}\left[\sqrt{{\displaystyle \frac{\mathrm{log}\left(\frac{1}{1\rho }\right)}{k_\mathrm{B}T/2\zeta }}}\mathrm{\Delta }u\right].`$
Thus $`\mathrm{\Phi }(\mathrm{\Delta }u)`$ is expected to obey also a Poisson distribution with the width $`\mathrm{\Omega }_{\mathrm{\Delta }u}=1/\sqrt{\mathrm{log}(\frac{1}{1\rho })/(k_\mathrm{B}T/2\zeta )}`$. The latter is proportional to $`\sqrt{k_\mathrm{B}T}`$ if the $`T`$-dependence of $`\rho `$ is neglected.

In Fig. 1, we show the normalized histograms $`r_n`$ and $`r_{\mathrm{\Delta }u}`$ in equilibrium for various $`T`$ which are obtained by the BD calculation. Here $`r_n=N_n/_nN_n`$ where $`N_n`$ is the number that the extended parts of a length $`n`$ appear in the chain configurations simulated, and $`r_{\mathrm{\Delta }u}`$ is similarly defined. In agreement with $`\mathrm{\Psi }(n)`$ of eqn. 8 $`r_n`$ obeys a Poisson distribution. The exponential tail of $`r_n`$ at $`k_\mathrm{B}T=1`$ is fitted to $`e^{n/\omega }`$ with $`\omega n5.7`$, which yields $`\rho 0.17`$ through the relation $`\omega =1/\mathrm{log}\left(\frac{1}{1\rho }\right)`$. In the simulation we have also directly counted the number of s-beads, which yields the same s-bead density $`\rho 0.17`$ as expected. The distribution $`r_{\mathrm{\Delta }u}`$ of the tangential movement of the extended parts, which we regard as that of the s-beads, also obeys a Poisson distribution as seen in the same figure. Its exponential tail, specified by $`\mathrm{\Omega }_{\mathrm{\Delta }u}`$ as $`e^{n/\mathrm{\Omega }_{\mathrm{\Delta }u}}`$, significantly depends on $`T`$. As shown in Fig. 2, its dependence on $`T`$ is consistent with eqn. 11. Even the pre-factor $`0.30`$ of the proportionality to $`\sqrt{k_\mathrm{B}T}`$ agrees with the predicted value $`1/\sqrt{2\mathrm{log}\frac{1}{1\rho }}0.37`$ within $`20\%`$.

The results obtained above confirm the role of the ‘defects’, or the s-beads, on the static and dynamic properties of a polymer in equilibrium. In particular, it is natural to regard $`cnl`$ as the ‘persistence length’ of the chain, where $`n`$ is the mean value of $`n`$ and $`c`$ is a numerical constant nearly equal to unity. In this context, we note $`\sqrt{2}>l>1`$, where $`l`$ is the mean distance between two adjacent beads. In fact we have already argued in I that this interpretation yields an appropriate conversion factor between the length scales of the BD simulation and the experiment. Here we have to emphasize that the value of $`n`$ is almost independent of $`T`$ in the range $`k_\mathrm{B}T<8`$, as well as of $`M`$ and $`a_{\mathrm{gel}}`$ we have examined. It is observed, however, to reduce about by a half when the excluded volume effect between beads is discarded. This enlargement of the ‘persistence length’ is one of the important consequences of the excluded volume effect on the present BD chain in a 3D space.

## IV ‘Defects’ Dynamics under CFGE

Under relatively strong fields DNA undergo the quasi-periodic behavior Oana et al. (1994): real time images of individual DNA exhibit elongated and contracted shapes alternatively. Corresponding behavior was observed in the previous 2D BD models Deutsch (1988); Matsumoto and Doi (1994); Masubuchi et al. (1995) likewise by our present BD simulation which has been discussed in details in I. We show in Fig. 3, a normalized curve of $`R_l(t)`$, i.e., $`\overline{R}_l(tt_{\mathrm{max}})/\overline{R}_l(t_{\mathrm{max}})`$ of chains with $`M=160`$ in a lattice of gel with the lattice distance of $`a_{\mathrm{gel}}=20`$ under a field $`EqE_\mathrm{b}=0.032`$. In order to get further insights into the mean elongation-contraction motion through the defect dynamics, we have examined histograms of length of the extended parts $`N_n(t)`$ and of their tangential movement $`N_{\mathrm{\Delta }u}(t)`$ at four representative time intervals, i.e., $`[t_{\mathrm{min}}0.1D,t_{\mathrm{min}}]`$, $`[t_{\mathrm{min}},t_{\mathrm{min}}+0.1D]`$, $`[t_{\mathrm{max}}0.1D,t_{\mathrm{max}}]`$, and $`[t_{\mathrm{max}},t_{\mathrm{max}}+0.1D]`$, as indicated in Fig. 3. Here $`t_{\mathrm{min}}\overline{t_{\mathrm{min}:\mathrm{l}}^nt_{\mathrm{max}}^n}`$ and $`D\overline{t_{\mathrm{min}:\mathrm{r}}^nt_{\mathrm{min}:\mathrm{l}}^n}`$, where $`t_{\mathrm{max}}^n`$ is the time of local maximum in the $`n`$-th peak in $`R_l(t)`$, $`t_{\mathrm{min}:\mathrm{r}}^n`$ and $`t_{\mathrm{min}:\mathrm{l}}^n`$ are respectively local minima just before and after the local maximum at $`t_{\mathrm{max}}^n`$, and the overline stands for the average over peaks in $`R_l(t)`$. As for the procedure how to pick up the peaks in $`R_l(t)`$, the reader may refer to I. We denote the four intervals introduced here as Interval 1, 2, 3 and 4, respectively. It is expected that Interval 1 corresponds to the time range where the chain is in a coiled shape, Interval 2 to that of a contracted shape when a U-shape conformation has just started to grow due to trapping by gel, Interval 3 to that of an elongated U-shape when the chain is going to get rid of the trap, and Interval 4 to that of an I-shape when it is going to retrieve another coiled shape. We cut out sequential $`2000`$ MD steps (mds) in each of the four intervals of each peak in $`R_l(t)`$ and accumulate data associated with them over 80 peaks whose $`D`$ are about $`1.5\times 10^5`$mds or larger.

The normalized histograms $`r_n(t)`$ $`(=N_n(t)/_mN_m(t))`$ and $`r_{\mathrm{\Delta }u}(t)`$ thus obtained are shown in Figs. 4 and 5, respectively. They exhibit peculiar time evolution within the ‘elongation-contraction’ period. One may at once notice from $`r_n(t)`$ that the chain is pulled by the field in Intervals 2 and 3. Particularly in Interval 3, $`r_n(t)`$ at $`40<n<120`$ becomes about $`10^3`$ whose integrated weight reaches about $`10\%`$, i.e., among 10 extended parts one of them is as large as several tens. The latter is quite elongated and is considered to constitute the longer arm of the U-shaped chain. Interestingly, however, we can see from Fig. 5 that $`r_{\mathrm{\Delta }u}(t)`$ at $`\mathrm{\Delta }u`$ larger than ten is always negligibly small. This means that sliding dynamics of the chain even in Interval 3 consists of a large number of defect movements whose $`\mathrm{\Delta }u`$ are less than several times of the mean bead distance. Thus the dynamics even in the elongated portion of the chain is considered to be quite stochastic when it is looked at through $`\mathrm{\Delta }u`$, i.e., changes in a unit of time of our BD simulation. In this context, we note here that up to several hundreds multi-scattering processes between beads are involved in one MD step of our BD method as explained in I.

Another interesting observation is that both $`r_n(t)`$ and $`r_{\mathrm{\Delta }u}(t)`$ in Interval 4 is the most closer, among the four intervals, to a Poisson distribution which is observed in equilibrium under $`E=0`$. Namely, although $`R_l(t)`$ is still relatively larger in this interval just after the chain is released at $`tt_{\mathrm{max}}`$ from trapping by an obstacle(s), many s-beads have invaded into the extended part(s) and the chain contracts dominantly by the entropic effect. When $`R_l(t)`$ comes closer to its local minimum value due to trapping by a new obstacle(s), i.e., in Interval 1, $`r_n(t)`$ and $`r_{\mathrm{\Delta }u}(t)`$ at larger $`n`$ and $`\mathrm{\Delta }u`$ become bigger than those of Interval 4, indicating that some parts of the chain are already elongated by the field.

Now let us turn to Figs. 6 and 7 where $`r_n`$ and $`r_{\mathrm{\Delta }u}`$ averaged over the whole time window of the BD run under different fields are shown. Up to the field strength of $`E0.008`$ neither $`r_n`$ nor $`r_{\mathrm{\Delta }u}`$ significantly deviates from a Poisson distribution, while in $`r_n`$ under $`E>0.016`$ a tail clearly shows up at large $`n`$. This clearly indicates the existence of a crossover between dynamical regimes with and without the elongation-contraction motion at around $`E0.008`$ with $`M=160`$ and $`a_{\mathrm{gel}}=20`$. Even in the regime with the elongation-contraction motion where $`r_n`$ exhibits a significant deviation from a Poisson distribution, $`\mathrm{\Delta }u`$ remains in an exponential distribution, though its width moderately increases with increasing $`E`$. This result again indicates that, even in the deterministic regime where the time evolution of averaged $`R_l(t)`$ or $`v_\mathrm{G}(t)`$ looks apparently deterministic, defect dynamics is still purely stochastic.

## V Extended Bond Fluctuation Method

The stochastic nature of ‘defect’ dynamics found by the BD method in the present work was just the basic assumption a priori introduced in the extended bond fluctuation method (EBFM) when we developed it. Azuma and Takayama (1999) The conventional BFM (CBFM) Carmesin and Kremer (1988) is a well-established MC method which simulate thermodynamic properties of polymers by a coarse-grained chain model on a lattice. It consists of a set of local movements which suffice the conditions necessary for the chain to evolve by self-avoiding random walk. In addition to this local updating process we introduced non-local movements of ‘defects’ in the EBFM. For this purpose ‘defects’, or ‘s-monomers’ are defined in almost the same way as s-beads are defined in the present work. We then let s-monomers hop non-locally within an extended portion of the chain that each s-monomer lies. The movements are incorporated in such a way that they fulfill the detailed valance and the self-avoiding conditions. In the EBFM this non-local movement of s-monomers and the local updating of the CBFM are tried alternatively. We note that the number of s-monomers does not change except at the chain ends by the former process, but it does change anywhere in the chain by the latter updating.

In Fig. 8 we show $`r_n`$ and $`r_{\mathrm{\Delta }u}`$ of a chain with $`M=80`$ obtained in the EBFM under $`\mathrm{\Theta }=0`$, where we denote the normalized field as $`\mathrm{\Theta }qE_\mathrm{b}/k_\mathrm{B}T`$ for the MC analysis. As is the case for the BD model, both $`r_n`$ and $`r_{\mathrm{\Delta }u}`$ obtained by the EBFM obey a Poisson distribution except for regions $`n,\mathrm{\Delta }u<5`$. Furthermore, in the EBFM, their widths coincide with each other at large $`n`$ and $`\mathrm{\Delta }u`$. In this context, it is noted that we have adopted equal probability of semi-local movement of an s-monomer by any distance $`k`$ among possible $`n`$ ones. The width of the distribution $`r_n`$ of the EBFM is about a half of that of the BD method. Correspondingly, ‘defect’ density $`\rho (0.33)`$ of the EBFM is about twice larger than that of the BD method. Although there exist such quantitative differences between the two methods, the qualitative agreement in the results demonstrated in Figs. 1 and 8 supports our EBFM approach in Ref. Azuma and Takayama, 1999.

The elongation-contraction motion under CFGE is also observed by the EBFM. The normalized histograms $`r_n(t)`$ and $`r_{\mathrm{\Delta }u}(t)`$ obtained in the four intervals introduced in Sec. IV are shown in Figs. 9 and 10. Their overall features qualitatively agree with those obtained by the BD calculation. By a detailed inspection of Fig. 9, however, we see that $`r_n(t)`$ in Intervals 2 and 3 by the EBFM still exhibit nearly exponential distributions at larger $`n`$ in contrast to those of the BD results which have an enhanced tail at $`n`$ as large as several tens. This difference is considered to originate from the following process of the EBFM: s-monomers are created or annihilated, and so length $`n`$ of extended parts are changed, only through local updating of monomers originally incorporated by the CBFM. By such local stochastic processes an extremely long extended part is hardly created.

The normalized histogram $`r_{\mathrm{\Delta }u}(t)`$ shown in Fig. 10, on the other hand, exhibit exponential behavior at $`\mathrm{\Delta }u>2`$ quite similar to the corresponding BD result shown in Fig. 5. This is also the case for the integrated distributions $`r_{\mathrm{\Delta }u}`$ under various field strengths. They are shown in Fig. 11 which one can compare with the corresponding BD result in Fig. 7. Quantitatively, $`r_{\mathrm{\Delta }u}(t)`$ and $`r_{\mathrm{\Delta }u}`$ by the EBFM have significant weight at relatively larger $`\mathrm{\Delta }u`$ as compared with those by the BD calculation. This difference is attributed to the process in the EBFM that, once a long extended part of length $`n`$ is created, a tangential move of the s-monomer by the distance $`n`$ is accepted by the probability proportional to $`1/n`$. Although there exist such minor differences between the dynamical processes involved in the EBFM and the BD method, the results demonstrated above are sufficient for us to conclude that the EBFM properly takes into account the stochastic dynamics of slack segments which give rise to the elongation-contraction motion of a polymer under gel electrophoresis.

## VI Conclusion

In equilibrium under a vanishing field it is found, as expected, that s-beads, or ‘defects’, are distributed randomly. Both distributions of extended parts $`r_n(t)`$ and of their tangential movements $`r_{\mathrm{\Delta }u}(t)`$ obey a Poisson distribution. In the relatively strong field region where a chain exhibits the elongation-contraction motions giving rise to quasi-periodic behavior in $`R_l(t)`$ and $`v_\mathrm{G}(t)`$, $`r_n(t)`$ is observed to show the corresponding time evolution. Especially, immediately before $`R_l(t)`$ reaches to a local maximum, chain conformations with extended parts of more than a few tens of beads, which are entropically unfavorable, are observed. As for $`r_{\mathrm{\Delta }u}(t)`$, such a large fluctuation seems not to occur. We therefore conclude that the s-bead dynamics at a semi-microscopic level is rather stochastic, though the coarse-grained quantities such as $`R_l(t)`$ look apparently deterministic during the elongation-contraction motion, or in time ranges which we have called deterministic ones. We have also shown that these features of $`r_n(t)`$, $`r_{\mathrm{\Delta }u}(t)`$ and $`R_l(t)`$ are similar to those observed by the MC analysis based on the EBFM which we recently proposed. Such coincidences between the two simulations strongly support our algorithm of EBFM which requires an order of magnitude less CPU times than MD calculations based on the present BD method do.

###### Acknowledgements.

The author wish to thank H. Takayama for useful discussions and suggestions. The computation in the present work has been done using the facilities of the Supercomputer Center, Institute for Solid State Physics, University of Tokyo, and those of the Computer Center of University of Tokyo.

## References

FIG. 1

FIG. 2

FIG. 3

FIG. 4

FIG. 5

FIG. 6

FIG. 7

FIG. 8

FIG. 9

FIG. 10

FIG. 11


format: text

text:
# Large limit-from𝐵-Fields and Noncommutative SolitonsThis work is dedicated to my parents Juan Moreno and Martha Soto.

## INTRODUCTION

Type IIB Superstring Theory allows two interesting ingredients: stable-BPS RR-charged $`Dp`$-branes ($`p`$ odd) and a massless antisymmetric $`B_{\mu \nu }`$ field. Although a $`D3D3`$ brane system is stable; a $`D3\overline{D3}`$ case is not, due to the presence of a tachyon in its spectrumkey-1 . The job of the $`B`$-field in the low energy limit is to introduce noncommutativitykey-4 . In the $`D3\overline{D3}`$ brane configuration, we may turn on a large $`B`$-field along two worldvolume spatial coordinates. The effect of this is that the complex tachyon allows a GMS solitonSeiberg . This object appears to be superconducting in the presence of fermions arising from the open string sectorMoreno .

## $`B`$-FIELDS AND NONCOMMUTATIVITY

Consider an open string attached to a $`D`$-brane. The OPE has the form

$$e^{ik_1X}e^{ik_2X}\left(\tau \tau ^{}\right)^{2\alpha ^{}g^{\mu \nu }k_{1\mu }k_{2\nu }}\times e^{i\left(k_1+k_2\right)X}+\mathrm{}.$$
(1)
However, the introduction of a large $`B`$-field alters this to

$$e^{ik_1X}e^{ik_2X}\left(\tau \tau ^{}\right)^{2\alpha ^{}G^{\mu \nu }k_{1\mu }k_{2\nu }}\times \left[e^{\frac{i}{2}\mathrm{\Theta }^{\mu \nu }k_{1\mu }k_{2\nu }}\right]\times e^{i\left(k_1+k_2\right)X}+\mathrm{},$$
(2)
where $`G^{\mu \nu }=\left(\frac{1}{g+2\pi \alpha ^{}B}g\frac{1}{g2\pi \alpha ^{}B}\right)^{\mu \nu }`$ is the effective metric seen by the open string modes, and $`\mathrm{\Theta }^{\mu \nu }=\left(2\pi \alpha ^{}\right)^2\left(\frac{1}{g+2\pi \alpha ^{}B}B\frac{1}{g2\pi \alpha ^{}B}\right)^{\mu \nu }`$ is the noncommutativity parameter matrixSeiberg . Likewise, the new term in brackets is known -in configuration space- as the Moyal $``$ product. Thus, in the low energy effective theory, fields get $``$-multiplied.

## NONCOMMUTATIVE SOLITONS

The idea of GMS solitons was cleverly used in Harvey  to construct real solutions to the tachyon in bosonic $`D`$-branes. In this work, complex tachyons in RR $`D`$-brane pairs are considered instead. For a $`D3\overline{D3}`$ brane pair in the presence of a large $`B`$-field along the $`xy`$ plane, the action of the tachyon is

$$S^{\left(\mathrm{\Sigma }_{3+1}\right)}=_{\mathrm{\Sigma }_{3+1}}𝑑x𝑑y𝑑z𝑑t\left[\overline{D^\mu T}D_\mu TV_{}(T,\overline{T})\right],$$
(3)
where $`[x,y]=i\theta ,`$ and $`z`$ and $`t`$ are commutative coordinates. Also, $`D_\mu T=_\mu TiA_\mu T+i\stackrel{~}{A_\mu }T,`$ where $`A_\mu `$ and $`\stackrel{~}{A_\mu }`$ are the respective gauge fields in each of the $`D`$-brane’s Chan-Paton $`U\left(1\right)`$ symmetry groups.

We make three assumptions:

* The potential is a polynomial:

$$V_{}(T,\overline{T})=\underset{k=1}{\overset{n}{}}a_k\left(\overline{T}T\right)^k.$$
(4)
* The gauge field $`R_\mu A_\mu \stackrel{~}{A_\mu }`$ has the form:

$$R_\mu =R_\mu (z,t).$$
(5)
* We’ll focus on solutions of the form:

$$T=T(x,y),\overline{T}=\overline{T}(x,y).$$
(6)

Eventually, it is shown that the solution is

$$T=t_{}T_0,\overline{T}=\overline{t_{}}\overline{T_0},$$
(7)
where $`t_{}`$ and $`\overline{t_{}}`$ solve the equations of motion in the commutative case, while $`T_0,\overline{T_0}`$ satisfy

$$\left(\overline{T_0}T_0\right)^k=\overline{T_0}T_0.$$
(8)
Since $`[x,y]=i\theta `$ is analogous to $`[\widehat{q},\widehat{p}]=i`$ in quantum mechanics, we identify $`T_0i|00|`$ and $`\overline{T_0}i|00|`$ in the Simple Harmonic Oscillator basis.

Applying the Weyl-Wigner-Moyal correspondence yields the following lowest energy solitonic solution:

$$T(x,y)=2it_{}e^{r^2},\overline{T}(x,y)=2i\overline{t_{}}e^{r^2},$$
(9)
where $`r^2=x^2+y^2`$Moreno .

## THE NONCOMMUTATIVE SUPERCONDUCTING STRING

Given a $`\left(3+1\right)`$ Dirac spinor $`\mathrm{\Psi }`$, with two-component entries $`\psi _R`$ and $`\psi _L`$ obeying $`\stackrel{}{\sigma }\widehat{p}\psi _R=\psi _R`$ and $`\stackrel{}{\sigma }\widehat{p}\psi _L=\psi _L`$, the action for the spinor coupled to the complex soliton has the form:

$$S^{\left(\mathrm{\Sigma }_{3+1}\right)}=_{\mathrm{\Sigma }_{3+1}}𝑑t𝑑z𝑑x𝑑y\left[f\left(\overline{T}\right)\overline{\mathrm{\Psi }}g\left(T\right)D/\mathrm{\Psi }\right],$$
(10)
where $`f`$ and $`g`$ are polynomials, and $`D/\mathrm{\Psi }=\gamma ^\mu \left(_\mu \mathrm{\Psi }iR_\mu \mathrm{\Psi }\right).`$

In terms of operators, we may express our spinors asPilo :

$$\widehat{\psi }_{L,R}\left(x^\mu \right)=\underset{m,n0}{}\psi _{mn}^{L,R}(z,t)|mn|.$$
(11)
In order to find the effective theory along the *string* (the $`z,t`$ coordinates,) we make two assumptions:

* $`\theta \mathrm{}`$, which means that the noncommutative kinetic part is negligible.
* $`\psi \psi ^L`$ (as in Witten’s superconducting string Witten ).

Therefore, the action for the *Noncommutative $`D`$-string* is

$$S^{\left(\mathrm{\Sigma }_{1+1}\right)}=2\pi i\theta f\left(\overline{t_{}}\right)g\left(t_{}\right)_{\mathrm{\Sigma }_{1+1}}\left[i\overline{\psi }^L\sigma ^aD_a\psi ^L\overline{\psi }^Lm\psi ^L\right],$$
(12)
where $`\psi ^L`$ denotes $`\psi _{00}^L`$ in (11), $`m`$ is a “mass” matrix ($`m(z,t)=\sigma ^\alpha R_\alpha `$)111$`a=z,t`$ and $`\alpha =x,y`$. .

## SUPERCONDUCTIVITY

In the massless case (after rescaling the action and getting rid of the unnecessary $`L`$ subscript)

$$S^{\left(\mathrm{\Sigma }_{3+1}\right)}=_{\mathrm{\Sigma }_{3+1}}𝑑z𝑑t\left(i\overline{\psi }\sigma ^aD_a\psi \right).$$
(13)
According to the bosonization technique Witten , in two dimensions we can equivalently describe the theory by either bosons or fermions.

This is done by introducing a scalar field $`\zeta (z,t)`$ such that

$$\overline{\psi }\sigma ^a\psi =\frac{1}{\sqrt{\pi }}\epsilon ^{ab}_b\zeta .$$
(14)
It can be shown that the kinetic term corresponds to

$$i\overline{\psi }\sigma ^aD_a\psi =\frac{1}{2}\left(_a\zeta \right)\left(^a\zeta \right)\frac{1}{\sqrt{\pi }}R_a\epsilon ^{ab}_b\zeta ,$$
(15)
which is associated to a conserved current

$$J^a=^a\zeta +\frac{1}{\sqrt{\pi }}\epsilon ^{ab}R_b.$$
(16)
However, this current may be expressed in terms of another scalar:

$$J^a=\epsilon ^{ab}_b\phi .$$
(17)
Thus, $`_b\phi =\epsilon _{ba}J^a`$ and

$$^b_b\phi =^b\epsilon _{ba}\left(^a\zeta +\frac{1}{\sqrt{\pi }}\epsilon ^{ab}R_b\right)=\frac{1}{\sqrt{\pi }}^bR_b.$$
(18)
This means that from $`_a\phi =\epsilon _{ac}J^c,`$ we get that $`\epsilon _{ac}^aJ^c=\frac{1}{\sqrt{\pi }}^aR_a.`$ In other words, for $`J^3`$ (the current along the string):

$$\frac{dJ^3(z,t)}{dt}=\frac{1}{\sqrt{\pi }}\frac{dR_0(z,t)}{dz},$$
(19)
which has the following solution:

$$J^3(z,t)=\frac{1}{\sqrt{\pi }}\left[R_0(z,\tau )R_0(z,\tau _i)\right].$$
(20)
This means that the current is nondecaying as long as the gauge fields has different values when it’s turn on ($`t=\tau _i`$) and turn off ($`t=\tau `$).

## ACKNOWLEDGMENTS

It is a pleasure to thank my advisor Dr. H. García-Compeán for useful discussions and enjoyable collaboration. I’m very grateful to J.M.L. Bauche A. for helping me write this manuscript. I also thank the organizers of the *VIII Mexican Workshop on Particles and Fields* for their hospitality. This work was supported in part by the CONACyT Fellowship No. 33951E.


text:
# Optical Holonomic Quantum Computation*footnote **footnote *ICQI 2001 - International Conference on Quantum Information, Rochester, June 10-13, 2001 (Contribution given by D. E)

## Abstract

Holonomic quantum computation (HQC) is materialized here with quantum optics components. Holonomies are the generalization of the Berry phases to unitary matrices with dimensionality the same as the degree of degeneracy of the system. In a nonlinear Kerr medium the degenerate states of laser beams are interpreted as qubits. Control manipulations with displacers, squeezers and two-mode interfering devices performed in a cyclic, adiabatic fashion produce holonomies. Here, they are employed as logical gates for our HQC proposal. The effects of errors from imperfect control of classical parameters, the looping variation of which builds up holonomic gates, are investigated.

today

©1999 Optical Society of America

OCIS codes:(999.9999) Quantum Computing

Quantum computers consist of a tensor product structure of two level quantum sub-systems (qubits) which are manipulated coherently with a large set of external operations (logical gates). Here an application of holonomic manipulations  (see also ) with optical components of laser Fock states is described. Holonomy is the generalization of the well known Berry phase to the case of a higher dimensional unitary matrix . It is produced by an adiabatic cyclic evolution of external control parameters that drive a degenerate Hamiltonian. This Hamiltonian in our case is produced by a Kerr medium through which the lasers are passing . Hence, a two fold degeneracy is created, where each qubit is encoded. As the first part of this presentation shows the logical gates are constructed by employing displacings and squeezings of each mode for producing one-qubit gates. Two-qubit gates are produced by two mode interfering devices like two-mode displacing or squeezing. Each laser encodes one qubit in its degenerate sub-space, spanned by the Fock states $`|0`$ and $`|1`$. By alternative manipulations with the displacing operator $`D(\lambda )=\mathrm{exp}(\lambda a^{}\overline{\lambda }a)`$ and the squeezing operator $`S(\mu )=\mathrm{exp}(\mu a_{}^{}{}_{}{}^{2}\overline{\mu }a^2)`$ in an adiabatic fashion it is possible to generate any one-qubit rotation. Adiabaticity along the cyclic evolution of the parameters $`\lambda =x+iy`$ and $`\mu =r_1\mathrm{exp}i\theta _1`$ guarantees that the states $`|0`$ and $`|1`$ will not mix with the rest of the Fock states of the laser. On the other hand as they have the same energy they are allowed to exchange population along this cyclic evolution. The exact evolution they undertake at the end of this manipulation is given by unitaries which depend on the geometrical characteristics of the cyclic procedure . The coordinates of the space where the loops belong are not abstract theoretical objects but rather the parameters the experimentalist is tuning in the lab, like the squeezing or the displacing amplitude. For a loop $`C_I`$ on the $`(x,r_1)_{\theta _1=0}`$ plane we obtain the holonomy

$$\mathrm{\Gamma }_A(C_I)=\mathrm{exp}(i\sigma _1\mathrm{\Sigma }_I),$$
(1)
where $`\mathrm{\Sigma }_I:=_{D(C_I)}𝑑x𝑑r_12e^{2r_1}`$ . Hence it represents the area on a hyperbolic surface enclosed by the loop $`C_I`$. For this procedure displacings and squeezings are alternated. Similarly we obtain for a loop $`C_{II}`$ on $`(x,r_1)_{\theta _1=\pi /2}`$ the holonomy

$$\mathrm{\Gamma }(C_{II})=\mathrm{exp}(i\sigma _2\mathrm{\Sigma }_{II})$$
(2)
with $`\mathrm{\Sigma }_{II}:=_{D(C_{II})}2e^{2r_1}`$ .

For the case of two-qubit gates we use two laser beams passing through a Kerr medium. This produces a 4-dimensional degenerate space with a tensor product structure. Two-mode interfering devices produce two-mode displacing transformation $`N(\xi )=\mathrm{exp}(\xi a_1^{}a_2\overline{\xi }a_1a_2^{})`$ or two-mode squeezing $`M(\zeta )=\mathrm{exp}(\zeta a_1^{}a_2^{}\overline{\zeta }a_1a_2)`$. Interchanging these procedures in a cyclic, adiabatic fashion produces two-qubit unitaries. In particular, for $`\zeta =r_2e^{i\theta _2}`$ and $`\xi =r_3e^{i\theta _3}`$ we have a four-parameter real manifold. Taking a loop $`C_{III}(r_2,r_3)_{\theta _2=\theta _3=0}`$ we obtain

$$\mathrm{\Gamma }(C_{III})=\mathrm{exp}(i\sigma ^{12}\mathrm{\Sigma }_{III})$$
(3)
with $`\mathrm{\Sigma }_{III}:=_{D(C_{III})}𝑑r_2𝑑r_32\mathrm{sinh}(2r_2)`$. This is a nontrivial unitary transformation and hence, together with the general one-qubit holonomic transformations they can produce any desired logical gate resulting to universality. Moreover, the ability to produce a big variety of gates by simply changing the experimental parameters and spanning different loops reduces the amount of resources needed to construct certain gates in the traditional way by the successive combination of basic gates.

In the previously proposed setup the degeneracy has to be present at all times. To achieve this along with the implementation of the adiabatic transformations with the control devices it is necessary to employ the “kicking” method. The laser beams are inserted into the nonlinear Kerr medium for time $`\mathrm{\Delta }t`$ and alternated with infinitesimal transformation $`D(\delta \lambda )`$, $`S(\delta \mu )`$, $`M(\delta \zeta )`$ or $`N(\delta \xi )`$ which take place in times $`\delta t\mathrm{\Delta }t`$. This ensures the presence of the degeneracy along the adiabatic evolution.

Holonomic quantum computation is an all geometric setup of quantum evolution. The parameters of the gates are areas of particular surfaces. Deformations of these contours due to statistical experimental errors does not change the value of the areas to the first order of approximation. Moreover, the weight factors in the parameters $`\mathrm{\Sigma }_I`$ and $`\mathrm{\Sigma }_{II}`$ of the corresponding gates are exponentially decreasing for large $`r_1`$. This makes the gates resilient to displacing control errors as they are suppressed exponentially due to squeezing . Such characteristics are important for error avoiding computation.

Next we aim to discuss the error avoiding features of the optical holonomic proposal in the presence of systematic errors. More specifically the parametric loops producing holonomies, when subject to imperfections while spanned, introduce an error in the final gates through their accordingly fluctuated parameters. This is systematically studied for the Hadamard (H) gate and a non trivial two qubit gate. In order to study the errors introduced by imperfect control of the external parameters we adopt an imperfectly spanned loop, $`C^{}`$ . This by its turn will produce an area as a gate parameter given by $`\mathrm{\Sigma }^{}=\mathrm{\Sigma }+ϵ`$. If the errors are statistical rather than systematical then the area spanned by this loop are, to the first order, zero. Let us consider how systematic errors in the area effect one and two qubit gates. The Hadamard gate is given by

$$\begin{array}{cc}U_H=& \left[\begin{array}{ccc}\mathrm{cos}\mathrm{\Sigma }& \mathrm{sin}\mathrm{\Sigma }& \\ \mathrm{sin}\mathrm{\Sigma }& \mathrm{cos}\mathrm{\Sigma }& \end{array}\right],\end{array}$$
(4)
for $`\mathrm{\Sigma }=\pi /4`$. Up to a corrective phase it may be produced by a loop $`C_{II}`$, with spanning area given by $`\mathrm{\Sigma }=_{D(C_{II})}𝑑x𝑑r_12e^{2r_1}`$ with $`D(C_{II})`$ taken to be a rectangular surface enclosed by $`\{0x\pi /4,0r_1\mathrm{ln}2\}`$. Introduce an error in this surface by translating the boarders of $`x`$ and $`r_1`$ by $`\alpha `$ and $`\beta `$ respectively, where $`\alpha ,\beta 1`$. This is a kind of systematic error. The imperfect Hadamard gate is given to the first order in $`ϵ`$ by $`U_H(ϵ)=U_H+ϵh`$, with $`\begin{array}{cc}h=& \frac{1}{\sqrt{2}}\left[\begin{array}{ccc}1& 1& \\ 1& 1& \end{array}\right]\end{array}`$. For large values of $`r_1`$ the dependence of $`ϵ`$ on $`\beta `$ is exponentially small, while the error of $`x`$ is introduced linearly, i.e. $`ϵ=\alpha `$. As the squeezing of the light beam is much harder to perform and to control in experiments compared to its displacing the above feature is very appealing. Theoretically this characteristic makes the passage from the geometrical to the topological quantum computation as $`\beta `$ errors in the spanning of the loop $`C^{}`$ are suppressed even if they are very large.

For a two qubit gate we perform a loop $`C_{III}`$ between the variables $`r_2`$ and $`r_3`$ with area $`\mathrm{\Sigma }_{III}=\frac{\pi }{4}`$ which eventually gives the gate

$$\begin{array}{cc}U=& \frac{1}{\sqrt{2}}\left[\begin{array}{cccc}\sqrt{2}& 0& 0& 0\\ 0& 1& 1& 0\\ 0& 1& 1& 0\\ 0& 0& 0& \sqrt{2}\end{array}\right]\end{array}$$
(5)
Allowing the area $`\mathrm{\Sigma }_{II}`$ to be enclosed by the following rectangular $`\{0r_2\text{arccosh}(2),0r_3\pi /8\}`$ then systematic small errors in the definition of the boarders of the rectangular of the form $`(\alpha ^{},\beta ^{})`$ gives an error in the gate parameter of the form $`\mathrm{\Sigma }^{}=\mathrm{\Sigma }+\delta `$ with $`\delta =1.7\alpha ^{}+\beta ^{}`$. The gate in this case is given by $`U(\delta )=U+\delta u`$ where

$$\begin{array}{cc}u=& \frac{1}{\sqrt{2}}\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 1& 1& 0\\ 0& 1& 1& 0\\ 0& 0& 0& 0\end{array}\right]\end{array}$$
(6)
We notice that with the ordering $`\{|00,|10,|11,|01\}`$, adopted for the basis vectors that span the four dimensional degenerate subspace of coding the quantum information, where the first qubit is the control and the second is the target, the $`U`$ gate is a control $`\pi /4`$-rotation. Introducing the control phase gate $`P_\varphi =diag(1,e^{i\varphi },1,1)`$, we can construct the control-not gate with a control sign flip operation as $`U_{CN}=P_\pi U^2`$. Then the holonomic control not gate subjected to systematic errors is $`U_{CN}(\delta )=U_{CN}+\delta P_\pi u`$.

By means of the imperfect one and two qubit quantum gates so obtained by the optical model of holonomic computation , we can now proceed along the lines of , to provide an all optical physical implementation of the holonomic teleportation  and the holonomic remote gate construction , that have been put forward in , in the more mathematical framework of the $`CP^n`$ model.

Finally we should mention that alternatives to the optical implementation of holonomies for quantum computation have been suggested recently , in a scheme that uses adiabatic manipulations of parameters determining the interaction between laser beams and trapped ions in order to generate abelian and non-abelian geometrical gates. In such a model in addition to the systematic control errors presented here one has to take in account other imperfection mechanisms such as interaction with the environment and the effect of the spontaneous atomic emission upon Berry connection, see e.g. .


text:
# Anomalous diffusion in non-equilibrium relativistic heavy-ion rapidity spectra

## 1 Introduction

In relativistic heavy-ion collisions nuclear matter is believed to be compressed to high baryon density and nuclei undergoing central collisions reduce strongly their original longitudinal momentum. This loss of rapidity, usually referred as baryon stopping, is an important characteristic to understand the reaction mechanism at high energy density.

Isotropic thermal models fail to reproduce the broad and flat rapidity distribution observed in the central Pb-Pb collisions at $`158`$ A GeV/c . Collective and anisotropic flows must be considered (with the introduction of new parameters) to take into account the experimental spectra . Complicate evolution codes take only partially into account of such distributions and the physical evolution of the system is not completely understood .

In this paper we study the evolution of the rapidity distribution from a macroscopic point of view by using a non-linear relativistic Fokker-Planck equation and we show that the observed broad rapidity shape could be a signal of non-equilibrium properties of the system. Let us note that a similar approach, within a linear Fokker-Planck equation, has been previously studied in Ref. by using a linear drift in the space of the rapidity and a free parameter diffusion coefficient. With this choice, the author found a strongly violation of the fluctuation-dissipation theorem. We will see that generalizing the Brownian motion to the relativistic kinetic variables, the standard Einstein relation is satisfied and Tsallis non-extensive statistics emerges in a natural way from the non-linearity of the Fokker-Planck equation.

## 2 Non-linear Fokker-Planck equation in relativistic dynamics

A class of anomalous diffusions are currently described through the non-linear Fokker-Planck equation (NLFPE)

$$\frac{}{t}[f(y,t)]^\mu =\frac{}{y}\left[J(y)[f(y,t)]^\mu +D\frac{}{y}[f(t,y)]^\nu \right],$$
(1)
where $`D`$ and $`J`$ are the diffusion and drift coefficients, respectively. Tsallis and Bukman  have shown that, for linear drift, the time dependence solution of the above equation is a Tsallis-like distribution with $`q=1+\mu \nu `$. The norm of the distribution is conserved for all times only if $`\mu =1`$, therefore we will limit the discussion to the case $`\nu =2q`$.

We are going to study the time evolution of the particle distribution by means the above NLFPE in the rapidity space. Basic assumption of our analysis is that the rapidity distribution $`y`$ is not appreciably influenced by transverse dynamics which is considered in thermal equilibrium. Such hypothesis is well confirmed by the experimental data and adopted in many theoretical works . A crucial rôle in the solution of the above NLFPE plays the choice of the diffusion and the drift coefficients. Such a choice influences the time evolution of the system and its equilibrium distribution.

Since the temperature at freeze-out exceeds $`100`$ MeV, Boltzmann approximation is usually adopted and the single particle equilibrium distribution is written as

$`E{\displaystyle \frac{d^3N}{d^3p}}E\mathrm{exp}(E/T)m_{}\mathrm{cosh}(y)\mathrm{exp}(m_{}\mathrm{cosh}(y)/T),`$ (2)
where $`y`$ is the rapidity, $`m_{}=\sqrt{m^2+p_{}^2}`$ is the transverse mass, $`T`$ is the temperature. Imposing the validity of the Einstein relation for Brownian particle, we can generalize to the relativistic case the standard expressions of diffusion and drift coefficients as follows

$$D=\gamma T,J(y)=\gamma m_{}\mathrm{sinh}(y)\gamma p_{},$$
(3)
where $`p_{}`$ is the longitudinal momentum and $`\tau `$ is a common constant. It is easy to see that the above coefficients give us the Boltzmann stationary distribution (2) in the linear Fokker-Planck equation ($`q=\nu =1`$)111This result cannot be obtained if one assumes a linear drift coefficient as in Ref.., while the stationary solution of the NLFPE (1) with $`\nu =2q`$ is a Tsallis-like distribution with the relativistic energy $`E=m_{}\mathrm{cosh}(y)`$. In the light of the above results, we generalize the results of Tsallis-Bukman by searching as a solution of the NLFPE (1) the following time dependent Tsallis-like distribution

$$f_q(y,m_{},t)=\left\{1(1q)\beta (t)m_{}\mathrm{cosh}[yy_m(t)]\right\}^{1/(1q)}.$$
(4)
The unknown functions $`\beta (t)=1/T(t)`$ and $`y_m(t)`$ have been derived by means numerical integration of Eq.(1) with initial $`\delta `$-function condition depending on the value of the experimental projectile rapidities. The rapidity distribution at fixed time is then obtained by numerical integration over the transverse mass $`m_{}`$ (or transverse momentum) as follows

$$\frac{dN}{dy}(y,t)=c_m^{\mathrm{}}m_{}^2\mathrm{cosh}(y)f_q(y,m_{},t)𝑑m_{},$$
(5)
where $`c`$ is the normalization constant fixed by the total number of the particles. Rapidity spectra calculated from (5) will ultimately depend on two parameters: the “interaction” time $`\tau =\gamma t`$ and the non-extensive parameter $`q`$. Therefore, no more free parameters are used in this analysis respect to previous theoretical studies .

## 3 Results and Conclusions

We have studied the time dependent rapidity distribution (5) in the case of central Pb+Pb collisions at 158 GeV/c and our results are compared with the net proton experimental NA49 data . The obtained spectra are normalized to 164 protons and the beam rapidity is fixed to $`y_{\mathrm{cm}}=2.9`$ (in the c.m. frame) . In Fig.1 we show the calculated rapidity spectra compared with the experimental data. The full line corresponds to the NLFPE solution (5) at $`\tau =0.82`$ and $`q=1.25`$; the dashed line corresponds to the solution of the linear case ($`q=1`$) at $`\tau =1.2`$. Only in the non-linear case ($`q1`$) exists a (finite) time for which the obtained rapidity spectra well reproduces the broad experimental shape. A value of $`q1`$ implies anomalous superdiffussion in the rapidity space, i.e., $`[y(t)y_M(t)]^2`$ scale like $`t^\alpha `$ with $`\alpha >1`$ .

We have studied the rapidity spectra from a macroscopic point of view. A good agreement with the experimental data is reached assuming a deviation from the equilibrium and a non-linear evolution which implies non-extensive Tsallis-like distribution. Non-extensive features in relativistic collisions are outlined in several works , even if a microscopic justification of these effects is still lacking. In Ref., the authors show that diquark breaking during inelastic collisions in dual parton model could be responsible for the observed baryon stopping. The possible connection between the nature of these inelastic collisions and the above observed non-extensive statistical effects is under investigations.

Acknowledgments. It is a pleasure to thank Prof. P. Quarati and Prof. G. Wilk for fruitful discussions and comments.


text:
# The gluon plasma at RHIC

## 1 Introduction

At present we have a very limited understanding of the properties of the partonic environment created in heavy-ion collisions at RHIC. Predictions for the density of the gluons produced vary by a factor of five depending on the model considered. Elliptic flow, $`v_2(p_{})=\mathrm{cos}(2\varphi )_p_{}`$, the differential second moment of the azimuthal momentum distribution, and the high-$`p_{}`$ suppression of the particle spectra have been the subject of increasing interest because they provide important constraints on the density and effective energy loss of partons.

The origin of the remarkable saturation of elliptic flow $`v_2(p_{})0.2`$ above $`p_{}2`$ GeV reported by STAR at Quark Matter 2001 is an open question for theory. Calculations based on inelastic parton energy loss do predict saturation or decreasing $`v_2`$ at high $`p_{}`$. These calculations are valid for high $`p_{}`$, where collective transverse flow from lower-$`p_{}`$ partons can be neglected and Eikonal dynamics is applicable. However, a constant spatial anisotropy was assumed throughout the evolution, while in reality, it decreases and probably even changes sign. This is likely to reduce the generated elliptic flow much below the preliminary data.

Ideal hydrodynamics, the simplest theoretical framework to study elliptic flow, agrees remarkably well with the measured elliptic flow data up to transverse momenta $`1.5`$ GeV$`/c`$. However, it fails to saturate at high $`p_{}>2`$ GeV as does the preliminary data.

A theoretical problem with ideal hydrodynamics is that it assumes local equilibrium throughout the whole evolution. This idealization is marginal for conditions encountered in heavy ion collisions. Covariant Boltzmann transport theory provides a convenient framework for nonequilibrium dynamics that depends on the local mean free path $`\lambda (x)1/\sigma n(x)`$.

Parton cascade simulations show on the other hand, that the initial parton density based on HIJING is too low to produce the observed elliptic flow unless the pQCD cross sections are artificially enhanced by a factor $`23`$. However, gluon saturation models predict up to five times higher initial densities, and these may be dense enough to generate the observed collective flow even with pQCD elastic cross sections.

In this study, we explore the dependence of elliptic flow and the high-$`p_{}`$ suppression of the particle spectra on the initial density and the elastic $`gg`$ cross section. Though parton cascades lack at present covariant inelastic energy loss, elastic energy loss alone may account for the observed high-$`p_{}`$ azimuthal flow pattern as long as the number of elastic collisions is large enough.

## 2 Covariant parton transport theory

We consider here, as in , the simplest nonlinear form of Lorentz-covariant Boltzmann transport theory in which the on-shell phase space density, evolves with an elastic $`22`$ rate. We solve the transport equation via the MPC algorithm, which maintains Lorentz covariance using the parton subdivision technique. See Ref.  and references therein for details.

For a given nuclear geometry and formation time, the solutions of the nonlinear transport equation has been shown to depend mainly on the transport opacity $`\chi 𝑑z\sigma _t\rho _g=N\mathrm{sin}^2\theta _{cm}`$ and the impact parameter $`b`$. Here $`\sigma _t`$ is the elastic transport cross section, $`N`$ is the average number of collisions per parton during the whole evolution, while $`\theta _{cm}`$ is the collision deflection angle in the c.m. frame. To good accuracy the transport opacity factorizes as $`\chi =C(b)\sigma _t(T_0)dN_g/d\eta `$.

We label our results by the transport opacity $`\chi `$ and impact parameter $`b`$. Furthermore, we quote our transport opacities relative to that for the pQCD minijet gluons predicted by HIJING \[$`dN/d\eta =210`$, $`\sigma _{gggg}3`$ mb $`\sigma _t1`$ mb, $`\chi _{b=0}0.3`$\].

## 3 Numerical results

The initial condition was a longitudinally boost invariant thermal Bjorken tube at proper time $`\tau _0=0.1`$ fm/$`c`$ with uniform pseudo-rapidity distribution between $`|\eta |<5`$ and transverse density distribution proportional to the binary collision distribution for the two gold nuclei. Based on HIJING, the pQCD jet cross section was normalized to yield $`dN_g/d\eta =210`$ in central collisions and the initial temperature was chosen to be $`T_0=700`$ MeV.

Two different hadronization schemes were applied. One is based on local parton-hadron duality, where each gluon is assumed to convert to a pion. The other hadronization prescription is independent fragmentation, where we considered only the $`g\pi ^\pm `$ channel. See Ref.  for details.

Fig. 1 shows the impact-parameter-averaged elliptic flow as a function of $`p_{}`$. With increasing $`p_{}`$, elliptic flow increases until $`p_{}1.52`$ GeV, where it saturates. To reproduce the preliminary STAR data, $`80`$ times more opaque initial gluon plasma is needed than the pQCD prediction from HIJING. Surprisingly, the results show no sensitivity to the applied hadronization prescription.

Fig. 2 shows that such large transport opacities are also consistent with the preliminary charged hadron spectra measured by STAR. Hadronization via parton-hadron duality yields too little suppression at high $`p_{}`$ because it only incorporates quenching due to elastic energy loss. However, with the additional quenching due to independent fragmentation, the parton cascade results approach the preliminary STAR data.

## Acknowledgements

We acknowledge the Parallel Distributed Systems Facility at the National Energy Research Scientific Computing Center for providing computing resources.

This work was supported by the Director, Office of Energy Research, Division of Nuclear Physics of the Office of High Energy and Nuclear Physics of the U.S. Department of Energy under contract No. DE-FG-02-93ER-40764.


text:
# Fractal extra dimension in Kaluza-Klein theory

\[

## Abstract

Kaluza-Klein theory in which the geometry of an additional dimension is fractal has been considered. In such a theory the mass of an elementary electric charge appears to be many orders of magnitude smaller than the Planck mass, and the ”tower” of masses which correspond to higher integer charges becomes aperiodic.

PACS no.: 04.50.+h, 11.10.Kk

\]

In modern Kaluza-Klein theories  the extra N-4 space-time dimensions are considered to be compact and small (with characteristic size on the order of the Planck length). The symmetries of this internal space are chosen to be the gauge symmetries of some gauge theory , so a unified theory would contain gravity together with the other observed fields. In the original form of the theory a five-dimensional space-time was introduced where the four dimensions $`x^1,\mathrm{},x^4`$ were identified with the observed space-time. The associated 10 components of the metric tensor $`g_{\alpha \beta }`$ were used to describe gravity. After a compactified fifth dimension $`x^5`$ with a small circumference $`L`$ was added, the extra four metric components $`g_{\alpha 5}`$ connecting $`x^5`$ to $`x^1,\mathrm{},x^4`$ gave four extra degrees of freedom which were interpreted as the electromagnetic potential (here and everywhere else in the paper we use the following convention for greek and latin indices: $`\alpha =1,\mathrm{},4`$; $`i=1,\mathrm{},5`$). An additional scalar field $`g_{55}`$ or dilaton may be either set to a constant, or allowed to vary.

When a quantum field $`\psi `$ coupled to this metric via an equation

$$\mathrm{}_5\psi +a\psi =0$$
(1)
is considered, where $`\mathrm{}_5`$ is the covariant five-dimensional d’Alembert operator, the solutions for the field $`\psi `$ must be periodic in the $`x^5`$ coordinate. This leads to the appearance of an infinite ”tower” of solutions with quantized $`x^5`$-component of the momentum:

$$q_n^5=2\pi n/L$$
(2)
where n is an integer. In our four-dimensional space-time on a large scale such solutions with $`n0`$ interact with the electromagnetic potential $`g_{\alpha 5}`$ as charged particles with an electric charge $`e_n`$ and mass $`m_n`$:

$$e_n=\mathrm{}q_n(16\pi G)^{1/2}/c$$
(3)
$$m_n=\mathrm{}(q_n^2a)^{1/2}/c$$
(4)
where $`G`$ is the gravitational constant (see for example the derivation in ). In this theory the conservation of the electric charge is a simple consequence of the conservation of the $`x^5`$-component of the momentum. Although a unified description of gravity and electromagnetism has been achieved, the result obtained for the mass of an elementary electric charge is unsatisfactory large if $`L`$ is of the order of the Planck length, and the constant $`a`$ is small. Fine tuning of $`a`$ to 20 decimal places may bring the mass to a reasonably low value, but such an arbitrary fine tuning is hardly a satisfactory solution to the problem.

Similar problems remain in all the modern versions of Kaluza-Klein theory, which introduce up to a dozen additional space-time dimensions and sometimes assume a complicated topology of the ”internal space”. Despite this theoretical diversity, a common feature of all the Kaluza-Klein theories introduced so far is the integer number of additional dimensions used. But what if the topological properties of the ”internal space” are substantially different from the properties of our familiar four-dimensional space-time on a large scale? After all, we make a huge leap of faith in assuming that space-time on the Planck scale remains ”smooth” and differentiable. It is reasonable at least to consider the possibility that on the Planck scale space-time experiences such substantial quantum fluctuations that it may be better described by fractal geometry with some non-integer dimension (some attempts to consider fractal space-time geometries may be found in the literature, see for example  and the references therein). What would be the consequences of the fractal extra dimensions in the Kaluza-Klein theory? The goal of this paper is to answer this question qualitatively in the most simple situation of a fractal extra dimension $`D`$ within the range $`0<D<2`$.

Fractal objects and fractal dimensions are very useful mathematical concepts. A typical problem where the fractal dimension arises naturally is an attempt to measure the perimeter of an island in the ocean. The result would depend on the resolution used in the measurements. The value of the perimeter measured on the large scale from an aerial map would be much smaller than the value obtained by walking along the beach with a ruler, when every tiny curve of the beach is measured. The fractal dimension $`D`$ is defined from the variation with resolution of the main fractal variable (a length $`L`$ of a fractal curve, an area of a fractal surface, etc.). If $`D_T`$ is the topological dimension ($`D_T=1`$ for a curve, $`D_T=2`$ for a surface), the scale dimension $`\delta =DD_T`$ is defined as

$$\delta =\frac{d(lnL)}{d(ln(l/\lambda ))},$$
(5)
where $`\lambda `$ is the resolution of the measurements. If $`\delta `$ is constant we obtain

$$L=L_0(l/\lambda )^\delta ,$$
(6)
where the length $`L_0`$ is measured when $`\lambda =l`$. Let us assume that the circumference of the internal $`x^5`$ coordinate obeys expression (6) on a sufficiently small resolution scale instead of been constant. As before, the internal space must look small from our four-dimensional macroscopic point of view, so we may assume $`L_0`$ to be of the order of the Planck length when measured on a scale compatible with the electron’s Compton length $`l_c=2\pi \mathrm{}/(m_ec)`$ or larger (here we are not concerned with the behavior of $`\delta `$ on macroscopic scales, so we may assume it to be constant or $`\delta 0`$ on the large scale). An immediate consequence of (6) in the Kaluza-Klein theory is the drastic reduction of mass of an elementary electric charge for $`\delta >0`$. Qualitatively this is clear from the following simple arguments.

If we start from the circumference of the internal space $`L_0`$ measured at large scale and small energies, we would obtain the expression (4) for the elementary charge mass which corresponds to $`q=2\pi /L_0`$. A particle with such a large mass would have the Compton wavelength of the order of the Planck length. It should ”see” quite a different circumference of the internal space. For $`\delta >0`$ the internal space circumference $`L`$ measured on such a small scale will be much larger than $`L_0`$. Thus, $`qL>>1`$ and such a solution can not be the ground state. The elementary electric charge solution must be obtained in a self-consistent manner, and will correspond to a $`x^5`$-component of the momentum, which is much smaller than $`\mathrm{}/L_0`$.

In order to find the ground state self-consistently we must find solutions of the equation $`q_nL(q_n)=2\pi n`$, which describes field solutions periodic in the $`x^5`$-coordinate, where $`L(q_n)`$ satisfies (6). We may assume that the measurement scale $`\lambda `$ in (6) corresponds to $`2\pi /q_n`$, since any constant factor may be included in the value of $`l`$. Thus, we obtain

$$q_nL_0(\frac{lq_n}{2\pi })^\delta =2\pi n$$
(7)
Assuming $`a=0`$ the spectrum of mass looks like

$$m_n=\frac{2\pi \mathrm{}}{c}(\frac{n}{L_0l^\delta })^{1/(1+\delta )}$$
(8)
This ”tower” of solutions is quite different from the ”tower” (4) obtained in the regular five-dimensional Kaluza-Klein theory. It is no longer periodic at $`a=0`$, and for $`\delta >0`$ the elementary charge mass $`m_1`$ is much smaller than the Planck mass: $`m_1<<\frac{2\pi \mathrm{}}{cL_0}`$. For example, in the case of $`\delta =1`$ we obtain

$$m_1=\frac{2\pi \mathrm{}}{c(L_0l)^{1/2}}.$$
(9)
Both these developments are good, since they bring the theoretical picture closer to physical reality.

The qualitative result obtained above finds support in the following more precise theoretical consideration. Since it is difficult to define a length interval along a fractal $`x^5`$-direction, let us introduce an angular coordinate $`\varphi ^5`$ varying within an interval from 0 to $`2\pi `$. Despite the fractal character of $`x^5`$ we may unambiguously define $`\varphi ^5`$ for each point $`A`$ as

$$\varphi ^5=2\pi \frac{x_0^5}{L_0},$$
(10)
where $`x_0^5`$ is the distance along the $`x^5`$-coordinate from point $`A`$ to some point $`A_0`$ (designated as a zero-point) measured on a large scale, and $`L_0`$ is the circumference along the $`x^5`$-direction measured on a large scale. Although all the coordinates are supposed to be treated equally on the small scale, we leave the question of the possible fractal nature of the other four coordinates of our common space-time out of consideration, and write the metric as

$$ds^2=g_{\alpha \beta }dx^\alpha dx^\beta +2g_{\alpha 5}dx^\alpha d\varphi ^5+g_{55}d\varphi ^5d\varphi ^5,$$
(11)
where $`g_{\alpha 5}`$ depend explicitly on $`x^\alpha `$ and the scale of measurements, and $`g_{55}`$ also depends explicitly on the measurements scale. This dependence on the scale accounts for the fractal nature of the $`x^5`$-coordinate. Here we are not interested in possible spatial dependence of $`g_{55}`$ and consider all $`g_{i5}`$ components to be independent of $`\varphi ^5`$. Scale dependence of the metric will be addressed later in the discussion. In any event, we consider all $`g_{i5}`$ to be small at all scales, so $`ds^2`$ remains reasonably well defined in our four-dimensional space-time.

An analog of equation (1) with $`a=0`$ for the quantum field $`\psi `$ in this metric is

$`{\displaystyle \frac{}{x^\alpha }}(g^{\alpha \beta }{\displaystyle \frac{\psi }{x^\beta }})+{\displaystyle \frac{}{x^\alpha }}(g^{\alpha 5}{\displaystyle \frac{\psi }{\varphi ^5}})+`$ (12)
$`{\displaystyle \frac{}{\varphi ^5}}(g^{5\alpha }{\displaystyle \frac{\psi }{x^\alpha }})+{\displaystyle \frac{}{\varphi ^5}}(g^{55}{\displaystyle \frac{\psi }{\varphi ^5}})=0`$ (13)
Since we assume that the $`g_{i5}`$ do not depend explicitly on $`\varphi ^5`$, we should not address the meaning of terms like $`g^{i5}/\varphi ^5`$, which would be ambiguous for a fractal $`x^5`$-coordinate. Thus, equation (12) remains well defined, and we may search for its solution in the usual form as $`\psi =\mathrm{\Psi }(x^\alpha )e^{iq\varphi ^5}`$, where periodicity in $`\varphi ^5`$ requires $`q_n=n`$. As a result, we obtain

$$\mathrm{}\psi q_n^2\frac{1g_{\alpha 5}g^{\alpha 5}}{g_{55}}\psi +2iq_ng^{\alpha 5}\frac{\psi }{x^\alpha }+iq_n\frac{g^{\alpha 5}}{x^\alpha }\psi =0$$
(14)
This is the same as the Klein-Gordon equation in the presence of an electromagnetic field: in four-dimensional space-time it describes a particle of mass

$$m=\frac{2\pi \mathrm{}q_n}{cg_{55}^{1/2}(q_n)}$$
(15)
which interacts with a vector field $`g^{\alpha 5}`$ through a quantized charge $`e_nq_n`$. If we identify $`g_{55}^{1/2}(q)`$ as the scale dependent circumference $`L(q)`$ of the internal space, equations (14) and (8) will be equivalent to each other, and we arrive at the same result for the elementary charge mass and the aperiodic ”tower” of solutions described qualitatively earlier.

Thus far the question of scale dependence of $`g_{i5}`$ components in (11) has been left without detailed consideration. In order for $`ds^2`$ to remain a four-dimensional scalar these components must be functions of other four-dimensional scalars. Thus, a natural choice of the scalar scale would be the Compton wavelength $`l_c=2\pi \mathrm{}/(mc)`$, where $`m`$ is the effective four-dimensional mass, obtained similar to (14). This choice is consistent with the choice of scale that led to (7). Thus, despite the limited and heuristic nature of our approach (we left the question of possible fractal nature of the other four space-time coordinates without any consideration) some self-consistency has been achieved.

The numerical value of mass of an elementary charge obtained from (8) or (9) remains quite large. If we assume $`l`$ to be equal to the electron’s Compton length equation (9) gives $`m_110^{11}m_e`$. This estimate may be reduced by a few orders of magnitude by selecting somewhat larger values of $`L_0`$ and $`l`$, or by increasing the scale dimension $`\delta `$, so the effect of fractal extra dimension may in principle show up at an energy scale of hundreds of TeV. Overall, addition of the fractal extra dimension appears to be an alternative way of introducing large extra dimensions which has become very popular recently .

In conclusion, we have considered a Kaluza-Klein theory where the geometry of an extra dimension is fractal. In such a theory the mass of an elementary electric charge appears to be many orders of magnitude smaller than the Planck mass, and the ”tower” of masses which correspond to higher integer charges becomes aperiodic.


text:
# Numerical studies of left-handed materials and arrays of split ring resonators.

## I Introduction

Very recently, a new area of the research, called left-handed materials (LHM) has been experimentally demonstrated by Smith et al.  based on the work of Pendry et al. . LHM are by definition composites, whose properties are not determined by the fundamental physical properties of their constituents but by the shape and distribution of specific patterns included in them. Thus, for certain patterns and distribution, the measured effective permittivity $`ϵ_{\mathrm{eff}}`$ and the effective permeability $`\mu _{\mathrm{eff}}`$ can be made to be less than zero. In such materials, the phase and group velocity of an electro-magnetic wave propagate in opposite directions giving rise to a number of novel properties . This behavior has been called “left-handedness”, a term first used by Veselago  over thirty years ago, to describe the fact that the electric field, magnetic intensity and propagation vector are related by a left-handed rule.

By combining a 2D array of split-ring resonators (SRRs) with a 2D array of wires, Smith et al.  demonstrated for the first time the existence of left-handed materials. Pendry et al.  has suggested that an array of SRRs give an effective $`\mu _{\mathrm{eff}}`$, which can be negative close to its resonance frequency. It is also well known  that an array of metallic wires behaves like a high-pass filter, which means that the effective dielectric constant is negative at low frequencies. Recently, Shelby et al.  demonstrated experimentally that the index of refraction $`n`$ is negative for a LHM. Negative refraction index was obtained analytically and also from numerically simulated data . Also, Pendry  has suggested that a LHM with negative $`n`$ can make a perfect lens.

Specific properties of LHM makes them interesting for physical and technological applications. While experimental preparation of the LHM structures is rather difficult, especially when isotropic structures are required, numerical simulations could predict how the transmission properties depends on various structural parameters of the system. It will be extremely difficult, if not impossible, to predict the transmission properties of such materials analytically. Mutual electro-magnetic interaction of neighboring SRRs and wires makes the problem even more difficult. Numerical simulations of various configurations of SRRs and of LHMs could be therefore very useful in searching of the direction of the technological development.

In this paper, we present systematic numerical results for the transmission properties of LHMs and SRRs. An improved version of the transfer-matrix method (TMM) is used. Transfer matrix was applied to problems of the transmission of the electro-magnetic (EM) waves through non-homogeneous media many years ago . It was also used in numerical simulations of the photonic band gap materials (for references see Ref. 15). TMM enables us to find a transmission and a reflection matrices from which the transmission, reflection and absorption could be obtained. The original numerical algorithm was described in Ref. 13. In our program we use a different algorithm which was originally developed for the calculation of the electronic conductance of disordered solids .

The paper is organized as follows: In Section II we describe briefly the structure. We concentrate on the structure displayed in Figure 1. In Section III we present and discuss our results. The dependence of the transmission of the LHM and SRR on the electrical permittivity of the metallic components of our structure is given in Section III A. In Section III B we present the dependence of the transmission of the LHM on the size of the unit cell and size of the metallic wires. In Section III C we show the dependence of the resonance frequency of SRR on the parameters of the SRR. Section III D deals with the dependence of the resonance frequency on the permittivity of the board and embedding media. In Section IV we summarize our results and give some conclusions. Finally, in the Appendix A we give detailed description of the transfer matrix method.

## II Structure of the LHM meta-material

Both in the experiment and in the numerical simulations, the left-handed meta-materials consist from an array of unit cells, each containing one SRR and one wire. Figure 1a shows a realization of the unit cell that we have simulated. The size of the unit cell $`L_x\times L_y\times L_z`$ and the size of SRR itself are of order of mm. Waves propagate along the $`z`$-direction. The SRR lies in the $`yz`$ plane, and the wire is parallel to the $`y`$ axis.

As we are interested mostly in the transmission properties of the left-handed meta-material, the configuration as presented in Figure 1a, should be considered as one-dimensional. Indeed, such meta-materials possesses the left-handed properties only for the electro-magnetic wave incoming in the $`z`$-direction and even then only for a given polarization. Two-dimensional structures have been realized in experiments , in which two SRRs have been positioned in each unit cell in two perpendicular planes. For such structures, left-handed transmission properties have been observed for waves coming from any direction in the $`xz`$ plane. No three - dimensional structure has been realized so far.

Figure 1b shows a single square SRR of the type used for our simulations and also for experiments . The structure of the SRR is defined by the following parameters: the ring thickness $`c`$, the radial gap $`d`$, the azimuthal gap $`g`$ and the inner diameter $`r`$. The size of the SRR is

$$w=4c+2d+r.$$
(1)
Another parameter is the thickness of the SRR itself (in the $`x`$-direction). This thickness is very small in the experiments ($`0.02`$ mm). We can not simulate such thin structures yet. In numerical simulations, we divide the unit cell into $`N_x\times N_y\times N_z`$ mesh points. For homogeneous discretization, used throughout this paper, the discretization defines the minimum unit length $`\delta =L_x/N_x`$. All length parameters are then given as integer of $`\delta `$. This holds also for the thickness of the SRR. Generally, the thickness of SRR used in our simulations is 0.25-0.33 mm. Although we do not expect that the thickness will considerably influence the electro-magnetic properties of the SRR, it still could cause small quantitative difference between our data and the experimental results.

## III Structural parameters

### A Metallic permittivity

The existence of LHM has been experimentally demonstrated  for structures that have resonance frequencies in the GHz region. In this frequency region, we do not know the exact values of electrical permittivity $`ϵ_\mathrm{m}`$ of the metal. We know that Im $`ϵ_\mathrm{m}`$ is very large, and/or the Re $`ϵ_\mathrm{m}`$ is large but negative. In our previous studies  we have found that the resonance frequency $`\nu _0`$ of the LHM depends only on the absolute value of $`ϵ_\mathrm{m}`$. In fact $`\nu _0`$ reaches the saturated value provided that $`|ϵ_\mathrm{m}|>10^4`$. Since we do not know the exact values of the metallic permittivity, we have studied the transmission of the LHM with different values of $`ϵ_\mathrm{m}`$. In the results presented in Figure 2, we choose $`ϵ_\mathrm{m}=1+i\mathrm{Im}ϵ_\mathrm{m}`$ with different values of Im $`ϵ_\mathrm{m}`$.  The last is proportional to $`\sigma (\omega )/\omega `$ . For simplicity, we neglect the $`\omega `$-dependence of Im $`ϵ_\mathrm{m}`$ and consider Im $`ϵ_\mathrm{m}`$ = 8000, 18000 and 38000 for the three cases presented in Figure 2. For each case, we present results of transmission for different number (1 to 10) of unit cells. Notice that the higher imaginary part of the metal the higher is the transmission. Also the losses due to the absorption are smaller, as can be seen from the decrease of the transmission peak as the length of the system increases. This result is consistent with the formula presented by Pendry et al.  for the effective permeability of the system

$$\mu _{\mathrm{eff}}=1\frac{F\nu ^2}{\nu ^2\nu _0^2+i\gamma \nu }$$
(2)
with the damping factor

$$2\pi \gamma =\frac{2L_x\rho }{r\mu _0}$$
(3)
and the resonance frequency

$$(2\pi \nu _0)^2=\frac{3L_xc_{\mathrm{light}}^2}{\pi \mathrm{ln}\left(\frac{2c}{d}\right)r^3}.$$
(4)
where $`\rho `$ is the resistance of the metal, $`L_x`$ is the size of the system along the $`x`$ axis, $`c_{\mathrm{light}}`$ is a the velocity of light in vacuum and parameters $`r`$, $`c`$ and $`d`$ characterize the structure of SRR. They are defined in Figure 1b. Notice that the damping term $`\gamma 0`$ as $`\sigma \mathrm{}`$. Since the Im $`ϵ_\mathrm{m}`$ is proportional to $`\sigma `$, $`\gamma `$ is inversely proportional to Im $`ϵ_\mathrm{m}`$. Our numerical results suggest that it is reasonable to expect that the LHM effect will be more pronounced in systems with higher conductivity.

In Figure 3, we present the frequency dependence of the transmission for SRRs with the same parameters as those in Figure 2. Notice that the transmission is more pronounced as the length of the system is increased. Note also that the resonance gap becomes narrower when Im $`ϵ_\mathrm{m}`$ increases. This is in agreement with Eq. (2). The frequency interval, in which the effective permeability is negative, becomes narrower when the damping factor $`\gamma `$ decreases.

In Figure 4 we show the transmission through the LHM, in which the SRR are turned around their axis by 90 degrees. If we keep the same size of the unit cell as that of Figure 2, we do not obtain any LHM peak in the transmission, although there is a very well defined gap for the SRR alone. It seems that for this orientation of the SRR there is no overlap of the field of the wire with that of the SRR. The results shown in Figures 4 and 5 are therefore obtained with a reduced unit cell of $`3.66\times 3.66\times 3.66`$ mm (the size of SRR is still $`3\times 3`$ mm). The LHM transmission peak is located close to the lower edge of the SRR gap (shown in Figure 5). This is in contrast to the results presented in Figures 2 and 3, where the LHM transmission peak is always located close to the upper edge of the SRR gap. Finally, the gap shown for the “turned” SRR shown in the Figure 5, is deeper and broader than the gap for the “up” SRR.

In fact, for the SRR “up” structure, we found that the transmission in the gap is always of order of $`10^710^8`$. We can explain this effect by non-zero transmission from the $`p`$ to $`s`$ polarized wave (and back). If the transmission $`t(ps)`$ and $`t(sp)>0`$, then there is always the non-zero probability $`t(ps)t(sp)`$ for the $`p`$-polarized wave to switch into the $`s`$ state, at the beginning of the sample, move throughout the sample as the $`s`$ wave (for which neither wires nor SRR are interesting), and in the last unit cell to switch back into the $`p`$-polarized state. This process contributes to the transmission probability $`T(pp)`$ of the whole sample and determines the bottom level of the transmission gap for SRR. We indeed found that $`t(sp)10^4`$ for the “up” SRR. In the “turned” SRR case, both $`t(ps)`$ and $`t(sp)`$ should be zero due to the symmetry of the unit cell. Our data give $`t(ps)t(sp)10^6`$ for the “turned” SRR which determines the decrease of the transmission in the gap below $`10^{11}`$.

### B Dependence on the size of the unit cell and the width of the metallic wire.

As we discussed in Section III A, the transmission peak for the LHM with cuts in the SRR in the horizontal direction appears only when the size of the unit cell is really small. The effect of the size of the unit cell was demonstrated already in Figures 2 and 3, where we compared the transmission for the “up” SRR and LHM of different size of the unit cell. Both the transmission gap for an array of SRRs and the transmission peak for LHM are broader for smaller unit cell.

In Figure 6, we show the transmission for the LHM structure with q unit cell of size $`5\times 3.66\times 5`$ mm for the “turned” SRR. Evidently, there is no transmission peak for all the system lengths studied. The size of the unit cell must be reduced considerably to obtain a transmission peak.

Figure 7 presents the transmission peak for various sizes of the unit cell. Resonance frequency decreases as the distance between the SRRs in the $`x`$ direction decreases. This agree qualitatively (although not quantitatively) with theoretical formula given by Eqn. (4). We see also that an increase of the distance between SRR in the $`z`$ direction while keeping the $`L_x`$ constant causes sharp decrease and narrowness of the transmission peak.

### C Resonance frequency of SRR

In this section we study how the structure of the SRR influences the position of the resonance gap. In order to simulate various forms of the SRR, we need to have as many as possible mesh points in the $`yz`$ plane. Keeping in mind the increase of the computer time when the number of mesh points increases, we used a unit cell with $`L_x<L_y,L_z`$. The actual size of the unit cell in this section is $`L_x=2.63`$ mm and $`L_y=L_z=6.05`$ mm. and we use uniform discretization with $`N_x\times N_y\times N_z=10\times 23\times 23`$ mesh points. This discretization defines a minimum unit length $`\delta =0.263`$ mm. SRR with size of $`5\times 5`$ mm, is divided into $`19\times 19`$ mesh points.

The electrical permittivity of the metallic components is chosen to be $`ϵ_\mathrm{m}=1000+10.000i`$. We expect that larger value of Im $`ϵ_\mathrm{m}`$ will increase a little the position of the resonance gap . However, the dependence on the different structural parameters will remain the same. Higher values of $`|ϵ_\mathrm{m}|`$, however, will require more CPU time because of shorter interval between the normalization of the transmitted waves (see Appendix A for details).

We have considered 23 different SRR structures and studied how the resonance frequency $`\nu _0`$ depends on the structure parameters. The lowest $`\nu _0=3.75`$ GHz was found for SRR with c : d : r : g = 1 : 1 : 13 : 1. On the other hand, SRR with c : d : r : g = 2 : 3 : 5 : 3 exhibits $`\nu _0=6.86`$ GHz.

We present our results on the dependence of $`\nu _0`$ on the azimuthal gap $`c`$ (Figure 8), radial gap $`g`$ (Figure 9) and ring thickness $`d`$ (Figure 10). In all these cases $`\nu _0`$ increases as the different parameters increase. The dependence shown in Figures 8-10 agrees qualitatively with those done by a different numerical method . When compare our results with the analytical arguments presented by Pendry et al. , we have to keep in mind that that various assumptions about the structural parameters have been done in derivation of Eqn. (4), which are not fulfilled for our structure. Note also that the azimuthal gap does not enter the formula for the resonance frequency given by Eq. (4). Moreover, as the size of the SRR is constant, the structural parameters are not independent each form other. Thus, due to the Eq. (1), increase of the azimuthal gap causes decrease of the inner diameter and vice versa. When taking these restrictions into account, the agreement with analytical results is satisfactory.

### D Material parameters

In Figures 11 and 12 we show how the resonance frequency depends on the permittivity of the dielectric board and on the permittivity of the embedding media. As expected, the resonance frequency decreases considerably with the increase of the value of both permittivities.

## IV Conclusion

In summary, we have used the transfer matrix method to calculate the transmission properties of the left-handed materials and arrays of split ring resonators. The role of absorption of the metallic components of our SRR and LHM has been simulated. It is found that the LHM transmission peak depends on the imaginary part of the metallic permittivity $`ϵ_\mathrm{m}`$, the length of the system and the size of the unit cell. Higher conductivity of the metal guarantees better transmission properties of LHM.

For an array of SRR, the resonance frequency $`\nu _0`$ was computed and is found to agree with experimental data. The dependence of the resonance frequency $`\nu _0`$ on various structural parameters of the SRR were numerically obtained and compared with analytical estimates and also with other numerical techniques.

The main disadvantage of the presented transfer-matrix method is that it can not treat structures with smaller length scales than our discretization mesh. For example, the thickness of the SRR is an order of magnitude smaller in experiments than in our simulation. Also structural parameters of SRR can be changed only discontinuously as multiplies of the unit mesh length. This could be partially overcomed by generalizing the present code to a non-uniform mesh discretization. Nevertheless, already uniform discretization enables us to obtain credible data. Comparison of our results with those obtained by the commercial software MAFIA  confirmed that both methods find the same position of the resonant gap provided that they use the same mesh discretization.

Our numerical data agree qualitatively with the experimental results.  As we can not tune the exact parameters of SRR (as well as its circular shape), and when taken into account the strong dependence of the resonance frequency on the permittivity of the board, we do not expect to obtain very accurate quantitative agreement with experimental data.

Our studies demonstrate that the transfer matrix method can be reliable used to calculate the transmission and reflection properties of left-handed materials and split-ring resonators. Thus, numerical simulations could answer some practical questions about different proposed structures, which might be too complicated to be treated by analytical studies. The transfer matrix method can be used in the future for detailed studies of two-dimensional and even three-dimensional structures. These structures should contain more SRRs and wires per unit cell, which makes their analytical analysis extremely difficult. On the other hand, it is extremely important to find the best design and test the transmission properties of proposed meta-material even before their fabrication and experimental measurements start.

## A

Transfer matrix calculations are based on the scattering formalism. The sample is considered as the scatterer of an incoming wave. The wave normalized to the unit current is coming from the $`\mathrm{}`$, and is scattered by the sample. Scatterer is characterized by four parameters: transmission of the wave from the left to the right ($`t_+`$), from the right to the left ($`t_{}`$), and by reflection coefficient from the right to the right ($`r_+`$) and from the left to the left ($`r_{}`$). Corresponding scattering matrix $`S`$ reads:

$$S=\left(\begin{array}{cc}t_+\hfill & r_+\hfill \\ r_{}\hfill & t_{}\hfill \end{array}\right)$$
(A1)
which determines the amplitudes of the outgoing waves $`B,C`$ in terms of the amplitudes of the incoming waves $`A,D`$:

$$\left(\genfrac{}{}{0pt}{}{C}{B}\right)=S\left(\genfrac{}{}{0pt}{}{A}{D}\right)$$
(A2)
Relation (A2) can be re-written into the form

$$\left(\genfrac{}{}{0pt}{}{D}{C}\right)=𝒯\left(\genfrac{}{}{0pt}{}{B}{A}\right)$$
(A3)
$`𝒯`$ is the transfer matrix, which determines the fields on one side of the sample with the fields on the another side. Its explicit form reads

$$𝒯=\left(\begin{array}{cc}t_{}^1\hfill & t_{}^1r_{}\hfill \\ r_+t_{}^1\hfill & t_+r_+t_{}^1r_{}\hfill \end{array}\right).$$
(A4)
Transfer matrix (TM) fulfills the composition law. If the sample consists from two subsystems, then the transfer matrix $`𝒯_{12}`$ of the whole sample can be calculated form transfer matrices of its subsystems as

$$𝒯_{12}=𝒯_2𝒯_1.$$
(A5)
Resulting TM $`𝒯_{12}`$ has again the form (A4). This composition law enables us to calculate transmission of complicated structure from the transfer matrices of its parts (thin slices).

In numerical calculations, the total volume of the system is divided into small cells and fields in each cell are coupled to those in the neighboring cell. We discretize the Maxwell equations following the method described in Ref. 13. In each point of the lattice we have to calculate four components of the EM field: $`E_x`$, $`E_y`$, $`H_x`$, $`H_y`$.

We assume that our system is connected to two semi-infinite leads (with $`ϵ=1`$ and $`\mu =1`$). EM wave is coming from the right and is scattered by the sample. Resulting waves either continues to the left on the left side lead, or are traveling back to the right on the right side lead. Periodic boundary conditions in the directions perpendicular to the direction of the wave propagation are used.

We decompose the system into $`n`$ thin slices and define a TM for each of them. Explicit form of the TM for a thin slice is in Refs. 12 and 13. The EM field in the $`(k+1)`$th slice can be obtain from the $`k`$th slice as

$$\mathrm{\Phi }_{k+1}=𝒯_k\mathrm{\Phi }_k.$$
(A6)
with $`𝒯_k`$ being the transfer matrix corresponding to the $`k`$th slice. The transfer matrix $`𝒯`$ of the whole sample reads

$$𝒯=𝒯_n𝒯_{n1}\mathrm{}𝒯_2𝒯_1$$
(A7)
If there is $`N`$ mesh points in the slice, then the length of the vector $`\mathrm{\Phi }`$ is $`4N`$ (it contains 4 components of EM field in each point).

Note that we are able to find the explicit form of the TM only in the real space representation. To obtain the transmission, we have to transform the TM into the “wave” representation, which is defined by the eigenvectors of the TM in the leads. Therefore, in the first step we have to diagonalize the TM in the leads.

Each eigenvalue of the TM is two time degenerate because there are two polarizations $`p`$ and $`s`$ of the EM wave. Moreover, if $`\lambda `$ is an eigenvalue, then $`\lambda ^1`$ is also an eigenvalue corresponding to the wave traveling in the opposite direction. In general, the TM has some eigenvalues with modulus equal to 1: $`\lambda =\mathrm{exp}ik`$. The corresponding eigenvectors represent propagating waves. Others eigenvalues are of the form $`\lambda =\mathrm{exp}\pm \kappa `$. They correspond to the evanescent modes. For the frequency range which is interesting for the LHM studies, the TM has only one propagating mode.

As the TM is not Hermitian matrix, we have to calculate left and right eigenvectors separately. From the eigenvectors we construct three matrices: The $`2N\times 4N`$ matrix $`R_1`$ contains in its columns $`2N`$ right eigenvectors which correspond to the wave traveling to the left. Matrices $`L_1`$ and $`L_2`$ are $`4N\times 2N`$ matrices which contains in their rows the left eigenvectors corresponding to waves traveling to the left and to the right, respectively.

The general expression of the TM given by Eq. (A4) enables us to find the transmission matrix explicitly

$$t_{}^1=L_1𝒯R_1$$
(A8)
and the reflection matrix from the relation

$$r_+t_{}^1=L_2𝒯R_1.$$
(A9)
At this point we have to distinguish between the propagating and the evanescent modes. For a frequency range of interest, the TM in leads has only one propagating mode for each direction. We need therefore only $`2\times 2`$ sub-matrices $`t_{}(ij)`$ and $`r_+(ij)`$ with $`i,j=1`$ or 2 for the $`p`$ or $`s`$ polarized wave. The transmission and reflection are then

$$T_{ij}=t_{}(ij)t_{}^{}(ij)R_{ij}=r_+(ij)r_+^{}(ij)$$
(A10)
and absorption

$$\begin{array}{cc}A_p\hfill & =1T_{pp}T_{ps}R_{pp}R_{ps}\hfill \\ A_s\hfill & =1T_{ss}T_{sp}R_{ss}R_{sp}\hfill \end{array}$$
(A11)
It seems that relations (A8) and (A9) solve our problem completely. However, the above algorithm must be modified. The reason is that the elements of the matrix $`t_{}^1`$ are given by their larger eigenvalues. We are, however, interesting in the largest eigenvalues of the matrix $`t_{}`$. As the elements of the transfer matrix increase exponentially in the iteration procedure given by Eq. (A7), an information about the smallest eigenvalues of $`t_{}^1`$ will be quickly lost. We have therefore to introduce some re-normalization procedure. We use the procedure described in Ref. 16.

Relation (A8) can be written as

$$t_{}^1=L_1r^{(n)}$$
(A12)
where we have defined $`2N\times 4N`$ matrices $`r^{(k)}`$, $`k=0,1,\mathrm{}n`$ as

$$r^{(k)}=𝒯_kr^{(k1)},r^{(0)}=R_1$$
(A13)
Each matrix $`r`$ can be written as

$$r=\left(\genfrac{}{}{0pt}{}{r_1}{r_2}\right)$$
(A14)
with $`r_1`$, $`r_2`$ being the $`2N\times 2N`$ matrices. We transform $`r`$ as

$$r=r^{}r_1r^{}=\left(\genfrac{}{}{0pt}{}{1}{r_2r_1^1}\right)$$
(A15)
and define $`r^{(k)}=𝒯_k(r^{})^{(k1)}`$. In contrast to $`r_1`$ and $`r_2`$, all eigenvalues of the matrix $`r_2r_1^1`$ are of order of unity. Relation (A12) can be now re-written into the form

$$t_{}^1=L_1\left(\genfrac{}{}{0pt}{}{1}{r_2^{(n)}\left[r_1^n\right]^1}\right)r_1^{(n)}r_1^{(n1)}\mathrm{}r_1^{(1)}r_1^{(0)}$$
(A16)
from which we get that

$$\begin{array}{cc}t_{}=\left[r_1^{(0)}\right]^1\left[r_1^{(1)}\right]^1\mathrm{}\left[r_1^{(n)}\right]^1\left[L_1\left(\genfrac{}{}{0pt}{}{1}{r_2^{(n)}\left[r_1^{(n)}\right]^1}\right)\right]^1.\hfill & \end{array}$$
(A17)
From Eqn. (A9) we find

$$r_+=\left[L_2\left(\genfrac{}{}{0pt}{}{1}{r_2^{(n)}\left[r_1^{(n)}\right]^1}\right)\right]\times \left[L_1\left(\genfrac{}{}{0pt}{}{1}{r_2^{(n)}\left[r_1^{(n)}\right]^1}\right)\right]^1.$$
(A18)
The matrix inversion in the formulae (A15-A18) can obtained also by the soluton of a system linear equations. Indeed, matrix $`BA^1`$ equals to matrix $`X`$, which solves the system of linear equations $`B=XA`$. CPU time could be reduces considerably in this way, especially for large matrices.

All elements of the matrices on the rhs of Eqn. (A17) are of order of unity. The price we have to pay for this stability is an increase of the CPU time. Fortunately, if the elements of the transfer matrix are not too large (which is not the case in systems studied in this paper), then it is enough to perform described normalization procedure only after every 6-8 steps.

We thank D.R. Smith, M. Agio and D. Vier for fruitful discussions. Ames Laboratory is operated for the U.S.Department of Energy by Iowa State University under Contract No. W-7405-Eng-82. This work was supported by the Director of Energy Research, Office of Basic Science, DARPA and NATO grant PST.CLG.978088. P.M. thanks Ames Laboratory for its hospitality and support and Slovak Grant Agency for financial support.