text:
# Structured Mixed Phase is Favored in Neutron Stars

## Abstract

We give a general thermodynamical argument showing that in neutron stars, the Coulomb structured mixed phase is always favored for any first order phase transition involving systems in equilibrium with baryon number and electric charge as the two independent components. This finding is likely to have important consequences for many neutron star properties, e.g., glitch phenomena, transport and superfluid properties, r-mode instabilities, and the braking index.

A number of phase transitions may take place in neutron stars. Described in most detail is the first order nuclear liquid-gas transition in the inner crust of a neutron star, where nuclei with different shapes are in equilibrium with a neutron gas . Above the nuclear saturation density, phase transitions involving, for instance, quark matter, pion, and kaon condensation have been suggested and described, see e.g., Ref.  for an overview.

Before 1992, local charge neutrality was imposed on all the possible phase transitions above saturation density, leaving the systems with only one independent component, namely the baryon number. This treatment made any first order transition resemble, e.g., the water-ice transition. However, in 1992 it was pointed out by Glendenning  that because neutron star matter has two independent components, namely the baryon number and the electric charge, charge neutrality must be applied as a global and not a local condition. Gibbs general conditions for thermodynamic equilibrium for both chemical potentials (usually taken to be the neutron and electron chemical potentials) cannot be satisfied otherwise. As a consequence of this, the charge does not vanish identically; rather, there results a nonvanishing charge density of opposite sign for each bulk phase .

A first order transition with multiple independent components behaves qualitatively differently from a single component one as seen in Fig. 1. For a one-component system, the transition takes place at a constant pressure where the two phases, having different densities, for all proportions are in equilibrium. The equilibrium pressure is found through a Maxwell construction. Systems with two components have an additional degree of freedom, which means the proportions of the two phases in equilibrium vary as a function of the pressure. The phases still have different densities, but the densities are now also functions of the pressure . Mixtures of two miscible liquids in equilibrium with their vapor phase are well described examples from physical chemistry . In a neutron star this means that a two-component phase transition may have a considerable radial extent in contrast to a one-component transition, which takes place at a single radial point in the star, the point where the pressure of the two phases are equal. This region of two coexisting phases in equilibrium is usually referred to as the mixed phase.

In  surface and Coulomb effects were neglected in the description of the mixed phase, although it was pointed out that their inclusion would result in a Coulomb lattice, analogous to the lattice resulting from the nuclear liquid-gas transition in the inner crust of a neutron star. Such a structured mixed phase above the saturation density was first studied in  for the deconfinement transition. It was there concluded that the mixed phase may not be energetically favored, compared to the locally charge neutral, one-component system, if the surface and Coulomb energies are ‘sufficiently large’.

We show in the following that for any first order phase transition of a two-component system of nuclear matter in full thermodynamic equilibrium, the energy (including Coulomb and surface energies) of this system is always smaller than the energy of the same system where one component is frozen out . For a neutron star, the important implication is that for matter in equilibrium there will always be a structured mixed phase if any first order phase transition takes place in its dense interior.

Gibbs conditions for thermodynamical equilibrium between two phases (I and II) with $`n`$ independent components separated by an arbitrary boundary are

$`T^I`$ $`=`$ $`T^{II}`$ (1)
$`\mu _i^I`$ $`=`$ $`\mu _i^{II}i=1,2,\mathrm{},n`$ (2)
$`P^I(\mu _i^I,T^I)`$ $`=`$ $`P^{II}(\mu _i^{II},T^{II})+\sigma \left({\displaystyle \frac{1}{R_1}}+{\displaystyle \frac{1}{R_2}}\right)+\mathrm{},`$ (3)
where $`T`$ is the temperature, $`\mu _i`$ is the chemical potential of the $`i`$’th independent component, $`P`$ is the pressure, $`\sigma `$ is the surface tension, and $`R_1`$ and $`R_2`$ are the principal curvature radii (higher order corrections like curvature effects have not been explicitly included in the pressure equilibrium condition). These conditions describe thermal, chemical, and mechanical equilibrium, respectively. In addition the energy (as well as the free energies) is at a minimum for a system in equilibrium. All these conditions are direct consequences of the second law of thermodynamics. In the case of a neutron star, Gibbs conditions have to be solved consistently with the condition of charge neutrality, which as emphasized, cannot be imposed as a condition of local vanishing charge density else there are too many unknowns for the number of equations.

Cold nuclear matter in a single phase in equilibrium is of course electrically neutral. $`\beta `$-equilibrium ensures that the sum of the proton and electron chemical potentials equals the chemical potential of the neutron, $`\mu _n=\mu _p+\mu _e`$ (the neutrinos escape the system freely). This local charge neutrality requires the electron and proton densities to be equal in bulk. Earlier, local charge neutrality was likewise imposed on both phases of the system during a phase transition. This assumption has the effect of freezing out the electric charge component and leaving the baryon number as the single independent component. The Maxwell construction ensures pressure equilibrium and chemical equilibrium between the baryons across the phase boundary for such a phase transition. However, the electron chemical potential is generally discontinuous across the phase boundary.

In the current picture of the bulk phases, the local charge neutrality of both phases is relaxed to a global charge neutrality condition, which allows the system to freely explore the additional degree of freedom in a system with two independent components. Compliance with Gibbs phase equilibrium conditions for both chemical potentials naturally result in a nonvanishing charge density for each phase. The normal nuclear matter phase will be positively charged, and thus more isospin symmetric; whereas the high density phase (e.g., quark matter) will be negatively charged, and partly replace electrons (and muons) as global neutralizing agents. The bulk energy of the one-component system lies above the bulk energy of the mixed phase system if Coulomb and surface energies are neglected, since the charge neutrality constraint is relaxed in the latter approach .

In the current picture of the structured mixed phase as presented in Ref. , the mixed phase is energetically favored only if the sum of Coulomb and surface energies in the mixed phase is smaller than the gain in bulk energy between the two systems. A schematic plot of this comparison is shown in Fig. 2 (see also the figure in  or Fig. 30 in ). There are three possible scenarios, corresponding to the three curves of Coulomb and surface energies: $`a`$, the structured mixed phase is favored in the whole bulk mixed phase region, $`b`$, it is favored in some regions of the bulk mixed phase, $`c`$, it is not favored anywhere. The properties of the equations of state including surface and screening effects for the two phases determine which scenario is valid for the phase transition in question. Screening effects reduce the Coulomb energy and will play an increasingly important role with increasing size of the charged structures. Screening lengths for different components have been estimated, the ranges are from about 5 fm to 13 fm. Structures that are smaller than the screening lengths will have almost uniform particle densities within each phase and a nonvanishing local charge density. On the other hand, the larger the structures are compared to the screening lengths, the more the system will resemble two electrically neutral bulk phases .

The new picture we present is based on the thermodynamical equilibrium conditions described in Eq. (1) and below. These conditions are direct consequences of the second law of thermodynamics and therefore have general validity. The standard Maxwell construction can only ensure mechanical equilibrium and one chemical potential to be in equilibrium across a phase boundary. The other chemical potential will generally be discontinuous across the boundary. Thus, if one component is frozen out by imposing constraints on a system with two independent components, this system cannot be in full thermodynamical equilibrium. Therefore the (free) energy of such a system is not at a minimum, and there must be another system in thermodynamical equilibrium in which the energy is a global minimum. This statement is valid over the whole range of the mixed phase region and not only at the pressure $`P_t`$, where the actual discontinuity in one chemical potential is. To see this consider the system (I) in the neutron star with the global energy minimum at $`P_t`$ and assume that the single phase (II) with the imposed local charge neutrality is preferred at an infinitesimal smaller (or larger) pressure. There will be a phase boundary between these phases, but they will have different chemical potentials and therefore will only be in mechanical equilibrium. Thus there must be a system different from (II) at this pressure which is in thermodynamical equilibrium with (I) and being at the global energy minimum. This argument can be repeated until the whole mixed phase region has been spanned. Notice, that the above argument is valid irrespective of whether the central pressure of the neutron star is larger than $`P_t`$ of the Maxwell system or not. Therefore, a two-component system with one component frozen out cannot be at a global energy minimum since it is not in full thermodynamical equilibrium. The system with the global energy minimum, in the pressure region where the equations of state indicate a two-phase equilibrium is possible, must therefore be the mixed phase with the particular Coulomb lattice structure which has the lowest energy, since this is the only system which is in full compliance with Gibbs equilibrium conditions.

Screening and surface effects do not change the above conclusion. Phase boundaries are inevitable in two-phase systems in equilibrium, and likewise are screening effects an integral part of all charged systems in equilibrium, whether the systems consist of a single phase or two (or more) phases separated by a phase boundary. In the latter case screening and surface effects across the phase boundary become inseparable. Charged systems larger than the typical screening length can reduce their Coulomb energy by reducing the local charge density. But this does not imply that these large systems will generally resemble that of a one-component system with the electric charge frozen out. As long as the charge density does not vanish identically everywhere in the systems, the Wigner-Seitz cells will survive . On this basis the scenarios $`b`$ and $`c`$ in Fig. 2 can be excluded for any first order phase transition of such two-component systems in full thermodynamical equilibrium.

For a given phase transition a figure like Fig. 2 is very useful though. If the model describing the structured mixed phase has an energy corresponding to scenario $`b`$ or $`c`$, we can immediately conclude that the model does not give a correct description, since we already know the Maxwell-like, one-component approximation represents an excited state. This is a powerful test of whether the approximations in the treatment of Coulomb and surface energies are good enough. An important implication is that a rough upper limit on the surface tension in a given model can be estimated from this comparison. However, one should be careful not to use (reasonable) guesses for the value of the surface tension to draw conclusions about the physical possibilities of the system as it was done in . We note that in the description of the first order nuclear liquid-gas transition, the findings of Ref.  are that the structured mixed phase is indeed favored except close to the upper end of the mixed phase region. Thus, their model gives at least a physically sound description of the structures in most of the mixed phase region except close to its upper end where the details of the model, here becoming increasingly important, can be improved upon.

There are situations where scenario $`b`$, if it had been favored, could have led to buoyancy instabilities analogous to a situation where mercury is placed on top of liquid water. Assume for example that the structured mixed phase for the system shown in Fig. 1 is favored only if the pressure is less than 65 MeV/fm3 corresponding to baryon densities below 0.48 fm-3. Since hydrostatic equilibrium requires that the pressure increases monotonically from the surface to the center of the neutron star, there will be a region in the star where a denser structured mixed phase is located further from the center than a less dense Maxwell-like phase. This is a gravitationally unstable situation commonly seen in some stars. A similar situation arises within the structured mixed phase, but here the two phases in each Wigner-Seitz cell have a non-vanishing charge density, and the Coulomb force which is much stronger than gravity stabilizes the cell and secure that it remains intact.

In specific studies of possible phase transitions, reasonable models for the equations of state of the phases are crucial. For some choices of models (or parameter values within a model) it may happen that it is not possible to assure full compliance with Gibbs equilibrium conditions, even if a Maxwell construction, indicating the transition is first order, can be found. Examples relating to kaon condensation can be found in . Since the Maxwell construction represents an excited state in a two-component first order phase transition, the model must be unphysical since by construction it cannot describe the system in full thermodynamical equilibrium. We note that the authors of Ref.  are aware of this.

There is a single case where the transition can be first order without there being mixed phase. If it happens that also the electron chemical potential is continuous at the constant pressure transition of the Maxwell-like system, then there will be no mixed phase (or more precise, the radial extent of the mixed phase region has shrunk to nothing). But this is by no means related to screening effects or energy considerations between the Maxwell-like system and the mixed phase. It would be a completely coincidental property of the two phases and the equations of state describing them, and therefore it has vanishing probability.

We have, based on Gibbs equilibrium conditions, presented a new picture showing that in neutron stars, the structured mixed phase will always be energetically favored for arbitrary first order phase transitions involving systems in full thermodynamical equilibrium with baryon number and electric charge as the two independent components. Therefore, models that impose local charge neutrality on such two-component systems and thereby make the phase transition resemble that of one-component systems cannot be physically sound. Additionally, we have shown a way of checking whether the applied models are physically reasonable.

This conclusion is likely to have important consequences for a broad range of neutron star properties, e.g., glitch phenomena, transport and superfluid properties, r-mode instabilities, and the braking index, see  and references therein. Details concerning these phenomena depend of course on the specific properties of any first order phase transition taking place in the dense interior of neutron stars, especially the radial extent of the mixed phase region, the density at which it starts forming, and whether unsurmountable energy barriers prevent the neutron star matter from actually achieving its equilibrium configuration.

M.B.C. thanks Jes Madsen for useful discussions. M.B.C. was supported in part by the Carlsberg Foundation and Direktør Ib Henriksens Fond. N. K. G. was supported by the Director, Office of Science, Office of High Energy and Nuclear Physics, Division of Nuclear Physics, and by the Office of Basic Energy Sciences, Division of Nuclear Sciences, of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098.


format: text

text:
# 1 Introduction

## 1 Introduction

Recently some magnetic systems have been discovered in which simple quantum mechanical principles lead to some striking macroscopic phenomena. In particular, the observation of discrete steps in the magnetization of single crystals of the compound $`Mn_{12}Ac`$ in the presence of a time-dependent magnetic field has evoked considerable theoretical interest. This has been termed variously as quantum hysteresis and quantum or resonance tunneling of magnetization.

The basic underlying physics is easy to understand . The system consists of magnetic molecules which interact only weakly with each other. Each molecule is a cluster consisting of a core tetrahedron of four $`Mn^{4+}`$ ions each with a spin of $`\frac{3}{2}`$, and an outer crown consisting of eight $`Mn^{3+}`$ ions each with spin $`2`$. The intracore as well as intracrown magnetic interactions are ferromagnetic, while the interactions between the spins in the core and the crown is antiferromagnetic (see Fig. 1). Thus each molecule is a ferrimagnetic cluster with a ground state spin of $`10`$. In the cluster, the dipolar interaction between the spins leaves only the $`M_s=10`$ and $`10`$ states degenerate. The application of a magnetic field lifts this degeneracy, resulting in a nonzero magnetization. As the field is increased, different pairs of $`M_s`$ states become degenerate at certain values of the field. At those particular fields, the presence of matrix elements between the degenerate states, provided either by a weak transverse component in the applied magnetic field or by higher order spin-spin interactions, causes tunneling between the states. This results in a jump in the magnetization. At all other values of the field where there are no energy degeneracies, the plot of magnetization vs field shows plateaus or discrete steps provided the sweep rate of the magnetic field is not too low . This is because, according to the Landau-Zener theory , the tunneling amplitude to go from one magnetization state to another is very small unless the sweeping frequency is so low that it is comparable to the matrix element connecting the two states. Such steps in the magnetization are also seen in another magnetic cluster $`Fe_8`$ , although the effect is less dramatic there than in the $`Mn_{12}`$ cluster.

In recent years, there have been many model calculations which illustrate such steps in the $`M`$ vs $`H`$ curves . These models involve either the presence of a transverse magnetic field, or higher order spin couplings which leads to a term of the type $`c(\widehat{S}_x^4+\widehat{S}_y^4)`$ allowed by the symmetry of the cluster. However, most of these calculations have been restricted to the spin-$`10`$ manifold, i.e., a total of $`21`$ states. In contrast, we have carried out an explicit calculation of the low-lying states of a $`Mn_{12}Ac`$ cluster using a Heisenberg spin model; we then consider $`117`$ of these states which involve several different values of the total spin. We have studied quantum tunneling between these low-lying states by setting up a Hamiltonian in this subspace of states which includes, besides the multipolar spin-spin interactions and a transverse magnetic field, different gyromagnetic ratios for the core and crown spins. This last interaction, which is reasonable to introduce because of the different environments around the core and crown spins, has the effect of mixing up the different spin states of the cluster when a magnetic field is applied; it is therefore essential to keep states with different values of the total spin in the calculation, rather than restrict oneself to the spin-$`10`$ ground state manifold. We have then evolved an initial state, which is taken to be the ground state with a specific value of $`M_s=S_z`$ in the absence of the magnetic field, by using the time-dependent formulation of the problem in the restricted subspace. In the following section, we formulate the problem and present some numerical details. This will be followed by a section on the results of our time evolution studies and a discussion of the results.

## 2 Theoretical Model and Computational Details

The dimension of the Fock space spanned by the magnetic cluster $`Mn_{12}Ac`$ is exactly a hundred million. When specialized to the states of interest in the low energy sector, the dimensionalities reduce considerably even though they are still large. Table 1 gives the dimensionalities of the Hilbert spaces corresponding to the various $`M_s`$ values of interest. These states can be represented by $`32`$ bit integers by associating two bits with each of the core spins of $`s=\frac{3}{2}`$ (namely, $`M_s=\frac{3}{2}00,M_s=\frac{1}{2}01,M_s=+\frac{1}{2}10`$, and $`M_s=+\frac{3}{2}11`$), and three bits with each of the outer spins with $`s=2`$ ($`M_s=2000,M_s=1001,M_s=0010,M_s=+1011`$ and $`M_s=+2100`$); these states are generated in an ascending order of the $`32`$ bit integers that represent them.

The Hamiltonian describing the spin interactions in the cluster is given by

$$\widehat{H}=J\underset{<ij>_{core}}{}\widehat{S}_i\widehat{S}_jJ\underset{<ij>_{crown}}{}\widehat{S}_i\widehat{S}_j+J^{}\underset{<i_{core}j_{crown}>}{}\widehat{S}_i\widehat{S}_j,$$
(1)
where $`J`$ is ferromagnetic, $`J^{}`$ is antiferromagnetic, and the summations are taken over nearest neighbors. We have taken $`J=1.0`$ and $`J^{}=0.2`$ in our calculations. Every crown spin has two nearest neighbors in the crown while every core spin has three nearest neighbors in the core corresponding to the tetrahedral geometry of the core spin. Every core spin also has three nearest neighbor crown spins.

The Hamiltonian matrix is set-up in the desired $`M_s`$ space, which in our case is restricted to $`M_s=9,10`$ and $`11`$. In each subspace we have obtained a few low-lying states using the Davidson algorithm . We have also calculated the spin densities and the spin-spin correlation functions in each of the states. Using the spin-spin correlation functions, we have computed the expectation value of $`S_{total}^2`$ operator, from which we have identified the total spin of the state. We also compute the total spin density of the core and crown spins in each of these states. From the total spin and the spin density of the core and crown spins in the given $`M_s`$ state, we have used the spin ladder operators to obtain the spin densities in the core and crown of the cluster in all the allowed $`M_s`$ states. These are later used in computing the magnetization response of the system. We find that the lowest five states with the above values of $`M_s`$ belong to $`S_{total}=10,11,12,10`$ and $`13`$; their energies are given in Table 2. We therefore have a total of $`117`$ low-lying states.

To study quantum tunneling we have considered all these low-lying states and the following Hamiltonian,

$$\widehat{H}=E_sD\widehat{S}_{z,total}^2+c(\widehat{S}_{x,total}^4+\widehat{S}_{y,total}^4)g_{crown}h(t)\widehat{S}_{z,crown}g_{core}h(t)\widehat{S}_{z,core}.$$
(2)
Here $`D`$ is the quadratic anisotropy factor, $`g_{crown}`$ and $`g_{core}`$ are the Landé $`g`$-factors for the crown and core spin respectively, and $`h(t)`$ is the time-dependent magnetic field. We have chosen $`D=10^3`$ and $`c=10^3`$ (in units of J) in accordance with the experimental values . We take $`g_{crown}=1.85`$ and $`g_{core}=2.0`$ in order to keep the average $`g`$ value equal to the experimental value of $`1.9`$. The constants $`E_s`$ in (2) correspond to the five lowest energies obtained from Eq. (1) and given in Table 2. The fourth order anisotropy term allows transition between states with $`\mathrm{\Delta }M_s=\pm 4`$. To study the evolution of the magnetization as a function of the applied magnetic field, we start from an axial field equal to $`1`$ (in units of $`J/\mathrm{}`$) and then slowly increase it in steps till it equals $`+1`$; the exact procedure for sweeping the field is described below.

We have studied the time evolution of the system by solving the time-dependent Schr$`\ddot{o}`$dinger equation

$$i\mathrm{}\frac{d\psi }{dt}=\widehat{H}(t)\psi .$$
(3)
We assume the system to be in the state with $`S=10,M_s=10`$ (all-spins-down state) at time $`t=0`$. This is the initial state which is time evolved according to the equation

$$\psi (t+\delta t)=e^{i\widehat{H}(t+\frac{\delta t}{2})\delta t/\mathrm{}}\psi (t).$$
(4)
The evolution is carried out by explicit diagonalization of the Hamiltonian matrix $`𝐇(t+\frac{\delta t}{2})`$, and using the resulting eigenvalues and eigenvectors to evaluate the matrix of the time evolution operator $`e^{i\widehat{H}(t+\frac{\delta t}{2})\delta t/\mathrm{}}`$. Since the Hamiltonian matrix is in a truncated basis of $`117`$ eigenstates of the magnetic cluster, it is possible to repeatedly carry out the time evolution in small time steps of size $`\mathrm{\Delta }t`$.

## 3 Results and Discussion

We have carried out a systematic investigation of the dependence of the magnetization steps on the field sweep rate, the fourth order anisotropy term, and the presence of a transverse field. Before doing that, it is useful to have an idea of the energy levels of the system as a function of the magnetic field.

Fig. 2 shows the the energy levels of the Hamiltonian in Eq. (2) with a constant axial magnetic field $`h(t)=H_z`$. We can see from that figure that there are at least three values of the axial field where level crossings occur. In the case of interacting spins with exchange constant $`J`$ and field value $`H_z=0`$, levels with opposite magnetizations $`\pm M_s`$ are degenerate. But when we introduce quantum fluctuations, such as the fourth order anisotropy term, then the magnetization is no longer conserved. Crossings occur around $`H_z=\pm 0.25`$, $`\pm 0.4`$ and $`\pm 0.5`$ (in units of $`J/\mathrm{}`$). The energy spectrum is symmetric about $`H_z=0`$. We should expect to see jumps in the magnetization value at those values of the field where the crossings occur. Besides, as we will see, the occurrence and widths of plateaus in the magnetization are strongly dependent on the field sweep rate. The number of plateaus and their locations and widths depend on the probability of tunneling from one magnetization state to another. This probability increases when the time scale of sweeping matches with the time scale of tunneling. In that case, the probability of staying in the same eigenstate is small; the state is scattered into another eigenstate which produces a step in the magnetization plot.

Following the technique used by De Raedt et al , we now study the behavior of the magnetization as the field is changed with time. The magnetic field $`H_z`$ is increased from $`1`$ to $`1`$ in steps of $`0.02`$. At each value of the field, the state is evolved for $`300`$ time steps of size $`\mathrm{\Delta }t`$. We have considered two different time step values given by $`\mathrm{\Delta }t=0.01\mathrm{}/J`$ and $`0.1\mathrm{}/J`$; thus each value of the field is kept fixed for a time equal to $`3\mathrm{}/J`$ and $`30\mathrm{}/J`$ respectively. The field sweep rate is given by $`0.02/(300\mathrm{\Delta }t)`$; we therefore have two sweep rates differing by a factor of $`10`$. At each time step, the time evolved state is used to compute the magnetization $`M`$ given by

$$M(t)=\psi (t)|\widehat{S}_{z,total}|\psi (t).$$
(5)
The magnetization at each value of the field is then taken to be the average of the magnetization computed over all the time steps for which the field is held fixed.

In Fig. 3, we show the step behavior of the magnetization with the applied field. The upper curve is for a time step equal to $`0.1`$ (in units of $`\mathrm{}/J`$), while the lower curve is for a time step of $`0.01`$. We observe jumps and steps in the magnetization plot at field values of approximately $`H_z=0.25,0.4`$ and $`0.5`$. Before the first jump in the magnetization the magnetization value remains almost constant at $`10`$. The reason why we do not see any jumps at $`H_z=0.25,0.4`$ or $`0.5`$ is because of the significant energy difference between the spin-$`10`$ ground state and the higher spin excited states; this makes the scattering probability extremely small. We observe a remarkable thing that the magnetization value seems to saturate after a certain time evolution, but it never approaches the state with all spins up. We can argue that in our model the system can only gain or loose energy by interacting with the time-dependent field, and there is no interaction with the environment through, for example, spin-phonon or nuclear spin-electron spin interactions. So even for very large $`H_z`$ the magnetization does not reach the saturation value in a finite time. However, the magnetization does reach a higher value for large fields if the sweep rate is lower, since there is more time to tunnel to the lowest energy states in that case. The inset of Fig. 3 shows the result obtained when the field is held fixed for a longer time equal to $`300\mathrm{}/J`$ corresponding to $`3000`$ time steps of size $`0.1\mathrm{}/J`$ each. Note that the plateau in the inset occurs at a different value of the magnetization compared to the plateaus in the two curves in the main figure where the sweep rates were faster. This is due to tunneling to nearly degenerate states with different values of the magnetization.

Experimentally we see transitions between states for which $`\mathrm{\Delta }m=\pm 1`$ . This is because, even in a single crystal of $`Mn_{12}Ac`$, all the clusters do not have exactly the same orientation. This would imply that the magnetic field seen by a cluster would also have a transverse component due to a slight misalignment of the axial field. To account for this fact we have added a term $`gH_x\widehat{S}_x`$ to the Hamiltonian in Eq. (2). In Fig. 4, we show the magnetization plateaus in the presence of a fixed and small transverse field equal to $`0.01J/\mathrm{}`$ when the axial field is applied in the same manner as before. For the lower sweep rate, we do see one more plateau here than in Fig. 3 where there was only a fourth order term and no transverse field.

In Fig. 5, we show the magnetization when there is only a transverse field equal to $`0.05`$ and no fourth order term. On comparing this with Fig. 4, we can see that in the presence of a higher transverse field the magnetization is much larger for positive values of the axial field; in fact, the magnetization almost reaches the maximum possible value of $`10`$ for the lower sweep rate. This is because excited states start getting populated more easily when the magnitude of the leaking term is larger. We also note that there are no intermediate plateaus for the lower sweep rate in Fig. 5 because the system has more time to evolve to the lowest energy state at each field.

It is instructive to contrast the effects of the fourth order term vs. a transverse field on the dispersion of the low-lying energy levels. In Fig. 6, we show the energy levels for two different axial fields equal to $`0`$ and $`0.26`$. The two figures on the left have the fourth term but no transverse field, while the two figures on the right have a transverse field of $`0.05`$ but no fourth order term. (For ease of comparison, we have subtracted out appropriate constants to make the ground state energy equal to zero in all the four figures). It is clear that the energy levels have a much greater dispersion in the presence of the fourth order term; this means that the gaps are larger, the tunneling amplitudes will be smaller and the magnetization plateaus would be more prominent in that case. With a transverse field present, the energy levels almost form a continuum with very small gaps; hence the system can easily make transitions from one state to another and the plateaus would be smaller. This emphasizes that the fourth order term is essential for observing magnetization plateaus, and explains why there are more plateaus in Fig. 3 than in Fig. 5.

To study the populations of different $`M_s`$ states at the different plateau regions, we show histograms in Figs. 7 (a) and (b). $`|C|^2`$ is the sum of the squares of the modulus of the coefficients for each value of $`M_s`$ in the time evolved state $`\psi (t)`$ averaged over the time for which the field is held constant. We have followed the changes in $`|C|^2`$ for different $`M_s`$ states at three values of the axial fields where the plateaus begin in Fig. 3. In Fig. 7 (a), we show the change of weights of different $`M_s`$ values at the three onset values of the field. We can clearly see that the $`\widehat{S}_x^4+\widehat{S}_y^4`$ term connects states with $`\mathrm{\Delta }M_s=\pm 4`$. This plot also explains why the saturation magnetization value never reaches $`+10`$. The $`M_s=10`$ state has a considerable weight in all three plateaus although it decreases gradually at higher magnetic fields. In Fig. 7 (b), we have studied the case when the system is experiencing a transverse field in addition to the fourth order term. In the presence of the leaking term $`gH_x\widehat{S}_x`$, the population distribution changes completely. Here we see transitions between states with $`\mathrm{\Delta }m=\pm 1`$. The population of the $`M_s=10`$ state decreases drastically from the first plateau to the third plateau. We do not see any weight for $`M_s=10`$ state in this case also.

In Figs. 8 (a) - (c), we show the time evolution of the magnetization at a constant axial field equal to $`0.74J/\mathrm{}`$. The oscillations are most complex when we have both a transverse field and a fourth order term as in Fig. 8 (b). The oscillations are somewhat simpler when we have only the fourth order anisotropy term and no transverse field as in Fig. 8 (a); we see a beat structure superimposed on regular oscillations. Finally, when we only have a small transverse field equal to $`0.05`$ and no fourth order term, the magnetization evolves as a smooth cosine function with time as shown in Fig. 8 (c); the behavior is just like that of a two-level system. Note also that the ranges of the magnetization oscillations change considerably with the sweep rate in all the three figures.

Finally, in Fig. 9, we present the weights of the different $`M_s`$ states when there is a constant axial field equal to $`0.74`$ and a transverse field equal to $`0.05`$, but no fourth order term. The weights are dramatically different for the two sweep rates; for the lower sweep rate shown in the lower panel, the system has sufficient time to reach states with higher magnetization.

## 4 Summary and Outlook

We have studied the effects of the sweep rate, fourth order anisotropic term and a transverse magnetic field on the quantum tunneling of magnetization in $`Mn_{12}Ac`$. An effort to reproduce the experimental details will require the introduction of terms which take into account the other degrees of freedom necessary for the spin system to relax back to the ground state after a scattering event. It would also be worthwhile to extend this approach to finite temperatures to investigate the thermodynamic properties of the molecule and to study the phenomenon of thermally assisted tunneling . For models of the Ising kind, a nontrivial response of the magnetization is known for an alternating field . A study of the same in our model is currently in progress. In conclusion, we would like to say that incorporating the spin-phonon , hyperfine  and dipolar interactions  to quantitatively explain the experiments will be a challenge for the future.

Acknowledgments

We thank Prof. Indrani Bose for introducing us to this problem, and the Council of Scientific and Industrial Research, India for their grant No. 01(1595)/99/EMR-II.

| $`M_s`$ value | Dimension of Hilbert Space |
| --- | --- |
| $`\pm 13`$ | 139672 |
| $`\pm 12`$ | 269148 |
| $`\pm 11`$ | 484144 |
| $`\pm 10`$ | 817176 |
| $`\pm 9`$ | 1299632 |
| $`\pm 8`$ | 1954108 |
| $`\pm 7`$ | 2785384 |
| $`\pm 6`$ | 3772176 |
| $`\pm 5`$ | 4862352 |
| $`\pm 4`$ | 5974048 |
| $`\pm 3`$ | 7003944 |
| $`\pm 2`$ | 7842070 |
| $`\pm 1`$ | 8390440 |
| $`0`$ | 8581300 |

Table 1. Different $`M_s`$ values and the corresponding Hilbert space dimensions.

| $`S_{total}`$ | Energies |
| --- | --- |
| 10 | -53.251824 |
| 11 | -52.419914 |
| 12 | -51.512498 |
| 10 | -50.910154 |
| 13 | -50.529748 |

Table 2. Energies (in units of J) of the lowest five states and the corresponding $`S_{total}`$ values for the Hamiltonian in Eq. (1).

Figure Captions

1. A schematic diagram of the exchange interactions between the $`Mn`$ ions in the $`Mn_{12}Ac`$ molecule.

2. Energy spectrum of the Hamiltonian in Eq. (2) in the presence of a time-independent axial field. Only a few low-lying energy levels are shown.

3. Evolution of magnetization at two different sweep rates. The upper curve is for a time step equal to $`0.1`$ (in units of $`\mathrm{}/J`$), while the lower curve is for a time step of $`0.01`$. The inset shows the result obtained when the field is held fixed for a time equal to $`300\mathrm{}/J`$.

4. Magnetization vs. the axial field for two different sweep rates when there is both a fourth order term and a transverse field equal to $`0.01J/\mathrm{}`$.

5. Magnetization vs. the axial field for two different sweep rates when there is only a transverse field equal to $`0.05`$ and no fourth order term.

6. Energy level distribution of the system for two different axial fields. The two figures on the left have a fourth order term and no transverse field, while the two figures on the right have a transverse field and no fourth order term.

7. Weights of different $`M_s`$ states at the onset fields of three different plateaus is shown in plot (a), and in the presence of an additional small transverse field $`H_x=0.01`$ is shown in plot (b). The fourth order term is present in both cases. Note that $`\mathrm{\Delta }M_s`$ is equal to $`4`$ and $`1`$ in cases (a) and (b) respectively.

8. Evolution of magnetization with time is shown at a constant value of the axial field equal to $`0.74`$. In (a) there is only a fourth order anisotropy term; in (b) there is both a fourth order term and a transverse field equal to $`0.01`$; in (c) there is only a transverse field equal to $`0.05`$. Note that the ranges of the oscillations change considerably with the sweep rate.

9. Weights of different $`M_s`$ states in the presence of an axial field equal to $`0.74`$ and a transverse field equal to $`0.05`$ for two different sweep rates.


text:
# Spatial Distribution of Initial Interactions in High Energy Collisions of Heavy Nuclei

## 1 Introduction

In this note we discuss the spatial distribution of initial interactions in the high energy collisions of heavy nuclei, in order to gain some insight into geometrical effects in such collisions. We apply a simple, widely used formalism (e.g. ) that is often attributed to Glauber , to calculate the density of interactions projected onto the plane transverse to the beam axis for four commonly used weighting functions: number of wounded nucleons, number of binary collisions, colliding hard spheres, and colliding disks. The wounded nucleon and binary collision calculations incorporate the usual Woods-Saxon parametrization of nuclear density. We compute the mean and dispersion of the radius and the eccentricity of the interaction zone as a function of impact parameter for symmetric collisions of heavy nuclei. Using the same calculation we tabulate the number of participants for the wounded nucleon and hard sphere paramterizations, together with the number of binary collisions, as a function of impact parameter.

The distributions calculated are implicit in all model calculations of low or high sophistication that incorporate similar parameterizations of nuclear geometry, and in that sense nothing new is presented here. However, we have not found these distributions collected and discussed explicitly in the literature, and that is the the purpose of this note. A more detailed discussion of these and related results was presented previously in .

The formalism is described briefly in section 2. Section 3 discusses the transverse profile of the interaction volume as a function of impact parameter. Section 4 presents a parametrization of the eccentricity of the interaction region for non-central collisions. Section 5 tabulates the number of participants for the wounded nucleon and hard sphere paramterizations, together with the number of binary collisons, as a function of impact parameter.

## 2 Formalism

The coordinate system is defined in Figure 1. We utilize the standard nuclear thickness function $`\mathrm{T}_\mathrm{A}`$ ,

$$\mathrm{T}_\mathrm{A}(|\stackrel{}{s}|)=𝑑z\rho _\mathrm{A}(z,\stackrel{}{s}).$$
(1)
For the nuclear density we use a Woods-Saxon distribution,

$$\rho _\mathrm{A}(r)=\rho _0\frac{1}{1+e^{(r\mathrm{R}_\mathrm{A})/\mathrm{a}}},$$
(2)
where $`r=\sqrt{s^2+z^2}`$, $`\mathrm{R}_\mathrm{A}=1.12\mathrm{A}^{1/3}`$, and $`\rho _0=0.159`$ GeV/$`\mathrm{fm}^3`$ and a=0.535 fm for $`{}_{}{}^{197}\mathrm{Au}`$. $`\rho _\mathrm{A}`$ is normalized so that $`d^2s\mathrm{T}_\mathrm{A}(|\stackrel{}{s}|)=\mathrm{A}`$.

In Fig. 1, the nuclear impact parameter is denoted b, and the distance to a point in the transverse plane from the center of either nucleus is denoted $`\mathrm{b}_\mathrm{A}`$ and $`\mathrm{b}_\mathrm{B}`$. The vector from the origin to this point is written

$$\stackrel{}{s}=\stackrel{}{\mathrm{b}_\mathrm{A}}\frac{b}{2}\widehat{x}=\stackrel{}{\mathrm{b}_\mathrm{B}}+\frac{b}{2}\widehat{x}.$$
(3)
where $`\widehat{x}`$ is the unit vector in the $`x`$ direction.

To calculate the density of interactions as a function of impact parameter for a given process, we utilize the limiting cases of high and low nucleon-nucleon cross section for that process:

* Wounded Nucleon Scaling: We consider the number of nucleons at $`\stackrel{}{s}`$ that are struck at least once by the nucleons in the oncoming nucleus, where “struck” means inelastically excited with nucleon-nucleon collision cross section111$`\sigma _{\mathrm{NN}}`$$``$ 30 mb at the SPS ($`\sqrt{s_{\mathrm{NN}}}`$=20 GeV) and 40 mb at RHIC ($`\sqrt{s_{\mathrm{NN}}}`$=200 GeV). $`\sigma _{\mathrm{NN}}`$. In the transverse projection, the density of wounded nucleons per nuclear collision is given in units $`1/\mathrm{fm}^2`$ by:

$$\frac{\mathrm{d}^2\mathrm{N}_{\mathrm{WN}}}{\mathrm{ds}^2}=\mathrm{T}_\mathrm{A}(\mathrm{b}_\mathrm{A})(1e^{\mathrm{T}_\mathrm{B}(\mathrm{b}_\mathrm{B})\sigma _{\mathrm{NN}}})+\mathrm{T}_\mathrm{B}(\mathrm{b}_\mathrm{B})(1e^{\mathrm{T}_\mathrm{A}(\mathrm{b}_\mathrm{A})\sigma _{\mathrm{NN}}}),$$
(4)
where the dependence on nuclear impact parameter is via Eq. 3.
* Binary Collision Scaling: More generally, we mean those processes with sufficiently small nucleon-nucleon cross section $`\sigma _{\mathrm{NN}}^{\mathrm{hard}}`$ that $`(\mathrm{T}_\mathrm{A}\sigma _{\mathrm{NN}}^{\mathrm{hard}})1`$ (we apply the label “hard” because we are usually refering to high momentum transfer processes, or hard scattering). Nucleon-nucleon interaction probabilities can then be summed for the total interaction probability, and the number of hard scatterings per nuclear collision goes as the number of binary nucleon-nucleon collisions. The density in the transverse projection of hard processes per nuclear collision, in units $`1/\mathrm{fm}^2`$, is then

$$\frac{\mathrm{d}^2\mathrm{N}_{\mathrm{hard}}}{\mathrm{ds}^2}=\sigma _{\mathrm{NN}}^{\mathrm{hard}}\mathrm{T}_\mathrm{A}(\mathrm{b}_\mathrm{A})\mathrm{T}_\mathrm{B}(\mathrm{b}_\mathrm{B}).$$
(5)
Integrating over the transverse plane (and changing integration variables), the total number of hard processes per nuclear collision is given as a function of nuclear impact parameter b by:

$$\mathrm{N}_{\mathrm{hard}}(b)=\sigma _{\mathrm{NN}}^{\mathrm{hard}}d^2s\mathrm{T}_\mathrm{A}(|\stackrel{}{s}|)\mathrm{T}_\mathrm{B}(|\stackrel{}{b}\stackrel{}{s}|)\sigma _{\mathrm{NN}}^{\mathrm{hard}}\mathrm{T}_{\mathrm{AB}}(b)$$
(6)

## 3 Transverse Profile of the Interaction Volume for Non-central Collisions

The azimuthal anisotropy of momentum distributions in nuclear collisions can be related to the orientation of the event plane (the x-z plane in Figure 1) for noncentral collisions over a wide range of energies , and has been quantitatively studied up to the highest energy nuclear collisions at the SPS  and RHIC . The lowest order harmonics are referred to as directed and elliptic flow, and are characterized by coefficients $`\mathrm{v}_1`$ and $`\mathrm{v}_2`$ respectively in a Fourier expansion of the azimuthal angle or momentum distributions .

Elliptic flow at midrapidity has received particular attention recently because of a possible connection to the Equation of State (see discussion and references in ). To help distinguish dynamics from purely geometrical effects, it has been suggested  that the measured $`\mathrm{v}_2`$, the elliptic anisotropy, be scaled by the eccentricity of the reaction volume. This is defined to be

$$ϵ\frac{<y^2><x^2>}{<y^2>+<x^2>}\frac{b}{2\mathrm{R}_\mathrm{A}}$$
(7)
where $`<\mathrm{}>`$ indicates the spatial average over the transverse plane weighted by a density such as that of wounded nucleons (equation 4). The approximation is the ratio of the lengths of the axes of the overlap region in Figure 1, $`(y_{x=0})/(x_{y=0})`$, not weighted by nuclear density .

We calculate the transverse density distribution of the reaction volume as a function of impact parameter, utilizing four different weighting functions:

* Wounded Nucleons: The transverse density profile is calculated using a Woods-Saxon density distribution and Eq. 4, appropriate for the bulk of particle production. We use $`\sigma _{\mathrm{NN}}`$=30 mb, but also investigate the sensitivity of the computed quantities to this parameter.
* Binary Collisions: The transverse density profile is calculated using a Woods-Saxon density distribution and Eq. 5, appropriate for hard processes such as jet and $`\mathrm{J}/\mathrm{\Psi }`$ production.
* Hard Sphere: The transverse density profile is calculated for colliding sharp-edged spheres, with the density defined as $`\mathrm{T}_\mathrm{A}`$+$`\mathrm{T}_\mathrm{B}`$. This corresponds to the limit of the Wounded Nucleon density in which the Woods-Saxon parameter $`a`$ in Eq. 2 is small and $`\sigma _{\mathrm{NN}}`$ in Eq. 4 is large. The radius of the hard sphere corresponding to an Au nucleus is 7.24 fm, increased from the Woods-Saxon value of 6.52 fm in order to obtain the same total interaction cross section.
* Two Dimensional: The transverse density profile is simply the area of overlap region in Fig. 1, with uniform weighting (i.e. without taking into account the nuclear thickness in the z-direction and corresponding variation in density when projected onto the plane). The radius corresponding to a Au nucleus is 7.24 fm, as in the Hard Sphere calculation.

We first present the transverse density profiles graphically, and then calculate moments to compare the distributions quantitatively. Figures 2 to 5 show the transverse density profiles for Au+Au collisions as linear contour plots in one quandrant of the coordinate system defined in Fig. 1, for the four weighting functions at different impact parameters. The Wounded Nucleon and Binary Collision profiles have very similar shape, with the Binary Collision distribution falling off slightly faster from the origin. The less realistic Hard Sphere and Two Dimensional calculations exhibit larger aspect ratios for noncentral collisions and unphysical sharp declines in the density at the boundaries of the overlap region.

Figure 6 shows the mean transverse radius ($`<|\stackrel{}{\mathrm{s}}|>`$ for Fig. 1) and its dispersion ($`\sqrt{<|\stackrel{}{\mathrm{s}}|^2><|\stackrel{}{\mathrm{s}}|>^2}`$) for the four weighting functions, as a function of impact parameter for Au+Au collisions. As is seen in the density profiles themselves, the Two Dimensional and Hard Sphere functions give larger mean radii than the Wounded Nucleon and Binary Collision functions. The dispersion in the radius is similar for all four functions, and is a weak function of impact parameter. The values at large impact parameter are dominated by the treatment of the nuclear surface.

Note that in the case of Binary Collision weighting, the mean (and median) radius is about 3.5 fm for the most central collisions. In other words, about one half of all produced e.g. $`\mathrm{J}/\mathrm{\Psi }`$ or jets are generated farther than 3.5 fm. in the transverse plane from the center of the reaction zone in central collisions, with a distribution in radius which has a half-width of about 1.5 fm. Rather few of the $`\mathrm{J}/\mathrm{\Psi }`$ or jets are produced near the center of the reaction, due simply to the geometry of nucleus-nucleus collisions.

Figure 7 shows the eccentricity (Eq. 7) of the interaction region as a function of impact parameter for Au+Au collisions, for the four weighting functions. Also shown is the approximation b/(2$`\mathrm{R}_\mathrm{A}`$) from Eq. 7. The eccentricity of the Wounded Nucleon and Binary Collision models are similar and significantly smaller than those of the other models. Parametrization of $`ϵ`$ in the Wounded Nucleon Model and its dependence upon $`\sigma _{\mathrm{NN}}`$ is given in the next section.

## 4 Parametrization of $`ϵ`$ and $`S`$ in the Wounded Nucleon Model

In this section we give a more detailed discussion of the characterization of the geometry of the interaction region in the Wounded Nucleon Model. The elliptic anisotropy of the overlap region $`ϵ`$, defined in Eq. 7, is shown as a function of impact parameter in Fig. 7 for all four models for Au+Au collisions. $`ϵ`$ is parametrized for the Wounded Nucleon Model ($`\sigma _{\mathrm{NN}}`$=30 mb) as

$$ϵ(\mathrm{b})=0.0469x+2.754x^24.797x^3+4.852x^42.492x^5$$
(8)
where $`x=\mathrm{b}/(2\mathrm{R})`$ and $`2\mathrm{R}=15`$ fm. for Pb+Pb collisions and $`2\mathrm{R}=14.7`$ fm. for Au+Au collisions.

Increase of $`\sigma _{\mathrm{NN}}`$ to 40 mb decreases $`ϵ`$ by 13% at b=2.5 fm and 4% at b=10.5 fm, with the change monotonic in b. The decrease of $`ϵ`$ with increasing $`\sigma _{\mathrm{NN}}`$ can be understood by noting that the value of $`ϵ`$ is dominated by surface effects. Increasing $`\sigma _{\mathrm{NN}}`$ effectively extends the surface to larger radii, with the consequence that the shape of the interaction region in the transverse projection is less eccentric than for $`\sigma _{\mathrm{NN}}`$=30 mb at the same impact parameter.

For reference, we also include a parametrization in the same model of the impact parameter dependence of the quantity $`S=\pi R_xR_y`$ . Here, $`R_x=\sqrt{<x^2>}`$, $`R_y=\sqrt{<y^2>}`$, and $`<\mathrm{}>`$ denotes the spatial average used in Equation 7. In , $`S`$ is used to calculate the particle density in the overlap region. The parametrization is

$$S(\mathrm{b})/(\pi \mathrm{R}^2)=0.164+0.0141x0.684x^2+1.026x^30.763x^4+0.284x^5$$
(9)
where, again, $`x=\mathrm{b}/(2\mathrm{R})`$ and $`2\mathrm{R}=15`$ fm. for Pb+Pb collisions and $`2\mathrm{R}=14.7`$ fm. for Au+Au collisions.

## 5 Number of Participants and Binary Collisions

We conclude with a tabulation of the number of participants and binary collisions as a function of impact parameter. The distributions in Figures 2 to 4 can be integrated to calculate the total number of nucleon participants or binary collisions per nuclear collision. These are presented in Tables 1 and 2 for Au+Au and Pb+Pb collisions. Columns 2 and 3 give the number of Wounded Nucleons (Equation 4) with $`\sigma _{\mathrm{NN}}`$=30 and 40 mb. Column 4 gives the fraction of volume of the Hard Spheres that interact, normalized to the total number of incoming nucleons. Column 5 gives the number of binary collisions for an interaction cross section of 1 $`\mu \mathrm{barn}`$ (the rate for other small cross sections is obtained by linear scaling of this number).

Tables 1 and 2 show only a weak dependence of the number of participants on $`\sigma _{\mathrm{NN}}`$, and good agreement between the number of participants calculated with the Hard Sphere and Wounded Nucleon models. The latter point is in sharp contrast to the strong difference seen between the models for $`ϵ`$ (Figure 7), and is due to the fact that the number of participants is dominated by the bulk volume, whereas $`ϵ`$ is dominated by the surface overlap.

## 6 Acknowledgements

We thank Art Poskanzer, Raimond Snellings, and Sergei Voloshin for valuable discussions.


text:
# Public key cryptography and error correcting codes as Ising models

## Abstract

We employ the methods of statistical physics to study the performance of Gallager type error-correcting codes. In this approach, the transmitted codeword comprises Boolean sums of the original message bits selected by two randomly-constructed sparse matrices. We show that a broad range of these codes potentially saturate Shannon’s bound but are limited due to the decoding dynamics used. Other codes show sub-optimal performance but are not restricted by the decoding dynamics. We show how these codes may also be employed as a practical public-key cryptosystem and are of competitive performance to modern cyptographical methods.

Error-correcting codes are of significant practical importance as they provide mechanisms for retrieving the original message after possible corruption due to noise during transmission. They are being used extensively in most means of information transmission from satellite communication to the storage of information on hardware devices. The coding efficiency, measured in the fraction of informative transmitted bits, plays a crucial role in determining the speed of communication channels and the effective storage space on hard-disks. Rigorous bounds Shannon  have been derived for the maximal channel capacity for which codes, capable of achieving arbitrarily small error probability, can be found. However, most existing practical error-correcting codes are significantly far from saturating this bound and the quest for more efficient error-correcting codes has been going on ever since.

One family of codes, introduced originally by Gallager Gallager , and abandoned in favour of other codes due to the limited computing facilities of the time, has recently been re-introduced MacKay , showing excellent performance with respect to most existing codes. In fact, it has recently been discovered that irregular constructions of Gallager’s code result in better performance than any other method Richardson ; casdec  and nearly saturate Shannon’s bound for infinite message size. Gallager-type methods are generally based on the introduction of random sparse matrices for generating the transmitted codeword as well as for decoding the received corrupted codeword. Various decoding methods have been successfully employed; we will mainly focus here on the leading technique of Belief Propagation (BP) pearl .

Most studies of Gallager-type codes conducted so far have been carried out via numerical simulations. Some analytical results have been obtained via methods of information theory MacKay , setting bounds on the performance of certain code types, and by combinatorial/statistical methods Richardson . Here we analyze their typical performance for several parameter choices via the methods of statistical physics, and validate the analytical solutions against results obtained by the Thouless-Anderson-Palmer (TAP) approach TAP  to diluted systems and via numerical methods.

In a general scenario, the $`N`$ dimensional Boolean message $`𝝃`$ is encoded to the $`M`$ dimensional vector $`𝑱^\mathrm{𝟎}`$ which is then transmitted through a noisy channel with flipping probability $`p`$ per bit (other noise types may also be considered). The received message $`𝑱`$ is decoded to retrieve the original message.

One can identify several slightly different versions of Gallager-type codes. The one used here, termed the MN codeMacKay  is based on choosing two randomly-selected sparse matrices $`A`$ and $`B`$ of dimensionality $`M\times N`$ and $`M\times M`$ respectively; these are characterized by $`K`$ and $`L`$ non-zero unit elements per row and $`C`$ and $`L`$ per column respectively. The finite, usually small, numbers $`K`$, $`C`$ and $`L`$ define a particular code; both matrices are known to both sender and receiver. Encoding is carried out by constructing the modulo 2 inverse of $`B`$ and the matrix $`B^1A`$ (modulo 2); the vector $`𝑱^\mathrm{𝟎}=B^1A𝝃`$ (modulo 2, $`𝝃`$ in a Boolean representation) constitutes the codeword. Decoding is carried out by taking the product of the matrix $`B`$ and the received message $`𝑱=𝑱^\mathrm{𝟎}+𝜻`$ (modulo 2). The equation

$$A𝝃+B𝜻=A𝑺+B𝝉\text{(mod 2)},$$
(1)
is solved via the iterative methods of BP MacKay  to obtain the most probable Boolean vectors $`𝑺`$ and $`𝝉`$; BP methods in this context have recently been shown to be identical to a TAP based solution of a similar physical system us\_sourlas .

The similarity between error-correcting codes of this type and Ising spin systems was first pointed out by SourlasSourlas , who formulated the mapping of a simpler code onto an Ising spin system Hamiltonian. To facilitate the current investigation we first map the problem onto that of an Ising model with finite connectivity. We employ the binary representation $`(\pm 1)`$ of the dynamical variables $`𝑺`$ and $`𝝉`$ and of the vectors $`𝑱`$ and $`𝑱^\mathrm{𝟎}`$ rather than the Boolean $`(0,1)`$ one; the vector $`𝑱^\mathrm{𝟎}`$ is generated by taking products of the relevant binary message bits $`J_\mu ^0=_{i\mu }\xi _i`$, which correspond to the non-zero elements of $`B^1A`$, producing a binary version of $`𝑱^\mathrm{𝟎}`$. As we use statistical mechanics techniques, we consider the message and codeword dimensionality ($`N`$ and $`M`$ respectively) to be infinite, keeping the ratio $`R=N/M`$, which constitutes the code rate, finite. Using the thermodynamic limit is natural as Gallager-type codes are used to transmit long ($`10^410^5`$) messages, where finite size corrections are likely to be negligible. We examine the Hamiltonian

$$=\underset{\mu ,\sigma }{}𝒟_{\mu \sigma }\delta [1;𝒥_{\mu \sigma }\underset{i\mu }{}S_i\underset{j\sigma }{}\tau _j]\frac{F_s}{\beta }\underset{i=1}{\overset{N}{}}S_i\frac{F_\tau }{\beta }\underset{j=1}{\overset{M}{}}\tau _j.$$
(2)
The tensor product $`𝒟_{\mu \sigma }𝒥_{\mu \sigma }`$, where $`𝒥_{\mu \sigma }=_{i\mu }\xi _i_{j\sigma }\zeta _j`$ and $`\sigma =j_1,\mathrm{}j_L`$, is the binary equivalent of $`A𝝃+B𝜻`$, treating both signal ($`𝑺`$ and index $`i`$) and noise ($`𝝉`$ and index $`j`$) simultaneously. Elements of the sparse connectivity tensor $`𝒟_{\mu \sigma }`$ take the value 1 if the corresponding indices of both signal and noise are chosen (i.e., all corresponding indices of $`A`$ and $`B`$ are 1) and 0 otherwise; it has $`C`$ unit elements per $`i`$-index and $`L`$ per $`j`$-index, representing the system’s degree of connectivity. The $`\delta `$ function provides $`1`$ if the selected sites’ product $`_{i\mu }S_i_{j\sigma }\tau _j`$ is in disagreement with $`𝒥_{\mu \sigma }`$, recording an error, and $`0`$ otherwise. Note that this term is not frustrated, as there are $`M+N`$ degrees of freedom and $`M`$ constraints (1). The two terms on the right represent our prior knowledge in the case of biased messages $`F_s`$ and of the noise level $`F_\tau `$, and require assigning certain values to these additive fields. The choice of $`\beta \mathrm{}`$ imposes the restriction of Eq.(1), while the last two terms remain finite. Note that the noise dynamical variables $`𝝉`$ are irrelevant to measuring the retrieval success $`m=\frac{1}{N}_{i=1}^N\xi _i\text{sign}S_i_\beta _\xi .`$ The latter monitors the normalized mean overlap between the Bayes-optimal retrieved message, corresponding to the alignment of $`S_i_\beta `$ to the nearest binary valueSourlas , and the original message; the subscript $`\beta `$ denotes thermal averaging. The selection of elements in $`𝒟`$ introduces disorder to the system; we calculate the partition function $`𝒵(𝒟,𝑱)=\text{Tr}_{\{𝑺,𝝉\}}\mathrm{exp}[\beta ]`$ averaged over $`𝒟`$, $`𝝃`$ and $`𝜻`$ using the replica method us\_sourlas . We employ the replica symmetry ansatz to obtain a set of saddle point equations with respect to the emerging continuous order parameters, representing local field probability distributions and the respective conjugate distributions us\_gallager .

For unbiased messages and either $`K3`$ $`(L2)`$ or $`L3`$ $`(K2)`$ we obtain both the ferromagnetic and paramagnetic solutions either by applying the TAP approach or by solving the saddle point equations numerically. The former was carried out at the values of $`F_\tau `$ and $`F_s=0`$) which correspond to the true noise and input bias levels (for unbiased messages $`F_s=0`$) and thus to Nishimori’s conditionNishimori . This is equivalent to having the correct prior within the Bayesian framework Sourlas .

The most interesting quantity to examine is the maximal code rate, for a given corruption process, for which messages can be perfectly retrieved. This is defined in the case of $`K,L3`$ by the value of $`R=K/C=N/M`$ for which the free energy of the ferromagnetic solution becomes smaller than that of the paramagnetic solution, constituting a first order phase transition. The critical code rate obtained $`R_c=1H_2(p)=1+\left(p\mathrm{log}_2p+(1p)\mathrm{log}_2(1p)\right),`$ coincides with Shannon’s capacity.

The MN code for $`K,L3`$ seems to offer optimal performance. However, the main drawback is rooted in the co-existence of the stable $`m=1,0`$ solutions, which implies that from most initial conditions the system will converge to the undesired paramagnetic solution. Studying the ferromagnetic solution numerically shows a highly limited basin of attraction, which becomes smaller as $`K`$ and $`L`$ increase, while the paramagnetic solution at $`m=0`$ always enjoys a wide basin of attraction.

Studying the case of $`K=L=2`$ , indicates the existence of paramagnetic and ferromagnetic solutions depicted in the inset of Fig.1. For corruption probabilities $`p>p_s`$ one obtains either a dominant paramagnetic solution or a mixture of ferromagnetic ($`m=\pm 1`$) and paramagnetic ($`m=0`$) solutions. Reliable decoding may only be obtained for $`p<p_s`$, which corresponds to a spinodal point, where a unique ferromagnetic solution emerges at $`m=1`$ (plus a mirror solution at $`m=1`$). Initial conditions for both simulations (TAP/BP) and the numerical solutions were chosen almost randomly, with a slight bias of $`𝒪(10^{12})`$, in the initial magnetization. The results obtained point to the existence of a unique pair of global solutions to which the system converges (below $`p_s`$) from all initial conditions.

The main question that emerges is the possibility of producing more complicated constructions for which the spinodal point is closer to Shannon’s critical flip rate. This has been based mainly on the introduction of irregular constructions Richardson ; casdec , where the number of unit elements per row/column in the matrices $`A`$ and $`B`$ is not fixed. Analytical investigation irregular  aimed at optimising the construction are yet to provide a principled method for carrying out the optimisation.

The study of parity check codes and the insight gained from the analysis led us to suggest the potential use of a similar system as a public-key cryptosystem crypto .

Public-key cryptography is based on a distribution of a public key which may be used to encrypt messages in a manner that can only be decrypted, in practical time scales, by the service provider. Several quite safe and efficient cryptosystems are currently in use such as RSA, elliptic-curves and the McEliece cryptosystem book , most of which are based on number theory methods. Public-key cryptography plays an important role in many aspects of modern information transmission, for instance, in the areas of electronic commerce and internet-based communication. It enables the service provider to distribute a public key which may be used to encrypt messages in a manner that can only be decrypted by the service provider.

In the suggested cryptosystem, a plaintext represented by an $`N`$ dimensional Boolean vector $`𝝃(0,1)^N`$ is encrypted to the $`M`$ dimensional Boolean ciphertext $`𝑱`$ using a predetermined Boolean matrix $`G`$, of dimensionality $`M\times N`$, and a corrupting $`M`$ dimensional vector $`𝜻`$, whose elements are 1 with probability $`p`$ and 0 otherwise, in the following manner $`𝑱=G𝝃+𝜻,`$ where all operations are (mod 2). The matrix $`G`$ and the probability $`p`$ constitute the public key; the corrupting vector $`𝜻`$ is chosen at the transmitting end. The matrix $`G`$, which is at the heart of the encryption/decryption process is constructed by choosing two randomly-selected sparse matrices $`A`$ ($`M\times N`$) and $`B`$ ($`M\times M`$), and a dense matrix $`D`$ ($`N\times N`$), defining $`G=B^1AD\text{(mod 2)}.`$ The matrices $`A`$ and $`B`$ are generally characterised by $`K`$ and $`L`$ non-zero unit elements per row and $`C`$ and $`L`$ per column respectively; all other elements are set to zero. The finite, usually small, numbers $`K`$, $`C`$ and $`L`$ define a particular cryptosystem. The dense invertible Boolean matrix $`D`$ is arbitrary and is added for improving the system’s security. It may be constructed as $`D=TP`$, where $`T`$ and $`P`$ are $`N\times N`$ triangular and random permutation matrices respectively, for minimising the computational costs. All matrices are known only to the authorised receiver. Suitable choices of probability $`p`$ will depend on the maximal achievable rate for the particular cryptosystem as discussed below.

The authorised user may decrypt the ciphertext $`𝑱`$ by taking the (mod 2) product $`B𝑱=A(D𝝃)+B𝜻`$, and finding the most probable solution to Eq.(1) using the methods of BP; obtaining the estimate of $`𝝃`$ is obtained by taking the product of the $`(D𝝃)`$ estimate and $`D^1`$. Studying the case of $`K=L=2`$ and $`p<p_s`$ we learned that iterative BP decoding converges to the ferromagnetic solution from any initial conditions. Cryptosystems with other $`K,L`$ values generally suffer from a decreasingly small basin of attraction as $`K`$ and $`L`$ increase, although specific matrices with higher $`K`$ and $`L`$ values (such as in casdec ) may still be used successfully.

The cryptosystem offers a guaranteed convergence to the plaintext solution, in the thermodynamic limit $`N\mathrm{}`$, as long as $`p<p_s`$. The main consequence of finite plaintexts would be a decrease in the allowed corruption rate. Experimental results with systems size as small as $`N=1000`$ still show good performance.

The unauthorised receiver, on the other hand, faces the task of decrypting the ciphertext $`𝑱`$ knowing only $`G`$ and $`p`$. The straightforward attempt to try all possible $`𝜻`$ constructions is clearly doomed, provided that $`p`$ is not vanishingly small, giving rise to only a few corrupted bits. We study the problem by exploiting the similarity between the task at hand and the error-correcting model suggested by SourlasSourlas . In this case, the matrix $`G`$ generated in the case of $`K=L=2`$ is dense and has a certain distribution of unit elements per row. The fraction of rows with a low (finite, not of $`𝒪(N)`$) number of unit elements vanishes as $`N`$ increases, allowing one to approximate this scenario by the diluted random energy model studied in us\_sourlas .

To investigate the typical properties of this (frustrated) model, we calculate again the partition function and the free energy by averaging over the randomness in choosing the plaintext, the corrupting vector and the choice of the random matrix $`G`$ (being generated by a product of two sparse random matrices). To assess the likelihood of obtaining spin-glass/ferromagnetic solutions, we calculated the free energy landscape (per plaintext bit - $`f`$) as a function of overlap $`m`$. This can be carried out straightforwardly using the analysis of us\_gallager , and provides the golf-course-like energy landscape with a relatively flat area around the one-step replica symmetry breaking (frozen) spin-glass solution and a very deep but extremely narrow area, of $`𝒪(1/N)`$, around the ferromagnetic solution crypto .

It is worthwhile mentioning that this free-energy landscape may be related directly to the marginal posterior $`P(S_i=1|𝑱)1iN`$ and is therefore indicative of the difficulties in obtaining ferromagnetic solutions when the starting point for the search is not infinitesimally close to the original plaintext (which is clearly highly unlikely). Numerical studies of similar energy landscapes show that the time required increases exponentially with the system size Parisi .

Most attacks on this cryptosystems, by an unauthorised user, will face the same difficulty: without explicit knowledge of the current plaintext and/or the decomposition of $`G`$ to the matrices $`A`$, $`B`$ and $`D`$ it will require an exponentially long time to decipher a specific ciphertext. We investigated attacks of several types, some of which appear in crypto , concluding that the suggested system is secure.

A brief comparison of our method and the leading technique of RSA book  shows that: 1) RSA decryption takes $`𝒪(N^3)`$ operations while our method only requires $`𝒪(N\mathrm{log}N)`$ operations. 2) Encryption costs are of $`𝒪(N^2)`$ (as in RSA); inverting the matrices $`B`$ and $`D`$ is carried out only once and is of $`𝒪(N^3)`$. Two drawbacks of our method: 1) The public key is a dense matrix of dimensionality $`M\times N`$. However, as public key transmission is carried out only once we do not expect it to be of great significance. 2) The ciphertext/plaintext bit ratio is greater than one (as is the case in RSA). Choosing the $`N/M`$ ratio is in the hands of the user and is related to the security level required. In addition, the increased transmission time is compensated by a very fast decryption and the added robustness against noise.

We discussed the relation between Ising models, certain error-correcting codes and public-key cryptosystems. Important aspects that are yet to be investigated include the relation between our results and the bounds obtained in the information theory literature for error-correcting codes, finite size effects and methods for alleviating the drawbacks of the new cryptosystem.

Support by JSPS-RFTF (YK), The Royal Society and EPSRC-GR/N00562 (DS) is acknowledged.


text:
# Angular structure of diffraction dissociation at LHC

## Introduction

During recent years CERN, DESY and FNAL have been producing interesting results on diffractive production in hadron and deep-inelastic processes . Discovery of hard diffraction at CERN S$`\overline{p}p`$S and diffractive events in the deep-inelastic scattering at HERA were among the most surprising results obtained recently. Significant fraction of high-$`t`$ events among the diffractive events in deep-inelastic scattering and hadron-hadron interactions were also observed at HERA and Tevatron respectively. These experimental discoveries renewed interest in the experimental and theoretical studies of the diffractive production processes.

The understanding of the diffractive interactions plays fundamental role in the studies of a high-energy limit of the modern strong interaction theory — QCD. This field is most interesting one since there are no firmly established computational methods within QCD. However, the progress in understanding of high–energy basic features and the underlying dynamics of diffraction mechanism can be traced in the recent theoretical papers (cf. review ).

It is of high importance to have the experimental data at highest possible energy and there is no doubt that soft and hard diffractive processes should be measured at the LHC. Besides elastic scattering (total and elastic cross-sections) it is important to measure the inelastic diffractive final states, jet production, hadron transverse momentum distribution, strange and charmed particle productions. The single diffraction dissociation is a most simple inelastic diffractive process and studies of soft and hard final states in this process could become a next step after elastic measurements. Such measurements are crucial for understanding of the microscopic nature of driving mechanism called Pomeron, its possible parton structure, high-energy limit of strong interactions and approaching the asymptotia.

At the currently available energies there were obtained several sets of the experimental data for single diffraction production process

$$h_1+h_2h_1+h_2^{}$$
(1)
when the hadron $`h_2`$ is excited to the state $`h_2^{}`$ with invariant mass $`M`$ and the same quantum numbers. Its subsequent decay results in the multiparticle final state. The inclusive differential cross–section shows a simple dependence on the invariant mass $`M`$:

$$\frac{d\sigma _{diff}}{dM^2}\frac{1}{M^2}.$$
(2)
Meanwhile, energy dependence of the diffractive production cross–section $`\sigma _{diff}(s)`$ is not so evident from the experimental data. The data obtained at FNAL seem to demonstrate a growing diffractive cross–section but some data obtained at CERN may indicate a falling diffractive cross–section. Resolution of this problem is important in particular for the study of the role of unitarity in the inelastic diffraction and its asymptotical properties.

The angular dependence of the diffractive dissociation is even less clear. Similarity between elastic and inelastic diffraction in the $`t`$-channel approach suggests that the latter one would have similar to elastic scattering behavior of the differential cross-section. However, it cannot be taken for granted and e.g. transverse momentum distribution of diffractive events in the deep-inelastic scattering at HERA shows a power-like behavior with no apparent dips . Similar behavior was observed also in the hadronic diffraction dissociation process at CERN  where also no dip and bump structure was observed. Angular dependence of diffraction dissociation together with the measurements of the differential cross–section in elastic scattering would allow to determine the geometrical properties of elastic and inelastic diffraction, their similar and distinctive features and origin.

The experimental regularities of diffractive production can be described in the framework of different approaches, e.g. models for inelastic diffraction based on $`s`$-channel absorptive unitarity were considered in , optical approach was used in , dipole Pomeron and effective Pomeron with energy-dependent intercept were applied in Refs.  and , respectively. $`1/M^2`$ dependence is naturally described by the triple–pomeron diagrams in the framework of Regge–model . The proposed in Ref.  similarity between the pomeron and photon exchanges allowed to calculate diffractive dissociation cross–section in terms of structure function $`\nu W_2`$ measured in deep inelastic lepton scattering. The explanation of $`M^2`$–dependence in the framework of optical model considering diffractive dissociation as a bremsstrahlung where virtual quanta are released from a strong field was made in Ref. . Explanation of the smallness of a large mass diffraction as a result of the strong non-perturbative interaction of gluons was given in . The list of the above references is certainly far from being complete and many theoretical papers devoted to the various aspects of diffractive production were not mentioned.

In Ref.  for the description of single diffractive processes we used geometrical notions on quark scattering in approach based on unitarity for the scattering amplitude and chiral quark model for hadron structure. Motivated by the importance of this process for the study of long distance dynamics we have shown how the energy and $`M^2`$–dependencies can be obtained in the approach to hadron interactions with account for the spontaneous breaking of chiral symmetry and presence of a quark condensate inside a hadron.

In this paper we consider transverse momentum distribution in the framework of the above approach and show that the diffraction cone would be suppressed and disappears at LHC energy in the production process $`p+pp+X`$.

## 1 Amplitude of the diffractive production

To construct the amplitude of the diffractive production process let us remind that the unitary equation for the scattering amplitude

$$\text{Im}f(s,b)=|f(s,b)|^2+\eta (s,b)$$
(3)
allows one to express the inelastic channel contribution

$$\eta (s,b)=\underset{n}{}\sigma _n(s,b)$$
(4)
through the $`U`$-matrix

$$\eta (s,b)=\text{Im}U(s,b)|1iU(s,b)|^2,$$
(5)
and to get respectively the total inelastic cross-section

$$\sigma _{inel}(s)=8\pi _0^{\mathrm{}}b𝑑b\frac{\text{Im}U(s,b)}{|1iU(s,b)|^2}.$$
(6)
The quantity $`\text{Im}U(s,b)`$ can be represented in the following form

$$\text{Im}U(s,b)=\underset{n}{}\overline{U}_n(s,b),$$
(7)
where

$$\overline{U}_n(s,b)=𝑑\mathrm{\Gamma }_n|U_n(s,b,\{\xi _n\})$$
(8)
and $`d\mathrm{\Gamma }_n`$ is an element of $`n`$-particle phase space volume and $`\{\xi _n\}`$ is the set of kinematical variables related to the $`n`$–particle final state. Sum in the right hand side of this equation runs over all inelastic final states which include as well diffractive as non-diffractive ones. To obtain the cross-section of the diffractive dissociation process (1) we should single out in this sum the final states corresponding to the process (1) . Let for simplicity consider the case of pure imaginary $`U`$-matrix, i.e. $`UiU`$. Then we can represent $`d\sigma _{diff}/dM^2`$ in the form

$$\frac{d\sigma _{diff}}{dM^2}=8\pi _0^{\mathrm{}}b𝑑b\frac{U_{diff}(s,b,M)}{[1+U(s,b)]^2}$$
(9)
where $`U_{diff}(s,b,M)`$ includes contributions from all the final states $`|n_{diff}`$ which result from the decay of the excited hadron $`h_2^{}`$ of mass $`M`$: $`h_2^{}|n_{diff}`$. The corresponding impact parameter amplitude $`F_{diff}(s,b,M)`$ can be written in this pure imaginary case as a square root of the cross-section, i.e.

$$F_{diff}(s,b,M)=\sqrt{U_{diff}(s,b,M)}/[1+U(s,b)]$$
(10)
and the amplitude $`F_{diff}(s,t,M)`$ is

$$F_{diff}(s,t,M)=\frac{is}{\pi ^2}_0^{\mathrm{}}b𝑑bJ_0(b\sqrt{t})\sqrt{U_{diff}(s,b,M)}/[1+U(s,b)].$$
(11)
In Refs.  we used the notions of effective chiral quark model for the description of elastic scattering and diffractive production. Different aspects of hadron dynamics were accounted in the framework of effective Lagrangian approach. The picture used for the hadron structure implies that overlapping and interaction of peripheral condensates at hadron collision occurs at the first stage. In the overlapping region the condensates interact and as a result the massive quarks appear. Those quarks are transient ones: they are transformed back into the condensates of the final hadrons in elastic scattering and the diffraction dissociation processes. In elastic scattering each of the constituent valence quarks located in the central part of the hadron is supposed in the model to scatter in a quasi-independent way by the produced virtual quark pairs at given impact parameter and by the other valence quarks.

For consideration of the diffractive production at the quark level we have extended the picture of hadron interaction in case of elastic scattering . Since the constituent quark is an extended object there is a non–zero probability of its excitation at the first stage of hadron interaction when peripheral condensates interact. Therefore it seems rather natural to assume that the origin of diffractive production process is an excitation of one of the valence quarks in colliding hadron: $`QQ^{}`$, its subsequent scattering and then decay into the final state. The excited constituent quark is scattered similar to other valence quarks in a quasi-independent way. The function $`U_{diff}(s,b,M)`$ can be represented then as a product

$$U_{diff}(s,b,M)=\underset{Q=1}{\overset{N1}{}}f_Q(s,b)f_Q^{}(s,b,M_Q^{}),$$
(12)
where $`M_Q^{}`$ is the mass of excited constituent quark. This mass is taken to be proportional to the mass $`M`$ of excited hadron $`h_2^{}`$. This assumption is based on the additivity of constituent quark masses in the hadron and the absence of the diffractive radiation from the other constituent quarks.

In the model the $`b`$–dependence of the amplitudes $`f_Q`$ and $`f_Q^{}`$ is related to formfactor of the constituent quark and excited constituent quark respectively. The strong interaction radius of constituent quark is determined by its mass. We suppose that the same is valid for the size of excited quark, i.e. $`r_Q^{}=\xi /M_Q^{}`$. The expression for $`U_{diff}(s,b,M)`$ can be rewritten then in the following form:

$$U_{diff}(s,b,M)=g^{}U(s,b)\mathrm{exp}[(M_Qm_Q)b/\xi ],$$
(13)
where constant $`g^{}`$ is proportional to the relative probability of excitation of the constituent quark. The value of $`g^{}`$ is a non-zero one, however, $`g^{}<1`$ since we expect that the excitation of any constituent quark has lower probability compared to the probability for this quark to stay unexcited. The excited quark is unstable and its subsequent decay is associated with the decay of hadron $`h_2^{}`$ into the multiparticle final state $`|n_{diff}`$. The expression for $`U(s,b)`$ is the following :

$$U(s,b)=g(s)\mathrm{exp}(\stackrel{~}{M}b/\xi )\stackrel{~}{G}\left[1+\alpha \frac{\sqrt{s}}{m_Q}\right]^N\mathrm{exp}(\stackrel{~}{M}b/\xi ),$$
(14)
where $`\stackrel{~}{M}=_{Q=1}^Nm_Q`$ and $`N`$ is the total number of the constituent quarks in colliding hadrons.

## 2 Total and differential cross-sections of the diffractive production

The cross-section of diffractive dissociation process is given by Eq. (9) and has the following $`s`$ and $`M^2`$ dependence

$$\frac{d\sigma _{diff}}{dM^2}\frac{8\pi g^{}\xi ^2}{(M_Q^{}m_Q^2)^2}\eta (s,0)\frac{8\pi g^{}\xi ^2}{M^2}\eta (s,0)$$
(15)
Thus, we have a familiar $`1/M^2`$ dependence of the diffraction cross-section which is related in our model to the geometrical size of excited constituent quark. The energy dependence of single diffractive cross-section has the form

$$\sigma _{diff}(s)=8\pi g^{}\xi ^2\eta (s,0)_{M_0^2}^{M_1^2}\frac{dM^2}{M^2}=8\pi g^{}\xi ^2\eta (s,0)\mathrm{ln}\frac{s(1x_1)}{M_0^2},$$
(16)
where $`x_1`$ is the lower limit of the relative momentum of hadron $`h_1`$ and corresponds to the experimental constraint on diffractive process $`x_10.80.9`$. Eq. (16) shows that the total cross-section of diffractive dissociation has a non-trivial energy dependence which is determined by the contribution of inelastic channels into unitarity equation at zero value of impact parameter. The dependence of $`\eta (s,0)`$ is determined by Eq. (5), where $`U(s,b)`$ is given by Eq. (14).

At $`ss_0`$, ($`s_0`$ is determined by equation $`|U(s_0,0)|=1`$) $`\eta (s,0)`$ increases with energy. This increase as it follows from Eq. (14) and from the experimental data  is rather slow one. However at $`ss_0`$, $`\eta (s,0)`$ reaches its maximum value $`\eta (s,0)=1/4`$ and at $`s>s_0`$, the function $`\eta (s,0)`$ decreases with energy and at asymptotical energies the inelastic diffraction cross section drops to zero. But at the LHC energy $`\sqrt{s}=14`$ $`GeV`$ the single diffractive inelastic cross-sections can reach the value of $`2.4`$ mb  and therefore might be quite significant.

Hence, it worth to consider the structure of the corresponding angular distribution. The corresponding amplitude $`F_{diff}(s,t,M)`$ can be calculated analytically. To do so we continue the amplitudes $`F_{diff}(s,\beta ,M),\beta =b^2`$, to the complex $`\beta `$–plane and transform the Fourier–Bessel integral over impact parameter into the integral in the complex $`\beta `$ – plane over the contour $`C`$ which goes around the positive semiaxis. Then for the amplitude $`F_{diff}(s,t,M)`$ the following representation takes place:

$$F_{diff}(s,t,M)=\frac{is}{2\pi ^3}_C𝑑\beta F_{diff}(s,\beta ,M)K_0(\sqrt{t\beta })$$
(17)
where $`K_0(x)`$ is the modified Bessel function. The amplitude $`F_{diff}(s,\beta ,M)`$ has the poles in the $`\beta `$–plane determined by equation

$$1+U(s,\beta )=0.$$
(18)
The solutions of this equation can be written as

$$\beta _n(s)=\frac{\xi ^2}{\stackrel{~}{M}^2}\left\{\mathrm{ln}g(s)+i\pi n\right\},n=\pm 1,\pm 3,\mathrm{}$$
(19)
The amplitude $`F_{diff}(s,\beta ,M)`$ besides the poles has a branching point at $`\beta =0`$.

Therefore the amplitude $`F_{diff}(s,t,M)`$ can be represented as a sum of the pole contribution and the contribution of the cut:

$$F_{diff}(s,t,M)=F_{diff,p}(s,t,M)+F_{diff,c}(s,t,M)$$
(20)
Up to this point the calculation is similar to the case of elastic scattering. For elastic scattering amplitude $`F(s,t)`$ the pole and cut contributions are decoupled dynamically when $`g(s)\mathrm{}`$ at $`s\mathrm{}`$ . Contribution of the poles determines the elastic amplitude in the region $`|t|/s1(t0)`$. The amplitude in this region can be represented in a form of series over the parameter $`\tau (\sqrt{t})`$:

$$F(s,t)=s\underset{k=1}{\overset{\mathrm{}}{}}\tau ^k(\sqrt{t})\phi _k[R(s),\sqrt{t}],$$
(21)
where $`\phi _k[R(s),\sqrt{t}]`$ are the oscillating functions of the variable $`\sqrt{t}`$. The parameter $`\tau `$ decreases exponentially with $`\sqrt{t}`$:

$$\tau (\sqrt{t})=\mathrm{exp}(\frac{2\pi \xi }{\stackrel{~}{M}}\sqrt{t}).$$

This series reproduces diffraction peak and familiar dip-bump structure of the differential cross-section in elastic scattering. In the region of moderate $`t`$ it is sufficient to take into account few or even one of the terms of series Eq. 21. The differential cross-section in this region has familiar Orear behavior.

However, the situation is different in the case of diffraction production. Instead of dynamical separation of the pole and cut contribution discussed above when

$$F_p=O(s\mathrm{ln}^{1/2}g(s)),F_c=O(s[g(s)]^1),$$
(22)
we have a suppression of the pole contribution at high energies since at fixed $`t`$

$$F_{diff,p}=O(s[g(s)]^{\frac{M}{2\stackrel{~}{M}}}\mathrm{ln}^{1/2}g(s)),F_{diff,c}=O(s[g(s)]^{\frac{1}{2}}),$$
(23)
i.e. the pole contribution is suppressed at high energies where $`g(s)>1`$ since $`M>\stackrel{~}{M}`$. Therefore, at all $`t`$ values we will have

$$F_{diff}(s,t,M)F_{diff,c}(s,t,M),$$
(24)
where

$$F_{diff,c}(s,t,M)ig^{}g^{1/2}(s)(1\frac{t}{\overline{M}^2})^{3/2},$$
(25)
where $`\overline{M}=(M\stackrel{~}{M}1)/2\xi `$. This means that the differential cross-section of the diffraction production will have smooth dependence on $`t`$ with no apparent dips and bumps

$$\frac{d\sigma _{diff}}{dtdM^2}(1\frac{t}{\overline{M}^2})^3.$$
(26)
It is interesting to note that at large values of $`M\stackrel{~}{M}`$ the normalized differential cross-section $`\frac{1}{\sigma _0}\frac{d\sigma }{dtdM^2}`$ ($`\sigma _0`$ is the value of cross-section at $`t=0`$) will exhibit scaling behavior

$$\frac{1}{\sigma _0}\frac{d\sigma }{dtdM^2}=f(t/M^2),$$
(27)
and explicit form of the function $`f(t/M^2)`$ is the following

$$f(t/M^2)=(14\xi ^2t/M^2)^3.$$
(28)
The numerical parameter $`\xi `$ in the model is about 2. The function Eq. (27) against the variable $`t/M^2`$ is depicted in the Fig. 1.

The above scaling has been obtained in the model approach, however it might have a more general meaning.

The angular structure of diffraction dissociation processes given by Eq. (26) takes place at high energies where $`g(s)>1`$ while at moderate energies where $`g(s)1`$ the both contributions from poles and cut are significant. At low energies where $`g(s)1`$ the situation is similar to the elastic scattering, i.e. there diffraction cone and possible dip-bump structure should be present in the region of small values of $`t`$. In this region

$$F_{diff}(s,t,M^2)=s\underset{k=1}{\overset{\mathrm{}}{}}\tau ^k(\sqrt{t})\phi _{diff,k}[R(s),\sqrt{t},M^2],$$
(29)
where the parameter $`\tau (\sqrt{t})`$ is the same as in the elastic case and

$$\phi _{diff,k}=O(s[g(s)]^{\frac{M}{2\stackrel{~}{M}}}\mathrm{ln}^{1/2}g(s)).$$
(30)
## Conclusion

We considered behavior of the differential cross-section for single inelastic diffraction. At low and moderate energies behavior of the differential cross-section will be rather complicated and incorporate diffraction cone, Orear type and power-like dependencies.

However, at high energies a simple power-like dependence on $`t`$ is predicted. It was shown that the normalized differential cross-section has a scaling form and only depends on the ratio $`t/M^2`$ at large values of $`M^2`$.

In fact, our particular comparative analysis of the poles and cut contributions has very little with the model form of the $`U`$–matrix. The approach is based on the combination of unitarity with the diffraction production mechanism via excitation and subsequent decay of a constituent quark. This is why it may have a more general meaning.

At the LHC energy the diffractive events with the masses as large as 3 TeV could be studied. It would be interesting to check this prediction at the LHC where the scaling and simple power-like behavior of diffraction dissociation differential cross-section should be observed. Observation of such behavior would confirm the diffraction mechanism based on excitation of the complex hadronlike colorless object - constituent quark. This mechanism can in principle explain angular structure of diffraction in the deep - inelastic scattering at HERA where smooth angular dependence on the thrust transverse momentum was observed . If it is the case, then diffraction in DIS at lower energies should manifest typical soft diffractive behavior with exponential peak at small $`t`$ as it does in hadronic reactions. It would be interesting to check it at lower energies, e. g. at Jefferson Lab.

There could be envisaged various experimental configurations at the LHC; e.g. soft diffraction goes well to the interest of the TOTEM experiment, while hard diffractive final states can be measured by CMS detector and possible correlations between the features of the soft and hard diffractive processes can be obtained using combined measurements of TOTEM and CMS .

Measurements of the angular dependence of diffraction dissociation simultaneously with the measurements of the differential cross–section of elastic scattering are important for the determination of the global geometrical properties of elastic and inelastic diffraction, their similar and distinctive features and the origin and nature of the driving asymptotic mechanism. In general these studies could be an essential tool to probe non-perturbative QCD.

## Acknowledgements

We are grateful to V. A. Petrov for the interesting discussions and comments. This work was supported by RFBR (Grant No. 99-02-17995).