text:
# Surface plasmon toy-models of black holes and wormholes.

\[

## Abstract

Surface plasmons at metal interfaces are collective excitations of the conduction electrons and the electromagnetic field. They exist in ”curved three-dimensional space-times” defined by the shape of the metal surface and the spatial distribution of the dielectric constant near the surface. Here we show that surface plasmon toy models of many non-trivial space-time metrics, such as wormholes and black holes, can be easily built and studied in experiments. For example, a droplet of dielectric on the metal surface behaves as a black hole for surface plasmons within a substantial frequency range. On the other hand, a nanohole in a thin metal membrane may be treated as a wormhole connecting two ”flat” surface plasmon worlds located on the opposite interfaces of the membrane.

PACS no.: 78.67.-n, 04.70.Bw

\]

The realization that solid-state toy models may help in an understanding of electromagnetic phenomena in curved space-time has led to considerable recent effort in developing toy models of electromagnetic  and sonic  black holes. In the case of media electrodynamics this is possible because of an analogy between the propagation of light in matter and in curved space-times: it is well known that Maxwell s equations in a general curved space-time background $`g_{ik}(x,t)`$ are equivalent to the phenomenological Maxwell equations in the presence of a matter background with nontrivial electric and magnetic permeability tensors $`ϵ_{ij}(x,t)`$ and $`\mu _{ij}(x,t)`$ . In this analogy, the event horizon corresponds to a surface of singular electric and magnetic permeabilities, so that the speed of light goes to zero, and light is ”frozen” near such a surface. In the absence of established quantum gravitation theory the toy models are helpful in understanding electromagnetic phenomena in curved space-times, such as Hawking radiation  and the Unruh effect . Very recently it was also realized  that many results obtained in the optics of random media may be applicable to the case of light propagation in stochastic space-time metrics. Thus, an exchange of ideas between media electrodynamics and gravitation theory proves to be useful for both fields. Unfortunately, up to now all the suggested electromagnetic black hole toy-models were very difficult to realize and study experimentally, so virtually no experimental work has been done in this field. For example, the electromagnetic toy black hole suggested in  requires a peculiar spatial distribution of temperature inside a ferroelectric and ferromagnetic crystal, which would be very difficult to create and maintain over any reasonable period of time.

In this Letter we propose novel electromagnetic toy models of wormholes and black holes, which are based on surface plasmons and are very easy to build and study. We are going to consider a close analogy between the optics of surface plasmons propagating along metal interfaces, which are curved and/or covered with dielectric films, and the field theory in a curved space-time background. Surface plasmons are collective excitations of the conductive electrons and the electromagnetic field . They exist in ”curved three-dimensional space-times” defined by the shape of the metal-dielectric interface. Since in many experimental geometries surface plasmons are weakly coupled to the outside world (to free-space photons, phonons, single-electron excitations, etc.) it is reasonable to treat the physics of surface plasmons separately from the rest of the surface and bulk excitations, so that a field-theory of surface plasmons in a curved space-time background may be considered. For example, a nanohole in a thin metal membrane may be treated as a ”wormhole” connecting two ”flat” surface plasmon worlds located on the opposite interfaces of the membrane . This nanohole may be filled with a material with some nontrivial electric and magnetic permeability tensors $`ϵ_{ij}(x,t)`$ and $`\mu _{ij}(x,t)`$, so that nontrivial space-time metrics can be emulated. This point of view allowed us to justify and explore an analogy between the nonlinear optics of cylindrical surface plasmons of nanowires and nanoholes and the lower-dimension Kaluza-Klein theories .

On the other hand, near the plasmon resonance (which is defined by the condition that $`ϵ_m(\omega )=ϵ_d`$ at the metal-dielectric interface, where $`ϵ_m(\omega )`$ and $`ϵ_d`$ are the dielectric constants of metal and dielectric, respectively ) the surface plasmon velocity vanishes, so that the surface plasmon ”stops” on the metal surface, and the surface charge and the normal component of the electric field diverge. For a given frequency of light, the spatial boundary of the plasmon resonance (”the event horizon” of our toy model) may be defined at will using the geometry $`ϵ_d(x,y)`$ of the absorbed layer of dielectric on the metal surface. Thus, the plasmon resonance becomes a natural candidate to emulate the event horizon of a black hole. As a result, toy two-dimensional surface plasmon black holes can be easily produced and studied. In what follows we are going to show that consideration of such surface plasmon toy wormholes and black holes may be useful both for the fields of nanooptics and for gravitation theory.

Let us consider in detail the dispersion law of a surface plasmon (SP), which propagates along the metal-dielectric interface. The SP field decays exponentially both inside the metal and the dielectric. Inside the dielectric the decay exponent is roughly equal to the SP wave vector. As a first step let us assume that both metal and dielectric completely fill the respective $`z<0`$ and $`z>0`$ half-spaces. In such a case the dispersion law can be written as

$$k^2=\frac{\omega ^2}{c^2}\frac{ϵ_dϵ_m(\omega )}{ϵ_d+ϵ_m(\omega )},$$
(1)
where we will assume that $`ϵ_m=1\omega _p^2/\omega ^2`$ according to the Drude model, and $`\omega _p`$ is the plasma frequency of the metal. This dispersion law is shown in Fig.1(b) for the cases of metal-vacuum and metal-dielectric interfaces. It starts as a ”light line” in the respective dielectric at low frequencies and approaches asymptotically $`\omega =\omega _p/(1+ϵ_d)^{1/2}`$ at very large wave vectors. The latter frequency corresponds to the so-called surface plasmon resonance. Under the surface plasmon resonance conditions both phase and group velocity of the SPs is zero, and the surface charge and the normal component of the electric field diverge. Since at every wavevector the SP dispersion law is located to the right of the ”light line”, the SPs of the plane metal-dielectric interface are decoupled from the free-space photons due to the momentum conservation law.

If a droplet of dielectric (Fig.1(a)) is placed on the metal surface, the SP dispersion law will be a function of the local thickness of the droplet. Deep inside the droplet far from its edges the SP dispersion law will look similar to the case of a metal-dielectric interface, whereas near the edges (where the dielectric is thinner) it will approach the SP dispersion law for the metal-vacuum interface. As a result, for every frequency between $`\omega _p/(1+ϵ_d)^{1/2}`$ and $`\omega _p/2^{1/2}`$ there will be a closed linear boundary inside the droplet for which the surface plasmon resonance conditions are satisfied. Let us show that such a droplet of dielectric on the metal interface behaves as a ”surface plasmon black hole” in the frequency range between $`\omega _p/(1+ϵ_d)^{1/2}`$ and $`\omega _p/2^{1/2}`$, and that the described boundary of the surface plasmon resonance behaves as an ”event horizon” of such a black hole.

Let us consider a SP within this frequency range, which is trapped near its respective ”event horizon”, and which is trying to leave a large droplet of dielectric (see Fig.2). The fact that the droplet is large means that ray optics may be used. Since the component of the SP momentum parallel to the droplet boundary has to be conserved, such a SP will be totally internally reflected by the surface plasmon resonance boundary back inside the droplet at any non-zero angle of incidence. This is a simple consequence of the fact that near the ”event horizon” the effective refractive index of the droplet for surface plasmons is infinite (according to eq.(1), both phase and group velocity of surface plasmons is zero at surface plasmon resonance). On the other hand, even if the angle of incidence is zero, it will take the SP infinite time to leave the resonance boundary. Thus, the droplet behaves as a black hole for surface plasmons, and the line near the droplet boundary where the surface plasmon resonance conditions are satisfied plays the role of the event horizon for surface plasmons.

The toy black holes described above are extremely easy to make and observe. In our experiments a small droplet of glycerin was placed on the gold film surface and further smeared over the surface using lens paper, so that a large number of glycerin microdroplets were formed on the surface (Fig.3(a)). These microdroplets were illuminated with white light through the glass prism (Fig.1(a)) in the so-called Kretschman geometry . The Kretschman geometry allows for efficient SP excitation on the gold-vacuum interface due to phase matching between the SPs and photons in the glass prism. As a result, SPs were launched into the gold film area around the droplet. Photograph taken under a microscope of one of such microdroplets is shown in Fig.3(b). The white rim of light near the edge of the droplet is clearly seen. It corresponds to the effective SP event horizon described above. Near this toy event horizon SPs are stopped or reflected back inside the droplet. In addition, a small portion of the SP field may be scattered out of the two-dimensional surface plasmon world into normal three-dimensional photons. These photons produced the image in Fig.3(b). If we are to continue our analogy with real black holes in the worlds with other than three spatial dimensions, we should state that near the event horizon the surface plasmons are scattered from their two-dimensional world into ”extra dimensions” (which in this case is our normal three-dimensional space).

Unfortunately, the described toy SP black hole model does not work outside the frequency range between $`\omega _p/(1+ϵ_d)^{1/2}`$ and $`\omega _p/2^{1/2}`$. On the other hand, this is a common feature of every electromagnetic toy black hole model suggested so far. All such toy models necessarily work only within a limited frequency range. However, the ease of making and observing such toy SP black holes makes them a very promising research object. If we forget about the language of ”black holes” and ”event horizons” for a moment, the SP optics phenomenon represented in Fig.3(b) remains a potentially very interesting effect in surface plasmon optics. Namely, this photo shows the existence of a two-dimensional SP analog of whispering gallery modes, which are well-known in the optics of light in droplets and other spherical dielectric particles. Whispering gallery modes in liquid microdroplets are known to substantially enhance nonlinear optical phenomena due to cavity quantum electrodynamic effects . One may expect even higher enhancement of nonlinear optical mixing in liquid droplets on the metal surfaces due to enhancement of surface electromagnetic field inherent to surface plasmon excitation, and in addition, due to accumulation of SP energy near the surface plasmon event horizons at the droplet boundaries. This strong enhancement of nonlinear optical effects in liquid droplets may be very useful in chemical and biological sensing applications.

Looking back at the consequences of our model for gravitation theory, we may anticipate that such effects as recently predicted optical second harmonic generation near the black hole event horizon  can be easily observed in experiments with toy SP black holes. One can even look forward to modeling of more exotic (and maybe more realistic) gravitation theory situations where a wormhole exists inside a black hole. In the more prosaic language of surface plasmon optics this situation corresponds to nanoholes drilled in a free standing metal membrane and covered with droplets of dielectric (Fig.4(a,c)). In fact, our recent experiments on single-photon tunneling  and optical control of photon tunneling  were performed with exactly these kinds of samples. Theoretical analysis of these experiments performed in  indicated some usefulness of the gravitational theory analogy: parallels can be drawn between the nonlinear optics of surface plasmons in cylindrical nanoholes (cylindrical surface plasmons) and areas of gravitation theory such as Kaluza-Klein theories .

This analogy stems from the way in which electric charges are introduced in the original five-dimensional Kaluza-Klein theory (see for example  and Fig.4(b)). In this theory the electric charges are introduced as chiral (nonzero angular momentum) modes of a massless quantum field, which is quantized over the cyclic compactified fifth dimension. Electromagnetic forces between the electric charges appear as nonlinear coupling of these chiral modes, so that four-dimensional electrodynamics described by the Maxwell equations may be understood as nonlinear optics of these modes. Similar Kaluza-Klein theories may be formulated in lower-dimensional space-times. Such theories reproduce electrodynamics of electric charges in worlds with less than three spatial dimensions. Since nonlinear interactions of surface plasmons of a cylindrical nanowire or a nanohole look as if they occur in a space-time which besides an extended z-coordinate has a small ”compactified” angular $`\varphi `$-dimension along the circumference of the cylinder, the theory of cylindrical surface plasmon (CSP) mode propagation and interaction appears to be similar to the three-dimensional Kaluza-Klein theory: solutions of the nonlinear Maxwell equations for interacting cylindrical surface plasmons may be found , which behave as interacting Kaluza-Klein charges. According to these solutions, higher $`(n>0)`$ CSP modes posses quantized effective chiral charges proportional to their angular momenta $`n`$. In a metal nanowire these slow moving effective charges exhibit long-range interaction via exchange of fast massless CSPs with zero angular momentum. These zero angular momentum CSPs may be considered as massless quanta of the gyration field (the field of the gyration vector $`\stackrel{}{g}`$), which relates the $`\stackrel{}{D}`$ and $`\stackrel{}{E}`$ fields in an optically active medium :

$$\stackrel{}{D}=ϵ\stackrel{}{E}+i\stackrel{}{E}\times \stackrel{}{g},$$
(2)
In order for such a physical picture to be valid the medium in or around the nanohole or nanowire should be chiral or optically active, and exhibit the magneto-optical effect. This condition is fulfilled automatically in the case of any metal , since all metals are optically active in the presence of a magnetic field. Thus, optical experiments performed on nanowires and nanoholes may become the proving ground for experimental testing of many ideas in theoretical physics, such as compactified extra dimensions and wormholes.

Detailed illustration of how the nonlinear optics of cylindrical surface plasmons may be formulated in a way that is similar to Kaluza-Klein theories can be found in . The value of the Kaluza-Klein and black hole models for the field of nanooptics is quite apparent: these models provide very useful intuitive guides for finding solutions of nonlinear Maxwell equations in the situations when nonlinear interactions are strong and can not be considered as small perturbations. In such situations, often one can only guess the general form of the solutions, based on the comparison of the nonlinear system of interest with other better understood nonlinear systems described by similar equations. We can also ask if experimental and theoretical studies of nanoholes and microdroplets can be useful for gravitation theory. Some arguments can be made that this might be the case. For example, the solutions for the electromagnetic eigenmodes of the wormholes and their interaction are often sought for very narrow wormholes, which are understood as strings connected to D-branes (D-dimensional membranes) . The multi-component Kaluza-Klein charges of the wormholes in such configurations have some resemblance to the multi-component chiral charges of metal nanowires and nanoholes described in . Thus, solid state models may provide some useful insights. It is customary in theoretical description of wormholes in higher-dimensional Kaluza-Klein theories to distinguish between electric and magnetic Kaluza-Klein charges and fields . According to this classification, we should identify the chiral charges of the cylindrical surface plasmons as magnetic charges, while (quite naturally) the conductive electrons in the metal nanowire perform the role of the electric Kaluza-Klein charges . Another important question, which follows from the analogy between the nanoholes and the wormholes arises from the fact that many optical waveguides (from the electrodynamics point of view a wormhole is a waveguide) have cut-off frequencies for the propagating electromagnetic modes. For example, unlike metal nanowires, the nanoholes do not support propagating electromagnetic modes with small frequencies. This may mean that in gravitation theory an attempt to describe a wormhole using a frequency-independent space-time metric may fail in some cases. Excitations with small energies below the wormhole cut-off may not see the wormhole throat at all.

In conclusion, we have introduced and observed experimentally surface plasmon analogues of such nontrivial space-time topologies as black holes and wormholes. We anticipate that further consideration of such analogues will produce a mutually beneficial exchange of ideas between nanooptics and gravitation theory.

Figure captions.

Fig.1 (a) Experimental geometry of surface plasmon toy black hole observation. (b) Surface plasmon dispersion law for the cases of metal-vacuum and metal-dielectric interfaces

Fig.2 A surface plasmon trapped inside a droplet near the effective ”event horizon”: The projection of surface plasmon momentum parallel to the droplet edge must be conserved. Due to effectively infinite refractive index near the droplet edge surface plasmons experience total internal reflection at any angle of incidence.

Fig.3 Microscopic photographs of a toy surface plasmon black hole. (a) Droplet of glycerin on a gold film surface (illuminated from the top). The droplet diameter is approximately 15 micrometers. (b) The same droplet illuminated in the Kretschman geometry, which provides efficient coupling of light to surface plasmons on the gold-vacuum interface (Fig.1(a)). The white rim around the droplet boundary corresponds to the effective surface plasmon ”event horizon”.

Fig.4 (a) A nanohole in a metal membrane may be treated as a wormhole for surface plasmons which exist on the flat top and bottom interfaces of the membrane. Alternatively, this figure may represent a nanowire connecting two metal-vacuum interfaces. (b) Similarity between Kaluza-Klein electric charges and cylindrical surface plasmons. (c) Array of 20 nm diameter nanoholes drilled in a free standing gold membrane, which has been studied experimentally . Such nanoholes may be treated as surface plasmon wormholes.


format: text

text:
# Noncommutative field theory and violation of translation invariance

## 1 Introduction

Noncommutativity of coordinates has been intensively studied in the literature as it arises in the context of string theory , but also because it has interesting properties and implications in field theory . In many treatments of noncommutative field theory the noncommutative parameters are not regarded as Lorentz tensors, but instead a set of numbers that do not transform covariantly which implies naturally in the breaking of Lorentz invariance down to the stability subgroup of the noncommutative parameter . Furthermore, noncommutative structures of the Lie-type imply in the violation of the energy-momentum conservation as the coordinate commutation relations break translational invariance. Alternatively, one could consider instead the noncommutative parameter as a Lorentz tensor covariant under Lorentz boosts. This approach has been studied earlier and in this framework we have shown that a noncommmutative scalar field coupled to gravity admits a covariant formulation (where associativity111The relation between a nonassociative star product on D-branes and noncommutative theories on curved spaces has been discussed in . is mantained only at perturbative level) which is compatible with a homogeneous and isotropic space-time . Further attempts along these lines include work on noncommutative scalar field theory in three-dimensions  and on QED . We mention that in the latter, gauge invariance is implemented via the introduction of a noncommutative gauge field, a “star” gauge invariance and “star” commutators.

In this work we shall analyse classical and quantum field theory features of models where the noncommutativity of the coordinates has the following form

$$[x^\mu ,x^\nu ]=i\mathrm{\Lambda }_\omega ^{\mu \nu }x^\omega ,$$
(1)
with the condition of nilpotency as specified below. This leads to a violation of the symmetry under translations and, consequently, requires a reformulation of energy-momentum conservation. This reformulation is proven to be Poincaré-invariant and reduces to the usual momentum conservation in the commutative limit. The formalism we develop follows the study of ref. , where the Baker-Hausdorff formula is related to the Kontsevich (see  and ealier references therein) noncommutative product, for Lie algebras of the form (1).

Before closing our introduction, let us point out that our nilpotency condition (cf. eq. (14) below) excludes noncommutative structures of the Lie-type such as the semisimple Lie algebras ($`SU(2)`$ for the fuzzy sphere) and the $`\kappa `$-deformed Minkowsky space that are often discussed in the literature (see ref.  for an extensive review). Furthermore, we mention that recent work by Robbins and Sethi  is closely related with ours, even though it considers examples that are more directly inspired by string theory.

## 2 Mathematical formulation

### 2.1 Noncommutative algebra

A noncommutative associative product may be defined through the Lie-algebra commutator eq. (1), where $`\mathrm{\Lambda }^{\mu \nu \omega }`$ is a real tensor with units of $`\text{mass}^1`$ and $`\mathrm{\Lambda }^{\mu \nu \omega }=\mathrm{\Lambda }^{\nu \mu \omega }`$. On its hand, associativity implies the Jacobi identity

$$\mathrm{\Lambda }_\omega ^{\mu \nu }\mathrm{\Lambda }_\beta ^{\omega \alpha }+\mathrm{\Lambda }_\omega ^{\nu \alpha }\mathrm{\Lambda }_\beta ^{\omega \mu }+\mathrm{\Lambda }_\omega ^{\alpha \mu }\mathrm{\Lambda }_\beta ^{\omega \nu }=0.$$
(2)
A noncommutative Fourier mode is defined by

$$e_{}^{ikx}=\underset{n=0}{\overset{\mathrm{}}{}}\frac{i^n}{n!}\stackrel{n\text{ factors}}{\stackrel{}{\left(kx\right)\mathrm{}\left(kx\right)}}=\underset{n=0}{\overset{\mathrm{}}{}}\frac{i^n}{n!}\left(kx\right)_{}^n,$$
(3)
and we study the functional space spanned by these Fourier modes, with elements of the form

$$f\left(x\right)=\frac{d^nk}{\left(2\pi \right)^n}\stackrel{~}{f}\left(k\right)e_{}^{ikx}$$
(4)
which, in the commutative limit, reduces to the usual Hilbert space. Notice that in eq. (3) the star product acts only on the configuration variables and not on the momentum ones.

The product of two generic functions is then given by

$$fg=\frac{d^nk}{\left(2\pi \right)^n}\frac{d^nq}{\left(2\pi \right)^n}\stackrel{~}{f}\left(k\right)\stackrel{~}{g}\left(q\right)e_{}^{ikx}e_{}^{iqx},$$
(5)
where we have expressed the functions in terms of their noncomutative Fourier expansion. This product is completely determined if the product of two Fourier modes $`e_{}^{ikx}e_{}^{iqx}`$ can be evaluated. This can be achieved by making use of the Baker-Hausdorff formula

$$e_{}^{ikx}e_{}^{iqx}=\mathrm{exp}_{}\left\{i\left(k+q\right)x+\frac{1}{2}[ikx,iqx]+\mathrm{}\right\},$$
(6)
where the dots stand for higher order commutators. Since the commutators obey

$$[x^{\mu _1},[x^{\mu _2},\mathrm{},[x^{\mu _n},x^\nu ]]\mathrm{}]i^nx^\omega ,$$
(7)
the product of two Fourier modes is a Fourier mode

$$e_{}^{ikx}e_{}^{iqx}=e_{}^{i\left[k+q+V(k,q)\right]x}$$
(8)
with $`V`$ determined by the Baker-Hausdorff expansion:

$$V_\omega (k,q)=k_\mu q_\nu \mathrm{\Lambda }_\lambda ^{\mu \nu }\left[\frac{1}{2}\delta _\omega ^\lambda +\frac{k_\alpha q_\alpha }{12}\mathrm{\Lambda }_\omega ^{\alpha \lambda }\right]+O(\mathrm{\Lambda }^3).$$
(9)
### 2.2 Quadratic actions

In order to build actions, a star-integration must be defined. In the functional space whose elements are of the form (4), any function can be integrated if the integral of a Fourier mode is known. Hence, we introduce the following star-integration

$$_{}d^nxe_{}^{irx}=\left(2\pi \right)^n\delta \left(r\right),$$
(10)
which yields the usual integration in the commutative limit.

Consider now the star-integral

$$I=_{}d^nxfg.$$
(11)
which, in Fourier space, is written as

$$I=\frac{d^nk}{\left(2\pi \right)^n}\frac{d^nq}{\left(2\pi \right)^n}\stackrel{~}{f}\left(k\right)\stackrel{~}{g}\left(qk\right)_{}d^nxe_{}^{i\left[q+V(k,qk)\right]x},$$
(12)
and implies

$$I=\frac{d^nk}{\left(2\pi \right)^n}d^nq\stackrel{~}{f}\left(k\right)\stackrel{~}{g}\left(qk\right)\delta \left(q+V(k,qk)\right).$$
(13)
If the structure constants are nilpotent, that is, for $`n>n_{}`$

$$\mathrm{\Lambda }_{\omega _1}^{\mu _1\nu }\mathrm{\Lambda }_{\omega _2}^{\mu _2\omega _1}\mathrm{}\mathrm{\Lambda }_{\omega _n}^{\mu _n\omega _{n1}}=0,$$
(14)
then

$$\delta \left(q+V(k,qk)\right)=\frac{\delta \left(q\right)}{\left|det\left(\delta _\nu ^\mu \frac{V_\nu }{q_\mu }\right)\right|}=\delta \left(q\right)$$
(15)
since $`det\left(1+M\right)=1`$ if $`M^n=0`$, which holds if $`\mathrm{\Lambda }`$ is nilpotent. Thus

$$I=\frac{d^nk}{\left(2\pi \right)^n}\stackrel{~}{f}\left(k\right)\stackrel{~}{g}\left(k\right)=d^nxf_C\left(x\right)g_C\left(x\right)$$
(16)
where $`f_C`$, $`g_C`$ are inverse Fourier transforms *using commutative Fourier modes*

$$f_C\left(x\right)=\frac{d^nk}{\left(2\pi \right)^n}\stackrel{~}{f}\left(k\right)e^{ikx}.$$
(17)
Equation (16) states that, in momentum space, quadratic terms in the lagrangian are the same as their commutative counterparts. In particular, free propagators will remain unchanged.

## 3 Violation of momentum conservation

We have concluded that the quadratic part of a lagrangian is not changed and, hence, the free theory is the same as the commutative one. In particular, the free Green function is equal to the commutative case and the dispersion relation $`ϵ^2=p^2+m^2`$ is unchanged, since it is given by the poles of the free propagator. Yet, we shall see that interactions are altered by non-commutativity.

Consider a non-commutative field theory, with generic fields $`A_i`$ and an interaction term

$$S_I=_{}d^nxM_{i_1\mathrm{}i_m}A_{i_1}\mathrm{}A_{i_m},$$
(18)
where $`M_{i_1\mathrm{}i_m}`$ are constants.

Writing the fields in momentum space we get

$$S_I=\left[\underset{i=1}{\overset{m}{}}\frac{d^nk_i}{\left(2\pi \right)^n}\right]\stackrel{~}{M}_{i_1\mathrm{}i_m}\left(\underset{¯}{k}_m\right)\stackrel{~}{A}_{i_1}\left(k_1\right)\mathrm{}\stackrel{~}{A}_{i_m}\left(k_m\right)$$
(19)
where we use the notation $`\underset{¯}{k}_m=(k_1,\mathrm{},k_m)`$. The interaction in momentum space is given by

$$\stackrel{~}{M}_{i_1\mathrm{}i_m}\left(\underset{¯}{k}_m\right)=M_{i_1\mathrm{}i_m}_{}d^nxe_{}^{ik_1x}\mathrm{}e_{}^{ik_mx}.$$
(20)
In eq. (19) the variables $`k_i`$ are mute, so we can sum over all $`\pi `$ permutations of the indices $`i_m`$:

$$S_I=\left[\underset{i=1}{\overset{m}{}}\frac{d^nk_i}{\left(2\pi \right)^n}\right]\stackrel{~}{M}_{i_1\mathrm{}i_m}^{\mathrm{symm}}\left(\underset{¯}{k}_m\right)\stackrel{~}{A}_{i_1}\left(k_1\right)\mathrm{}\stackrel{~}{A}_{i_m}\left(k_m\right),$$
(21)
where

$$\stackrel{~}{M}_{i_1\mathrm{}i_m}^{\mathrm{symm}}\left(\underset{¯}{k}_m\right)=\frac{1}{m!}\underset{\pi \mathrm{perm}.}{}()^{N\left(\pi \right)}\stackrel{~}{M}_{i_{\pi \left(1\right)}\mathrm{}i_{\pi \left(m\right)}}\left(\underset{¯}{k}_{\pi \left(m\right)}\right).$$
(22)
To evaluate eq. (20), we use the expression

$$e_{}^{ik_1x}\mathrm{}e_{}^{ik_mx}=\mathrm{exp}_{}\left\{i\underset{j=1}{\overset{m}{}}k_jx+iV^m\left(\underset{¯}{k}_m\right)x\right\}$$
(23)
where

$$V^m\left(\underset{¯}{k}_m\right)=V^{m1}\left(\underset{¯}{k}_{m1}\right)+V(\underset{i=1}{\overset{m1}{}}k_i+V^{m1}\left(\underset{¯}{k}_{m1}\right),k_m)$$
(24)
with $`V^2\left(\underset{¯}{k}_2\right)=V(k_1,k_2)`$. This yields both the noncommutative energy-momentum law and the noncommutative vertex

$$\stackrel{~}{M}_{i_1\mathrm{}i_m}\left(\underset{¯}{k}_m\right)=\left(2\pi \right)^n\delta \left(\underset{i=1}{\overset{m}{}}k_i+V^m\left(\underset{¯}{k}_m\right)\right)M_{i_1\mathrm{}i_m}.$$
(25)
Hence, the new energy-momentum law for the vertex reads

$$\underset{i=1}{\overset{m}{}}k_i+V^m\left(\underset{¯}{k}_m\right)=0.$$
(26)
The full theory involves $`\stackrel{~}{M}_{i_1\mathrm{}i_m}^{symm}`$ , which will have contributions whenever

$$\underset{i=1}{\overset{m}{}}k_i+V^m\left(\underset{¯}{k}_{\pi \left(m\right)}\right)=0$$
(27)
for all $`m!`$ permutations of indices, $`\pi `$.

Thus, we see that the energy-momentum conservation is violated as the theory is not invariant under translations. In fact, in a translation $`x^\mu x^\mu +b^\mu `$, the commutator of the coordinates is changed by

$$[x^\mu ,x^\nu ]i\mathrm{\Lambda }_\omega ^{\mu \nu }x^\omega +i\theta ^{\mu \nu },$$
(28)
that is, a constant term $`\theta ^{\mu \nu }=\mathrm{\Lambda }_\omega ^{\mu \nu }b^\omega `$ is added to the commutator of the coordinates. So, the interaction vertex becomes

$$\stackrel{~}{M}_{i_1\mathrm{}i_m}\left(\underset{¯}{k}_m\right)\left(2\pi \right)^n\delta \left(\underset{i=1}{\overset{m}{}}k_i+V^m\left(\underset{¯}{k}_{\pi \left(m\right)}\right)\right)M_{i_1\mathrm{}i_m}\mathrm{exp}\left\{i\theta ^m\left(\underset{¯}{k}_m\right)\right\}$$
(29)
where

$$\theta ^m\left(\underset{¯}{k}_m\right)=\theta ^{m1}\left(\underset{¯}{k}_{m1}\right)+\theta (\underset{i=1}{\overset{m1}{}}k_i+V^{m1}\left(\underset{¯}{k}_{m1}\right),k_m)$$
(30)
and $`\theta ^2\left(\underset{¯}{k}_2\right)=k_{1\mu }\theta ^{\mu \nu }k_{2\nu }`$. Hence, the interaction vertex is altered by an overall oscillating momentum-dependent factor and, thus, invariance under translations is broken. This example shows that translations give always rise to a constant term in the noncomutative tensor. However, the new energy- momentum law is unchanged, so it is a Poincaré-invariant expression, even though the theory is not.

## 4 Kinematical applications

### 4.1 Preliminaries

The first non-trivial behaviour arising from the new interaction vertex occurs with three particles. The energy-momentum equation is found to be

$$k_1+k_2+k_3+V(k_1,k_2)+V(k_1+k_2+V(k_1,k_2),k_3)=0$$
(31)
and similar expressions for all permutations of the indices.

On physical grounds, eq. (31) represents three particles interacting. Formally, these equations can be treated in the context of the usual momentum conservation if one thinks in terms of four interacting particles, with the fourth particle’s energy-momentum vector being given by a nonlinear function of the others. This reasoning may be extended to the $`m`$-particle case, eq. (27).

For the time being, let us consider a simple model with

$$\mathrm{\Lambda }_{\omega _1}^{\mu _1\nu }\mathrm{\Lambda }_{\omega _2}^{\mu _2\omega _1}=0$$
(32)
which complies with the Jacobi identity eq. (2).

The energy-momentum equation becomes

$$k_1+k_2+k_3+V(k_1,k_2)=0,$$
(33)
where

$$V_\omega (k_1,k_2)=\frac{1}{2}k_{1\mu }k_{2\nu }\mathrm{\Lambda }_\omega ^{\mu \nu }.$$
(34)
There are nontrivial covariant solutions to eq. (32). For instance, consider a constant antisymmetric tensor $`\mathrm{\Lambda }^{\mu \nu }=\mathrm{\Lambda }^{\nu \mu }`$ with nontrivial kernel, that is, $`det\mathrm{\Lambda }=0`$, and a non-vanishing vector $`r^\omega `$ belonging to this kernel. Hence a solution is given by

$$\mathrm{\Lambda }^{\mu \nu \omega }=\mathrm{\Lambda }^{\mu \nu }r^\omega .$$
(35)
In four dimensions we can parametrize $`\mathrm{\Lambda }^{\mu \nu }`$ with two spatial vectors $`\stackrel{}{E}`$ and $`\stackrel{}{B}`$

$$\mathrm{\Lambda }^{\mu \nu }=\left(\begin{array}{cccc}0& E_x& E_y& E_z\\ E_x& 0& B_z& B_y\\ E_y& B_z& 0& B_x\\ E_z& B_y& B_x& 0\end{array}\right),r_\nu =\left(\begin{array}{c}r_0\\ r_x\\ r_y\\ r_z\end{array}\right).$$
(36)
Condition eq. (35) implies that

$$r^2=\left|\stackrel{}{r}\right|^2\left[\left(\frac{B}{E}\mathrm{sin}\delta \right)^21\right],$$
(37)
with $`\delta `$ being the angle between $`\stackrel{}{B}`$ and $`\stackrel{}{r}`$. The massless, massive and tachyon regimes of $`V`$ are readily identifiable. Since we assume that $`\mathrm{\Lambda }^{\mu \nu }`$ is a Lorentz tensor, there are always inertial frames where $`\stackrel{}{E}`$ is non-vanishing, and the above expression holds only for such frames. If $`B<E`$ (a Lorentz-invariant inequality) then $`r^\omega `$ behaves like a tachyon; otherwise, the behaviour of $`r^\omega `$ will depend on $`\delta `$.

From the momentum conservation, eq. (33), we get the following result

$$\mathrm{\Lambda }^{\mu \nu }\left(k_1+k_2+k_3\right)_\nu =0,$$
(38)
which states that the vector sum of the momenta belongs to the (nontrivial) kernel of the noncommutative tensor. We also have the following expressions

$$\mathrm{\Lambda }^{\mu \nu }k_\nu =\left(\begin{array}{c}\stackrel{}{E}\stackrel{}{k}\\ k_0\stackrel{}{E}+\stackrel{}{B}\times \stackrel{}{k}\end{array}\right)$$
(39)
and

$$q_\mu \mathrm{\Lambda }^{\mu \nu }k_\nu =\stackrel{}{E}\left(q_0\stackrel{}{k}k_0\stackrel{}{q}\right)+\stackrel{}{B}\left(\stackrel{}{k}\times \stackrel{}{q}\right).$$
(40)
Eqs. (38) and (39) imply that the three-momentum is conserved along the direction of $`\stackrel{}{E}`$. Energy is conserved if the total three-momentum $`\stackrel{}{k_i}`$ is along the direction of $`\stackrel{}{B}`$. Also,

$$k_1\mathrm{\Lambda }k_2=k_1\mathrm{\Lambda }k_3=k_2\mathrm{\Lambda }k_3$$
(41)
and one is required only to study eq. (33) with $`V(k_1,k_2)`$ and $`V(k_1,k_2)`$. Note that the second case is obtained by performing $`\stackrel{}{E},\stackrel{}{B}\stackrel{}{E},\stackrel{}{B}`$. Thus, once computations have been performed in the first case, the results in the second one are obtained by performing this substitution.

While performing calculations, the following dimensionless combinations of the masses and the noncommutative parameters arise:

$$x_i=\frac{1}{2}E\left|\stackrel{}{r}\right|m_i,y_i=\frac{1}{2}Bs_\delta \left|\stackrel{}{r}\right|m_i.$$
(42)
where the notation $`s_\omega =\mathrm{sin}\omega `$ is used.

### 4.2 Massive particle decay

Consider now the decay of a massive particle $`\mathrm{\Phi }_3`$ into two particles $`\mathrm{\Phi }_1`$ and $`\mathrm{\Phi }_2`$, that is

$$\mathrm{\Phi }_3\mathrm{\Phi }_1+\mathrm{\Phi }_2.$$
(43)
Let $`m_i`$ be the mass of particle $`\mathrm{\Phi }_i`$ and $`m_1m_2`$. In the rest-frame of $`\mathrm{\Phi }_3`$ the angle $`\alpha `$ between particles $`\mathrm{\Phi }_1`$ and $`\mathrm{\Phi }_2`$ is given by

$$\mathrm{cos}\alpha =\frac{1+x_3c_\theta c_\phi }{\sqrt{(1+x_3c_\theta c_\phi )^2+(x_3c_\theta s_\phi )^2}}$$
(44)
where $`\theta `$ is the angle between $`\stackrel{}{p}_1`$ and $`\stackrel{}{E}`$, and $`\phi `$ the angle between $`\stackrel{}{p}_1`$ and $`\stackrel{}{r}`$. Also, the notation $`c_\omega =\mathrm{cos}\omega `$ is used. The absolute value of the right-hand side of this equation is always smaller than one, meaning that this decay is always possible. In the high energy regime ($`x_31`$) the variable $`x_3`$ decouples and one obtains $`\alpha \pm \phi `$ or $`\alpha \pi \pm \phi `$. In the low-energy regime ($`x_31`$) one finds the first-order correction $`\alpha \pi \pm x_3s_\phi c_\theta `$.

The new equation for the energy is

$$m_3=ϵ_1\left(1y_3c_\theta v_1\right)+ϵ_2,$$
(45)
where we have used $`ϵ_i,v_i`$ as the energy and velocity of particle $`i`$.

If $`\left|y_i\right|<1`$ then spontaneous decay will occur if $`m_3`$ satisfies the condition

$$m_3>\frac{m_1+m_2}{1+y_1c_\theta v_1}.$$
(46)
If the denominator is zero or negative, then the decay is impossible. We can see that if the velocity of particle 1 is high and mainly along the direction of $`\stackrel{}{E}`$ then the decay can occur for a value of $`m_3`$ smaller than the sum of the masses $`m_1+m_2`$

$$m_3>\frac{m_1+m_2}{1+|y_1|}.$$
(47)
### 4.3 Massless particle decay

The decay of a massless particle into two massive particles is kinematically forbidden. However, in the present model, all massless particles become unstable and can decay into two massive particles. If the particle $`\mathrm{\Phi }_3`$ is massless we have the following condition for the energy of the photon in the rest frame of particle $`\mathrm{\Phi }_1`$

$$\omega >\frac{2}{|Bs_\delta Ec_\phi ||\stackrel{}{r}||c_\theta |}\frac{\omega _0(\phi )}{|c_\theta |},$$
(48)
where we have defined $`\theta `$ as the angle between $`\stackrel{}{E}`$ and $`\stackrel{}{\omega }`$ and also $`\phi `$ as the angle between $`\stackrel{}{r}`$ and $`\stackrel{}{\omega }`$. If the denominator is zero, the decay is impossible. Note that this limit is independent of the mass of the decaying particles. This result is only valid in the low-energy limit where $`\left|x_1\right|<1`$, that is, when the decay produces particles with low mass.

The above limit for angle $`\theta `$ implies that

$$c_\theta ^2>\left(\frac{\omega _0}{\omega }\right)^2,$$
(49)
which states that, as the energy of the photon grows larger, the decay is possible for a wider range of $`\theta `$ angles. If $`\theta =0`$, the decay is impossible.

### 4.4 The GZK cutoff

The Greisen-Zatsepin-Kuz’min (GZK) cutoff mechanism asserts that ultra-high-energy (UHE) protons with energies $`ϵ_p>4\times 10^{19}eV`$ from sources beyond $`50100Mpc`$ should not be observed, due to their interaction with Cosmic Microwave Background (CMB) photons. It has been proposed (for brief review see ref. ) that Lorentz-violating terms in the kinematics of hadronic reactions may be the answer to this puzzle. The GZK cutoff has the following dominant resonance

$$p+\gamma _{CMB}\mathrm{\Delta }_{1232}.$$
(50)
It is easily shown that the model eq. (35) does not account for a violation of the GZK cutoff. In fact, in the case of head-on collision, the new equation for the energy yields

$$ϵ_p[1y_\omega c_\theta (1+v_p)]+\omega =ϵ_\mathrm{\Delta },$$
(51)
where $`ϵ_i,v_i`$ denote the energy and velocity of particle $`i`$ and $`\theta `$ is the angle between $`\stackrel{}{E}`$ and $`\stackrel{}{\omega }`$. In order to occur any appreciable deviation that renders this reaction impossible (for instance $`ϵ_\mathrm{\Delta }<m_\mathrm{\Delta }`$), one should have $`y_\omega 1`$, which, given the low energy of the $`CMB`$ photon, would yield a very small mass for the noncommutative parameters.

Nevertheless, the violation of the GZK limit may be explained in the context of the model

$$\mathrm{\Lambda }_{\omega _1}^{\mu _1\nu }\mathrm{\Lambda }_{\omega _2}^{\mu _2\omega _1}\mathrm{\Lambda }_{\omega _3}^{\mu _3\omega _2}=0,$$
(52)
with $`\mathrm{\Lambda }_{\omega _1}^{\mu _1\nu }\mathrm{\Lambda }_{\omega _2}^{\mu _2\omega _1}0`$. This cannot be implemented by model eq. (35), which complicates the analysis. The equation for the momentum is given by

$$\left(k_1+k_2+k_3\right)_\omega =k_{1\mu }k_{2\nu }\mathrm{\Lambda }_\lambda ^{\mu \nu }\left[\frac{1}{2}\delta _\omega ^\lambda +\frac{\left(k_1k_2\right)_\alpha }{12}\mathrm{\Lambda }_\omega ^{\alpha \lambda }\right]$$
(53)
where we have used eq. (31) recursively and the fact that cubic terms in $`\mathrm{\Lambda }`$ vanish.

This condition can be modeled by a simpler one, more suitable for phenomenological considerations which, however, breaks Lorentz invariance. As we have seen, the quadratic term in the momentum does not account for the violation of the GZK cutoff, so it will be dropped. Taking into account that the proton has the highest energy and the $`\mathrm{\Delta }`$ the second highest energy, we can write the new momentum equation for the reaction (50) as

$$\left(k_p+k_\gamma \right)^\mu =k_\mathrm{\Delta }^\mu s^\mu \frac{ϵ_p^2}{M^2}ϵ_\mathrm{\Delta }$$
(54)
where adimensional vector $`s^\mu `$ is of the order of unity and $`M`$ is the typical noncommutative mass scale. In this case, the process is impossible if $`s^0>0`$ and $`ϵ_p>M`$, which sets the scale of noncommutativity.

Note that one must consider all permutations of the indices in eq. (53). Due to the low energy of CMB photons, the permutations that lead to a term of the type $`ϵ_\gamma ^2`$ will not violate the GZK cutoff. Since there are six permutations of the indices and only two lead to this type of term, we can estimate that $`2/3`$ of the events leading to the ressonance (50) will violate the GZK cutoff.

It is generally believed  that a cubic term in the equations of dispersion will explain the violation of this cutoff. In fact, eq. (54) can be obtained by assuming the usual momentum conservation and postulating a new equation of dispersion by the substitution

$$k^\mu k^\mu +s^\mu \frac{ϵ^2}{M^2}\lambda $$
(55)
where $`\lambda `$ represents the typical energy of the product of the reaction. This will lead to the following dispersion relation

$$m^2=ϵ^2p^2+2s^\mu v_\mu \frac{\lambda }{M^2}ϵ^3$$
(56)
where $`v^\mu =(1,\stackrel{}{v})`$ is the four-vector velocity, which we assume to be nearly light-speed. Only the lower order terms of the correction were kept.

Thus, it is as if a cubic term is added to the dispersion relation. However, this model differs from the one of ref.  in the sense that eq. (56) is sensible to the typical energy of the product of the reaction, that is, to the process in question. Also, there is a dependence on the geometry of the propagation of the particle, through the term $`s^\mu v_\mu `$. In addition, since the free theory is not altered by our approach, this effective dispersion equation may only be used to study particle reactions and not classical free-particle propagation.

## 5 Quantum field theory

To study the quantum aspects of the noncommutative model discussed in sections II and III we consider the $`\lambda \mathrm{\Phi }^4`$ theory. The action is given by

$$S=_{}d^nx\left[\frac{1}{2}_\mu \mathrm{\Phi }^\mu \mathrm{\Phi }+\frac{m^2}{2}\mathrm{\Phi }\mathrm{\Phi }+\frac{\lambda }{4!}\mathrm{\Phi }\mathrm{\Phi }\mathrm{\Phi }\mathrm{\Phi }\right]$$
(57)
and can be evaluated in Fourier space to yield the corresponding Feynman rules. The vertices are already calculated in eq. (25) and the propagators are the same as the commutative ones. We shall consider the euclidean formulation.

The free propagator obeys the usual energy-momentum conservation. However, since interactions do not, it is expected that the quantum corrections to the propagator due to the self-interaction term in the action will not obey the usual momentum conservation. Therefore, it is of interest to compute the first order correction to the two-point function, which is proportional to

$$\mathrm{\Gamma }_1(p_i,p_f)=\frac{1}{2}\frac{\lambda }{4!}\frac{1}{p_i^2+m^2}\frac{1}{p_f^2+m^2}I(p_i,p_f)$$
(58)
with

$$I(p_i,p_f)=\frac{d^nq}{\left(2\pi \right)^n}\frac{\underset{\pi }{}\delta \left[p_i+p_f+V^4(p_i,p_f,q,q)\right]}{q^2+m^2},$$
(59)
where the sum is computed over all $`\pi `$ permutations of the arguments of $`V^4`$.

We proceed by evaluating this integral in the simple model eq. (32). From the several contributions to the above integral, there are two which differ from the commutative case. They are formally identical, and the relevant integrals are given by

$$J(p_i,p_f)=\frac{d^nq}{\left(2\pi \right)^n}\frac{1}{q^2+m^2}\delta [p_i+p_f+2V(q,k)]$$
(60)
where $`k=p_i,p_f`$ for each case.

Using a Schwinger parametrization and the usual Fourier representation for the delta-function, the integrals are gaussian and yield

$$J(p_i,p_f)=\frac{1}{\left(4\pi \right)^n\sqrt{detN}}\frac{1}{\left(p_i+p_f\right)N^1\left(p_i+p_f\right)+m^2},$$
(61)
where

$$N_{\omega \lambda }\left(k\right)=k_\nu \mathrm{\Lambda }_\omega ^{\mu \nu }k_\beta \mathrm{\Lambda }_{\mu \lambda }^\beta .$$
(62)
Notice that this matrix is singular for model eq. (35), but not in general. The fact that $`\mathrm{\Lambda }`$ is nilpotent is not important, as $`\mathrm{\Lambda }_\alpha ^{\mu \nu }`$ is nilpotent only regarding the indices $`\mu `$ and $`\alpha `$, and $`N`$ involves only $`\alpha `$ type indices. Also, matrix $`N`$ is singular if $`k=p_i=0`$ or $`k=p_f=0`$ and, hence, there is an IR divergence, which is usual in noncommutative quantum field theories.

The integral

$$J_C(p_i,p_f)=\delta \left(p_i+p_f\right)\frac{d^nq}{\left(2\pi \right)^n}\frac{1}{q^2+m^2},$$
(63)
which is ultraviolet divergent for $`n=4`$, also arises from eq. (59). Thus one concludes that its regularization is still required, although this is not necessary in eq. (60). Hence, the UV renormalizability properties at the one-loop approximation are not altered.

From the correction to the two-point function in eq. (61) two interesting features arise. First, the conservation of momentum is lost, since the delta function $`\delta \left(p_i+p_f\right)`$ is no longer present. Second, the correction to the dispersion relation, which is given by the pole of eq. (61), manifests itself in a quite specific way, involving the Lorentz algebra of the momentum vectors and matrix $`N^1`$. In fact, the poles of eq. (61) suggest the particle is subjected to a momentum-dependent metric given by $`N^1`$.

## 6 Discussion and conclusions

In this work we have presented a noncommutative field theory where the coordinates have a Lie-algebra commutator as eq. (1) with nilpotent structure constants. This breaks Lorentz as well as translational invariance. Free theory is unchanged so the propagators and the dispersion relations are not altered. The vertices show a new energy-momentum law, which steems from the breaking of translational invariance. The kinematical studies of such law where established in particle decay physics and shown how it can be applied as a possible explanation for the violation of the GZK cutoff, setting the noncommutative mass scale at $`M4\times 10^{19}eV`$. A link between these kinematics and Lorentz-violating theories was established, using a simplified model. However, there are well definite differences between our approach and the ones usual discussed in the literature (see ref. ), the most important one being that an effective dispersion law always depends on the energy and geometry of the processes in question.

It is tempting to speculate that our results have a bearing on the other known astrophysical puzzles, namely the observation of high energy photons, $`ϵ20TeV`$, from far away sources and the pion stability in extensive air showers (see ref.  and references therein). Indeed, since both phenomena can be understood via a cubic deformation in the relativistic dispersion relation, so at pair creation through the process $`\gamma +\gamma _{IRbacground}e^++e^{}`$ cannot occur and pion decay into photons has a smaller width, it is plausible to assume that these paradoxes can be explained in our model as well.

In the context of quantum field theory, it was shown that it is possible to carry out explicit calculations regarding the first-order correction to the two-point function in $`\lambda \mathrm{\Phi }^4`$ theory. The interaction terms violate momentum conservation and this is expressed in the two-point function, where the usual delta function structure $`\delta \left(p_i+p_f\right)`$ is lost. Even though the noncommutative contributions are UV finite, usual commutative integrals are still present and are UV divergent. Thus, the UV renormalization properties of one-loop calculations are not altered. Also, in a strict sense, it is shown that free theory is unchanged and so the propagators and the dispersion relations are not altered. New IR divergences arise in the noncommutative corrections, a feature which is shared with constant commutator noncommutative field theories, known as UV/IR mixing. The poles of the noncommutative terms indicate that there is a correction to the dispertion relation, through the Lorentz algebra of matrix $`N`$, eq. (62), which seems to indicate that the particle satifies a dispersion relation arising from a momentum-dependent metric, eq. (62).

###### Acknowledgments.

The authors would like thank Justin Conroy, Daniel Robbins, Naoki Sasakura, Ricardo Schiappa and Edward Witten for their comments and constructive remarks.


text:
# Reaching large lengths and long times in polymer dynamics simulations

## I Introduction

Lattice models are widely used as a theoretic tool to study polymer solutions. Most notably, Flory Flory (1942, 1944) and Huggins Huggins (1942a, b) introduced a mean-field theory of a particular simple model which has proven powerful enough to explain solubility properties and liquid-liquid phase separation of polymer mixtures. An extension to the model, describing the dynamics of such systems, has been introduced by Cahn and Hilliard Cahn and Hilliard (1958); Cahn (1959) and Hillert Hillert (1961). Some of the earliest computer simulations of lattice polymers were performed in 1962 by Verdier and Stockmayer Verdier and Stockmayer (1962), who studied the dynamics of polymers, modeled by self-avoiding walks on square and cubic lattices. Kolinski et al. Kolinski et al. (1987a, b) added two-monomer moves to the dynamics and simulated a homopolymeric melt.

To describe the polymer dynamics in a melt with more detail than the random-walk approach, Carmesin introduced in 1988 the bond fluctuation model Carmesin and Kremer (1988). In the two-dimensional version of this model, each monomer occupies four ($`2\times 2`$) lattice sites of a square lattice. Multiple occupation of lattice sites is not permitted. Monomers adjacent in the chain are connected by bonds with lengths between 2 and $`\sqrt{13}`$. The dynamics of the polymer chain consists of the displacement of single monomers to nearest-neighbor lattice sites, restricted by the constraints on bond length and excluded volume. The bond-fluctuation model was extended to three dimensions in 1991, by Deutsch and Binder Deutsch and Binder (1991). This model has been used by several groups Paul et al. (1991); Kreer et al. (2001) to simulate dense polymer melts.

In 1971, de Gennes de Gennes (1971) proposed that the main mechanism of polymer dynamics in gels and dense polymer mixtures is reptation, i.e., diffusion of stored length. To verify the theoretical predictions of de Gennes, Evans and Edwards Evans and Edwards (1981) introduced the cage model. in this model, a single polymer is simulated as a random walk on a square or cubic lattice. In the limit of tight gels, the dynamics reduces to the re-orientation of pairs of connected segments pointing in opposite direction (“kinks”). The cage model has been used for the study of diffusion and relaxation times of a single polymer in a gel Evans and Edwards (1981); Deutsch and Madden (1989a, b); Barkema and Krenzlin (1998), as well as for star polymers Barkema and Baumgärtner (1999). Extending the model to include excluded-volume effects or interactions between polymers is not easy.

An alternative model for a single reptating polymer was proposed by Rubinstein Rubinstein (1987). In this model, known as the repton model, the polymer is represented as a chain of monomers, residing on the sites of a square or cubic lattice. Monomers adjacent along the chain reside in either nearest-neighbor lattice sites, or in the same site. The latter case corresponds to the presence of stored length, a central concept in de Gennes’ description of entangled polymer dynamics. The dynamics of the polymer is strictly limited to reptation: the only kind of moves allowed to a monomer in the interior of the chain is that it can move to an adjacent lattice site, provided that one of its adjacent neighbors is already on the site to which it is moving, and the other is in the site it leaves. Per unit of time, each monomer attempts to move in each direction statistically once. The repton model is illustrated in Figure 1.

The repton model has been studied extensively by one of us Newman and Barkema (1997); Barkema and Newman (1997).

The repton model can be extended to simulate polymers with self-avoiding walk statistics, as well as many interacting polymers; this is the route followed in this manuscript. In the next section, we discuss the extensions to the repton model which are required to describe equilibrium and dynamical properties of dense polymer solutions and melts. Since the dynamics of polymeric systems occurs on long time-scales and large length-scales, computational efficiency is of paramount importance. We therefore discuss in detail how the resulting extended repton model can be implemented with high computational efficiency, exploiting a technique known as multispin coding. At the end of this manuscript, we show the usefulness of our work by presenting simulation results on the phase separation of a dense binary polymer mixture.

## II Extensions of the repton model to simulate polymer melts

The repton model as proposed by Rubinstein simulates a single polymer obeying random-walk statistics, with dynamics limited to the diffusion of stored length. Two extensions to the repton model are needed to use it for the simulation of entangled chains: the chains have to obey volume exclusion and on the other hand, the constraint that disallows sideways motion has to be lifted.

The repton model describes “stored length” as two or more monomers on the same lattice site. Therefore, we cannot limit the polymers to self-avoiding-walk statistics by disallowing two monomers to occupy the same lattice site, without destroying the reptation dynamics that so successfully described polymer-diffusion in a gel. As reptation is also expected to play an important role in the melt, we should keep this kind of dynamics. The solution to this problem is to limit only the contours of the polymers to self-avoiding-walk statistics. Thus, multiple occupation of a lattice site is allowed only for two or more monomers which are adjacent along the chain. One convenient side-effect of this choice is that since the reptation moves in the interior of the chain do not affect the tube, they never cause a violation of volume exclusion.

Sideways motion is implemented by also allowing monomers to move by a single lattice spacing if this changes the tube of the polymer. Single-monomer moves work well on lattice structures that contain loops of three sites, like the triangular and face-centered-cubic (FCC) lattice; on lattices without such three-site loops, like the square and cubic lattices, the sideways movement can lead to lots of unwanted artifacts or, alternatively, must be implemented by allowing bonds to extend to next-nearest neighbors.

The repton model is often studied in its projected version, in which the state of the polymer is characterized by the bonds $`s_i`$, $`i=1\mathrm{}N1`$, between the monomers $`i=1\mathrm{}N`$, where $`s_ix_{i+1}x_i`$ is equal to zero or $`\pm 1`$. Usually all allowed moves are attempted with the same (unit) rate, with the consequence that the three values for $`s_i`$ occur with equal probability. The density of stored length (probability that $`s_i`$ is zero) is then 1/3. To obtain the same density of stored length in higher-dimensional repton models (before projection), the rates for moves which decrease the stored length should equal $`2/z`$ times those of the reverse moves, in which $`z`$ is the lattice coordination number: 4, 6 and 12 for the square, cubic or triangular, and FCC lattices, respectively. To conserve equilibrium properties, sideways moves which decrease the density of stored length should also be attempted with rates which are $`2/z`$ times the reverse moves.

Figure 2

shows model-polymers on a triangular lattice. In the upper polymer in the figure, reptation moves are possible for monomers 2, 4, 6, 9, 10, and 11; sideways moves can be made by monomers 2, 4, 6, 8, 9, 10, and 11. Monomers 3, 5, and 7 cannot move. The end-monomers 1 and 12 can move to any empty nearest-neighbor site. In the lower polymer, reptation moves are possible for monomers 3, 5, 6, 10, and 11; sideways moves can be made by monomers 3, 4, 5, 6, 7, 10, and 11. Monomers 2, 8, and 9 cannot move. The end-monomer 1 can move to any empty nearest-neighbor site, while end-monomer 12 can only move to the site occupied by monomer 11.

## III multispin coding

Multispin coding is a programming technique which makes use of the low-level bit manipulation functions provided by the CPU of our computer, effectively doing many simple computations in parallel. These functions include logical operations like AND ($``$), OR ($``$) exclusive-or (XOR, $``$), negation (NOT, $`\neg `$), as well as bitwise shifts ($``$, integer multiplication by a power of 2; $``$, integer division by a power of 2); and arithmetic operations like addition and multiplication. For a general introduction to multispin coding we refer to Ref. Newman and Barkema (1999), Chapter 15. Here, we will explain this technique analogous to Ref. Barkema et al. (1994), as applied to the projected repton model. Next we will discuss the implementation of this technique in our model for polymer melts.

### III.1 repton model

In the projected repton model, the state of the chain with monomers $`i=1\mathrm{}N`$ is specified by the set of their coordinates $`\{x_1,\mathrm{},x_N\}`$. Monomers which are adjacent along the chain have to reside in either the same lattice site, or nearest neighbors, with the consequence that the step $`s_ix_{i+1}x_i`$ can only take the values $`\pm 1`$ or 0. An alternative way to store the polymer state is therefore to store the coordinate $`x_1`$ of the first monomer, and to store for step $`i=1\mathrm{}N1`$ two arrays of bits $`(A_i,B_i)`$ that can take the combination $`(0,0)`$ if $`s_i=0`$, $`(0,1)`$ if $`s_i=1`$, and $`(1,0)`$ if $`s_i=1`$.

In one unit of time, each monomer attempts to move in each direction statistically once. Thus, $`2N`$ elementary moves, i.e. one monomer attempting to move in one direction, should be attempted per unit of time.

It is a property of the repton model, that if monomer $`i`$ in the interior of the chain can join its neighbor $`i+1`$, it cannot join its neighbor $`i1`$, and vice versa. We exploit this by attempting these two elementary moves simultaneously. The first and last monomer can move to both directions, if they are located in the same site as their neighbor along the chain, so this trick does not apply for these two monomers; they are therefore selected with a double probability.

In detail, each monomer in the interior of the chain is selected with a probability $`1/(N+2)`$. If the selected monomer $`i`$ is located in the same site as exactly one of its two adjacent neighbors, it will join the other neighbor. The bits that contain information on whether a move is possible or not, are $`(A_i,B_i)`$ and $`(A_j,B_j)`$ with $`j=i1`$. A move is possible if one of the two pairs is $`(0,0)`$ and the other is not. If the move is carried out, $`(A_i,B_i)`$ and $`(A_j,B_j)`$ are exchanged. A sequence of logical operations that achieves this is:

$`\begin{array}{ccc}\hfill y& =& (A_iB_i)(A_jB_j)\hfill \\ \hfill m_A& =& (A_iA_j)y\hfill \\ \hfill m_B& =& (B_iB_j)y\hfill \\ \hfill A_i^{}& =& A_im_A\hfill \\ \hfill B_i^{}& =& B_im_B\hfill \\ \hfill A_j^{}& =& A_jm_A\hfill \\ \hfill B_j^{}& =& B_jm_B\hfill \end{array}`$ (8)
The first and last monomers are selected, each, with a probability $`2/(N+2)`$; twice the probability to select a specific monomer in the interior of the chain. After the selection of an end monomer, the intended direction is also randomly selected, with 50% probability to be in the positive or negative direction. The rates at which elementary moves are attempted is thus equal for interior monomers and those at the ends. If the first or last monomer is selected, similar statements suffice to update $`(A_1,B_1)`$ or $`(A_{N1},B_{N1})`$, respectively. The following statements can be used to update the end-monomer 1, trying to displace it in the negative or positive direction based on the value of $`r_1`$.

$`\begin{array}{ccc}\hfill A_1^{}& =& \neg (B_1r_1)\hfill \\ \hfill B_1^{}& =& (\neg A_1)r_1\hfill \end{array}`$ (11)
However, in case the first monomer is selected, also its coordinate $`x_1`$ needs to be updated. Since this coordinate is an integer number that can take a wide range of values, its update is not implemented in a multispin coding fashion. Luckily, the probability to select the first monomer decreases with increasing polymer length.

The motivation for multispin coding lies in its efficiency. The simulation of the dynamics of 64 polymers involves $`128N`$ elementary moves per time unit. With multispin coding, it involves that the above sequence of logical operations for interior monomers should be performed statistically $`N2`$ times, which involves only 11 logical operations, 4 loads, and 4 stores. On a fast workstation, this takes 41 ns of cpu-time, or 0.32 ns per elementary move. For updating the last monomers, even less logical operations suffice. For the first monomers, however, the update of the 64 values of $`x_1`$ cannot be achieved as efficiently, and a loop over all 64 polymers is inevitable, with the consequence that the simulation of the first monomer requires 5.8 ns cpu-time per elementary move.

The multispin implementation outlined above will thus perform 64 simulations in parallel. These simulations are however correlated, since they share the sequence of selected monomers. In fact, if at some point two simulations are in identical polymer configurations, they will stay in identical configurations ever after; in long simulations of small systems, this “locking” will inevitably happen. Complete locking is avoided by an uncorrelated choice in the directions in which end monomers attempt to move, using a 64-bit random bit pattern rather than a binary choice between all-up or all-down. If desired, more de-locking can be obtained at the expense of a lower efficiency, by rejecting a fraction of the allowed moves, also using a 64-bit random bit pattern. An important remark is, however, that as long as each of the 64 simulations are correct in themselves, one obtains 64 unbiased results; one should just be careful to assign significance to the spread in those 64 simulations.

### III.2 model for polymer melts

While monomers in the projected repton model live on a one-dimensional lattice, the monomers in the model that we propose for polymer melts live on a FCC lattice. It is helpful to note that the three-dimensional hyperplane, located in a four-dimensional hypercubic space through the origin and with perpendicular vector $`(1,1,1,1)`$, is such an FCC lattice. Stated differently, the set of points $`\stackrel{}{x}=(a,b,c,d)`$ with integer-valued coordinates, constrained to $`a+b+c+d=0`$, forms a FCC lattice. The twelve vectors pointing to nearest-neighbor sites are $`\widehat{t}=(1,1,0,0)`$, $`\widehat{u}=(0,1,1,0)`$, $`\widehat{v}=(0,0,1,1)`$, $`\widehat{w}=(\widehat{t}+\widehat{u}+\widehat{v})=(1,0,0,1)`$, and some of their combinations, as listed in Table 1.

The vector $`(0,0,0,0)`$ is used as the representation for a zero-length bond.

As in the projected repton model, the state of a polymer on our FCC lattice can be specified by the location $`\stackrel{}{x}_1`$ of the first monomer, plus additionally the direction in which adjacent neighbors are located. We choose for the polymer melt model to store the four directional bits in a single word: bits $`k`$, $`k+16,k+32`$ and $`k+48`$ of the 64-bit word $`D_i`$ indicate the vector pointing from monomer $`i`$ to monomer $`i+1`$. Thus, sixteen polymers are updated simultaneously.

The coordinates of the $`i^{\mathrm{th}}`$ monomer in polymer $`k`$ can then be retrieved by summing over all words $`D_i`$ bits $`k`$, $`k+16`$, $`k+32`$ and $`k+48`$, yielding respectively the numbers $`a`$, $`b`$, $`c`$ and $`d`$; the monomer position is then $`\stackrel{}{x}_i=a\widehat{t}+b\widehat{u}+c\widehat{v}+d\widehat{w}`$. Note that these summations require only $`3i`$ operations (a right shift of $`k`$ bits, masking the direction bits, and adding to the sum), since the summation in the different bits can be done in a single operation. Since in the polymer melt model we need the coordinates frequently (every time we attempt a possible sideways move), and since our polymers are often several hundred monomers long, we do not keep track explicitly of only the position $`\stackrel{}{x}_1`$ of the first monomer, but also of the last monomer as well a few other monomers along the chain, such that the distance along the chain to a monomer with known position is always less than 15. Of course, to retrieve the position $`\stackrel{}{x}_i`$, we start from the nearest monomer with known position in either direction along the polymer.

The implementation proceeds analogous to the projected repton model. Also here, if an interior monomer can move in one direction, its move in the other direction is blocked; this can be exploited as in the repton model, by combining two elementary moves. The precise operations are:

$`\begin{array}{ccc}\hfill y_0& =& D_i(D_i32)\hfill \\ \hfill y_1& =& y_0(y_016)\hfill \\ \hfill z_0& =& D_j(D_j32)\hfill \\ \hfill z_1& =& z_0(z_016)\hfill \\ \hfill m_0& =& (y_1z_1)(2^{16}1)\hfill \\ \hfill m_1& =& (D_iD_j)(Mm_0)\hfill \\ \hfill D_i^{}& =& D_im_1\hfill \\ \hfill D_j^{}& =& D_jm_1\hfill \end{array}`$ (20)
Here, $`Ak`$ stands for shifting the word $`A`$ over $`k`$ bits to the left, and $`M=2^0+2^{16}+2^{32}+2^{48}`$ is a constant, used to duplicate the low 16 bits in the higher bits of the word. Thus, with only 15 operations, 2 loads, and 2 stores, we have performed 32 elementary moves. On a fast workstation, the above implementation requires 1.25 ns cpu-time per elementary move.

As in the projected repton model, if the first monomers are displaced, whose positions are tracked, these positions have to be updated in a loop over the 16 polymers. We succeeded in implementing moves of the first monomers in 82 ns cpu-time per elementary move. Other monomers whose position is tracked require roughly the same computational effort. Displacement of the first and last monomers is attempted at twice the rate of the other monomers, for the same reason as in the projected repton model.

Besides the reptation moves, the dynamics consists of sideways moves. If a sideways move is attempted on monomer $`i`$, the relevant bit patterns are those indicating the direction from monomer $`i1`$ to $`i`$ and from $`i`$ to $`i+1`$. These bit patterns, as listed in Table 1, are all numbers in the range 0 to 15, except for the values 5, 10, and 15 that do not occur; the number 0 denotes stored length, while the other 12 numbers denote bonds to the twelve nearest-neighbor sites. For every combination $`(D_i,D_j)`$ of those bit patterns, we have pre-computed the lists of all bit patterns $`D_i^{}`$ and $`D_j^{}`$ after a sideways move. Depending on the combination $`(D_i,D_j)`$, at most four different sideways moves can be proposed. We therefore have precomputed four such lists for $`D_i^{}`$ and for $`D_j^{}`$. In one step, we first select randomly the monomer number $`i`$ and the list number $`k`$; then we attempt a sideways move of monomer $`i`$ to the position determined by the $`k^{\mathrm{th}}`$ list; finally, if this move does not lead to overlapping monomers, it is accepted. This check for overlap requires computing the position of monomers, which requires computing the distance to the nearest tracked monomer. In our implementation, in which this distance is at most 15 monomers, the total cpu-time required per such step equals 90 ns.

Not all moves are attempted with the same frequency. Since the long-time dynamics is determined by reptation, as argued by de Gennes, the time scale is set such that reptation moves in the interior of the polymers are attempted with unit rate. Sideways moves which do not increase the amount of stored length are attempted with some rate $`r_s`$; the most natural choice for this rate would be unit rate once more, but because of the much higher computational cost for these moves, we often chose some value of $`r_s<1`$. As discussed above, consistency with the density of stored length of the projected repton model demands that moves in which the amount of stored length is increased are attempted with a rate of $`2r_s/z`$. Since the most mobile monomers hop away to other sites with a total rate of $`2r_s`$, and since 16 monomers are potentially moved in each step, $`r_s/8`$ steps as described above should be performed per monomer and per unit of time.

Moves in which the first or last monomer joins its neighbor along the chain can be viewed partly as reptation moves, partly as sideways movement. On these grounds, we have decided to attempt these moves with rate $`1+r_s`$. Consistency with respect to the density of stored length requires that moves in which the first or last monomer leaves its neighbor along the chain in a specific direction are attempted with a rate of $`2(1+r_s)/z`$.

## IV Application: phase separation of a binary polymer mixture

To illustrate the efficiency of the above computational approach, we performed a simulation of the phase separation of a binary polymer mixture with polymer types $`A`$ and $`B`$. The $`A`$ and $`B`$ polymers interact with a short-range repulsion, described by the Hamiltonian

$$H=J\underset{\stackrel{}{r},\stackrel{}{r}^{}}{}\delta (\sigma _\stackrel{}{r},A)\delta (\sigma _\stackrel{}{r}^{},B)+\delta (\sigma _\stackrel{}{r},B)\delta (\sigma _\stackrel{}{r}^{},A)\text{,}$$
(21)
where the summation runs over all pairs of nearest-neighbor sites, and $`\sigma (\stackrel{}{r})`$ is $`A`$, $`B`$, or 0 if site $`\stackrel{}{r}`$ is occupied by a polymer of type $`A`$, $`B`$ or empty, respectively. The repulsion between $`A`$ and $`B`$ polymers provides the driving force for the phase separation.

We simulated a system containing in total $`46,080`$ polymers of length $`L=100`$ on a lattice of $`N=13,824,000`$ sites, at inverse temperature $`\beta J=0.1`$. The system evolves in time through reptation moves, at unit rate, in combination with sideways moves at a rate of $`r_s=1/30`$. Figure 3 shows two-dimensional slices of the three-dimensional system at times $`t=0`$, $`t=4.710^5`$, and $`t=3.710^6`$ MC time units.

At various times, we determine the two-point correlation function

$$g_{AB}(\stackrel{}{r})=\frac{\frac{1}{N}_\stackrel{}{r}^{}\delta (\sigma _\stackrel{}{r}^{},A)\delta (\sigma _{\stackrel{}{r}^{}+\stackrel{}{r}},B)}{\left(\frac{1}{N}_\stackrel{}{r}^{}\delta (\sigma _\stackrel{}{r}^{},A)\right)\left(\frac{1}{N}_\stackrel{}{r}^{}\delta (\sigma _\stackrel{}{r}^{},B)\right)}\text{.}$$
(22)
From this function we determine the spherically-averaged radial distribution function, defined as $`RDF(r)=1g(\stackrel{}{r})`$. This function is 1 at $`r=0`$, then decreases, and eventually approaches 0 for large $`r`$. After some time, the conserved dynamics gives rise to damped oscillations in the RDF. The frequency of these oscillations can be determined from the shortest distance $`r_0`$ at which the RDF equals zero. The typical domain size $`d(t)`$ is then obtained as twice this distance.

Figure 4

shows the domain size of the phase-separating mixture at different times. The domain size $`d(t)`$, as defined above, is plotted as a function of the cubic root of the simulation time in Monte Carlo sweeps. The straight line shows that the domain size grows with $`d(t)t^{1/3}`$. This is to be expected in a system with conserved order dynamics and without hydrodynamics, as in “Model B” Hohenberg and Halfein (1977).

## V Summary and outlook

In summary, we have presented a lattice model to study the dynamics of polymeric systems. Compared to other lattice polymer models in current use, this model lends itself for highly efficient simulation: employing multispin coding, the computer time required per elementary move is reduced by more than three orders of magnitude, to a few nanoseconds.

To demonstrate the strength of our model, we simulated the phase separation of a binary mixture of polymers. We performed simulations of a system with about $`50,000`$ polymers, each containing 100 monomers, located on a lattice with about 14 million sites. We verified that over the final two decades in time, the domain size $`d(t)`$ grows according to $`d(t)t^{1/3}`$, as expected for a system with overdamped dynamics and a local conservation law. This simulation involves $`3.410^{13}`$ elementary moves, but could nevertheless be carried out on a single-processor workstation in about 20 days.

Recently, we have used this lattice polymer model to study fractionation in quasi-binary (polydisperse) polymer mixtures van Heukelum et al. (2003). Currently, we are using the model presented here to study the dynamics of absorption of a polydisperse polymer mixture on a surface, as well as the sieving process of polydisperse polymer mixtures through nanopores.


text:
# Transition matrix Monte Carlo and flat-histogram algorithm

## 1 Introduction

In traditional Monte Carlo sampling method, the computation of a thermodynamic quantity $`Q`$ is usually through a simple arithmetic average:

$$Q=\frac{1}{M}\underset{i=1}{\overset{M}{}}Q(\sigma _i),$$
(1)
where the configurations $`\sigma _i`$ are generated according to a specified distribution, such as the Boltzmann distribution. However, it is possible to collect other information in the same simulation, from which we can obtain better statistics, or estimates of quantities other than that at simulation parameters. The histogram method  and multi-histogram method  collect energy histogram at a given temperature, from which the quantity at nearby temperature can be inferred. The key observation here is that the histogram of energy is related to density of states through $`H(E)n(E)\mathrm{exp}(E/kT)`$ (in canonical ensemble). From the histogram, we can determine the density of states $`n(E)`$. Once the density of states is known, we can compute most of the thermodynamic quantities at any temperature.

Histogram method has been found to be an excellent tool for study critical phenomena. Further improvement can be made by collecting ‘high-order’ statistics, i.e., the transition matrix . With histogram method, each configuration provides just an ‘1’ to a histogram entry, while in transition matrix method, each configuration gives several numbers of magnitude about $`N`$ to the transition matrix elements, thus variance reduction is expected. One of the most appealing feature of transition matrix Monte Carlo is an easy and straightforward way to combine several simulations. Additionally, we can use any valid sampling algorithm in a generalized ensemble which realizes the microcanonical property that states with same energy have the same probability. The flat-histogram algorithm  is such an algorithm that realizes multicanonical ensemble in which the energy histogram distribution is a constant. In the following, we present the transition matrix Monte Carlo method, introduce the flat-histogram and other related algorithms. We discuss the performance of algorithms with examples from Monte Carlo simulation results of the Ising models. We summarize in the last section.

## 2 Transition Matrix Monte Carlo Method

First, we give the definition for the transition matrix. Let $`W(\sigma \sigma ^{})`$ be the transition probability of the states from $`\sigma `$ to $`\sigma ^{}`$ of a Markov chain. To be definite, we consider a single-spin-flip dynamics with a canonical distribution, but the formalism is general. The transition matrix in the space of energy from $`E`$ to $`E^{}`$ is

$$T(EE^{})=\frac{1}{n(E)}\underset{E(\sigma )=E}{}\underset{E(\sigma ^{})=E^{}}{}W(\sigma \sigma ^{}),$$
(2)
where the summations are over all initial states $`\sigma `$ with energy $`E`$ and all final states $`\sigma ^{}`$ with energy $`E^{}`$. Estimates of the transition matrix can be obtained during a Monte Carlo sampling, where the summation over $`E`$ divided by $`n(E)`$ is interpreted as a microcanonical average of the state-space transition probabilities of the Markov chain, i.e.,

$$T(EE^{})=\underset{E(\sigma ^{})=E^{}}{}W(\sigma \sigma ^{})_E,$$
(3)
The expression can be further simplified if we consider single-spin-flip dynamics, with a spin choosing at random. In this case, $`_{E(\sigma ^{})=E^{}}W(\sigma \sigma ^{})=\frac{1}{N}N(\sigma ,\mathrm{\Delta }E)a(EE^{})`$, where $`N(\sigma ,\mathrm{\Delta }E)`$ is the number of sites such that a spin-flip causes the energy increasing by $`\mathrm{\Delta }E=E^{}E`$ in the current state $`\sigma `$. It is also the number of possible moves that one can make to change the energy by $`\mathrm{\Delta }E`$. Note that $`_{\mathrm{\Delta }E}N(\sigma ,\mathrm{\Delta }E)=N`$, where $`N`$ is the number of sites. A common choice of the single-spin-flip rate $`a(EE^{})`$ is the Metropolis rate $`\mathrm{min}(1,\mathrm{exp}(\mathrm{\Delta }E/kT))`$. Since this factor is a function of $`E`$ and $`E^{}`$, the microcanonical average $`\mathrm{}_E`$ is performed over $`N(\sigma ,\mathrm{\Delta }E)`$ only. We have

$$T(EE^{})=\frac{1}{N}N(\sigma ,\mathrm{\Delta }E)_Ea(EE^{})=T_{\mathrm{}}(EE^{})a(EE^{}),$$
(4)
where we have defined a normalized $`N(\sigma ,\mathrm{\Delta }E)`$ as the infinite temperature transition matrix.

The eigenvector corresponding to the eigenvalue 1 of the transition matrix is the probability of finding states with energy $`E`$. It is also proportional to the histogram $`H(E)`$. To determine the density of states, a numerically better choice is from the detailed balance. This gives us the relationship between histogram and transition matrix:

$$H(E)T(EE^{})=H(E^{})T(E^{}E).$$
(5)
If we use the fact that $`H(E)n(E)\mathrm{exp}(E/kT)`$ and Eq. (4), we obtain the so-called broad-histogram equation

$$n(E)N(\sigma ,E^{}E)_E=n(E^{})N(\sigma ^{},EE^{})_E^{}.$$
(6)
This is one of the basic equation for determining the density of states, as well as for the flat-histogram algorithm below.

## 3 Flat-Histogram Algorithm

Any sampling algorithm that can realize microcanonical property, i.e., the distribution of the states is a function of energy only, can be used to collect statistics for $`N(\sigma ,\mathrm{\Delta }E)_E`$. Using a canonical ensemble simulation, we need dozen temperatures in order to cover all the relevant energies. However, comparing to multi-histogram methods, the combination of data at different temperatures is very easy, we simply add up the matrix elements and then properly normalize.

Multicanonical ensemble  is a particularly good choice for the collection of transition matrix elements, since it reaches all energy levels with equal probability. Multicanonical ensemble is defined to be $`H(E)=\mathrm{const}`$, or the probability of configuration $`P(\sigma )1/n\left(E(\sigma )\right)`$. It is purely an artificial ensemble designed for computational efficiency. To realize the multicanonical ensemble, we can perform a single-spin flip with a flip rate of $`\mathrm{min}(1,n(E)/n(E^{}))`$. However, since the density of states $`n(E)`$ is not known beforehand, we have proposed to use the count number $`N(\sigma ,\mathrm{\Delta }E)`$. From the broad-histogram equation, Eq. (6), the ratio of $`n()`$ is related to the ratio of $`N()`$, we have

$$a(\sigma \sigma ^{})=\mathrm{min}(1,\frac{N(\sigma ^{},EE^{})_E^{}}{N(\sigma ,E^{}E)_E}).$$
(7)
This is our flat-histogram flip rate. Although the microcanonical average $`N(\sigma ,\mathrm{\Delta }E)_E`$ is also not available before the simulation, it can be obtained approximately during a simulation. We use the instantaneous value and running average to replace the exact microcanonical average. Numerical tests have shown that this procedure converges to the correct ensemble for sufficiently long runs. For realizing truly a Markov chain, it is sufficient for a two-pass algorithm. The first pass is as before. In the second pass, we use a multicanonical sampling rate, using the density of states determined from the first pass.

A variation of the algorithm is equal-hit algorithm which combines the N-fold way method  with a flip rate that gives an extended ensemble that is uniform in probability of visiting each new energy. Reference  gives more extensive discussions, as well as comparison with Wang-Landau method .

## 4 Some Results

As there are more detailed balance equations among the transitions of different energies than the number of energy levels, we determine the density of states from the transition matrix by solving a least-squares problem, or more generally, a nonlinear optimization problem. The optimization can be done either in the density of states $`n(E)`$, or in the transition matrix elements $`T(EE^{})`$. There are a number of constraints that the transition matrix must satisfy. The trivial one is the normalization, $`_{\mathrm{\Delta }E}T(EE+\mathrm{\Delta }E)=1`$. There exists a rather interesting constraint, known as $`TTT`$ identity, as well:

$$T(EE^{})T(E^{}E^{\prime \prime })T(E^{\prime \prime }E)=T(EE^{\prime \prime })T(E^{\prime \prime }E^{})T(E^{}E).$$
(8)
These constraints complicate the optimization algorithms.

While any of those extended ensemble methods reduce their efficiency as the system size increases, the accuracy of a two-pass flat-histogram/multicanonical simulation is rather good for a given fixed amount of CPU times . The method can also give excellent result for large systems, such as a $`256\times 256`$ Ising lattice , using a parallelized version of the program. The method is also applied to a lattice protein model, the HP model, with good performance .

A possible measure of computational efficiency is through the tunneling times. The tunneling time is defined to be the number of Monte Carlo steps in units of a lattice sweep ($`N`$ basic moves), for system making a pass from the highest energy level to lowest level, or vice versa. For the two-dimensional Ising model, this tunneling time diverges with system sizes according to $`L^{2.8}`$ which is worse than standard random walk. On the other hand, for spin-glasses with complicated low-temperature free-energy landscape, the tunneling time is much larger. It is about $`L^{4.7}`$  in two dimensions and $`L^{7.9}`$  in three dimensions. Another measure for spin glasses is given by the average first-passage times. It is defined as the average number of sweeps needed to reach a ground state. It is found in ref.  that the first-passage time diverges exponentially rather than according to a power. In any case, the equal-hit algorithm performs comparable to ‘extremal optimization’  which is an optimization algorithm inspired from self-organized criticality.

## 5 Conclusion

By collection the transition matrix, more information is obtained about the system, giving more accurate results. The effect of using transition matrix is more dramatic for small systems. Although the transition matrix analysis of data can be used with any simulation algorithms, extended-ensemble-based algorithms, such as flat histogram algorithm, are excellent choices. The efficiency of the flat-histogram related algorithms has be studied.

## Acknowledgments

The author thanks J. Gubernatis for invitation to this conference “The Monte Carlo method in the physical sciences: celebrating the 50th anniversary of the Metropolis algorithm” at Los Alamos. He would also like to thank Robert H. Swendsen, Lik Wee Lee, Zhifang Zhan, and Yutaka Okabe for collaborations on topics discussed here.