text:
# Remote information concentration using a bound entangled state

## Abstract

Remote information concentration, the reverse process of quantum telecloning, is presented. In this scheme, quantum information originally from a single qubit, but now distributed into three spatially separated qubits, is remotely concentrated back to a single qubit via an initially shared entangled state without performing any global operations. This entangled state is an unlockable bound entangled state and we analyze its properties.

Quantum entanglement has generated a great deal of interest recently as the most important resource in quantum information processing. The protocols of super dense coding , quantum teleportation  and telecloning  cannot be performed without some form of entanglement between the parties involved: classical correlations alone can never achieve the quantum efficiency arising from entanglement. Given that entanglement is a resource, it is important to be able to quantify it in order to deduce how effectively we can process information. There have been a number of suggestions for quantifying entanglement , but the most fruitful method comes from a procedure known as entanglement distillation . In this procedure, two distant users, Alice and Bob, share a certain number of entangled pairs all in the same state $`\rho `$. They then are allowed to perform local operations and communicate classically with each other (LOCC). The question is how many maximally entangled pairs they can obtain in this way. The limit of distillation in the infinite number of initial copies of $`\rho `$ is known as the entanglement of distillation . A natural question to ask is: which states $`\rho `$ can be distilled to maximally entangled states? Separable states $`\rho =_i\rho _A^i\rho _B^i`$ are clearly non-distillable. Surprisingly, however, a recent important discovery by the Horodecki family showed that there are also some entangled states which cannot be distilled . These states have appropriately been called bound entangled. They are peculiar as entanglement has to be invested in creating them by LOCC, but this invested entanglement cannot then be recovered by LOCC. Bound entanglement has been studied extensively in the last two years , nevertheless no information processing protocol has been found where bound entangled states perform better than just classically correlated states. Therefore it has seemed that they are useless for quantum information processing and that we always need to use some form of “free” (unbound) entanglement to achieve greater-than-classical efficiency. However, as we show in this letter, this is not the case.

We present an important protocol where bound entanglement can be utilized effectively and performs better than any classically correlated states. This protocol is remote information concentration, the inverse of telecloning . Quantum telecloning, as its name suggests, combines teleportation and cloning in such a way that a sender teleports an unknown qubit state $`|\varphi =\alpha |0+\beta |1`$ to a number of spatially separated receivers simultaneously. These teleported qubits cannot, of course, be exact replicas of the original qubit due to the linear laws of quantum evolution (“no-cloning theorem”) . However it has been shown that fidelities as high as allowed by the non-exact cloning (known as optimal cloning ) can be achieved. The optimal cloning state for $`|\varphi `$ is represented by a three qubit state

$`|\psi _c`$ $`=`$ $`\alpha \sqrt{{\displaystyle \frac{2}{3}}}\left\{|0|00+{\displaystyle \frac{1}{2}}|1\left(|01+|10\right)\right\}`$ (1)
$`+`$ $`\beta \sqrt{{\displaystyle \frac{2}{3}}}\left\{|1|11+{\displaystyle \frac{1}{2}}|0\left(|01+|10\right)\right\}.`$ (2)
where the first qubit is an ancilla and the last two qubits are two optimal clones. Now the question we ask is: once a state has been telecloned to spatially separated parties, can it then be recreated using only LOCC? The answer is yes and surprisingly involves a recently constructed unlockable bound entangled state .

The four particle unlockable bound entangled state presented by Smolin  is

$`\rho _{ub}={\displaystyle \frac{1}{4}}{\displaystyle \underset{i=0}{\overset{3}{}}}|\mathrm{\Phi }^i\mathrm{\Phi }^i||\mathrm{\Phi }^i\mathrm{\Phi }^i|`$ (3)
where $`|\mathrm{\Phi }^i`$ represents the four Bell states, $`|\mathrm{\Phi }^0=\left(|00+|11\right)/\sqrt{2}`$, $`|\mathrm{\Phi }^1=\left(|00|11\right)/\sqrt{2}`$, $`|\mathrm{\Phi }^2=\left(|01+|10\right)/\sqrt{2}`$ and $`|\mathrm{\Phi }^3=\left(|01|10\right)/\sqrt{2}`$. This state is not distillable if we do not allow joint quantum operations (i.e. if all four parties only operate locally), and is therefore, a bound entangle state. However, if we allow a two qubit joint operation, i.e. Bell joint measurement on any two qubits, we can obtain a maximally entangled state for the other two qubits via LOCC. Thus this state is unlockable. The unlocking mechanism is based on a joint operation for two out of four qubits.

Before we explain remote information concentration, we briefly summarize the forward process, telecloning. We focus on the $`1`$ to $`2`$ telecloning and its reverse in this letter. Generalizations to more qubits are possible and will be investigated elsewhere. The telecloning scheme  allows direct distribution of optimal clones from an single original qubit state $`|\varphi `$ to spatially separated parties using LOCC. In the telecloning scheme, we use an initially shared entangled channel (telecloning state)

$`|\xi _{tc}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}\left\{\right|00|00+|11|11`$ (4)
$`+`$ $`{\displaystyle \frac{1}{2}}\left(\right|01+|01)\left(\right|01+|01)\}.`$ (5)
where the first qubit is an input port of the distributor, the second qubit is an output port for the ancilla, and the third and forth qubits are output ports for the optimal clones. The telecloning protocol  is similar to teleportation; the distributor performs a Bell joint measurement between the unknown state and the input port qubit, and then the receivers, who hold output port qubits, perform a single qubit operations depending on the distributor’s measurement result.

Now we present our remote information concentration scheme. From the distributed optimal cloning qubits shared by the spatially separated parties (Alice who holds the ancilla qubit and Bob and Charlie who each hold a clone qubit), the original single qubit state is recreated at the location of an receiver, David in our scheme: $`|\psi _c_{ABC}|\varphi _D`$. We employ the unlockable bound entangled state (Eq.(3)) as an entangled channel for this scheme. The four qubits of the unlockable bound entangled state are initially distributed to Alice, Bob and Charlie (input port qubits) and David (output port qubit). The three senders, Alice, Bob and Charlie, individually perform Bell joint measurements between their qubits of the optimal cloning state and their input port qubits. We stress that no global operation is allowed between qubits belonging to different parties. One of the four possible outcomes $`\{\mathrm{\Phi }^i\}`$ is obtained by the measurement of a party. All three senders classically communicate their measurement results with David. ($`2\times 3=6`$ bits of classical information is communicated in total.) Each Bell measurement result $`\{\mathrm{\Phi }^i\}`$ is associated with the corresponding Pauli operators $`\{\sigma _i\}`$, where $`\sigma _0\mathrm{𝟏}`$, $`\sigma _1\sigma _z`$, $`\sigma _2\sigma _x`$, and $`\sigma _3\sigma _z\sigma _x`$. David performs a Pauli operation $`\sigma _j`$, which is the product (up to a global phase factor) of the three Pauli operators associated with the three Bell measurements, on his output port qubit. The output port of David is now in the original state $`|\varphi `$. A schematic picture of this protocol is shown in Fig.1. Since we do not allow joint operations on spatially separated qubits, the information channel in our scheme is indeed bound entangled. It is surprising that a bound entangled state can actually be useful for “transmitting” quantum information. In the following, we analyze this feature from two points of view: remote quantum operation and entanglement structure.

Remote quantum operation is performance of (global) unitary operations on remote qubits: A unitary operation $`U`$ is implemented by an initially shared entangled channel, Bell measurement, classical communication and (simple) single qubit operations, instead of directly running a quantum circuit. This is a generalization of quantum teleportation . Telecloning  and quantum information distribution via entanglement  are examples of remote quantum operation for $`1N`$ quantum optimal cloning with $`d`$-level particles, which requires one input port and $`2N1`$ output ports. More general cases requiring more than one input port have been studied by Gottesman and Chuang  in the context of “quantum computation using teleportation” and “quantum software”. The initially shared entanglement in the remote quantum operation scheme functions as quantum software. According to their result, unitary operations which belong to the Clifford group  can be implemented remotely, if we restrict the single qubit operations to be the Pauli operations. The compounding qubits of the shared entanglement need not be in the same location. In this case, the share entanglement functions as a transmission channel as well as quantum software. We consider this most restricted case of all-separated qubits.

To implement a unitary operation $`U`$ on a state of three input qubits ($`|\psi `$), the entangled channel consists of three input port qubits and three output port qubits. For unitary operations that can be decomposed into CNOT (controlled NOT) and Hadamard gates, which are members of the Clifford group , the entangled channel state is given by

$`|\xi ={\displaystyle \underset{\stackrel{~}{k}=\stackrel{~}{0}}{\overset{\stackrel{~}{7}}{}}}|\stackrel{~}{k}U|\stackrel{~}{k},`$ (6)
where $`\stackrel{~}{k}`$ is a 3-bit binary number, for example, $`\stackrel{~}{0}=000`$, $`\stackrel{~}{1}=001`$, …, $`\stackrel{~}{7}=111`$. The first three qubits are the input ports and the last three qubits are the output ports. All the qubits of this channel are spatially separated from each other. We assume that Alice, Bob, Charlie, David, Elizabeth, Fred each hold one qubit of the channel (in this order). Alice, Bob and Charlie individually perform Bell joint measurements on their input qubits (in the state $`|\psi _{ABC}`$, the qubits to be processed) and the input port qubits. David, Elizabeth and Fred perform an appropriate Pauli operation depending on the measurement results of Alice, Bob and Charlie. The mapping between measurement results and Pauli operations is initially agreed. The final state of David, Elizabeth and Fred is $`U|\psi _{DEF}`$.

Now we return to reverse optimal cloning. We define a reverse cloning unitary operator $`U_r`$

$`U_r|\psi _c=|\varphi \sqrt{{\displaystyle \frac{2}{3}}}\left(|00+{\displaystyle \frac{|01+|10}{2}}\right),`$ (7)
where the last two qubits are ancillas that are disentangled from the first qubit which holds the concentrated single qubit information. Note that $`U_r`$ does not initialize the ancilla qubits after the operation into the conventional ancilla state $`|00`$. $`U_r`$ can be decomposed into just CNOT gates as shown in Fig.2. Thus the reverse cloning operation is in the Clifford group and can be performed by remote quantum operation. $`U_r`$ is explicitly given in the computational basis by

$`U_r=\left(\begin{array}{cccccccc}1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0\end{array}\right).`$ (8)
Inserting $`U_r`$ in Eq.(6), we obtain the channel for the remote reverse cloning:

$`|\xi _{rc}`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{2}}}\{\left(\right|0000+|1111)|00+\left(\right|0101+|1010)|01`$ (9)
$`+`$ $`\left(\right|0011+|1100)|10+\left(\right|0110+|1001)|11\}.`$ (10)
In this expression, Alice, Bob, Charlie and David hold the first, second, third and fourth qubits, respectively. The last two qubits are ancillas and may be separated from the other qubits (their location is irrelevant). From Eq.(10), we obtain the unlockable bound entangled state $`\rho _{ub}`$, if we trace out the ancilla variables. Since the operations performed on ancillas do not effect the output port qubit, we can trace out the ancilla variables from the beginning. Then this remote quantum operation is equivalent to remote information concentration. This is why the unlockable bound entangled state actually functions as a channel for remote information concentration.

Next, we analyze the remote information concentration scheme from the view point of entanglement structure by considering the amount of entanglement in various separation cuts. First, we consider the pure state representation of the channel including the ancillas (Eq.(10)). For simplicity, we denote the six qubits as A, B, C, D, E, F. From Eq.(6), we see that the channel can be considered as a pair of maximally entangled 8-level particles, if we separate DEF from ABC. This means that there are 3-ebits of entanglement across the cut ABC:DEF. This is why the information of three input qubits can be processed and faithfully transmitted via teleportation from the three senders to the three receivers. We now consider the entanglement structure of the unlockable bound entangled state. In our scheme, at least 1-ebit of entanglement is required across the ABC:D cut for the faithful transmission of a single qubit quantum information. We calculate the relative entropy of entanglement  across the ABC:D cut. Relative entropy of entanglement for a mixed state is defined by $`E_{RE}=\mathrm{min}_{\sigma D}S(\rho ||\sigma )`$, where $`S(\rho ||\sigma )=\mathrm{Tr}\rho \mathrm{log}\rho \rho \mathrm{log}\sigma `$ is the quantum relative entropy, and the minimum is taken over $`D`$, the set of separable states. We can prove that in our case $`E_{RE}=1`$. The amount of entanglement across the ABC:D cut is indeed 1-ebit. However, as we have described before, there is no distillable pairwise entanglement in the unlockable bound entangled state, if no joint operations are allowed for qubits in different locations. Therefore the 1-ebit of entanglement across ABC:D does not explain the successful concentration of information. How can the information be processed and faithfully transmitted only by LOCC?

The answer lies in the optimal cloning state. The three qubits in the optimal cloning state are actually entangled each other. The Peres-Horodecki criterion , which is the smallest eigenvalue of partially transposed reduced density matrix of an ancilla qubit and a clone qubit $`\rho _{12}\mathrm{tr}_3|\psi _c\psi _c|`$, is $`(1\sqrt{17})/120.26`$. This negativity shows that each of two clone qubits is entangled (although not maximally entangled) with the ancilla qubits and this entanglement is not bound entanglement. We may conjecture how information processing and transmission have been achieved using only a bound entangled state and LOCC in our scheme as follows: the Bell joint measurement “combines” the optimal cloning state (the input state) and the bound entangled state (the channel state), which is initially “closed” for transmission. The unbound entanglement of the input state provides quantum correlation among the qubits of the bound entangled state. The quantum correlation “opens” the channel for transmitting concentrated single qubit information from distributed in three qubits of the input state. Entanglement of the input optimal cloning state and the unlockable bound entangled channel state function in a complementary fashion. This result explains the importance of the ancilla qubit in the optimal cloning state, since the ancilla qubit is necessary for holding entanglement.

Another interesting observation is that the unlockable bound entangled state is also valid for remotely concentrating information from the spatially separate 3-qubit error correction state: $`|\psi _e_{ABC}=\alpha |000+\beta |111|\varphi _D`$. The procedure is similar to the case of optimal cloning. The only difference is a modification to the mapping to Pauli operations. David performs an additional $`\sigma _2`$ if the measurement results from Bob or Charlie, but not both, belong to the set $`\{|\mathrm{\Phi }^0,|\mathrm{\Phi }^1\}`$. In this case, we may again consider that the (unbound) entanglement of the input state $`|\psi _e`$ “opens” the bound entangled channel for transmitting concentrated single qubit information from the input state. If we consider the quantum state $`|\varphi `$ as a quantum key , remote information concentration together with information distribution  may allow more secure distribution of the quantum key to David via spatially-separated, branched repeaters Alice, Bob and Charlie.

Finally, we show that no classically correlated state can achieve the same task (c.f. ). In optimal cloning scheme, due to the linearity of quantum transformations, mixed states can be cloned as well as pure. The same of course holds for telecloning. We consider the case when the qubit to be telecloned is maximally entangled with another qubit of George. After telecloning the qubit state into the qubits of Alice, Bob and Charlie, we perform the reverse process and remotely concentrate information at the location of David. Consequently, the qubits of David and George become entangled. If a shared state with only classical correlation could perform this remote information concentration, the procedure would create entanglement between David and George. This, however, is not possible: entanglement cannot be increased by LOCC. Therefore no classically correlated state can perform remote information concentration.

In this letter we have presented remote information concentration, the reverse process of quantum telecloning. It was shown that, surprisingly, the state needed for this operation is a bound entangled state. We have analyzed the remote information concentration scheme from two points of view, considering remote quantum operations and analyzing the entanglement structure of the bound state and the input state. We have shown that the unlockable bound entangled state is a reduced density matrix for the entanglement channel of remote reverse cloning, if we trace out the ancilla qubits of the output state. From our entanglement structure analysis, we have found that the functions of the entanglement of the optimal cloning state and the unlockable bound entangled state are complementary. We have also shown that the unlockable bound entangled state can be used for remotely concentrating information from a distributed 3-qubit error correction state, which may be useful for secure transmission of a quantum key. Furthermore, we showed that no purely classically correlated state can achieve this task. We hope that our work would stimulate more research into the nature of entanglement and its general usefulness in quantum information processing.

The authors are grateful to J. Watson for his help during the preparation of this manuscript. M.M. is supported by the Special Postdoctoral Researchers Program of RIKEN, and the Grant-in-Aid for Encouragement of Young Scientists (Grant No. 12740253) by Japan Society of the Promotion of Science. V.V. gratefully acknowledges funding by the European Union project EQUIP (Contract No. IST-1999-11053) and Hewlett Packard. Parts of this work were completed during the Banasque Center for Science program “Quantum Information Processing” (2000).


format: text

text:
# 1 Introduction

## 1 Introduction

Coarsening processes have been extensively studied both experimentally and theoretically over the past decades . Most of these studies deal with the way a system approaches its thermal equilibrium state. For example, when a liquid is quenched to temperatures below the liquid-gas transition point a coarsening process takes place as the system evolves towards the equilibrium phase separated state. In these systems the order parameter (the density) is conserved and the evolution proceeds by either nucleation growth or by spinodal decomposition . Similar processes take place in magnetic systems where a system evolves towards its magnetically ordered state when quenched below its critical point. In that case, however, the order parameter (the magnetization) need not be conserved by the dynamics and the details of the coarsening mechanism may differ from that of conserving dynamics.

Typically coarsening processes are characterized by a single growing length scale $`\mathrm{}(t)`$, for example the average domain size in the system. In many cases, at late times the system reaches a scaling regime where $`\mathrm{}(t)t^n`$. The value of the exponent $`n`$ usually depends on the symmetry of the system and its conservation laws. It has been found that for a scalar order parameter that is not conserved under the dynamics (as is the case in many magnetic systems) the growth exponent is $`n=1/2`$ . This has been demonstrated theoretically by numerous studies of Ising-type models with Glauber dynamics . On the other hand, when the order parameter is conserved (such as in liquid gas transitions and in phase separation in binary mixtures) the coarsening process is slower and the growth exponent was found to be $`n=1/3`$. This was first demonstrated by Lifshitz, Slyozov and Wagner  and has been confirmed by many studies of Ising models with Kawasaki dynamics .

More recently attention has focused on systems far from thermal equilibrium, driven by an external field . In many cases such systems reach a steady state in which, unlike the equilibrium case, detailed balance is not obeyed. The lack of detailed balance allows for many phenomena which do not occur in thermal equilibrium, such as phase separation and symmetry breaking in one dimensional systems . Most of the recent studies of these systems have been focused on the properties of the steady state itself rather than on the evolution towards it. There are however indications that coarsening processes in these systems may be rather different from those of systems evolving toward their equilibrium state. For example a study of the evolution of a driven Ising model with conserving dynamics has shown that in one dimension the average domain size grows as $`t^{1/2}`$  rather than the usual $`t^{1/3}`$ expected for non-driven systems.

The presence of the drive introduces a preferred direction in space making the systems inherently anisotropic. In many of these systems this results in striped structures. Typical examples are the stepped structures which occur in surface growth  and wind ripples formed in sand . Many models of driven systems have been introduced and studied in recent years. The striped structures that naturally emerge may be oriented either parallel or perpendicular to the driving field, depending on the details of the dynamics . For example in a driven Ising model introduced by Katz et al  stripes parallel to the driving field of alternating up and down spins are found at low temperatures. When the system is quenched to the low temperature phase at high values of the driving field narrow stripes are formed on a short time scale. On a longer time scale a coarsening process takes place in which the typical width of these stripes grows in time as the system evolves towards a fully phase separated state. This coarsening process is rather different from the one taking place in a system evolving towards its equilibrium state and is not well understood . On the other hand, stripes perpendicular to the drive direction are observed in related spin-1 type models where two types of oppositely charged particles move in the presence of neutral vacancies .

Recently, a different class of driven models, where phase separation takes place even in one dimension, has been introduced . In this class of models three or more types of particle are driven in a preferred direction under local dynamics. The important feature in these models is that the local dynamics results in a phase separated state even in one dimension ($`D=1`$). It has been demonstrated that in $`D=1`$ the coarsening process which accompanies the phase separation is slow, by which it is meant that the average domain size grows only logarithmically with time. An extension of one of the models belonging to this class to two and higher dimensions  showed that in dimensions greater than one stripes of alternating types of particles are formed perpendicular to the driving field. The width of the stripes along the direction of the driving field is found to grow as $`\mathrm{log}(t)`$, which is the same growth law as the one dimensional case. This is related to the fact that the interfaces which separate adjacent stripes are macroscopically smooth.

Given the common occurrence of striped structures in different classes of driven models, it is of interest to explore the possible coarsening phenomena within a wider range of models exhibiting these structures. For example, in the preceding paragraph we described a driven system with conserving dynamics exhibiting slow coarsening. This behavior can be attributed to the smoothness of the interfaces separating neighboring stripes. A natural question is as to how the coarsening of stripes is altered when one considers more general scenarios such as, say, non-conserving dynamics or dynamics which leads to rough interfaces separating the stripes.

To investigate this issue we study in this paper a simple model for the evolution of driven striped structures. In the model stripes emerge oriented perpendicular to the driving field. The microscopic dynamics is non-conserving and rough interfaces between the stripes arise. We consider a system of infinite extent in the direction parallel to the drive and of finite size $`L`$ in the other directions (see Fig. 1(a)). A detailed definition of the microscopic model and its dynamics is given in Section 2. Here we just describe the salient features:

* A microscopic state of the system is given by the configurations of the interfaces separating adjacent stripes (Fig. 1(a)). These interfaces are assumed to be single valued, that is no bubbles or overhangs are present.
* Each interface evolves under local driven dynamics. For example one can consider growth dynamics belonging to the KPZ universality class. All interfaces evolve under the same dynamical rules, and thus, in particular, they all move in the same direction and with the same average velocity.
* When two interfaces meet at a point, they locally merge to form a single interface. The evolution is such that after some time the entire two interfaces coalesce (Fig. 1(b)).

As a result of this dynamical process the number of interfaces keeps decreasing and the average width of a stripe, $`\mathrm{}(t)`$, increases with time. One is interested in studying the details of this coarsening process.

One can think of the model in analogy to a $`q\mathrm{}`$ state Potts system, where each site is assigned spin variable $`S=1,2\mathrm{}q`$. In the initial configuration the spins in each stripe are given the same state and the state associated with each stripe increases in an ordered sequence $`1,2,3\mathrm{}`$ (see Fig. 1(a)). Each state $`S=i`$ propagates into states $`S=j`$ where $`j>i`$. Clearly, the dynamics is such that the order parameter (the density of particles of type $`S`$) is not conserved since the number of particles of a given state changes in time. Related models have previously been considered in the context of cyclic food chain  and the evolution of the spatial mosaic of single-species domains was studied.

As a general motivation for the type of model we study, one may think of the spin state as representing the height of a $`D`$ dimensional terraced surface. Thus the striped structure corresponds to a sequence of terraces of increasing height; an interface between two stripes corresponds to a step. The surface evolves through particles being adsorbed or evaporated from the steps. In this picture two coalescing interfaces correspond to a terrace of a given height being eliminated from the system. However, the model to be studied in this paper should be not viewed as a microscopic description of this particular growth process. Rather it provides a very simple dynamics that leads to coarsening of stripes with no bubbles or overhangs in the interfaces between stripes. In order to construct a realistic microscopic model for a particular system (such as the terraced surface) one would have to add other features that would complicate analysis of the model.

The dynamical process described above is a generalization of the reaction diffusion process $`A+AA`$, in which diffusing $`A`$ particles (either with or without a drive) undergo a merging reaction as they collide. In the present generalization the coalescing objects are not particles but extended objects, i.e. the interfaces between stripes, which are manifolds in $`d=D1`$ dimensions. The $`A+AA`$ reaction diffusion process has been studied in detail over the years . It has been found that for both biased and unbiased diffusion the density (or equivalently the average distance between particles) decays as $`t^{D/2}`$ for $`D<2`$ and $`t^1`$ for $`D>2`$ . It is of interest to investigate how this behavior is changed for the generalization to coalescing manifolds. In the present context the manifolds themselves evolve in time, due to a roughening process, which could also affect the behavior.

The scaling properties of driven manifolds are different from those of stationary ones. It is of interest to consider the evolution of striped structures when the interfaces are non-driven. In this case the dynamics can be related to an energy function. As we shall see, the coarsening process is affected by whether or not the interfaces are driven, unlike the case of coalescing particles described in the previous paragraph.

In studying the model we find that the coarsening process may be well understood by considering the scaling properties of an isolated interface. We find that the evolution is different at early and late times. In both regimes the average domain size $`\mathrm{}(t)`$ grows algebraically in time. The different regimes, however, are characterized by different exponents. These exponents are determined by the universal properties of the isolated interfaces and can therefore be deduced using scaling arguments. Moreover the exponents depend on whether the interfaces are driven or not.

The paper is organized as follows: In Section 2 we define the model. In Section 3 a scaling argument which is based on the properties of an isolated interface is used to deduce the coarsening behavior. Monte-Carlo simulations which verify the scaling argument are described in Section 4. In Section 5 the special case where the interfaces are non-driven and energy function exists is studied using similar scaling arguments and Monte-Carlo simulations. We conclude in Section 6 with a discussion.

## 2 The coalescing interfaces model

Our model is defined on a $`D`$-dimensional hyper-cubic lattice with periodic boundary conditions. The lattice is infinite along one of its axes, which we label $`x`$, and is of size $`L`$ along the $`d`$ other axes labeled $`𝐲=(y_1,y_2,\mathrm{}y_d)`$. Each site $`𝐫=(x,𝐲)`$ can be occupied by one of $`q\mathrm{}`$ species of particles $`S_𝐫=1,2,\mathrm{}q`$. In the initial configuration particles are arranged in stripes of single species perpendicular to the $`x`$ direction. The stripes are in increasing order of particle type, their average extent along the $`x`$ direction is $`\mathrm{}_0`$ and the interfaces between them are flat. We are interested in constructing a model where species $`i`$ propagates into species $`j`$ when $`j>i`$. In view of the initial configuration this can be simply achieved through the dynamics we now describe.

The model evolves under random sequential dynamics according to the following rules: at each time step two neighboring sites along the $`x`$ direction $`𝐫=(x,𝐲),𝐫^{}=(x+1,𝐲)`$ are chosen randomly. If the particles at sites $`𝐫,𝐫^{}`$ are not of the same species, then we proceed in one of the following ways: with probability $`1/(1+p)`$, where $`p<1`$, the update $`S_𝐫^{}S_𝐫`$ is made, but only if the particles at $`𝐫`$ and its neighbors along the y directions are of the same species; with probability $`p/(1+p)`$ the update $`S_𝐫S_𝐫^{}`$ is made, but only if the particles at $`𝐫^{}`$ and its neighbors along the y directions are of the same species. For example, in $`D=2`$ the net effect of this algorithm is for the following moves to occur with the relative rates indicated above the arrows

$`\begin{array}{cc}& \\ S\hfill & \\ & \\ S\hfill & S^{}\hfill \\ & \\ S\hfill & \end{array}`$ $`\stackrel{1}{}`$ $`\begin{array}{cc}& \\ S\hfill & \\ & \\ S\hfill & S\hfill \\ & \\ S\hfill & \end{array}`$ (7)
$`\begin{array}{cc}& \\ & S^{}\hfill \\ & \\ S\hfill & S^{}\hfill \\ & \\ & S^{}\hfill \end{array}`$ $`\stackrel{p}{}`$ $`\begin{array}{cc}& \\ & S^{}\hfill \\ & \\ S^{}\hfill & S^{}\hfill \\ & \\ & S^{}\hfill \end{array}`$ (15)
where an empty square indicates that the move is attempted irrespective of the type of particle at that site.

Given the initial conditions, this dynamics ensures that the species index always increases to the right, along the $`x`$ direction. This is equivalent to stating that the interfaces are single valued and that the order among the interfaces is preserved. Furthermore, when $`p<1`$ the dynamics is such that species $`i`$ is driven preferentially to the right into neighboring species $`j`$ where $`j>i`$. The case $`p=1`$ is special in that interfaces are not driven. In the following we consider $`p<1`$ and we discuss the special case $`p=1`$ towards the end of the paper. Finally note that in this dynamics the densities of the various species are clearly not conserved.

Under this dynamics an isolated interface, evolves as in a restricted solid on solid (RSOS) model, by which it is meant that changes in the $`x`$ position of the interface as one moves in $`𝐲`$ directions are at most of magnitude 1. Moving interfaces of this kind belong to the Kardar-Parisi-Zhang (KPZ) universality class.

The coarsening process for a system of many interfaces can be understood in terms of the scaling properties of an isolated interface. These properties will be discussed in detail in the following section. However, at this point it is useful to outline the general picture that emerges. Starting from a dense collection of flat interfaces, at early times interfaces meet and coalesce due to the evolution of their width, $`W`$, which measures the lateral extent of the interface. Thus $`\mathrm{}`$, the average distance between neighboring interfaces, scales with the width of an individual interface. So for early times, the dependence of $`\mathrm{}`$ on time is governed by the time evolution of the width In this regime the growth is independent of $`L`$ and is determined by an exponent $`\beta `$

$$\mathrm{}(t)Wt^\beta \text{at early times.}$$
(16)
Since we are considering systems of finite size $`L`$ in the directions transverse to the drive, the width saturates to some final value $`W_{\text{sat}}(L)`$. Therefore at late times, when the width of the interfaces has already saturated, interfaces meet due to fluctuations, $`\mathrm{\Delta }h`$, in their center of mass position $`h(t)`$. That is,

$$\mathrm{}(t)\mathrm{\Delta }h\frac{t^\gamma }{L^\varphi }\text{at late times.}$$
(17)
where $`\gamma `$ and $`\varphi `$ are exponents to be determined. Notice the algebraic evolution of $`\mathrm{}`$ in the two regimes which contrasts with the logarithmic evolution of $`\mathrm{}`$ found in the model studied in .

The proposed picture is valid only at time scales much larger than that required for the process of coalescence of two interfaces. We find that at short times, the coalescence process takes in a finite time independent of $`L`$. This is because neighboring interfaces touch at a finite density of points. On the other hand in the long time regime the time scale for coalescence is found numerically to be proportional to $`L`$. This has to be borne in mind when testing the above scaling picture.

The exponents $`\beta `$, $`\gamma `$ and $`\varphi `$ are determined by the single interface behavior. In the next section we introduce the scaling analysis in a more quantitative way to obtain the exponents $`\beta `$, $`\gamma `$ and $`\varphi `$. Specifically, for $`p<1`$ the interfaces we consider belong to the KPZ universality class, so the two exponents are given in terms of known KPZ exponents. At the special point $`p=1`$ the interfaces belong to the Edwards-Wilkinson (EW) universality class, and the two exponents are changed accordingly.

## 3 Scaling analysis

As argued in the previous section the evolution of the system is governed by the width of an isolated interface at short times. On the other hand, at long times the width saturates due to the finite lateral extent of the system whereas the fluctuation of the center of mass (average $`x`$ position) keeps increasing. Thus at long times the center of mass fluctuations govern the coarsening behavior. We now define these quantities more precisely and quantify their scaling behavior.

The width of an interface, $`W`$, is defined by

$$W^2=\frac{1}{L^d}\underset{𝐲}{}(x_𝐲h)^2,$$
(18)
where $`x_𝐲`$ is the location of the interface along the $`x`$ direction and $`h`$ is the center of mass (average location) of the interface $`h=_𝐲x_𝐲/L^d`$. The angular brackets denote an average over the dynamics, starting from the given initial condition. The fluctuation in the average location of the interface $`\mathrm{\Delta }h`$ is defined through

$$\mathrm{\Delta }h^2=\left(hh\right)^2,$$
(19)
where as before angular brackets denote an average over the dynamics.

In general for a fluctuating interface the scaling behavior of the width is given by (see for example )

$$W=t^\beta f\left(\frac{t}{L^z}\right),$$
(20)
where $`z`$ is the dynamic exponent. The scaling function $`f(x)`$ is constant for $`x1`$ while $`f(x)x^\beta `$ for $`x1`$. This implies,

$$W\{\begin{array}{cc}t^\beta ,\hfill & \mathrm{for}t/L^z1\hfill \\ & \\ L^\alpha ,\hfill & \mathrm{for}t/L^z1.\hfill \end{array}$$
(21)
where $`\alpha =\beta z`$.

The fluctuation in the average location of the interface, $`\mathrm{\Delta }h`$, is dimensionally equivalent to the width so we may write, in a similar fashion to (20),

$$\mathrm{\Delta }h=t^\beta g\left(\frac{t}{L^z}\right).$$
(22)
For long times we expect the center of mass of the interface to diffuse, namely we expect $`\mathrm{\Delta }ht^{1/2}`$, so that $`g(x)x^{1/2\beta }`$ for $`x1`$. Thus

$$\mathrm{when}tL^z,\mathrm{\Delta }h\frac{t^{1/2}}{L^\phi }$$
(23)
where $`\phi =z/2\alpha `$.

For short times one needs to consider the specific universality class to which the interfaces belong. We start by restricting ourselves to the case $`p<1`$ (the case of $`p=1`$ will be discussed in Section 5). In this case the dynamics is biased along the $`x`$ direction, and the interfaces clearly belong to the KPZ universality class.

For KPZ interfaces at short times $`tL^z`$ it has been argued  that there are two regimes depending on the ratio $`(d+4)/z`$. That is,

$$\mathrm{when}tL^z,\mathrm{\Delta }h\{\begin{array}{cc}\frac{t^\theta }{L^{d/2}},\hfill & \mathrm{for}d+4<4z\hfill \\ & \\ \frac{t}{L^{2(z1)}}\hfill & \mathrm{for}d+4>4z\hfill \end{array}$$
(24)
where $`\theta =(d+4)/2z1`$.

We now use the width of the interface, $`W`$, and the fluctuation in its average position, $`\mathrm{\Delta }h`$, to deduce the average distance between the interfaces as a function of time $`\mathrm{}(t)`$. As stated above we will show that $`W`$ controls the early time behavior while $`\mathrm{\Delta }h`$ controls the late time behavior of the system. We assume that initially the average distance between the interfaces $`\mathrm{}_0`$ is much smaller than the width of a fully developed isolated interface $`L^\alpha `$. It will be shown that two distinct coarsening behaviors are expected for $`tL^z`$ and $`tL^z`$.

Consider first the long time ($`tL^z`$) behavior of the system. From Eqs. 21 and 24 one can see that $`\mathrm{\Delta }hW`$ and therefore $`\mathrm{\Delta }h`$ is expected to control the time scale on which interfaces coalesce. This implies that the average distance between interfaces should behave as

$$\mathrm{}\frac{t^{1/2}}{L^\phi }tL^z.$$
(25)
Next, consider the short time behavior ($`tL^z`$). In this case it is straightforward to verify, using the hyper-scaling relation $`\alpha +z=2`$, that $`W\mathrm{\Delta }h`$ for both regimes in Eq. 24. Thus at short times the coarsening behavior is controlled by the width of the interface and is expected to behave as

$$\mathrm{}t^\beta tL^z.$$
(26)
We have thus shown that $`\mathrm{\Delta }h`$ controls the long time behavior and $`W`$ the short time behavior. The short time behavior is expected to be lost if the initial spacing between the interfaces is larger than $`W_{\text{sat}}L^\alpha `$. Also, in models for which the interfaces are smooth ($`\alpha 0`$) the short time behavior is not expected to be seen.

In the following we verify the predictions of the scaling analysis by numerical studies of the models in $`D=2`$ and $`D=3`$ dimensions.

## 4 Monte Carlo simulations

We first consider the model in $`D=2`$ dimensions. The model is simulated on a lattice of size $`L_x\times L_y`$ where $`L_y`$ gives the lateral size denoted by $`L`$ in the previous sections. Periodic boundary conditions are used both in the $`x`$ and $`y`$ directions. Particles of different types are initially ordered in stripes parallel to the $`y`$ axis. The positions of the interfaces between the stripes are chosen randomly, such that the mean distance between them is $`\mathrm{}_0`$. At each step two neighboring sites in the $`x`$ direction, which are occupied by different species are chosen randomly, and a move is made according to Eq. LABEL:ratesKPZ where we have used $`p=0`$ in our simulations. In this way we maximize the speed of the simulation. After each such move, time is advanced by $`\tau =1/ϵ_x`$, $`ϵ_x`$ being the total number of nearest neighbor pairs in the $`x`$ direction occupied by different species. The algorithm would be equivalent to a usual Monte Carlo simulation if $`\tau `$ were drawn from a Poisson distribution with mean $`1/ϵ_x`$. Here we make the approximation $`\tau =1/ϵ_x`$ which is valid as long as $`ϵ_x`$ is large.

In Fig. 2 typical snapshots of the system at early, intermediate and late times are presented for a system of size $`L_y=120`$ and $`L_x=2000`$. The initial spacing between the interfaces is $`\mathrm{}_0=4/3`$. In Fig. 2(a) the early time behavior is shown. Here interfaces meet each other at several points indicating that their width is of the same order as the spacing between them. Note that in order to see such configurations the initial spacing between the interfaces must be much smaller than the final width $`W_{\text{sat}}`$ of an isolated interface. At the short time regime the mean domain size $`\mathrm{}`$ is expected to scale according to Eq. 26 which for $`d=1`$ reads $`\mathrm{}t^{1/3}`$.

At late times (Fig. 2(c)) the distance between interfaces is much larger than their width, and interfaces meet due to fluctuations in their center of mass location. In this regime $`\mathrm{}`$ is expected to scale according to Eq.(25), which in the $`d=1`$ case is just $`\mathrm{}t^{1/2}/L^{1/4}`$.

In the intermediate regime $`\mathrm{\Delta }h`$ is comparable with $`W`$ and we expect a crossover from the early to late time behavior. Fig. 2(b) shows a typical configuration at this stage of the dynamics. We see from the figure that in this regime there is a significant probability of three interfaces meeting at a point. This results in some arrested local configurations where an interface is temporarily frozen in locally triangular forms and the apparent width of the interface becomes large. These configurations are then released when a fourth interface approaches from the left and sweeps through the triangular forms. This effect is a direct consequence of the RSOS condition of the model. We have explored variations of the model in which the RSOS condition is relaxed, and found a variety of model specific effects in this intermediate regime. Nevertheless the early and late time behaviors are not affected by the details of the model.

The overall picture we have outlined is quantified in Fig. 3. In this figure the average distance between stripes $`\mathrm{}(t)`$ is plotted as a function of time and is compared with the width $`W`$ and fluctuation in the average interface location $`\mathrm{\Delta }h`$ of a single isolated interface. We calculate $`\mathrm{}`$ through $`\mathrm{}=L_xL_y/ϵ_x`$ where $`ϵ_x`$ is defined as above, namely, the total length in the $`y`$ direction of the interfaces. Similar results were obtained for various system sizes.

In Fig. 3 one can see that at early time indeed the behavior of $`\mathrm{}`$ follows closely that of $`2W`$, while at late times it follows closely that of $`2\mathrm{\Delta }h`$. The factors of 2 are because the growth of $`\mathrm{}`$ is actually controlled by the sum of the fluctuations of two neighboring interfaces. In the intermediate regime where $`\mathrm{\Delta }h`$ and $`W`$ are comparable $`\mathrm{}`$ deviates from both curves. This reflects the triangular configurations observed in this regime as discussed before.

We now make more precise tests of the scaling predictions. At late times ($`tL^z`$) the scaling argument predicts a specific dependence of $`\mathrm{}`$ on $`t`$ and $`L`$, $`\mathrm{}t^{1/2}/L^{1/4}`$. To test this, in Fig. 4 we plot $`\mathrm{}L^{1/4}/t^{1/2}`$ against $`t`$ for various system sizes. One can see that at late time indeed the data collapses confirming the predicted behavior.

To get a more quantitative test of the predictions of the scaling argument at early times ($`tL^z`$) is more difficult. The reason is that one has to access a regime in which $`\mathrm{\Delta }hW`$ while, at the same time, $`W`$ displays its early time growth behavior. To obtain a small $`\mathrm{\Delta }h`$, and to increase the time window in which $`W`$ keeps growing, we require a large lateral size $`L`$. Also note that at short times a coalescence process of two interfaces occurs in a finite time which does not grow with $`L`$. This is because in this regime the neighboring interfaces touch at a finite density of points.

In Fig. 5 the average distance between the domains $`\mathrm{}`$ is plotted as a function of $`t`$ for a system with $`L=10^4`$ and $`\mathrm{}_0=10`$. After a transient time the domain size $`\mathrm{}t^\beta `$ over more than a decade, where $`\beta =0.33\pm 0.01`$ as expected.

Next we present results obtained from simulation of a three dimensional system ($`D=3`$). In Fig. 6 the average distance between interfaces $`\mathrm{}`$ is plotted along with $`W`$ and $`\mathrm{\Delta }h`$ for a system with $`L_y=L_z=20`$. Again one can see that at early times $`\mathrm{}`$ follows closely $`W`$ while at late times it follows $`\mathrm{\Delta }h`$. While these simulations confirm the general picture outlined above, quantitative estimates of the exponents would require a more elaborate study.

## 5 Unbiased dynamics

We turn now to study the case $`p=1`$. In this case one can associate with the dynamics a well-defined local energy function. Since at each update the number of $`x`$ neighbor pairs of different species is reduced, the energy is Potts-like in the $`x`$ direction with unlike neighbor pairs costing energy. Unlike pairs in the $`y`$ direction do not cost energy provided that the RSOS condition is satisfied. Any step of the dynamics carried out, serves to lower the energy. In this way the dynamics can be thought of as a zero temperature system approaching an equilibrium state.

Single interfaces in this special case belong to the EW universality class, rather than to the KPZ class. Therefore the scaling analysis of Section 3 has to be applied to EW type interfaces. This modifies the values of the exponents and also the short time behavior of $`\mathrm{\Delta }h`$.

For interfaces belonging to the EW universality class $`\mathrm{\Delta }ht^{1/2}`$ for all times. Therefore the scaling form (22) implies

$$\mathrm{\Delta }h\frac{t^{1/2}}{L^\varphi }$$
(27)
for all times, with $`\varphi =z/2\alpha `$ as before. We use the exact values of the EW exponents: $`z=2`$, $`\alpha =1d/2`$, $`\beta =1/2d/4`$ to find

$$\mathrm{\Delta }h\frac{t^{1/2}}{L^{d/2}}.$$
(28)
As in the previous case one can now determine the behavior of $`\mathrm{}(t)`$. For long times ($`tL^2`$) using (21), (27) and the values of the EW exponents one finds that $`\mathrm{\Delta }hW`$ and $`\mathrm{\Delta }h`$ is expected to control the mean distance between interfaces. Therefore

$$\mathrm{}\frac{t^{1/2}}{L^{d/2}}tL^2.$$
(29)
For short times ($`tL^2`$) we have $`Wt^\beta \mathrm{\Delta }h`$. Therefore the coarsening is governed by the width of the interfaces and

$$\mathrm{}t^{1/2d/4}tL^2.$$
(30)
To summarize, in this case, just as in the $`p1`$ case, $`\mathrm{\Delta }h`$ controls the long time behavior and $`W`$ the short time behavior. The difference between the two cases is due to the different values of the exponents $`\alpha `$, $`\beta `$ and $`z`$ and is thus a direct consequence of the universality class.

Monte-Carlo simulations in $`D=2`$ with no bias support these findings. The simulations show that indeed at early times $`\mathrm{}t^{1/4}`$ and at late times $`\mathrm{}t^{1/2}/L^{1/2}`$ in agreement with the scaling argument and the known EW exponents. The results for the evolution of $`\mathrm{}`$ are qualitatively similar to those of the KPZ case and we do not display them here.

## 6 Discussion

In this work we have considered the evolution of a class of striped structures. A coarsening process takes place in which interfaces separating the stripes meet and coalesce. The dynamics may thus be viewed as a generalization of the particle reaction diffusion process $`A+AA`$ to extended objects, the interfaces. In this generalization the coalescing objects have internal structure that in turn affects the evolution. Thus one may have different scaling laws depending on the properties of the coalescing objects.

We studied a microscopic model defined through the dynamics of the interfaces that could realize both driven and undriven interfaces. Both cases lead to rough interfaces along with the average stripe width growing as a power law in time. However, the two cases lead to different scaling exponents. In the driven case an isolated interface exhibits KPZ behavior, while in the undriven case EW behavior is obtained. It is this scaling behavior of an isolated interface that determines both short and long time coarsening regimes of the stripes. We have shown that at early time the coarsening dynamics is determined by the width of the interfaces, while at short time it is determined by the fluctuations in the locations of the interfaces. The behavior is different from that of other driven models where the interfaces are smooth and the coarsening is logarithmically slow in time .

The analysis performed in this paper should apply also to other classes of interfaces. For example, the scaling behavior of moving fronts or interfaces has been studied for the case of ‘pulled’ fronts, moving into an unstable phase. Recently, it has been suggested that such an interface is not of the KPZ type . This indicates that there might be other possible types of scaling behavior than those we have studied in this paper. It would be of interest to explore such possible dynamics and study the resulting coarsening.

Acknowledgments: MRE and DM thank the Einstein Center and the University of Edinburgh for support and hospitality during mutual visits; the support of the Israeli Science Foundation is gratefully acknowledged.


text:
# Abstract

## Abstract

The direct $`\gamma `$-spectrometry method is used to measure more than 350 residual product nuclide yields from 0.1, 0.2, 0.8, and 2.6 GeV proton irradiated natHgO targets. The $`\gamma `$-spectrometer resolution is of 1.8 keV at the 1332 keV $`\gamma `$-line. The resultant $`\gamma `$-spectra are processed by the GENIE2000 code. The $`\gamma `$-lines are identified, and the cross-sections calculated, by the ITEP-designed SIGMA code using the PCNUDAT radioactive database. The 27Al(p,x)22Na reaction is used as monitor. The experimental results are compared with calculations by the LAHET, CEM95, CEM2k, INUCL, CASCADE, and YIELDX codes.

## Introduction

Mercury is planned to be used as a target material in all the present-day designs of the Spallation Neutron Source (SNS) facilities , thus necessitating that the proton-Hg interaction characteristics should be studied in a broad energy range from a few MeV to 2-3 GeV. Among the characteristics, the yields of residual product nuclei are of particular importance. They are independent nuclear constants to be used in practical calculations as well as to verify the codes for calculating the SNS facility design parameters, such as radioactivity (both current and residual), deterioration of resistance to corrosion, yields of gaseous products, poisoning, etc.

## Experiment

The experimental samples are 10.5-mm diameter discs manufactured by pressing fine-dispersed natHgO powder. The weight contents of impurities in the samples do not exceed 0.16%, of which 0.01% Si, 0.03% Cl, 0.02% Ca, 0.04% Ti, 0.03% Fe, and 0.01% Ba. The total content of the rest 60 elements, found by the spark mass-spectrometry, is below 0.02%.

The measurements were made by the relative method, using the 27Al(p,x)22Na reaction to monitor the process. The monitors are 10.5-mm Al foils with chemical impurities below 0.001%. Two independent proton beams from the ITEP U-2 synchrotron are used to irradiate the samples, namely, the low-energy (70-200 MeV) and high-energy (800-2600 MeV) beams.

Table 1 presents the characteristics of the experimental samples and monitors together with the main irradiation parameters.

The techniques for irradiating the samples and for processing the $`\gamma `$-spectra are presented in  together with formulas used to determine the fragment nuclide yields.

## Results

More than 350 yields of residual nuclei (from 22Na 111The 22Na and 24Na yields have been determined disregarding the contributions from the Al monitor samples. to 203Hg) from 0.1, 0.2, 0.8, and 2.6 GeV proton-irradiated natHg have been measured. Figs. 1-9 show the products that were measured at all the four energies.

The experimental data were compared with the LAHET, CEM95, CEM2k, INUCL, CASCADE, and YIELDX code-simulated yields. The comparison method and a short description of all codes together with further references may be found in . It should be noted that all these codes do not calculate the independent and cumulative yields individually for the ground and metastable states of the produced radionuclides, whereas the yields of either ground or metastable states alone are often measured. Therefore, those particular measurement results were excluded from the comparison procedure. The only exclusion is the case of measuring both states, so the total yields can be compared with the simulation results. Table 2 and Figs. 1-4 present the results of a detailed nuclide-by-nuclide comparison.

In Table 2, NT is the total number of the measured yields; NG is the number of the measured yields selected to be used in comparison with calculations; NS is the number of the products whose yields were simulated by a particular code; NC1,3 is the number of the comparison events when the simulation-experiment difference does not exceed 30%; NC2.0 is the number of the comparison events when the simulation-experiment difference does not exceed a factor of 2.

Figs. 1-4 show the results of the nuclide-by-nuclide comparison of the experimental data with the LAHET, CEM95, CEM2k, CASCADE, INUCL, and the YIELDX code-simulated results. One can see that all codes (except YIELDX) adequately predict the A$`>`$170 product yields for 100, 200, and 800 MeV protons and the A$`>`$120 product yields for 2600 MeV protons. The yields simulated by all the codes in the remaining ranges of masses, i.e., in the fission and fragmentation regions, are very different from experiment, with the most significant differences observed in the 80$`<`$A$`<`$103 range for 100 and 200 MeV protons, in the 48$`<`$A$`<`$130 range for 800 MeV protons, and in the 28$`<`$A$`<`$100 range for 2600 MeV protons. It should be noted that CEM95 and CEM2k do not calculate the process of fission itself, and do not provide fission fragments and a further possible evaporation of particles from them. When, during a Monte Carlo simulation of a compound stage of a reaction using the evaporation and fission widths these codes have to simulate a fission, they simply remember this event (that permits them to calculate fission cross section and fissility) and finish the calculation of this event without a real subsequent calculation of fission fragments. Therefore, the results from CEM95 and CEM2k shown here reflect the contribution to the total yields of the nuclides only from deep spallation processes of successive emission of particles from the target, but do not contain fission products. To be able to describe nuclide production in the fission region, these codes have to be extended by incorporating a model of high energy fission (e.g., in the transport code MCNPX, where these code are used, they are complemented by the RAL fission model).

Figs. 5-8 show the simulated mass yield of the reaction products. The experimental cumulative yields, which are often equal to the mass yields within measurement errors, are also shown for comparison. The following conclusions may be drawn from the comparison between the experimental and simulated mass yields:

* in the case of 0.1 GeV protons and A$`>`$190, all codes predict the mass curve shape quite adequately,
* in the case of 0.2 GeV protons and A$`>`$180, all code-simulated yields are in a good agreement with the data,
* in the case of 0.8 GeV protons, the best agreement with experiment is obtained by YIELDX, for A$`>`$130 and by INUCL, for A$`<`$130,
* in the case of 2.6 GeV protons and A$`>`$100, the CEM95 and CASCADE results agree with the data, while the LAHET calculations are underestimated and the YIELDX and INUCL yields represent the mass curve shape erroneously. None of the codes can describe well the experimental curve shape at A$`<`$100.

Fig. 9 illustrates the dependence of a part of measured yields on the proton energy, i.e., excitation functions. These data will be analyzed further after a final release of all the experimental results.

## Conclusion

Our experiment-to-simulation comparison has shown that, on the whole, the simulation codes can but poorly predict the experimental results. In some cases, the differences reach an order of magnitude or even more. The last version of the improved cascade-exciton model code, CEM2k , shows the best agreement with the data in the spallation region at all energies except 2.6 GeV, where the model overestimates the expected experimental fission cross section of about 75 mb  by a factor of 4. This overestimation of the fission cross section causes an underestimation of the yield of nuclei which are most likely to fission at the evaporation stage of a reaction, after the cascade and preequilibrium stages , i.e., for 170 $`<`$ A $`<`$ 185. (Similar disagreement with the data one can see as well for LAHET and CEM95, that is also related with an overestimation of the fission cross section at 2.6 GeV). The code CEM2k is still under development, its problem with the overestimation of fission cross sections at energies above 1 GeV has yet to be solved, and it has to be complemented with a model of fission fragment production as mentioned above, to be able to describe as well fission products.

This means that the nuclear data must be accumulated persistently in the above-mentioned ranges of energies and masses, to help improve hadron-nucleus interaction models released in codes used in applications.

## Acknowledgement

The work was carried out under the ISTC Project #839 supported by the European Community, Japan (JAERI), Norway and, partially, by the U. S. Department of Energy.


text:
# Theory of extraordinary optical transmission through subwavelength hole arrays

\[

## Abstract

We present a fully three-dimensional theoretical study of the extraordinary transmission of light through subwavelength hole arrays in optically thick metal films. Good agreement is obtained with experimental data. An analytical minimal model is also developed, which conclusively shows that the enhancement of transmission is due to tunneling through surface plasmons formed on each metal-dielectric interfaces. Different regimes of tunneling (resonant through a ”surface plasmon molecule”, or sequential through two isolated surface plasmons) are found depending on the geometrical parameters defining the system.

\]

In the last few years, and mainly due to advances in nanotechnology, there has been a renewed interest in exploiting the dielectric response of metals to make photonic materials. For instance, the photonic insulating properties of metals can be used to trap incident radiation, focusing light in very small volumes. Very recently, another interesting effect of light interacting with structured metals has been discovered: the transmission of light through subwavelength hole arrays made in a metal film can be orders of magnitude larger than expected from standard aperture theory. Apart from its fundamental interest, this extraordinary transmission effect has potential applications in subwavelength photolithography, near-field microscopy, wavelength-tunable filters, optical modulators, and flat-panel displays, amongst others. While the wavelength at which some transmission features appeared suggested that surface plasmons (SP) were involved in the process, the physical mechanism for the huge enhancement has not yet been elucidated. Some calculations have been performed for a simpler geometry: an array of slits, where high transmission was also predicted. However, although interesting in their own right, these results do not apply to the experimental situation.

In this letter, we present the first fully three-dimensional theoretical study of transmission through hole arrays, obtaining an excellent agreement with experimental data. Moreover, we develop a simplified version of the model that clearly captures the physics involved.

Fig. 1 shows the experimental “zero-order” transmittance of light at normal incidence ($`T_{00}`$), through an array of holes in a free-standing metal film. The free-standing metal film, of which the fabrication is described elsewhere  consisted of a $`220nm`$ thick Ni core, perforated with a square array of holes by focused-ion beam milling. The film was subsequently overcoated with $`50nm`$ of Ag on both sides by sputter deposition which resulted in a coating of the walls of the holes as well as the in-plane surfaces of the film. The total thickness of the film was $`h=320nm`$ and the lattice constant of the hole array was $`L=750nm`$. After coating the holes had an average diameter of $`280nm`$. It has been shown  that such a ”sandwich” structure has the same transmission properties as an equivalent perforated film made of silver throughout. The advantage of a free-standing metal film is that it is possible to mill much better defined holes than in a film on a substrate. Moreover, in a free-standing film the dielectric constant is the same in all non-metallic regions. In Fig. 1, $`T_{00}`$ shows the well known Rayleigh minima of Wood’s anomaly  which appear in any diffractive array roughly when an order of diffraction emerges tangent to the array. The data also shows the extraordinary transmittance effect: at $`\lambda 800nm`$, $`T_{00}`$ of the order of 15% was found. As the area covered by holes is only 11%, the normalised-to-area transmittance of light is 130%. Standard aperture theory for a single hole predicts transmission efficiencies of order 1%.

In our theoretical model, we consider a metal film of thickness $`h`$, perforated with square holes of side length $`d`$, periodically distributed, as in the experiment, in a square array of lattice parameter $`L`$. In the long-wavelength limit, experimental results have shown that the transmission coefficient depends on hole area, but does not appreciably depend on hole shape. The choice of square holes is, therefore, inessential, and was motivated purely for analytical convenience. When comparing with experimental results we take the square side as the square root of experimental hole area. We consider the general case in which light impinges from a uniform medium 1, into the perforated metal film (region 2) which lies on a substrate (medium 3), where the light in collected. The dielectric constants are, $`ϵ_1`$ in medium 1, $`ϵ_2`$ inside the holes, and $`ϵ_3`$ in medium 3. The metal is characterised by a frequency-dependent complex dielectric constant $`ϵ_m(\omega )`$.

We treat the EM fields in the metal through the surface impedance boundary condition, a procedure fully justified, as the frequencies considered are well below the metal plasma frequency. We expand the EM fields in the eigenmodes of the different regions. That is, the EM fields are a linear expansion of S and P plane waves in regions 1 and 3, and Bloch waves combining (evanescent and, if there is any, propagating) $`TE`$ and $`TM`$ waveguide modes inside the holes. Although, in principle, an infinite number of modes in each region should be taken into account, results quickly converge with the number of parallel wavevector components considered. The probability amplitude, $`t_{if}`$, for an arbitrary incident (from medium 1) EM plane wave, $`i`$, to be transmitted to a outgoing plane wave $`f`$, in medium 3, was calculated in a multiple scattering formalism. By appropriately matching the EM fields in the 1-2 and 2-3 interfaces, $`t_{if}`$ can be expressed in terms of transmission and reflection amplitudes for a single interface as:

$$t_{if}=\underset{\alpha ,\beta ,\gamma }{}\tau _{i\alpha }^{12}e_\alpha \left(\delta _{\alpha \beta }\rho _{\alpha \gamma }^Re_\gamma \rho _{\gamma \beta }^Le_\beta \right)^1\tau _{\beta f}^{23}$$
(1)
Latin indexes refer to plane waves, either in medium 1 or 3, and Greek indexes to waveguide modes inside the holes. $`\tau ^{12},\tau ^{23},\rho ^R`$ and $`\rho ^L`$ are components of the scattering matrix for a single interface, either 1-2 or 2-3. More precisely, $`\tau ^{12}`$ is a matrix giving the transmission amplitude between modes in medium 1 to modes in medium 2. $`\rho _{\alpha \gamma }^R`$ is the reflection amplitude for the waveguide mode $`\alpha `$, traveling towards medium 1, to be reflected into the mode $`\gamma `$ traveling away from medium 1, after scattering with the 2-1 interface. $`\rho _{\alpha \gamma }^L`$ denotes the same but for the interface 2-3. $`\tau _{\beta f}^{23}`$ connects mode $`\beta `$ with outgoing mode $`f`$, through interface 2-3. Finally, $`e_\alpha =e^{ıq_{z\alpha }h}`$, $`z`$ is the direction perpendicular to the interfaces, and $`q_{z\alpha }`$ is the wavevector component along the z-direction for the waveguide mode $`\alpha `$.

Figure 2 shows the computed results for $`T_{00}`$. Calculations were done for the nominal experimental values for the different parameters defining the structure. The metal dielectric constant has been taken as the measured one for Ag. It must be stressed that in our calculation there are no fitting parameters. Also shown are results for slightly different values for hole size, $`d`$, in order to give an idea of the sensitivity of the results to non-uniformities of hole area. As shown in Fig. 2, our calculation also shows peaks, with $`T_{00}`$ which exceeds by 2-3 orders of magnitude the predictions by a theory based on independent holes. Moreover, the wavelenghts at which peaks occur closely match the experimental values. Clearly the mechanism responsible for the extraordinary transmission is present in our model. Discrepancies between calculated and experimental values for linewidths and maximum transmittances, can well be due to small variations in the hole diameter throughout the sample, as supported by the sensitivity of the transmittance to $`d`$ shown in Fig. 2.

Although the calculation also shows extraordinary transmission through periodic subwavelength holes in a metal, it is not obvious what mechanism is involved. However, as we show in the following, a highly simplified version of the model, which can be analytically worked out, also shows the same behaviour in the regime of interest: $`\lambda Ld`$. In this minimal model, a strong truncation in parallel wavevectors is made: in regions 1 and 3, only first order diffraction is considered, i.e., the possible wavevectors in the direction of the incoming electric field, $`x`$, are $`k_{x0}=0`$ and $`k_{x\pm 1}=\pm \frac{2\pi }{L}`$; in region 2, where all modes are evanescent for $`L>2d`$, only the most slowly decaying evanescent mode, the $`TE_{01}`$ mode, is taken into account. Figure 3 shows the comparison of the fully converged calculation (solid line) with the results obtained with the minimal model (dashed line), for the nominal parameters of the experimental set up, both for (3a) the experimental dielectric constant for Ag and (3b) for an hypothetical case of non-absorbing Ag ($`Im\left[ϵ(\omega )\right]=0`$). The minimal model captures well the extraordinary transmission phenomena: the main effect of truncating the number of modes is a small unimportant shift in the position of the transmission peaks. In order to further simplify the discussion, in the following we consider the ideal non-absorbing case. Finite absorption merely reduces the height of the peaks (for zero absorption the transmittance can be as large as 100 %) without altering the physical picture.

For the symmetric case, $`ϵ_1=ϵ_2=ϵ_3=1`$, $`\rho _{11}^R=\rho _{11}^L\rho `$ and $`t_{00}`$ simplifies to

$$t_{00}=\frac{\tau _{01}^{12}e^{|q_{z1}|h}\tau _{10}^{23}}{1\rho ^2e^{2|q_{z1}|h}}$$
(2)
where $`\tau _{01}^{12}=2S_0/(G_2+G_1),\tau _{10}^{23}=2G_2/(G_2+G_1)`$ and $`\rho =(G_2G_1)(G_2+G_1)`$. $`G_2`$ and $`G_1`$ are “conductances” (inverse of total impedances) in region 2 and in region 1, respectively. They have somewhat involved analytical expressions that, for the case of a hole array in a perfect metal, simplify to $`G_2=q_{z1}/g,G_1=S_0^2+2S_1^2g/k_{z1}`$, with $`g=\omega /c,S_0=d/L,S_1=S_0\mathrm{sin}\left[k_{x1}d/2\right]/(k_{x1}d/2)`$, and $`k_{z1}=\sqrt{g^2k_{x1}^2}`$.

Wood’s anomaly is present in this model, associated to a zero in either $`\tau _{01}^{12}`$ or $`\tau _{10}^{23}`$, i.e., when the conductance in medium 1 or 3, respectively, is infinite.

With respect to the extraordinary transmission, transmittance peaks $``$ 1 seem to point to the presence of resonant phenomena. This idea is reinforced by the presence of multiple scattering denominators in Eq.(2) that could be close to zero. Naively, one would expect that this possibility never occurs: $`e^{2|q_{z1}|h}1`$ in the regime considered and $`\rho `$ is the amplitude for reflection, which can be expected to be always less than 1. Actually, this is not true, as the condition $`|\rho |1`$ applies only for propagating modes, where it corresponds to current conservation. In the regime we are considering, all modes inside the hole are evanescent and current conservation only restricts $`Im[(\rho ]0`$, with no restrictions on the real part or the modulus of $`\rho `$ . This opens up the possibility of resonant denominators for $`|\rho |1`$ . In fact, it can be shown analytically that $`|\rho |`$ has a peak at $`\lambda `$ slightly larger than $`L`$, with a maximum value that scales as $`(\frac{L}{d})^3`$. Figure 4a shows $`|\rho |`$ as a function of $`\lambda `$ for the parameters considered in Fig. 2. It is illuminating at this point to consider the dependence of the transmittance on $`h`$. Figure 4b illustrates graphically that the peaks in $`T_{00}`$ occur at the $`\lambda `$ for which the distance between $`|\rho |`$ and $`e^{|q_{z1}|h}`$ is minimal. For zero absorption, a $`T_{00}`$ close to unity occurs whenever this distance is zero, i.e. when the curves intersect. For large enough $`h`$, the curves do not intersect and $`T_{00}`$ maxima decay exponentially with $`h`$.

All this is mathematics, showing that the effect of extraordinary transmission has a resonant nature. But, what is the physical origin? Analytically, we find that the frequency at which the maximum of $`|\rho |`$ appears coincides with the SP frequency of a periodically (air-metal with holes) isolated interface at parallel momentum $`\frac{2\pi }{L}`$ (what we are going to call first-order SP). This is a localised EM mode, confined to the interface. It is not the only EM mode at that frequency, as there also exists the zero-parallel-momentum propagating mode, into which the SP can decay. This decaying channel provides a linewidth and a maximum value to $`\rho (\lambda `$), which otherwise would be a delta function. Furthermore, when we calculate the SP of the periodically perforated metal slab, we find that the single interface first-order SPs on the two metal surfaces combine to form a ”SP molecule”, in much the same way as electronic states of isolated atoms combine to form molecular levels. The frequencies of the first-order SP molecule are given precisely by the same condition for resonant denominators in the expression for $`t_{00}`$, if the coupling to the radiative modes is neglected. As in the case of the isolated interface, this SP molecule coexists with a continuum and, rather than an exact eigenstate, is a (more or less long lived) resonance.

This point on the formation of the SP molecule can be more precisely discussed by considering typical times in the transmission process. The width of the peak of $`\rho (\omega )`$ is related to the lifetime of the SP, $`t_{rad}`$, the time taken for a SP of an isolated surface to decay into the radiative mode. We obtain that, approximately, $`t_{rad}(c|k_{z1}|)^1`$. From the frequency splitting between molecule levels, taken as uncoupled from the radiation modes, the time needed for each of the molecular modes to form, $`t_{res}`$, can be estimated as $`t_{res}t_{rad}(d^3e^{|q_{z1}|h}/L^3)`$. So, in the small $`h`$ regime, where the intersection of $`|\rho |`$ with the exponential occurs at the bottom of the $`|\rho |`$ peak, $`t_{res}t_{rad}`$ and the molecule levels are fully formed. The photon then goes back and forth several times inside the hole, building up coherent constructive interference in the forward direction much as would occur in electron resonant tunneling or in a Fabry-Perot interferometer (but with evanescent waves instead of propagating ones). As we consider larger $`h`$, $`t_{res}`$ increases, while $`t_{rad}`$, being a property of an isolated interface, essentially remains constant. So, the photon can make fewer round trips inside the hole before being radiated to infinity, and the concept of plasmon molecule becomes less well-defined. The condition $`t_{res}=t_{rad}`$ marks the value of $`h`$ at which the transmittance maximum ceases to be one. For even larger $`h`$, $`t_{res}t_{rad}`$ and the process is more like sequential tunneling, where the incoming photon gets trapped in a SP, tunnels to the SP at the other interface and then couples to the outgoing radiative mode and exits. The transmittance is enhanced in this case, as the photon can use two intermediate states (the first-order SP at both interfaces) to cross the metal film, but the enhancement is not as efficient as when the hopping occurs back and forth several times building up a constructive interference.

Absorption introduces another time into the problem: $`t_{abs}`$, the typical time it takes for a photon to get absorbed. Provided this time is smaller than both $`t_{res}`$, and $`t_{rad}`$, the physical picture in not greatly altered. This is the case for the experimental parameters: absorption reduces the $`T_{00}`$ maxima by a factor around 2-3, without altering either the position or the width of the peaks.

The above discussion is for the special case of matching dielectrics ($`ϵ_1=ϵ_3`$) but the picture still applies in the non-symmetric case ($`ϵ_1ϵ_3`$), the situation in which extraordinary transmittance was initially discovered. Here, SP on the two different interfaces of the metal film have different frequencies and, instead of a SP diatomic homopolar molecule forming in the presence of a continuum, we have instead a diatomic heteropolar molecule. From our model we find that there still are resonant peaks, but the transmittance maxima is reduced with respect to the symmetrical case. A detailed comparison between experimental and theoretical results in this case will be presented elsewhere.

In summary, we have presented a detailed theoretical study of the extraordinary transmission of light through subwavelength hole arrays in metal films, obtaining good agreement with experimental data on free-standing films. An analytical minimal model, based on the numerical calculations, conclusively shows that the the enhanced transmission is due to tunneling through surface plasmons. Different regimes in the enhanced transmission are found: for small film thickness the tunneling is resonant through ”plasmon molecule” levels while, for large thicknesses, the photon hops from surface plasmon to surface plasmon, but exits the structure before the molecule level is developed, the tunneling being then sequential. Undoubtedly, this found ability of surface plasmons to transmit and focus light very efficiently is not restricted to this particular geometry and will be exploited for controlling light in new photonic devices.