text:
# Baryon Masses in the 1/𝑁_𝑐 Expansion

## 1 INTRODUCTION

Large-$`N_c`$ spin-flavor symmetry has proven to be a useful symmetry for studying the properties of baryons in QCD . Baryon observables can be analyzed in an expansion involving a complete operator basis with each operator transforming according to a definite representation under the unbroken spin $``$ flavor subgroup and suppressed by a definite power of $`1/N_c`$. From the $`1/N_c`$ expansion, it is possible to determine the symmetry structure of $`1/N_c`$ corrections which account for deviations from the large-$`N_c`$ spin-flavor limit. By including additional parameters to account for explicit flavor symmetry breaking, one obtains a combined operator expansion in $`1/N_c`$ and flavor symmetry breaking. The symmetry relations which occur when one neglects individual operators in the combined $`1/N_c`$ and flavor-symmetry breaking expansion are predicted to hold at a definite order in $`1/N_c`$ and $`SU(3)`$ flavor symmetry breaking. The hierarchy of symmetry relations predicted by this procedure results in a nontrivial pattern for baryons in QCD since $`1/N_c=1/3`$ and the $`SU(3)`$ flavor symmetry breaking $`ϵm_s/\mathrm{\Lambda }_\chi `$ are comparable in magnitude, so the pattern of relations is not dominated by either $`1/N_c`$ or $`ϵ`$, and neither parameter can be neglected relative to the other.

The hierarchy of baryon mass relations predicted by the combined expansion in $`1/N_c`$ and flavor symmetry breaking has been obtained for the lowest-lying baryon spin-flavor multiplet, which is the $`\mathrm{𝟓𝟔}`$ of $`SU(6)`$ for $`N_c=3`$ colors and $`N_f=3`$ light quark flavors. This spin-flavor multiplet contains the spin-$`1/2`$ octet and spin-$`3/2`$ decuplet baryon spin $``$ flavor representations. The mass hierarchy predicted for these ground-state baryons is particularly predictive since the expansion extends to relative order $`1/N_c^3`$ in the $`1/N_c`$ expansion and to third order in $`SU(3)`$ flavor symmetry breaking. The mass spectrum of the spin-$`1/2`$ octet and the spin-$`3/2`$ decuplet baryons is in striking accord with the $`1/N_c`$ hierarchy. Suppression factors of $`1/N_c=1/3`$ are clearly evident in the data, which provides the strongest single piece of evidence for the relevance of the $`1/N_c`$ expansion for QCD baryons where $`1/N_c=1/3`$. Even more importantly, the $`1/N_c`$ expansion resolves the long-standing puzzle of why certain famous mass relations of the spin-$`1/2`$ octet and spin-$`3/2`$ decuplet baryons, such as the Gell-Mann–Okubo formula, Gell-Mann’s Equal Spacing Rule, the Okubo formula and the Coleman-Glashow relation, work so extraordinarily well, which is to say, better than one would have expected based on flavor symmetry breaking suppression factors alone. The point is that the operators in the baryon mass expansion which are highly suppressed in $`SU(3)`$ flavor symmetry breaking parameters also are suppressed by powers of $`1/N_c`$, so that, generically, mass relations which are highly suppressed in flavor are highly suppressed in $`1/N_c`$ as well.

This proceedings provides a brief description of progress that has been made in understanding the mass spectrum of QCD baryons in an $`\mathrm{𝑎𝑛𝑎𝑙𝑦𝑡𝑖𝑐}`$ expansion in $`1/N_c`$ and flavor symmetry breaking. A brief description of the $`1/N_c`$ operator expansion for baryons is provided, with the example of the operator expansion for the $`I=0`$ masses of the ground state baryons presented in detail. Evidence of the $`1/N_c`$ mass hierarchy of the baryon octet and decuplet is obvious for the $`I=0`$ baryon mass combinations, and recent experimental evidence supporting the predicted $`1/N_c`$ hierarchy in the $`I=1`$ sector is noted for the isospin-breaking Coleman-Glashow mass splitting. Finally, the extension of the $`1/N_c`$ expansion to baryons containing a single heavy quark in Heavy Quark Effective Theory (HQET) is summarized briefly, and then applied to the charm and bottom baryon mass spectra. Using the $`1/N_c`$ expansion, a prediction was made for the $`\mathrm{\Xi }_c^{}`$ mass valid to a precision of a couple MeV, which has recently been confirmed experimentally.

## 2 $`1/N_c`$ EXPANSION

A contracted $`SU(2N_f)`$ spin-flavor algebra for baryons arises in the large-$`N_c`$ limit. For the simplest case of $`N_f=2`$ light quark flavors, the spin-flavor generators are given by baryon spin $`J^i`$, isospin $`I^a`$ and commuting generators $`X^{ia}`$. The contracted large-$`N_c`$ spin-flavor algebra is given by

$`[X^{ia},X^{jb}]=0,`$
$`[J^i,X^{ja}]=iϵ^{ijk}X^{ka},[I^a,X^{ib}]=iϵ^{abc}X^{ic},`$ (1)
$`[J^i,J^j]=iϵ^{ijk}J^k,[I^a,I^b]=iϵ^{abc}I^c,`$
where the commuting generators $`X^{ia}`$ are obtained by rescaling the axial vector baryon currents $`A^{ia}`$, which grow as $`𝒪(N_c)`$, by $`1/N_c`$ and taking the large-$`N_c`$ limit, which yields a finite well-defined result:

$$X^{ia}\underset{N_c\mathrm{}}{lim}\frac{A^{ia}}{N_c}.$$
(2)
When working at finite $`N_c`$, it is useful to work with the operator

$$G^{ia}=\left(q^{}\frac{\sigma ^i\tau ^a}{4}q\right)$$
(3)
of the large-$`N_c`$ quark model, and to use the uncontracted $`SU(2N_f)`$ algebra of the quark model. In this case, $`G^{ia}/N_c`$ has finite baryon matrix elements equal to $`X^{ia}`$ in the large-$`N_c`$ limit.

The $`1/N_c`$ expansion of an arbitrary QCD $`1`$-body operator is given by

$$\overline{q}\mathrm{\Gamma }q=N_c\underset{n=0}{\overset{N_c}{}}\frac{1}{N_c^n}c_n\left(\frac{1}{N_c}\right)𝒪_n,$$
(4)
where $`𝒪_n`$ is an $`n`$-body quark operator which equals an $`n^{th}`$ order polynomial in the spin-flavor generators $`J^i`$, $`T^a`$, $`G^{ia}`$, and the operator coefficients $`c_n(1/N_c)=𝒪\left(1\right)+\mathrm{}`$ are unknown, but have a $`1/N_c`$ expansion beginning at order unity. The baryon matrix elements of the operators $`𝒪_n`$ are known, since the baryon matrix elements of the spin and flavor generators are known, as are the matrix elements of $`G^{ia}`$. (Only the leading $`𝒪(1)`$ matrix elements of the spin-flavor generators $`X^{ia}`$ are determined; the matrix elements of the spin-flavor generators at subleading orders are not, but any reasonable choice will do. The large-$`N_c`$ quark model generator $`G^{ia}`$ or the Skyrme model collective coordinate $`X_0^{ia}`$ are two possible choices.) The complete set of independent $`n`$-body operator products are classified. There is a single $`0`$-body operator $`\mathrm{𝟏}`$. The independent $`1`$-body operators are the $`SU(6)`$ spin-flavor generators $`J^i`$, $`T^a`$ and $`G^{ia}`$, and $`n`$-body operators are products of these $`1`$-body operators. Since not all $`n`$-body operator products are independent, operator identities are needed to reduce the operator products to an independent set. The complete set of operator identities is known.

## 3 BARYON MASS HIERARCHY

The baryon mass operator transforms as a spin singlet. $`SU(3)`$ flavor symmetry breaking proportional to the parameter $`ϵ`$ transforms as the eighth component of an $`SU(3)`$ octet, so the baryon mass operator expansion also contains $`SU(3)`$ flavor-symmetry breaking operators which extend to third order in $`SU(3)`$ breaking and transform as the $`\mathrm{𝟖}`$, $`\mathrm{𝟐𝟕}`$ and $`\mathrm{𝟔𝟒}`$ representations of $`SU(3)`$. The baryon $`I=0`$ masses can be organized according to their $`SU(3)`$ representation,

$$M=M^\mathrm{𝟏}+M^\mathrm{𝟖}+M^{\mathrm{𝟐𝟕}}+M^{\mathrm{𝟔𝟒}},$$
(5)
with $`1/N_c`$ operator expansions

$`M^\mathrm{𝟏}=N_c\mathrm{𝟏}+{\displaystyle \frac{1}{N_c}}J^2,`$
$`M^\mathrm{𝟖}=T^8+{\displaystyle \frac{1}{N_c}}\{J^i,G^{i8}\}+{\displaystyle \frac{1}{N_c^2}}\{J^2,T^8\},`$ (6)
$`M^{\mathrm{𝟐𝟕}}={\displaystyle \frac{1}{N_c}}\{T^8,T^8\}+{\displaystyle \frac{1}{N_c^2}}\{T^8,\{J^i,G^{i8}\}\},`$
$`M^{\mathrm{𝟔𝟒}}={\displaystyle \frac{1}{N_c^2}}\{T^8,\{T^8,T^8\}\},`$
where each operator in the $`1/N_c`$ expansion is understood to be accompanied by a coefficient, which has been suppressed in Eq. (3) for simplicity. The above operator expansion consists of eight independent operators corresponding to the eight $`I=0`$ baryon masses $`N`$, $`\mathrm{\Lambda }`$, $`\mathrm{\Sigma }`$, $`\mathrm{\Xi }`$, $`\mathrm{\Delta }`$, $`\mathrm{\Sigma }^{}`$, $`\mathrm{\Xi }^{}`$ and $`\mathrm{\Omega }`$. Each mass operator occurs at a definite order in $`1/N_c`$ and $`ϵ`$ and corresponds to one specific mass splitting, so the combined expansion in $`1/N_c`$ and $`ϵ`$ predicts a hierarchy of $`I=0`$ masses . A similar analysis can be performed for the $`I=1`$ and $`I=2`$ mass splittings of the octet and decuplet baryons. The analysis introduces two new parameters $`ϵ^{}`$ and $`ϵ^{\prime \prime }`$ which are $`I=1`$ and $`I=2`$ flavor symmetry breaking parameters, respectively. The two parameters are expected to be similar in magnitude.

Figure 1 plots sixteen independent mass splittings of the spin-$`1/2`$ octet and the spin-$`3/2`$ decuplet. The plot shows, relative to the leading $`𝒪(N_c)`$ baryon mass, $`I=0`$ mass splittings of order $`\frac{1}{N_c^2}`$, $`\frac{ϵ}{N_c}`$, $`\frac{ϵ}{N_c^2}`$, $`\frac{ϵ}{N_c^3}`$, $`\frac{ϵ^2}{N_c^2}`$, $`\frac{ϵ^2}{N_c^3}`$ and $`\frac{ϵ^3}{N_c^3}`$; $`I=1`$ mass splittings of order $`\frac{ϵ^{}}{N_c}`$, $`\frac{ϵ^{}}{N_c}`$, $`\frac{ϵ^{}}{N_c^3}`$ or $`\frac{ϵ^{}ϵ}{N_c^2}`$, $`\frac{ϵ^{}ϵ}{N_c^2}`$, $`\frac{ϵ^{}ϵ}{N_c^2}`$, and $`\frac{ϵ^{}ϵ^2}{N_c^3}`$; and $`I=2`$ mass splittings of order $`\frac{ϵ^{\prime \prime }}{N_c^2}`$, $`\frac{ϵ^{\prime \prime }}{N_c^2}`$, and $`\frac{ϵ^{\prime \prime }ϵ}{N_c^3}`$, respectively.

The first seven points in Figure 1, the $`I=0`$ mass splittings, are in striking agreement with the hierarchy predicted by the $`1/N_c`$ and $`ϵ`$ expansions. The singlet mass splitting of relative order $`1/N_c^2`$ is comparable in magnitude to the octet mass splitting of relative order $`ϵ/N_c`$. The three flavor-octet mass splittings proportional to one power of $`ϵ`$ are in accord with the predicted $`1/N_c`$ hierarchy of $`1/N_c`$, $`1/N_c^2`$, and $`1/N_c^3`$. If one did not know about the $`1/N_c`$ expansion and arbitrarily looked at three different flavor-octet mass

splittings, they would generically all contain some portion of the $`1/N_c`$ mass combination and therefore be of this relative magnitude. Only when examining the linear combinations picked out by the $`1/N_c`$ expansion does the $`1/N_c`$ hierarchy become apparent. The two flavor-$`\mathrm{𝟐𝟕}`$ mass splittings are proportional to two powers of $`ϵ`$, and occur at relative orders $`1/N_c^2`$ and $`1/N_c^3`$ in the $`1/N_c`$ expansion. These two flavor-$`\mathrm{𝟐𝟕}`$ mass splittings of the $`1/N_c`$ expansion are linear combinations of the Gell-Mann–Okubo flavor-$`\mathrm{𝟐𝟕}`$ mass splitting of the spin-$`1/2`$ baryon octet

$$\frac{1}{4}\left(2N\mathrm{\Sigma }3\mathrm{\Lambda }+2\mathrm{\Xi }\right)$$
(7)
and the flavor-$`\mathrm{𝟐𝟕}`$ Equal Spacing Rule mass splitting of the spin-$`3/2`$ baryon decuplet

$$\frac{1}{7}\left(4\mathrm{\Delta }5\mathrm{\Sigma }^{}2\mathrm{\Xi }^{}+3\mathrm{\Omega }\right).$$
(8)
Since these two famous mass splittings are linear combinations of the two $`1/N_c`$ mass splittings, they are both suppressed by a factor of $`1/N_c^2`$ in addition to the flavor suppression factor $`ϵ^2`$. Thus, the $`1/N_c`$ expansion predicts that these flavor-$`\mathrm{𝟐𝟕}`$ mass splittings are about an order of magnitude smaller than expected from an analysis of flavor-symmetry breaking alone. The single flavor-$`\mathrm{𝟔𝟒}`$ mass splitting

$$\frac{1}{4}\left(\mathrm{\Delta }3\mathrm{\Sigma }^{}+3\mathrm{\Xi }^{}\mathrm{\Omega }\right)$$
(9)
is third order in $`SU(3)`$ flavor symmetry breaking $`ϵ^3`$ and relative order $`1/N_c^3`$ in the $`1/N_c`$ expansion. The experimental data clearly show that the mass splitting is suppressed by a greater factor than expected from flavor-symmetry breaking alone, and the observed suppression is consistent with the $`1/N_c`$ hierarchy. A better measurement of this splitting would reduce the experimental error bar and test the $`1/N_c^3`$ prediction of the $`1/N_c`$ expansion.

There also is clear evidence for the $`1/N_c`$ hierarchy in the $`I=1`$ mass splittings . A new, very precise, measurement of $`\mathrm{\Xi }^0=1314.82\pm 0.06\pm 0.2\mathrm{MeV}`$  significantly reduces the errors on the $`I=1`$ mass combinations plotted in Fig. 1. For the first time, the Coleman-Glashow mass splitting

$$\left[\left(pn\right)\left(\mathrm{\Sigma }^+\mathrm{\Sigma }^{}\right)+\left(\mathrm{\Xi }^0\mathrm{\Xi }^{}\right)\right]$$
(10)
is measured to be non-zero, though only at the $`1\sigma `$ level. The Coleman-Glashow mass combination, the eleventh point in Fig. 1, is predicted to be of relative order $`\frac{ϵϵ^{}}{N_c^2}`$, or suppressed by a factor of $`\frac{ϵ}{N_c}`$ relative to the leading $`I=1`$ mass splittings, the eighth and ninth points, which are of relative order $`\frac{ϵ^{}}{N_c}`$. A similar suppression factor is predicted and observed for the $`I=1`$ mass splitting that is plotted as the tenth point in Fig. 1. The leading $`I=2`$ mass splitting, point fourteen, is also in line with the prediction of the $`1/N_c`$ expansion. Remaining points have large error bars. Improved measurements of a number of isospin mass splittings, particularly of the $`\mathrm{\Delta }`$ baryons, would reduce these errors and provide a further test of the $`1/N_c`$ hierarchy.

## 4 HEAVY QUARK BARYONS

Large-$`N_c`$ baryons containing a single heavy quark respect a large-$`N_c`$ $`SU(6)_{\mathrm{}}\times SU(4)_h`$ spin-flavor symmetry  with light-quark generators $`J_{\mathrm{}}^i`$, $`T_{\mathrm{}}^a`$ and $`G_{\mathrm{}}^{ia}`$ and heavy-quark generators

$`J_h^i=Q^{}{\displaystyle \frac{\sigma ^i}{2}}Q=J_c^i+J_b^i,`$
$`I_h^a=Q^{}{\displaystyle \frac{\tau ^a}{2}}Q,`$ (11)
$`G_h^{ia}=Q^{}{\displaystyle \frac{\sigma ^i\tau ^a}{4}}Q.`$
Light-quark spin-flavor symmetry is valid in the $`N_c\mathrm{}`$ limit, whereas the heavy-quark spin-flavor symmetry is valid in the large-$`N_c`$ limit and the heavy-quark limit of HQET.

The $`1/N_c`$ operator expansion for a baryon with $`Qqq`$ flavor quantum numbers involves operator products which are at most $`1`$-body in the heavy-quark generators and at most $`2`$-body in the light-quark generators. Light-quark generators are accompanied by a factor $`1/N_c`$, and heavy-quark generators are accompanied by a factor of $`1/N_c`$ and $`\mathrm{\Lambda }/m_Q`$ since the heavy-quark spin-flavor generators separately violate large-$`N_c`$ heavy-quark spin-flavor symmetry and heavy-quark spin-flavor symmetry of HQET. Heavy-quark flavor symmetry breaking for the $`Q=c`$ and $`Q=b`$ baryons enters through the generators

$`I_h^3={\displaystyle \frac{1}{2}}\left(N_{\mathrm{charm}}N_{\mathrm{bottom}}\right),`$
$`G_h^{i3}={\displaystyle \frac{1}{2}}\left(J_c^iJ_b^i\right),`$
whereas light-quark $`SU(3)`$ flavor symmetry breaking enters through light-quark generators which transform as the eighth or third component of an $`SU(3)`$ octet. The $`1/N_c`$ mass hierarchy can be analyzed separately for charm- and bottom-quark baryons using only heavy-quark spin symmetry, or together using heavy-quark spin-flavor symmetry. The joint expansion leads to relations between $`Q=b`$ and $`Q=c`$ mass splittings such as

$$\frac{1}{3}\left(\mathrm{\Sigma }_b+2\mathrm{\Sigma }_b^{}\right)\mathrm{\Lambda }_b=\frac{1}{3}\left(\mathrm{\Sigma }_c+2\mathrm{\Sigma }_c^{}\right)\mathrm{\Lambda }_c.$$
(12)
Heavy quark baryons and baryons containing no heavy quarks also can be analyzed in a joint expansion which involves the heavy-quark number operator

$$N_h=Q^{}Q=N_{\mathrm{charm}}+N_{\mathrm{bottom}}.$$
(13)
This analysis leads to relations between heavy-quark baryon mass splittings and mass splittings of the spin-$`1/2`$ octet and the spin-$`3/2`$ decuplet, such as

$$\frac{1}{3}\left(\mathrm{\Sigma }_Q+2\mathrm{\Sigma }_Q^{}\right)\mathrm{\Lambda }_Q=\frac{2}{3}\left(\mathrm{\Delta }N\right).$$
(14)
The $`1/N_c`$ expansion for heavy-quark baryon masses successfully predicted the $`\mathrm{\Xi }_c^{}`$ mass: the predicted value of $`\mathrm{\Xi }_c^{}=2580.8\pm 2.1\mathrm{MeV}`$  is in excellent agreement with the subsequent measured value of $`\mathrm{\Xi }_c^{}=2576.5\pm 2.3\mathrm{MeV}`$. With this measurement, all of the lowest-lying charm baryon masses are measured except for $`\mathrm{\Omega }_c^{}`$. The $`\mathrm{\Omega }_c^{}`$ mass can be determined in terms of the other measured charm baryon masses using the mass relation

$$\frac{1}{4}\left[\left(\mathrm{\Sigma }_c^{}\mathrm{\Sigma }_c\right)2\left(\mathrm{\Xi }_c^{}\mathrm{\Xi }_c^{}\right)+\left(\mathrm{\Omega }_c^{}\mathrm{\Omega }_c\right)\right],$$
(15)
which is predicted to hold at the sub-MeV level using the $`1/N_c`$ expansion, and therefore is extremely accurate and can be taken to be exact for most purposes. Presently, the extraction of the $`\mathrm{\Omega }_c^{}`$ mass using this mass relation is dominated by experimental uncertainties in the $`\mathrm{\Omega }_c`$ mass, whose PDG value of $`\mathrm{\Omega }_c=2704\pm 4\mathrm{MeV}`$ is $`2.5\sigma `$ away from the new CLEO measurement of $`\mathrm{\Omega }_c=2694.6\pm 3.5\mathrm{MeV}`$. These two experimental values yield central values for the $`\mathrm{\Omega }_c^{}`$ of 2768 and 2777 MeV, respectively. Finally, all of the $`Q=b`$ baryons (given the measured $`\mathrm{\Lambda }_b`$ mass) can be predicted from the charm baryon masses to a precision of about $`10`$ MeV, which is set by the accuracy of the worst $`Q=c`$ and $`Q=b`$ mass relation. Observation of some of the $`Q=b`$ baryons will allow a determination of the remaining unobserved $`Q=b`$ baryon masses to greater precision.

## 5 CONCLUSIONS

The $`1/N_c`$ expansion of QCD has provided an understanding of the mass spectrum of the spin-$`1/2`$ octet and spin-$`3/2`$ decuplet baryons. The mass hierarchy of baryons is well-described by an expansion in $`1/N_c`$ about the large-$`N_c`$ spin-flavor limit and an expansion in $`SU(3)`$ flavor symmetry breaking about the $`SU(3)`$ flavor symmetry limit. $`1/N_c`$ suppression factors are necessary to understand the relative magnitudes of the various baryon mass splittings of the baryon octet and decuplet. In the case of baryons containing a single heavy quark in HQET, there is an additional expansion parameter $`\mathrm{\Lambda }/m_Q`$ which describes corrections to the heavy-quark spin-flavor symmetry limit. The hierarchy of heavy-quark baryon masses predicted from the violation of large-$`N_c`$ light- and heavy-quark spin-flavor symmetries is sufficiently restrictive that it can be used to predict heavy-quark baryon masses at a level which is interesting experimentally. The combined $`1/N_c`$ and flavor-symmetry breaking expansions successfully predicted the $`\mathrm{\Xi }_c^{}`$ mass to an accuracy of several MeV. Additional predictions for charm and bottom baryon masses can and will be tested in the future.


format: text

text:
# SOME EXACT SOLUTIONS OF THE DIRAC EQUATION

## Acknowledgments

One of the authors (A.S.C.) would like thank FAPESP for financial support.

## References


text:
# ”Cold Melting” of Invar Alloys

## I Introduction

Among various approaches to the theory of crystal melting (see review ) the dislocation model (see  and references therein) occupies the leading position, describing the relevant physical phenomena in the most adequate way. According to the simplest version of the model , one should consider the change of the crystal free energy $`\mathrm{\Delta }F=WT\mathrm{\Delta }S`$, caused by a single dislocation, calculated per its unit length. Here $`W`$ is the internal energy change, and $`\mathrm{\Delta }S`$ is the corresponding entropy change. At elevated temperatures the free energy $`\mathrm{\Delta }F`$ may become negative, so that a spontaneous multiple generation of dislocations becomes thermodynamically favorable. This generation destroys the long range order of the crystal and leads to its amorphization and melting. A reasonable assumption is that both the energy $`W`$ and the entropy $`\mathrm{\Delta }S`$ are related only to the elastic strains induced by the dislocation and, hence, they hardly depend on temperature. It means that the melting temperature can be determined by the equation

$$T_m^{hot}=\beta W_{el}a.$$
(1)
with $`a`$ being the interatomic spacing. Here $`\beta `$ is a proportionality coefficient, which is $`\beta =1/(a\mathrm{\Delta }S)`$ in the simplest case . Generally, the coefficient $`\beta `$ may have a more complicated shape, since it should additionally incorporates such effects, as interactions between the dislocations, entropy due to various dislocation configurations and so on (see discussion in  and in the most recent paper ). The notation $`W_{el}`$ in equation (1) emphasizes connection of the dislocation self energy to the elastic strains.

When considering ferromagnetic crystals, whose internal energy incorporates also the exchange interaction energy, one may think about an additional contribution due to the variation of this exchange interaction in the field of the elastic strains around the dislocations. In typical ferromagnetics, say, Fe or Ni, the exchange interaction contribution makes at best several tenths of one percent of the elastic energy and can be disregarded. However, a very strong volume magnetostriction typical of the Invar alloys, as demonstrated, e.g., in , may result in an enormous enhancement of the dislocation exchange self energy and become comparable with the elastic one. Our recent paper  discusses the influence of this enhancement on the interaction between dislocations and solute atoms in Invar alloys with the aim to explain some features of plastic deformation of Invars.

The exchange self energy of a dislocation is negative and its absolute value grows in the ferromagnetic phase with decreasing temperature. It will be demonstrated in this paper that at low enough temperatures ($`T<T_C`$, $`T_C`$ is the Curie temperature) the total dislocation energy, which is now the sum of the elastic and exchange contributions, decreases significantly with decreasing temperature, and in some Invar alloys may even change its sign. It means that the criterion (1) can be met not only at elevated temperatures but also at a low temperature. It may lead to a situation when spontaneous multiple generation of dislocations and crystal amorphization may become thermodynamically favorable below a certain critical temperature.

This phenomenon may be called ”cold melting”. We plan to discuss here the conditions of the cold melting and find the temperature below which it becomes possible. We have no intentions to address here in detail the problem of what is the state of the crystal below the transition temperature and what the order of the transition is. We may only mention that we do not currently see reasons why the analysis of, say, Edwards and Warner  is not applicable in our case as well. They have demonstrated for the usual hot melting that the dislocation mechanism results in a first order phase transition. We expect that a similar first order transition may lead to an amorphized glassy-like state at low temperatures.

In order to find the temperature of the cold melting we need to calculate the dislocation self energy. First, the elastic energy of a unit length of an edge dislocation is determined by the equation

$$W_{el}=\frac{\mu b^2}{4\pi (1\nu )}ln\frac{R}{r_0}.$$
(2)
Here $`b`$ is the magnitude of the dislocation Burgers vector, $`\mu `$ is the shear modulus, $`\nu `$ is the Poisson coefficient. When calculating the energy (2), the integral of the dislocation induced strain field is cut both at small distances $`r_0`$, of the order of the interatomic spacing, and at large distances $`R`$, of the order of the average distance between the dislocations.

The volume magnetostriction results in a variation of the exchange energy of a ferromagnetic in the strain field in the vicinity of a dislocation. In order to calculate the corresponding energy change we consider the density of the exchange energy of the ferromagnetic in the molecular field approximation (see, e.g.,)

$$𝒲_{ex}=\frac{\omega M^2}{2}$$
(3)
where $`M`$ is the local magnetization, and

$$\omega =\frac{3k_BT_C}{np_{eff}^2\mu _B^2}$$
(4)
is the constant of the molecular field. Here $`n`$ is the atomic density, $`p_{eff}`$ is the effective number of the Bohr magnetons $`\mu _B`$ per atom. In principle, one should also consider a term proportional to $`(M)^2`$ in the energy (3). Its contribution to the exchange energy is, however, two orders of magnitude smaller than the leading term in (3) and, hence, it has been neglected.

The hydrostatic pressure $`p`$ created by an edge dislocation causes a local change of the magnetization

$$M=\overline{M}(1+\alpha p)$$
(5)
where $`\overline{M}`$ is the spontaneous magnetization of the ferromagnetic in the absence of the dislocation, $`\alpha `$ is a proportionality coefficient known empirically. It is usually rather small but may take anomalously high values in Invar alloys .

The hydrostatic pressure in the vicinity of an edge dislocation is

$$p(\rho ,\theta )=\frac{\mu b}{3\pi }\frac{1+\nu }{1\nu }\frac{\mathrm{sin}\theta }{\rho }$$
(6)
with $`\rho `$ and $`\theta `$ being the cylindrical coordinates. Substituting (5) and (6) into the energy (3) and integrating, using the same cuts as when calculating the elastic energy, one finds the exchange self energy of a dislocation per its unit length

$$W_{ex}=\frac{\omega \overline{M}^2\alpha ^2b^2\mu ^2}{18\pi }\left(\frac{1+\nu }{1\nu }\right)^2\mathrm{ln}\frac{R}{r_0}.$$
(7)
Now adding the exchange energy (7) to the elastic one (2), the total dislocation self energy per its unit length takes the form

$$W=W_{el}+W_{ex}=f(T)W_{el}$$
(8)
where

$$f(T)=1\frac{\omega \overline{M}^2(T)\alpha ^2(T)E(T)}{9}\frac{1+\nu (T)}{1\nu (T)}.$$
(9)
$`E=2\mu (1+\nu )`$ is the Young modulus. It is emphasized in equation (9) that material parameters in Invars are temperature dependent.

A similar analysis can be carried out with respect to the interaction between the dislocations. It leads to the conclusion that equation (8) is, in fact, more general. A similar equation with the same factor (9) holds also for the total internal energy of the system of dislocations which includes both the internal self energies of the individual dislocation and the interactions between them.

The second term in equation (9) reflects the contribution of the dislocation exchange energy. It is usually very small, about 10-3, in conventional ferromagnetics and, hence, may be neglected. In Invars, however, the situation changes dramatically. It is connected with the fact that the coefficient $`\alpha `$, connected to the volume magnetostriction, may be several tens times larger than, say, in Fe or Ni . In the ferromagnetic phase this second term grows with lowering temperature and may become comparable to one, so that the function $`f(T)`$ may become very small or even change its sign.

Now we discuss how the exchange contribution to the dislocation energy influences the melting criterion (1) for the temperature at which the spontaneous generation of dislocations becomes possible. Direct analysis shows that now we should take the same condition (1) but substitute there the total self energy of the dislocation (8) instead of the elastic one. As for the entropy due the dislocation induced changes in the magnetic system, it can be still neglected. Then the melting temperature can be found as a solution of the equation

$$T_m=\beta f(T_m)W_{el}a.$$
(10)
This equation may have more than one solutions. One solution can be found at high temperatures and corresponds to the conventional hot melting temperature $`T_m^{(hot)}`$. Usually the melting temperature exceeds the Curie temperature where the local magnetization disappears, $`\overline{M}=0`$, hence, $`f(T)=1`$. Then equation (10) coincides with (1) and leads to the standard description of the hot melting. As for the cold melting temperature, it is connected with hot melting temperature by the condition

$$T_m^{(cold)}=T_m^{(hot)}f(T_m^{(cold)})$$
(11)
Equation (11) can be solved when using the low temperature values of the material parameters and Curie-Weiss equation

$$\overline{M}(T)=M_0\sqrt{1\frac{T}{T_C}}$$
(12)
for the spontaneous magnetization. Here $`M_0`$ is the spontaneous magnetization at zero temperature. Then

$$T_m^{(cold)}=T_C\frac{\gamma 1}{\gamma {\displaystyle \frac{T_C}{T_m^{(hot)}}}}$$
(13)
where

$$\gamma =\frac{\omega M_0^2\alpha ^2E}{9}\frac{1+\nu }{1\nu }$$
(14)
In the known ferromagnetics $`T_C<T_m^{(hot)}`$, hence, the temperature of the cold melting may be positive only if $`\gamma >1`$. Although rather rough approximations have been used (see discussion of Invar properties in ), when deriving equation (13), the condition $`\gamma >1`$ is in fact more general and reflects just the fact that the energy necessary for creation of a dislocation becomes negative. This condition can be used as a good indicator when looking for Invar alloys in which cold melting may be observed.

It can be worth considering Fe-Pt Invar alloys, with a content close to Fe0.72Pt0.28, which are characterized by high values of $`\alpha `$. The volume magnetostriction becomes especially large in Fe0.72Pt0.28 alloys with disordered distribution of Fe and Pt atoms over the lattice sites. It may reach the value $`\alpha =2.4\times 10^{11}`$(dyn/cm2)-1 at room temperatures . As for a detailed information on the temperature dependence of $`\alpha `$ in these alloys, unfortunately, it is not available. Nevertheless, we shall use this value in the estimates to be done below.

This alloy is characterized by the following parameters: $`T_C=371`$K, $`p_{eff}=2.13`$ , and $`n=7.6\times 10^{22}`$cm3. Then equation (4) leads to $`\omega =5190`$. At $`M_0=np_{eff}\mu _B=1.5`$kG, $`\nu =0.3`$, and $`E=1.2\times 10^{12}`$ dyn/cm2 one gets $`\gamma =1.66`$. Therefore, we see that the condition $`\gamma >1`$ is holds very well.

Then knowing the temperature of the hot melting, $`T_m^{(hot)}=1812`$K, for this alloy one could have found the temperature $`T_m^{(cold)}=168`$K for the cold melting which follows from equation (13). However, we prefer doing the same calculation by using the experimental data on the temperature dependence of the magnetization $`\overline{M}(T)/M_0`$ , and elastic constants  temperature. The resulting temperature is even higher $`T_m^{(cold)}=298`$K and is rather close to the room temperature. This also provides a justification for using the room temperature value of the parameter $`\alpha `$ in this estimate.

Therefore, the prediction of our model is that at temperatures below $`T_m^{(cold)}=298`$K a spontaneous multiple generation of dislocations Fe-Pt Invar alloys in may become possible. It is interesting to note that an intensive generation of dislocations is really observed in disordered Fe-Pt Invar alloys in the temperature range from 300 to 77K in several experiments (see, e.g.,  and references therein). Our estimate for the temperature of the cold melting lies close to the upper limit of this range. However, the experimental observed phenomenon is rather complicated. This generation goes hand in hand with the martensite transitions from FCC to BCC phase, also observed in the same temperature range, and the two effects can be hardly separated one from another.

There are reasons to believe that they are really tightly connected. On one hand, dislocations may be generated as a result of accumulation and discharge of elastic strains in the vicinity of the boundaries between the martensite and austenite (host) phases in the process of the martensite growth . A multiple generation of dislocations is observed only in disordered alloys with larger values of the parameter $`\alpha `$ which favors its connection with the mechanism of the cold melting discussed above. As for ordered alloys with a smaller value of $`\alpha `$ a dislocation generation is not observed in them.

On the other hand, we believe that the martensite transitions, which are much more intensive in disordered Fe-Pt alloys , may be themselves induced by the dislocation generation in the process of the cold melting. It is known that appearance of strain-induced martensites is strongly facilitated by defects consisting of properly arranged dislocations which play the role of active centers in the formation of the martensite . We assume that the dislocation generation leads to creation of such centers, and these centers induce martensite transitions.

A detailed experimental and theoretical study of the properties of Invar alloys at temperatures close and below the cold melting temperature present a special interest. One may expect highly unusual plastic properties of Invars in this region. As an example of what one may expect, we mention that at $`T<T_m^{(cold)}`$ dislocations of the same sign start attracting rather than repelling each other, meaning that the theory of the deformation hardening should be completely revisited.


text:
# 1 Introduction

## 1 Introduction

It is not yet possible to simulate the formation of proteins’ native structures on the computer in a controlled way. This goal has been achieved in the context of simple lattice and off-lattice models, where typically each amino acid is represented by a single interaction site corresponding to the Cα atom, and such studies have provided valuable insights into the physical principles of protein folding  and the statistical properties of functional protein sequences . However, these models have their obvious limitations. Therefore, the search for computationally feasible models with a more realistic chain geometry remains a highly relevant task.

In this paper, we discuss a model based on the well-known fact that the main degrees of freedom of the protein backbone are the Ramachandran torsional angles $`\varphi _i,\psi _i`$ . Each amino acid is represented by five or six atoms, which makes this model computationally slightly more demanding than Cα models. On the other hand, it also makes interactions such as hydrogen bonds easier to define. The formation of native structure is, in this model, driven by hydrogen-bond formation and effective hydrophobicity forces; hydrophobicity is widely held as the most important stability factor in proteins , and hydrogen bonds are essential to properly model the formation of secondary structure.

In this model, we study in particular a three-helix-bundle protein with 54 amino acids, which represents a truncated and simplified version of the four-helix-bundle protein de novo designed by Regan and DeGrado . This example was chosen partly because there have been earlier studies of similar-sized helical proteins using models at comparable levels of resolution . The behavior of small fast-folding proteins is a current topic in both theoretical and experimental research, and a three-helix-bundle protein that has been extensively studied both experimentally  and theoretically  is fragment B of staphylococcal protein A.

In addition to the three-helix protein, to study size dependence, we also look at the behavior of the corresponding one- and two-helix segments. By using the method of simulated tempering , a careful study of the thermodynamic properties of these different chains is performed.

Not unexpectedly, it turns out that the behavior of the model depends strongly on the relative strength of the hydrogen-bond and hydrophobicity terms. In fact, the situation is somewhat reminiscent of what has been found for homopolymers with stiffness , with hydrogen bonds playing the role of the stiffness term. Throughout this paper, we focus on one specific empirical choice of these parameters.

For this choice of parameters, we find that the three-helix-bundle protein has the following three properties. First, it does form a stable three-helix bundle (except for a 2-fold topological degeneracy). Second, its folding transition is abrupt, from an expanded state to the native three-helix-bundle state. Third, compared to the one- and two-helix segments, it forms a more stable secondary structure. It should be stressed that these properties are found without resorting to the popular Gō approximation , in which interactions that do not favor the desired structure are ignored.

## 2 The Model

The model we study is a reduced off-lattice model. The chain representation is illustrated in Fig. 1. As mentioned in the introduction, each amino acid is represented by five or six atoms. The three backbone atoms N, Cα and C are all included. Also included are the H and O atoms shown in Fig. 1, which we use to define hydrogen bonds. Finally, the side chain is represented by a single atom, Cβ, which can be hydrophobic, polar or absent. This gives us the following three types of amino acids: A with hydrophobic Cβ, B with polar Cβ, and G (glycine) without Cβ.

The H, O and Cβ atoms are all attached to the backbone in a rigid way. Furthermore, in the backbone, all bond lengths, bond angles and peptide torsional angles ($`180^{}`$) are held fixed. This leaves us with two degrees of freedom per amino acid, the Ramachandran torsional angles $`\varphi _i`$ and $`\psi _i`$ (see Fig. 1). The parameters held fixed can be found in Table 1.

Our energy function

$$E=E_{\text{loc}}+E_{\text{sa}}+E_{\text{hb}}+E_{\text{AA}}$$
(1)
is composed of four terms. The local potential $`E_{\text{loc}}`$ has a standard form with 3-fold symmetry,

$$E_{\text{loc}}=\frac{ϵ_\varphi }{2}\underset{i}{}(1+\mathrm{cos}3\varphi _i)+\frac{ϵ_\psi }{2}\underset{i}{}(1+\mathrm{cos}3\psi _i).$$
(2)
The self-avoidance term $`E_{\text{sa}}`$ is given by a hard-sphere potential of the form

$$E_{\text{sa}}=ϵ_{\text{sa}}\stackrel{}{}_{i<j}\left(\frac{\sigma _{ij}}{r_{ij}}\right)^{12},$$
(3)
where the sum runs over all possible atom pairs except those consisting of two hydrophobic Cβ. The hydrogen-bond term $`E_{\text{hb}}`$ is given by

$$E_{\text{hb}}=ϵ_{\text{hb}}\underset{ij}{}u(r_{ij})v(\alpha _{ij},\beta _{ij}),$$
(4)
where

$`u(r_{ij})`$ $`=`$ $`5\left({\displaystyle \frac{\sigma _{\text{hb}}}{r_{ij}}}\right)^{12}6\left({\displaystyle \frac{\sigma _{\text{hb}}}{r_{ij}}}\right)^{10}`$ (5)
$`v(\alpha _{ij},\beta _{ij})`$ $`=`$ $`\{\begin{array}{cc}\mathrm{cos}^2\alpha _{ij}\mathrm{cos}^2\beta _{ij}\hfill & \alpha _{ij},\beta _{ij}>90^{}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}`$ (8)
In Eq. 4 $`i`$ and $`j`$ represent H and O atoms, respectively, and $`r_{ij}`$ denotes the HO distance, $`\alpha _{ij}`$ the NHO angle, and $`\beta _{ij}`$ the HOC angle. Any HO pair can form a hydrogen bond. The last term in Eq. 1, the hydrophobicity term $`E_{\text{AA}}`$, has the form

$$E_{\text{AA}}=ϵ_{\text{AA}}\underset{i<j}{}\left[\left(\frac{\sigma _{\text{AA}}}{r_{ij}}\right)^{12}2\left(\frac{\sigma _{\text{AA}}}{r_{ij}}\right)^6\right],$$
(9)
where both $`i`$ and $`j`$ represent hydrophobic Cβ. To speed up the simulations, a cutoff radius $`r_c`$ is used,The cutoff procedure is $`f(r)\stackrel{~}{f}(r)`$ where $`\stackrel{~}{f}(r)=f(r)f(r_c)(rr_c)f^{}(r_c)`$ if $`r<r_c`$ and $`\stackrel{~}{f}(r)=0`$ otherwise. which is 4.5Å for $`E_{\text{sa}}`$ and $`E_{\text{hb}}`$, and 8Å for $`E_{\text{AA}}`$.

In this energy function, roughly speaking, the first two terms, $`E_{\text{loc}}`$ and $`E_{\text{sa}}`$, enforce steric constraints, whereas the last two terms, $`E_{\text{hb}}`$ and $`E_{\text{AA}}`$, are the ones responsible for stability. Force fields similar in spirit, emphasizing hydrogen bonding and hydrophobicity, have been used with some success to predict structures of peptides  and small helical proteins .

The parameters of our energy function were determined largely by trial and error. The final parameters are listed in Table 2. The parameters $`\sigma _{ij}`$ of Eq. 3 are given by

$$\sigma _{ij}=\sigma _i+\sigma _j+\mathrm{\Delta }\sigma _{ij},$$

where $`\sigma _i,\sigma _j`$ can be found in Table 2, and $`\mathrm{\Delta }\sigma _{ij}`$ is zero except for CβC, CβN and CβO pairs that are connected by three covalent bonds. In these three cases, we put $`\mathrm{\Delta }\sigma _{ij}=0.625`$Å. This could equivalently be described as a change of the local $`\varphi _i`$ and $`\psi _i`$ potentials. In Fig. 2, we show $`\varphi _i,\psi _i`$ scatter plots for nonglycine (A and B) and glycine for our final parameters, which are in good qualitative agreement with the $`\varphi _i,\psi _i`$ distributions of real proteins .

Finally, we determined the strengths of the hydrogen-bond and hydrophobicity terms on the basis of the resulting overall thermodynamic behavior of the three-helix sequence. For this purpose, we performed a set of trial runs for fixed values of the other parameters. An alternative would have been to use the method of Shea et al. . The result of our empirical determination of $`ϵ_{\text{hb}}`$ and $`ϵ_{\text{AA}}`$ does not seem unreasonable; at the folding temperature of the three-helix sequence (see below), we get $`ϵ_{\text{hb}}/kT4.3`$ and $`ϵ_{\text{AA}}/kT3.4`$.

In this model, we study the three sequences shown in Table 3, which contain 16, 35 and 54 amino acids, respectively. Following the strategy of Regan and DeGrado , the A and B amino acids are distributed along the sequence 1H in such a way that this segment can form a helix with all hydrophobic amino acids on the same side. The sequence 3H, consisting of three such stretches of As and Bs plus two GGG segments, is meant to form a three-helix bundle. This particular sequence was recently studied by Takada et al. , who used a more elaborate model with nonadditive forces.

## 3 Results

To study the thermodynamic behavior of the chains described in the previous section, we use the method of simulated tempering. This means that we first select a set of allowed temperatures and then perform simulations in which the temperature is a dynamical variable. This is done to speed up low-temperature simulations. In addition, it provides a convenient method for calculating free energies.

An example of a simulated-tempering run is given in Fig. 3, which shows the Monte Carlo evolution of the energy $`E`$ and radius of gyration $`R_\text{g}`$ (calculated over all backbone atoms) in a simulation of the three-helix sequence. Also shown, bottom panel, is how the system jumps between the different temperatures. Two distinct types of behavior can be seen. In one case, $`E`$ is high, fluctuations in size are large, and the temperatures visited are high. In the other case, $`E`$ is low, the size is small and almost frozen, and the temperatures visited are low. Interesting to note is that there is one temperature, the next-lowest one, which is visited in both cases. Apparently, both types of behavior are possible at this temperature.

In Fig. 4a we show the specific heat as a function of temperature for the one-, two- and three-helix sequences. A pronounced peak can be seen that gets stronger with increasing chain length. In fact, the increase in height is not inconsistent with a linear dependence on chain length, which is what one would have expected if it had been a conventional first-order phase transition with a latent heat.

Our results for the radius of gyration (not displayed) show that the specific heat maximum can be viewed as the collapse temperature. The specific heat maximum is also where hydrogen-bond formation occurs, as can be seen from Fig. 4b. Important to note in this figure is that the decrease in hydrogen-bond energy per amino acid with decreasing temperature is most rapid for the three-helix sequence, which implies that, compared to the shorter ones, this sequence forms more stable secondary structure. The results for the chain entropy shown in Fig. 4c provide further support for this; the entropy loss per amino acid with decreasing temperature is largest for the three-helix sequence.

It should be stressed that the character of the collapse transition depends strongly on the relative strength of the hydrogen-bond and hydrophobicity terms. Figure 4 shows that the transition is very abrupt or “first-order-like” for our choice $`(ϵ_{\text{hb}},ϵ_{\text{AA}})=(2.8,2.2)`$. A fairly small decrease of $`ϵ_{\text{hb}}/ϵ_{\text{AA}}`$ is sufficient to get a very different behavior with, for example, a much weaker peak in the specific heat. In this case, the chain collapses to a molten globule without specific structure rather than to a three-helix bundle. A substantially weakened transition was observed for $`ϵ_{\text{hb}}=ϵ_{\text{AA}}=2.5`$. If, on the other hand, $`ϵ_{\text{hb}}/ϵ_{\text{AA}}`$ is too large, then it is evident that the chain will form one long helix instead of a helical bundle.

We now turn to the three-dimensional structure of the three-helix sequence in the collapsed phase. It turns out that it does form a three-helix bundle. This bundle can have two distinct topologies: if we let the first two helices form a U, then the third helix can be either in front of or behind that U. The model is, not unexpectedly, unable to discriminate between these two possibilities. To characterize low-temperature conformations, we therefore determined two representative structures, one for each topology, which, following , are referred to as FU and BU, respectively. These structures are shown in Fig. 5. They were generated by quenching a large number of low-$`T`$ structures to zero temperature, and we feel convinced that they provide good approximations of the energy minima for the respective topologies. Given an arbitrary conformation, we then measure the root-mean-square distances $`\delta _i`$ ($`i=`$ FU, BU) to these two structures (calculated over all backbone atoms). These distances are converted into similarity parameters $`Q_i`$ by using

$$Q_i=\mathrm{exp}(\delta _i^2/100\text{Å}^2).$$
(10)
At temperatures above the specific heat maximum, both $`Q_i`$ tend to be small. At temperatures below this point, the system is found to spend most of its time close to one or the other of the representative structures; either $`Q_{\text{FU}}`$ or $`Q_{\text{BU}}`$ is close to 1. Finally, at the peak, all three of these regions in the $`Q_{\text{FU}},Q_{\text{BU}}`$ plane are populated, as can be seen from Fig. 6a. In particular, this implies that the folding transition coincides with the specific heat maximum.

The folding transition can be described in terms of a single “order parameter” by taking $`Q=\mathrm{max}(Q_{\text{FU}},Q_{\text{BU}})`$ as a measure of nativeness. Correspondingly, we put $`\delta =\mathrm{min}(\delta _{\text{FU}},\delta _{\text{BU}})`$. In Fig. 6b, we show the free-energy profile $`F(Q)`$ at the folding temperature. The free energy has a relatively sharp minimum at $`Q0.9`$, corresponding to $`\delta 3`$Å. This is followed by a weak barrier around $`Q=0.7`$, corresponding to $`\delta 6`$Å. Finally, there is a broad minimum at small $`Q`$, where $`Q=0.2`$ corresponds to $`\delta 13`$Å.

What does the nonnative population at the folding temperature correspond to in terms of $`R_\text{g}`$ and $`E_{\text{hb}}`$? This can be seen from the $`Q,R_\text{g}`$ and $`Q,E_{\text{hb}}`$ scatter plots in Fig. 7. These plots show that the low-$`Q`$ minimum of $`F(Q)`$ corresponds to expanded structures with a varying but not high secondary-structure content. Although a detailed kinetic study is beyond the scope of this paper, we furthermore note that the free-energy surfaces corresponding to the distributions in Fig. 7 are relatively smooth. Consistent with that, we found that standard fixed-temperature Monte Carlo simulations were able to reach the native state, starting from random coils.

Let us finally mention that we also performed simulations of some random sequences with the same length and composition as the three-helix sequence. The random sequences did not form stable structures and collapsed more slowly with decreasing temperature than the designed three-helix sequence.

## 4 Summary and Outlook

We have studied a reduced protein model where the formation of native structure is driven by a competition between hydrogen bonds and effective hydrophobicity forces. Using this force field, we find that the three-helix-bundle protein studied has the following properties:

* It does form a stable three-helix-bundle state, except for a 2-fold topological degeneracy.
* It undergoes an abrupt folding transition from an expanded state to the native state.
* It forms more stable secondary structure than the corresponding one- and two-helix segments.

An obvious question that remains to be addressed is what is needed to lift the topological degeneracy. Not obvious, however, is whether this question should be addressed at the present level of modeling, before including full side chains.

A first-order-like folding transition that takes the system directly from the unfolded state to the native one is what one expects for small fast-folding proteins. For the model to show this behavior, careful tuning of the relative strength of the hydrogen-bond and hydrophobicity terms, $`ϵ_{\text{hb}}/ϵ_{\text{AA}}`$, is required. This $`ϵ_{\text{hb}}/ϵ_{\text{AA}}`$ dependence may at first glance seem unwanted but is not physically unreasonable; $`ϵ_{\text{hb}}`$ can be thought of partly as a stiffness parameter, and chain stiffness has important implications for the phase structure, as shown by recent work on homopolymers . Note also that incorporation of full side chains makes the chains intrinsically stiffer, which might lead to a weaker $`ϵ_{\text{hb}}/ϵ_{\text{AA}}`$ dependence.

Our three-helix sequence has previously been studied by Takada et al. , who used a more elaborate force field. It was suggested that it is essential to use context-dependent hydrogen bonds for the three-helix-bundle protein to make more stable secondary structure than its one-helix fragments. Our model shows this behavior, although its hydrogen bonds are context-independent.

Let us finally stress that we find a first-order-like folding transition without using the Gō approximation. Evidence for first-order-like folding transitions has been found for proteins with similar lengths in some Cα models , but these studies use this approximation.

## Acknowledgements

This work was in part supported by the Swedish Foundation for Strategic Research.


text:
# Aggregation number distributions and mesoglobules in dilute solutions of diblock and triblock copolymers

## I Introduction

The phenomena of coil–to–globule transition for a single polymer molecule and phase separation in polymer solutions are foundations of the modern polymer science . In recent years there were numerous studies of the behaviour of water soluble polymers near the lower critical solubility temperature (LCST). Typical systems include poly-N-isopropylacrylamide (PNIPAM) homopolymer and block copolymers  of the poly(ethylene oxide)-poly(propylene oxide)-poly(ethylene oxide), or PEO-PPO systems briefly . These copolymers are widely used in pharmaceutical, agricultural and food industries due to their low toxicity and surfactant characteristics (EO is hydrophilic and PO is hydrophobic). Block and more complex random heteropolymers in solutions and melts at higher concentrations exhibit a variety of ordered micro–phase separated and disordered glassy phases . Block copolymers also play role of surfactants in ternary mixtures of two otherwise immiscible liquids such as water and oil, and these mixtures also produce micelles and lamellae .

There are two types of triblock copolymers having rather distinct properties. First, a PEO-PPO-PEO polymer with hydrophilic ends and a hydrophobic middle is the standard Pluronic , which forms stable micelles. Second, a polymer with the inverse structure would tend to network hydrophobic end ‘stickers’ by hydrophilic bridges. These are so–called telechelic associating polymers forming physical gels with nontrivial rheological properties at higher concentrations .

The aggregation number distribution during micellisation has been a matter of intensive studies for amphiphilic molecules . The classical situations, in which there is a clear phase separation transition from $`M_{cl}=1`$ to $`M_{cl}=M_{max}`$, or where one could introduce notions of CMC (critical micelle concentration), and possibly also CAC (critical aggregation concentration), are well known. Some cases of molecular self–assembly, in which relatively monodispersed distributions arise around CMC, have also been seen. It is a conventional wisdom that in solutions of amphiphilic polymers the notion of CMC does not really exist since single chain micelles already form in an infinitely dilute solution if the primary sequence of the polymer permits micro–phase separation within the collapsed globule.

The phase diagrams of diblock and in some cases triblock copolymers have been closely investigated . At low concentration these systems form micelles, but the question of their aggregation number is not so well understood as it appears to be quite sensitive on various parameters of the system, although it is generally believed that such a distribution should be fairly polydispersed. This situation is in a marked contrast with a two–phase separation onto the macro–aggregate and a gas of single globules in the simpler case of homopolymer solutions, which are well described by the mean–field Flory–Huggins type theories , and their extensions based on scaling arguments , self–consistent treatments  and the Lifshitz approach .

Ordinary micellisation and more complex aggregation number distributions in copolymer solutions are realised due to the complex interactions favouring micro–phase separation of hydrophobic and hydrophilic units. Since the CMC of these polymers is essentially zero as in the case of polysoaps , one would like to see how the aggregation number of the micelles would behave nearby, i.e. in a very dilute solution. Indeed, normally the character of such a distribution has a clearly distinct character around CMC vs other regimes.

In attempting to resolve this interesting question, unfortunately, both experiment and theory had considerable difficulties and little success. The aggregation phenomenon significantly obscures experimental observation of a single chain collapse by light scattering techniques, leading often to a diverse range of theoretical interpretations, particularly as it may be hard to separate purely equilibrium issues from the kinetic ones . Indeed, the rather skewed shape of the Flory–Huggins diagram requires extremely low concentrations to see single globules  thus making the scattered light intensity too low for a reliable observation.

In a number of recent experiments on PNIPAM polymers in aqueous solution, in which some of the units are either hydrophobically modified or replaced by ionomers , it has been observed that stable size monodispersed nanoparticles can be formed on heating above the LCST instead of chains simply aggregating. The mean size of these particles is thermodynamically controlled and it increases with increasing polymer concentration and increasing strength of inter–monomer attraction. Electron microscopy further confirms the existence of such spherical particles with mean size and distribution in agreement with the light scattering. We called such an unusual type of mesoscopic structures mesoglobules as they tend to have a compact spherical shape and a rather monodisperse size distribution. To be explicit then their precise definition is as follows. Mesoglobules are essentially equally sized globules obtained by association of more than one and less than all polymer chains in the system and corresponding to, at least, a local minimum of the free energy. Thus, the term mesoglobules refers to the globules possessing a well characterised size distribution, but not necessarily a well defined conformational structure such as e.g. in micelles.

In our previous paper  based on the Gaussian variational theory and lattice Monte Carlo simulation we have shown that in very dilute solutions the histograms of cluster distribution favour size monodispersed structures for a certain class of heteropolymer sequences and only in some narrow strips of the phase diagram. We also have seen that the mean size of the mesoglobules is related to a characteristic scale of the micro–phase separation. Thus, generally the larger this scale (e.g. block size) the larger are the mesoglobules produced.

Significantly, in the situation when the mesoglobules occur there is a quite narrow single peak distribution in the aggregation numbers and the mean squared radius of gyration, so that the populations of single globules or the macro–aggregate are absent. Moreover, we have demonstrated that the mesoglobules are truly thermodynamically stable, i. e. the clusters of exactly equal size possess the lowest free energy value. Clusters of slightly unequal size obtained from these ideal mesoglobules have a higher free energy value and the barriers separating the mesoglobular minimum from other local minima tend to grow quickly as more and more asymmetric divisions of the system into macromolecular clusters are considered. This situation is quite unusual and novel for spherical micelles and it occurs due to a rather delicate competition of energetic (micro–phase separation) and entropic (translational entropy) terms under the connectivity constraints. For the mesoglobules to form the degree of amphiphilicity (distinction in the interaction energies between ‘a’ and ‘b’ monomers with the solvent) and the mean strength of the monomers attractions should be in a certain relation. In Ref.  we have seen that this situation is possible for alternating monomers, repeating blocks (having smaller mesoglobules), as well as for a number of randomly generated sequences, but in all cases at rather low concentrations. Mesoglobules can form at infinitesimally small concentrations, corresponding to essentially zero CMC in polymeric surfactants. The mean size of the mesoglobules grows with increasing concentration, but when the latter reaches a certain value, which perhaps one may view as an analog of CAC, a polydisperse distribution of clusters with a large population of the macro–aggregate results, recovering a more conventional regime.

Our understanding of the diblock copolymers, however, was not complete in paper . We believed that a finite size effect was significant there due to the chosen chain lengths. Thus, in this Letter we shall revisit the diblocks and also consider the triblocks of both possible types.

## II Results

As in Ref.  we adopt the Metropolis technique in the lattice model of Ref. . This model, apart from the connectivity and excluded volume constraints, includes short–ranged pair–wise interactions between lattice sites. The system is completely characterised by three Flory interaction parameters, $`\chi _{aa}`$, $`\chi _{ab}`$ and $`\chi _{bb}`$, along with $`N`$, $`M`$ and linear lattice size $`L`$. In addition to local monomer moves , we include translational moves representing diffusion of chains. The latter moves are applied to all clusters of chains with a probability inversely proportional to the number of monomers within (Stokes law). We shall consider sequences consisting of strongly hydrophobic and slightly hydrophilic units. However, while in paper  the number of monomers in each chain were $`N=24`$ and the total number of chains was $`M=20`$, here we shall analyse the case of twice the number $`M=40`$ of half–long chains $`N=12`$ in the same simulation box of linear size $`L=60`$. Thus, the monomer concentrations $`cNM/L^3`$ are the same in both cases, and also we shall keep the numbers of the hydrophobic and hydrophilic monomers equal in each chain.

In Figs. 1 we present snapshots of two typical types of structures observed for triblock copolymers in the region of the phase diagram where mesoglobules are normally seen . Fig. 1a corresponds to triblocks with the hydrophilic (or P- for polar) ends and the hydrophobic (or H-) middle block, and we shall call these the P-ends triblocks briefly. These structures clearly have essentially equal size and thus could be considered as mesoglobules. Moreover, they possess a pronounced micellar character with a hydrophobic (black circles) core and a hydrophilic (white circles) corona. Fig. 1b corresponds to the inverted case of the H-ends triblocks. In this case typical structures often have two, three (as in this figure), or more large clusters inter–connected by the hydrophilic middle block bridges. Such dimers, trimers and multimers are indeed to be expected at very low concentrations as it is known that this type of associating polymers forms telechelic physical gels at higher concentrations .

To investigate the aggregation number and size distributions we have obtained a large ensemble of $`Q=1000`$ independent equilibrium states. In Figs. 2 we present the calculated histograms of the aggregation number (mass) and size probability densities for the diblock ($`a_6b_6`$) and triblock sequences (H-ends and P-ends). We have seen in Fig. 7b of Ref.  that the diblock sequence s4 had a large population of the macro–aggregate there. We have argued in Ref.  that this might be due to a finite–size effect since the total number of chains $`M=20`$ was comparable to the mean aggregation number of the mesoglobules. In Fig. 2a, however, one can see for the diblocks that their mean aggregation number in the mesoglobules is about $`12`$, whereas the total number of chains is $`M=40`$, so that there is no population of the macro–aggregate any more. Therefore, the diblock sequences do indeed produce mesoglobules in a certain region of the phase diagram in dilute solution and the corresponding size distribution in Fig. 2b is quite narrow.

Similarly to diblocks, the P-ends triblocks produce a clear single peak distribution for the aggregation number (see Fig. 2a) corresponding to the mesoglobules of a somewhat smaller size and the respective cluster size distribution is also rather narrow in Fig. 2b. On the contrary, the aggregation number distribution for the H-ends triblocks is rather broad and has a certain population of the macro–aggregate. The width of this distribution is explained mainly by formation of dimer and multimer clusters, so that one can no longer define a preferred cluster size. Fig. 2b shows that the size distribution of the H-ends has two peaks: a very broad multimeric peak (near the region $`6<R_{cl}<8`$) and a still reasonably narrow peak corresponding to the mesoglobules and their dimers (near $`R_{cl}4`$). It is interesting to note that the mean size of the latter objects for the H-ends triblocks is considerably smaller than for the diblock sequence of the same length, although the situation is reverse for the mean aggregation numbers. This could be easily understood by comparison of Figs. 6d in Ref.  and 1b. One could see that the P-corona in the former case is quite extended, whereas in the latter case it is more dense and compact, containing a considerable number of loops.

Finally, it is also instructive to compare the mesoglobules for diblock $`a_6b_6`$ with those of twice longer sequence $`(a_6b_6)_2`$ (see s3 in Figs. 7a,8a in Ref. ) formed by merging every two $`a_6b_6`$ chains at the same monomer concentration $`c`$. By a careful analysis of these data one can conclude that the mean cluster mass increases and the mean size decreases slightly for $`(a_6b_6)_2`$ as compared to $`a_6b_6`$. The effect of an increase in the mean cluster mass is rather weak and does not have a simple explanation, but the decrease of the mean size is clearly due to a more compact size of loops in mesoglobules of $`(a_6b_6)_2`$ than the size of free P-ends for $`a_6b_6`$, similarly to the discussion in the previous paragraph.

## III Conclusion

In this Letter we have addressed the question of sequence specificity for forming mesoglobules in dilute solutions. The aggregation number and size distributions for inter–molecular clusters in a poor solvent have been found to be monodisperse for diblocks and triblocks with outer hydrophilic ends, while rather polydisperse for triblocks with hydrophobic end ’stickers’. Comparative analysis of chains of varying length but at a fixed total monomer concentration has elucidated nontrivial dependences of the mesoglobules mass and size on the copolymer chain length and sequence and a special role of hydrophilic loops and bridges in determining the compactness of the resulting structures.

The main conclusion from these studies is that the mesoglobules formation thermodynamics is rather delicate and that such stable structures are only possible for certain heteropolymer sequences and in a narrow region of the phase diagram. Our results also suggest that multimeric clusters of a certain size may have a stabilisation mechanism similar to that of mesoglobules, and again only for certain heteropolymer sequences.

We believe that the observations of this Letter may shed some light on understanding the problems of competition aggregation vs folding in protein solutions and, possibly, on self–organisation of the quaternary structure in multimeric proteins. It appears that for the latter process to take place a considerable number of ‘sticky’ hydrophobic amino acid residues should be exposed on the exterior of each of folded subunits. The resulting association produces a well defined quaternary structure, which is biologically functional as in e.g. hemoglobin, rather than a disordered aggregate. Hopefully, work in this direction would also help unravel the mechanism of association of so–called Bence–Jones proteins believed to be related to amyloid fibril formation.

More importantly, detailed understanding of the mechanism of mesoglobules formation is important for being able to select polymers from which nanoparticles of required size and polydispersity can be prepared. The main advantage of this new approach is that nanoparticles can be formed reversibly in mild conditions from a prepurified polymer. It also seems quite reasonable that even a weak electrostatic repulsion may play a crucial role for further stabilisation of mesoglobules and improving their monodispersity . Nanoparticles carrying functional groups properly located in the mesoglobule can be obtained, which could be further modified. Such nanoparticles with well controlled size distribution can find a broad range of applications in pharmaceutical, biotechnological and cosmetic industries, where the technological control is one of the major problems.

## ACKNOWLEDGMENTS

The authors acknowledge interesting discussions with Professor H. Orland, Professor T. Garel and our colleague Dr A.V. Gorelov. This work was supported by grants SC/99/186 and FR/2000/019 from Enterprise Ireland.


text:
# Neutrino Unification

\[

hep-ph/0011150

FERMILAB-Pub-00/250-T

IFIC/00-68, IFT-00/26

## Abstract

Present neutrino data are consistent with neutrino masses arising from a common seed at some “neutrino unification” scale $`M_X`$. Such a simple theoretical ansatz naturally leads to quasi-degenerate neutrinos that could lie in the electron-volt range with neutrino mass splittings induced by renormalization effects associated with supersymmetric thresholds. In such a scheme the leptonic analogue of the Cabibbo angle $`\theta _{}`$ describing solar neutrino oscillations is nearly maximal. Its exact value is correlated with the smallness of $`\theta _{reactor}`$. These features agree both with latest data on the solar neutrino spectra and with the reactor neutrino data. The two leading mass-eigenstate neutrinos present in $`\nu _e`$ form a pseudo-Dirac neutrino, avoiding conflict with neutrinoless double beta decay.

\]

The standard model (SM) and its minimal supersymmetric extension (MSSM) fail only in accounting for the solar and atmospheric neutrino data  which strongly indicate the need for neutrino conversions, as would arise from neutrino-mass-induced oscillations . It is thus reasonable to investigate simple theoretically motivated ansatz that can account for these observations.

Inspired by the idea of unification of fundamental interactions  here we propose the idea of neutrino unification: namely that the neutrino mass and mixings observed at low energies take a very simple form at some high energy scale $`M_X`$. Thus, we add to the basic Lagrangean the dimension–five operator

$`{\displaystyle \frac{\lambda _0\delta _{ab}}{M_X}}(\varphi \mathrm{}_a)(\mathrm{}_b\varphi )+H.c.`$ (1)
where $`\mathrm{}_a`$ denote the three lepton doublets and $`\varphi `$ is the standard Higgs doublet. As a working hypothesis, we assume this operator to be characterized at the scale $`M_X`$, by a single real dimension–less parameter $`\lambda _0`$. The breaking of the electroweak symmetry due to a non-zero vacuum expectation value (VEV) $`\varphi `$ will generate, in addition to the known SM masses, the seesaw-type neutrino mass operator $`M_\nu `$

$$M_\nu =\frac{\varphi ^2}{M_X}\lambda _0$$
(2)
For $`M_XM_{Planck}`$ neutrino masses would be too small to account for the atmospheric neutrino anomaly. Thus we will adopt an intermediate value for the $`M_X`$ scale. In the basis in which the charged lepton Yukawa couplings are diagonal (weak basis), $`\lambda _0\delta _{ij}`$ gets transformed into

$$\mathrm{\Lambda }=\lambda _0\mathrm{\Omega }^T\mathrm{\Omega }$$
(3)
where $`\mathrm{\Omega }`$ is an arbitrary unitary matrix. After the electroweak symmetry breaking we can rotate back the neutrino fields, $`\nu =\mathrm{\Omega }^1\nu ^{}`$, and in this basis (the charged lepton mass matrix is still diagonal) we get

$`{\displaystyle \frac{\lambda _0v^2}{M_X}}\nu ^{}\nu ^{}{\displaystyle \frac{g}{\sqrt{2}}}\overline{e_L}\overline{\sigma }^\mu W_\mu ^{}U\nu ^{}+H.c.`$ (4)
where $`U=\mathrm{\Omega }^{}`$ and the fields $`\nu ^{}`$ are the neutrino mass eigenstates (in this paper we number the eigenvalues so that $`\mathrm{\Delta }m_{21}^2=m_2^2m_1^2\mathrm{\Delta }m_{}^2`$ corresponds to the solar neutrino oscillations and $`\mathrm{\Delta }m_{32}^2=m_3^2m_2^2\mathrm{\Delta }m_{atm}^2`$ corresponds to the atmospheric neutrino oscillations). In general, the mixing matrix $`U`$ is characterized by three mixing angles and three $`CP`$ violating phases, one Dirac plus two Majorana-type phases  and can be written as the product of $`V_{CKM}`$ by

$`\left(\begin{array}{ccc}1& 0& 0\\ 0& e^{i\alpha _1}& 0\\ 0& 0& e^{i\alpha _2}\end{array}\right)`$ (5)
where we have chosen to parametrize $`V_{CKM}`$ exactly as the Cabibbo-Kobayashi-Maskawa matrix describing quarks, with the familiar three angles and one phase $`\delta `$. It is clear that with our ansatz (1) we have an additional freedom of performing any arbitrary real $`3\times 3`$ rotation of the neutrinos. The general form of such a restricted matrix $`U`$ can be derived by applying a convenient parametrization  of a general unitary matrix to the matrix $`M\mathrm{\Omega }^T\mathrm{\Omega }`$ and by imposing the symmetry condition $`M=M^T`$. Using the freedom of re-phasing the charged leptons, one can show that the $`U`$ matrix cannot be eliminated by redefining the neutrinos. Thus, despite the mass degeneracy amongst the neutrinos, they exhibit non-trivial mixing effects . In contrast to the case investigated in ref.  where neutrinos were massless but non-orthonormal, here it is the fact that the mass eigenvalues have different phases which prevents one from rotating away the mixing amongst the neutrinos. Nevertheless, the structure of the mixing in our case is more restricted than the general case, since only two angles and one phase remain physical. The above construction has the virtue of explicitly listing the relevant mixing parameters for exactly degenerate Majorana neutrinos. Once the neutrino mass degeneracy is lifted by quantum corrections the mixing effects require the full matrix $`U`$ with, however, a non-trivial relationship amongst the three mixing angles.

Non-zero values of the Dirac phases mean $`CP`$ non-conservation, while Majorana phases of $`\pi /2`$ indicate different $`CP`$ parities of the neutrino mass eigenstates and do not imply $`CP`$ violation . Thus in our case $`CP`$ is conserved if the phase $`\delta =0`$ mod $`\pi `$ and $`\alpha _1,\alpha _2=0`$ mod $`\pi /2`$. Our main conclusions do not depend on whether or not $`CP`$ is conserved in the neutrino sector and from now on we will assume it is conserved. Different $`CP`$ parities of neutrino mass eigenstates can then be accounted for by different signs of neutrino masses (which for $`\alpha _1=\pi /2`$ and/or $`\alpha _2=\pi /2`$ can be achieved by simply absorbing these phases in the neutrino fields).

Now we turn to the renormalization effects. One-loop renormalization group equations (RGEs) for the $`\mathrm{\Lambda }`$ coefficients characterizing the dimension-5 non-renormalizable terms can easily be written down both for the SM and MSSM . The RGEs take a particularly useful form when written directly for mass eigenvalues and for the elements of the mixing matrix . The flavour independent corrections to the $`\mathrm{\Lambda }=𝒪\text{ }(1)`$ coefficients are irrelevant for our discussion. The flavour-dependent corrections are due to lepton Yukawa couplings and, to a good approximation, determined by the $`\tau `$ Yukawa coupling. Using e.g. the results of ref.  one can easily verify that, with the ansatz (1), such effects cannot explain $`|\mathrm{\Delta }m_{32}^2||\mathrm{\Delta }m_{21}^2|`$ with phenomenologically acceptable mixing angles.

It is remarkable that supersymmetry can induce flavour-dependent threshold corrections associated with slepton mass splittling  which can dominate over the $`\tau `$ Yukawa corrections. We show in this paper that such corrections can lead to the desirable low energy pattern with the ansatz (1). This is possible provided the soft supersymmetry breaking scalar mass terms deviate sufficiently from universality.

The quantum corrections to our ansatz (1) in the basis defined by eq. (4) are given by

$`m_\nu ^{ab}=m_a\delta _{ab}+m_a(U^TIU^{})_{ab}+m_b(U^{}IU)_{ab}`$ (6)
where in our case all $`|m_a|=m`$ up to a common flavour-independent renormalization $`𝒪\text{ }(1)`$ factor. Here the indices $`a,b`$ refer to our starting basis used in eq. (1). The correction $`I`$ (calculated in the electroweak charged lepton mass eigenstate basis) consists of two parts

$`I=I^{RG}+I^{TH}`$ (7)
where $`I^{RG}`$ is the renormalization group correction  and $`I^{TH}`$ are the electroweak scale threshold corrections . Assuming no lepton flavour violation in other sectors of the theory (e.g. in the case of supersymmetry) the matrix $`I`$ is diagonal, $`I_{AB}=I_A\delta _{AB}`$ (the index $`A`$ refers to the electroweak basis). Without loss of generality in our discussion we take $`I_\mu =0`$. We adopt the general form for the CP conserving mixing matrix $`U`$

$`\left(\begin{array}{ccc}c_{12}c_{13}& s_{12}c_{13}& s_{13}\\ s_{12}c_{23}c_{12}s_{23}s_{13}& c_{12}c_{23}s_{12}s_{23}s_{13}& s_{23}c_{13}\\ s_{12}s_{23}c_{12}c_{23}s_{13}& c_{12}s_{23}s_{12}c_{23}s_{13}& c_{23}c_{13}\end{array}\right)`$ (8)
in which $`\theta _{13}\theta _{reactor}`$, and allow for different CP parities of neutrino eigenstates. For degenerate masses $`|m_a|=m`$ the matrix (8) contains one redundant angle and this degree of freedom is used to re-diagonalize the mass matrix (6). As long as we assume CP conservation we can use the standard perturbation theory for degenerate zeroth order eigenvalues in order to determine the non-trivial relation amongst the three mixing angles in eq. (8). Taking $`|I_\tau ||I_e|`$ we easily recover the conclusions of ref. . In this letter we investigate the alternative possibility that $`|I_e||I_\tau |`$. We begin with $`m_1=m_2=m_3m`$ in which case

$`m_\nu ^{ab}=m\left(\begin{array}{ccc}1+2U_{A1}^2I_A& 0& 2U_{A1}U_{A3}I_A\\ 0& 12U_{A2}^2I_A& 0\\ 2U_{A1}U_{A3}I_A& 0& 1+2U_{A3}^2I_A\end{array}\right)`$ (9)
Since $`m_1=m_3`$, the ordinary perturbation calculus tells us that the neutrino mass basis should be chosen so that the off-diagonal $`13`$ entry of the perturbation is zero. We therefore require that $`U_{A1}U_{A3}I_A=0`$ which fixes the redundant angle in the matrix (8). With our assumptions about $`I_A`$’s ($`|I_e||I_\tau |`$) this gives

$`s_{13}={\displaystyle \frac{s_{12}}{c_{12}}}s_{23}c_{23}r+𝒪(r^2)`$ (10)
where $`rI_\tau /I_e`$. It is then trivial to compute the corrected mass eigenvalues and one finds

$`\mathrm{\Delta }m_{atm}^24m^2I_es_{12}^2,`$ (11)
$`\mathrm{\Delta }m_{}^24m^2I_e\left[C_{12}\left(1s_{23}^2r\right)+(1+2c_{12}^2)s_{13}^2\right]`$ (12)
where $`C_{12}c_{12}^2s_{12}^2`$. For maximal $`\theta _{12}=\pi /4`$ one obtains that $`\mathrm{\Delta }m_{atm}^2I_e`$ and $`\mathrm{\Delta }m_{}^2I_er^2`$, in agreement with the experimental requirement $`\mathrm{\Delta }m_{}^2\mathrm{\Delta }m_{atm}^2`$.

Eqs. (11) and (12) lead to the following relation between the neutrino oscillation parameters

$$\mathrm{\Delta }m_{}^22\mathrm{\Delta }m_{atm}^2\left[C_{12}\left(1s_{23}^2r\right)+2s_{13}^2\right]$$
(13)
which is the prediction of our model, independent of the value of $`I_e`$. Since we need $`c_{12}^2s_{12}^2`$ we conclude that small mixing solution to the solar neutrino problem is ruled out in our neutrino unification picture.

It is easy to check that with our assumptions about $`I_A`$’s the same mechanism works in the case $`m_1=m_2=m_3`$. However, in the case $`m_1=m_2=m_3`$ the condition $`U_{A1}U_{A2}I_A=0`$ means that $`s_{12}c_{12}c_{13}^2=𝒪(r)`$ which means that one cannot have simultaneously large solar neutrino mixing and remain consistent with the reactor experiments  (this is also incompatible with the SMA solution to the solar neutrino problem because for $`s_{12}c_{12}0`$ one gets $`\mathrm{\Delta }m_{}^2\mathrm{\Delta }m_{atm}^2`$).

Since $`\mathrm{\Delta }m_{atm}^23\times 10^3`$eV2 is rather fixed by the analysis of the data  and since, as we shall see, typically $`|I_e|{}_{}{}^{}{}_{}{}^{<}10^3`$, it follows that $`m{}_{}{}^{}{}_{}{}^{>}`$ $`𝒪`$ (1) eV. Thus this picture leads to neutrino masses in the range accessible to $`\beta `$ decay and hot dark matter searches. Different CP parities of the neutrinos are however sufficient to ensure a destructive interference in the amplitude for neutrinoless double beta decay which is suppressed by the fact that the two leading mass eigenstate neutrinos in the $`\nu _e`$ form a pseudo-Dirac neutrino .

How can one realize our neutrino unification scenario? First of all, it requires the existence of new states giving a dominant contribution to $`I_e`$. Rather than extending the gauge and/or multiplet structure of the SM, we prefer to ask whether in supersymmetry one can arrange the corrections $`I_A`$ to satisfy our requirements. Taking a diagonal slepton mass matrix (in the same basis in which the charged leptons mass matrix is diagonal) and taking into account only the contribution of pure Wino of mass $`\stackrel{~}{m}`$ one finds for $`I_AI(y_A)`$ the following form

$`I_A={\displaystyle \frac{g^2}{32\pi ^2}}\left\{{\displaystyle \frac{1}{y_A}}+{\displaystyle \frac{y_A^21}{y_A^2}}\mathrm{ln}(1y_A)(A\mu )\right\}`$ (14)
where $`y_A=1(M_A/\stackrel{~}{m})^2`$, $`M_A`$ is the $`A`$-th charged slepton mass. It is easy to see that for $`I_e10^3`$ and positive we need $`M_e1.7M_{\mu ,\tau }`$ for the charged slepton masses.

Such slepton mass patterns can arise in models with inverted hierarchy , although the proposed mechanism may have other realizations. Being an operator of dimension five, the nontrivial neutrino mass matrix does not affect the evolution of the lepton Yukawa and slepton mass squared matrices. Hence, if there exists a basis at $`M_X`$ in which both lepton Yukawa and slepton mass squared matrices are simultaneously diagonal (fermion-sfermion mass alignment  present e.g. in models with abelian $`U(1)`$ gauge symmetry ), this basis remains unchanged during the RG evolution to the electroweak scale and hence the slepton masses remain diagonal. Thus in this case lepton flavor-violating decays such as $`\mu e\gamma `$ are not induced.

In Fig. 1 we illustrate the neutrino unification idea by displaying the supersymmetric evolution of solar and atmospheric neutrino mass splittings as a function of energy, from the high energy scale $`M_X`$ where neutrino masses unify, down to the weak scale for tan$`\beta =2`$ and a starting neutrino mass at $`M_X`$ of 1.2 eV.

Above the supersymmetry breaking scale the neutrino mass-square differences are $`\mathrm{\Delta }m_{21}^2=\frac{1}{3}\mathrm{\Delta }m_{31}^2=\frac{1}{2}\mathrm{\Delta }m_{32}^2m^2I_\tau `$. After threshold corrections they change according to Eqs. (11,12,13). Depending on the chosen value of $`M_X`$ this can fit large mixing angle solutions to the solar neutrino problem. For example for $`M_X=10^{13}`$ GeV and $`I_e2\times 10^3`$, $`r=\frac{I_\tau }{I_e}=0.029`$ we obtain quasi-degenerate neutrino masses $`|m_i|=0.9`$ eV split by $`\mathrm{\Delta }m_{21}^2=1.52\times 10^5`$ eV2 and $`\mathrm{\Delta }m_{32}^2=3.13\times 10^3`$ eV2. Neutrino mixing is bi-maximal $`s_{12}s_{23}0.707`$ with $`s_{13}0.0143`$, nicely matching the LMA solution. In order to achieve the LOW solution one needs a smaller $`I_\tau `$ value. Barring a possible cancellation between the terms in the r.h.s of Eq. (13) this requires a low neutrino unification scale $`M_X=10^5`$ GeV. In this case we find that quasi-degenerate masses $`|m_i|=1.1`$ eV for $`r=\frac{I_\tau }{I_e}=0.0103`$ and $`I_e1.3\times 10^3`$ with $`\mathrm{\Delta }m_{21}^2=2.74\times 10^7`$ eV2 and $`\mathrm{\Delta }m_{32}^2=3.16\times 10^3`$ eV2. Mixing is nearly bi-maximal with $`s_{13}0.0051`$. Due to the observed flatness of the latest recoil electron energy spectrum, present solar neutrino data globally prefer large mixing angle MSW solutions over the small mixing solution (SMA) . This feature fits well in our neutrino unification scheme, where the SMA solution can not be realized. From ref.  we find that in our scheme the LMA solution is present at the 97$`\%`$ C.L. while the LOW solution appears already at 90$`\%`$ C.L.

Conclusions.

We have proposed a simple theoretical ansatz where neutrino masses arise from a common dimension-5 operator at some “neutrino unification” scale $`M_X`$. This scale must be lower than that which characterizes the unification of gauge and Yukawa couplings. The ansatz naturally implies quasi-degenerate neutrinos at the electron-volt range, thus accessible to future $`\beta `$ decay and hot dark matter searches. However neutrinoless double beta decay is naturally suppressed due to the intrinsic Majorana neutrino CP parities. We have shown that neutrino splittings can be induced by supersymmetric thresholds arising from non-universal soft breaking terms, without necessarily conflicting with limits on lepton flavour violating processes such as $`\mu e\gamma `$. A non-trivial consequence of the neutrino unification idea is that the solar and atmospheric data can only be accounted for in terms of large mixing angle-type MSW oscillations. Note also that the consistency of our scheme correlates the smallness of $`\theta _{13}`$ indicated by reactor data and, to a lesser extent, also by the atmospheric data , to the largeness of $`\theta _{12}`$ required by the solar neutrino data.

Acknowledgements: Work supported by Spanish DGICYT under grants PB98-0693 and SAB/1998-0136 (A.I.), by the European Commission TMR network HPRN-CT-2000-00148, by the European Science Foundation network grant N. 86 and by the Polish State Committee for Scientific Research grant KBN 2 P03B 060 18. P.H.Ch. thanks the CERN Theory group for hospitality during the completion of this paper.


text:
# On the Origin of Faint Intracluster Starlight in Coma

## 1 INTRODUCTION

Gregg & West (1998) have recently reported the discovery of intracluster starlight features in the Coma cluster of galaxies. The features are revealed in R-band images of the Coma cluster, and are also seen in the B, V and I bands, confirming their existence. The plume-like features are tens of kiloparsecs in width, and extend to radii of more than a hundred kiloparsecs. The mean R-band surface brightness of the plumes is low, about 26 mag/arcsec2, but their absolute magnitude $`19.0`$ mag is comparable to the absolute R-band magnitude of a normal elliptical galaxy. Possible explanations of the origin of the plumes invoke the idea of the disruption of galaxies by galaxy-galaxy interactions (Dubinski et al 1996, Moore et al 1996), or the tidal stripping of galaxies by the global cluster gravitational field (Merritt 1984). However, the stellar plumes observed in the Coma cluster can not be produced so easily. Using N-body simulations, we have studied the possibility of forming stellar plumes with surface brightnesses 6–7 magnitudes fainter than the surface brightness in the centers of the normal elliptical galaxies similar to those observed in the Coma cluster. We have found that the tidal interaction of galaxies with the gravitational field of the cluster core cannot be responsible for the origin of intracluster starlight features in Coma. A head-on collision of two normal elliptical galaxies with a relative velocity less than the cluster velocity dispersion, or tidal interactions with intracluster density inhomogeneities might explain the origin of stellar plumes with the observed properties.

## 2 MODELS AND NUMERICAL METHOD

To model galaxy disruptions by the cluster core or galaxy-galaxy collisions, we use the parameters of the elliptical galaxy NGC 3377. The disruptions of the disk galaxies might also lead to the formation of stellar plumes. Nearby clusters of galaxies, however, are enriched with elliptical and lenticular galaxies (Dressler 1980). The elliptical galaxy NGC 3377 was chosen by Gregg & West (1998) in their comparison between normal galaxies and the properties of the Coma stellar plumes. NGC 3377 is well studied observationally. Kormendy et al (1998) published an accurate rotation curve of this galaxy, along with a composite surface brightness profile along its major axis. The rotation curve of NGC 3377 shows a steep kinematic gradient near the center of the galaxy, but it remains flat for radii $`r>3^{\prime \prime }`$, with a circular velocity about 100 km/sec. This yields a total mass for NGC 3377 of about $`3\times 10^{10}M_{}`$, if the galactic radius is $`140^{\prime \prime }`$ or 15 kpc assuming a Hubble constant $`H=50\mathrm{k}\mathrm{m}\mathrm{s}^1\mathrm{Mpc}^1`$, as used in our further simulations.

Using the rotational curve near the center, the central density of NGC 3377 can be rather firmly evaluated to be $`2.2\times 10^{21}`$ g/cm3 within the central $`3^{\prime \prime }`$. Approximately half of the mass can be associated with the central black hole, and we adopt in our simulations the value $`10^{21}`$ g/cm3.

For the radial dependence of the galactic density, we use the Plummer model, and Hernquist’s (1990) density distributions with core radii of 0.8 kpc and 0.6 kpc respectively. We adopt an analytic King approximation to the isothermal sphere (Sarazin 1986) for the dark matter density distribution in the Coma cluster. Hughes (1998) re-examined mass estimates for the Coma cluster derived from X-ray spectral observations. At the 99% confidence level, he found that the mass of the Coma cluster within a radius of 1 Mpc is $`(6.2\pm 0.9)\times 10^{14}M_{}`$. The dark matter core radius for the adopted matter profile is about 270-370 kpc (Hughes 1998), which gives a central dark matter density of $`(1.02.9)\times 10^{15}M_{}`$ / Mpc3. In our simulations, we adopt “extreme” disruptive parameters of the Coma cluster which have the values $`r_c=270`$ kpc, and $`\rho _0=3\times 10^{15}M_{}`$ / Mpc3.

The simulations were performed with a TREE code using N=32768 particles, a tolerance parameter equal to 0.75, and a softening length of 0.03125.

## 3 RESULTS AND DISCUSSION

We considered different possibilities for the formation of stellar plumes in clusters of galaxies, namely the tidal disruption of an elliptical galaxy by the cluster core, direct galactic collisions, and the interaction of an elliptical galaxy with hypothetical cluster substructures.

Collisions or tidal interactions often produce plumes of particles in numerical simulations. We need, however, a quantitative estimate of whether such plumes could be responsible for the observed stellar structures seen in Coma. To do this, we calculated the theoretical surface brightness in the plume assuming constant mass-to-luminosity ratios for the plume and the galaxy. The characteristic surface brightness in the central regions of elliptical galaxies in the V-band is about 20 mag/arcsec2 (Binney & Merrifield 1998), which for the R-band, gives a value of about 19.5 mag/arcsec2 (Frei & Gunn 1994). We assumed that our model galaxy has the same surface brightness in central regions before an interaction, and that tidal interactions do not change the stellar content of the plume. These assumptions allow us to estimate surface brightness of resulting features in our numerical simulation by comparing corespondinc column densities of particles. The last assumption excludes from our consideration the possible role of star formation triggered in the plumes by tidal interactions. Although tidally induced star formation in the plumes is likely to take place for the collisions of gas rich disk galaxies, this is not the case for the stellar plumes observed in Coma. According to Gregg & West (1998), the colors in the plumes are consistent with those of evolved stellar polulations observed in intermediate brightness normal elliptical galaxies. The photometric uncertainties, and uncertainties of the removal of superimposed objects leave, however, open the possibility of significant star fromation in the plumes.

Direct galactic collisions. Figure 1 shows the result of a direct collision of two equal-mass Plummer spheres taken at time $`1.7\times 10^8`$ years after the beginning of the interaction. The initial separation of the galaxies and their initial relative velocity were taken to be 32 kpc and 408 km/sec. The collision forms a broad stellar plume connecting both galaxies, with the surface brightness in the plume about 26.5 mag/arcsec2(Figure 2). The direct collisions of equal-mass galaxies might therefore constitute a possible explanation of the origin of the observed stellar plumes in Coma. The high velocity dispersion of about 1000 km/sec in Coma, however, makes this mechanism rather ineffective. Gregg & West (1998) argue that the galactic subgroups in Coma, with their relatively low internal velocity dispersions, will increase the effectiveness of the disruption of galaxies, and our simulations concur with this conclusion. The characteristic velocity dispersion within groups, which are fairly uniformly distributed in space and which show no preference for rich clusters, is about 250 km/sec (Hickson 1997). This value is much smaller than the typical velocity dispersions in rich clusters of galaxies. If a cluster is not virialized, then strong local interactions can occur within in-falling subgroups of galaxies (Conselice & Gallaher 1998). Burns et al (1994) present convincing arguments that such infalling groups of galaxies do indeed exist. They found that the galactic group NGC 4839 has fallen into the Coma Cluster, passing through the cluster core about 2 Gyr ago. The NGC 4839 galactic group cannot be responsible, however, for the observed starlight features in Coma. All four stellar features are located within $`0.5`$ Mpc of the cluster center, and cannot survive for more than 1-2 Gyr (Gregg & West 1998), which is shorter than the epoch of passage of the NGC 4839 group through the central regions of Coma.

Core tidal interactions. The tidal interaction of galaxies with the gravitational potential of the cluster core is a likely explanation for the origin of stellar intracluster structures in Coma. We simulated an extreme case of such an interaction with a galaxy plunging into the cluster core. In these simulations, the cluster core radius adopted was 270 kpc, and the core density equal to $`3\times 10^{15}M_{}`$ / Mpc3 which are the extreme “disruptive” parameters of the mass distribution in Coma. The galaxy was placed 600 kpc away from the cluster center with zero initial relative velocity. Figure 3 shows surface density contour plots taken at time $`T=3.36\times 10^9`$ years after the beginning of the simulations, when the galaxy made one oscillation in the cluster core. The surface brightness in the patchy plume is 32-34 mag/arcsec2, and apparently cannot explain the origin of the bright starlight features in Coma. This result is also illustrated by Figure 4, which shows the fraction of escaped particles as a function of time. A total of 10 % of the mass is stripped during the interaction time. This stripped mass cannot be regarded as the main source of intracluster stellar material in Coma.

Tidal disruption by the “local core”. One more possibility exists for the formation of the bright intracluster structures seen in Coma. The X-ray surface brightness measurements of Coma cluster reveal the regions of the enhanced X-ray emission with characteristic angular size less than the core radius of the cluster ( Soltan & Fabricant 1990). The intracluster stellar features in Coma may result from tidal interactions with small scale density fluctuations associated with the enhanced X-ray emission. Vikhlinin et al (1994) found indeed such mass concentrations in the core region of the Coma cluster. Both of the discovered density enhancements are centered at the positions of the brightest elliptical galaxies in the cluster: NGC 4874 and NGC 4889. Assuming that X-ray surface brightness enhancements result from the density enhancements of the underlying dark matter, Vikhlinin et al (1994) estimated the masses of inhomogeneities within 80 kpc to be about $`10^{13}M_{}`$. Remarkably, the most prominent plume-like structure, the LSB object 1 discovered by Gregg and West, is adjacent to NGC 4874. The core radius of the density enhancement centered at the position of NGC 4874 is about 100 kpc (Vikhlinin et al 1994) which gives the density of the “local core”about $`5\times 10^{25}`$ g/cm3. This value is too low to disrupt noticeably a galaxy by tidal interaction.

The parameters of the density enhancement centered at the position of NGC 4889 are more favorable for the galactic disruption. This density enhancement has core radius about 24 kpc, and the central density about $`8\times 10^{24}`$ g/cm3. We simulated an interaction of a galaxy with the NGC 4889 density enhancement. A galaxy was placed 200 kpc away from the density enhancement with the initial velocity 500 km/s. The galaxy was allowed then to collide non-centrally with the density enhancement with the impact parameter equal to 24 kpc. Figure 5 shows the result of such high-speed tidal interaction. Interaction with NGC 4889 density inhomogeneity forms a plume with the surface brightness only 29-30 mag/arcsec2. This value is too low to account for the properties of the observed plumes. The density enhancements, however increase an effective cross-section of the direct collision with the giant cD galaxies nested in their central regions. This is a likely scenario of the formation of the most luminous plume LSB no 1 adjusted to NGC 4884.

Stars in the plume structures which are induced in direct galactic collisions have higher velocity dispersion compared to the plumes produced in tidal interactions. The collisional plumes live therefore shorter than tidally induced structures. Our numerical simulations show that the density inhomogeneities tidally disrupt an elliptical galaxy if the central density of a “local core” is about 10 percent of the central galactic density. Such value can be achieved if an inhomogeneity has total mass about one thousand of a mass of a normal elliptical galaxy, and a core radius $``$ 15 kpc. If such inhomogeneities of the intracluster matter do really exist, then the starlight features in clusters of galaxies give an opportunity to probe the cluster density fluctuations.

## 4 SUMMARY

Our results can be summarized as follows:

1. We found that the observed intracluster starlight structures in Coma cannot be produced by the tidal interaction of an elliptical galaxy with the gravitational potential of the cluster core.

2. The head-on collision of two normal ellipticals with relative velocities lower than the velocity dispersion of the Coma cluster might be a mechanism for the formation of the stellar plumes observed in Coma. The subgroups of galaxies within the Coma cluster have relatively low velocity dispersions, and therefore direct galactic collisions within the subgroups are currently the most plausible explanation for the origin of the stellar plumes.

3. The stellar plumes seen in Coma may result from the tidal interactions of galaxies with local inhomogeneities of the cluster density distribution. If the existence of such inhomogeneities is confirmed in future observations, then intracluster starlight features will provide an opportunity for studying their nature.

The authors thank W. van Altena, G. Laughlin and C. Boily for comments on the manuscript resulted in a much improved presentation. VK is grateful to the National Astronomical Observatory Japan for the financial support.