text:
# Photo-induced volume changes in selenium. Tight-binding molecular dynamics study.

## I Introduction

The photo-induced structural changes during illumination Shimakawa95 ; Morigaki99  in chalcogenide glasses can be classified into two groups, those which are irreversible and which can be reversed. An other possible sorting is the following: there are defect related metastability and there are structure related one. Several investigations have been carried out in amorphous selenium (a-Se) in order to provide an acceptable explanation of these phenomena Fritzsche96 ; Kolobov99 ; Zhang99 ; Ganjoo00 ; Palyok01 . Present capacity of fast computers allow us to simulate the photo-induced structural changes i.e. to follow the time development of these phenomena by molecular dynamics (MD). Our particular aim was to determine how the structural changes occur in a-Se due to illumination. We performed tight-binding molecular dynamics simulations to study the dynamics of a Se8 ring, a linear chain and glassy Se networks before, during and after the illumination.

## II Simulation method

We have developed a MD computer code (ATOMDEP program package) to simulate the real preparation procedures of amorphous structures (growth by atom-by-atom deposition on a substrate and rapid quenching). In our recent workKohary01 , the growth of amorphous carbon films were simulated by this MD method. This computer code is convenient to investigate photo-induced volume changes as well if the built-in atomic interaction can handle the photo-excitation. Standard velocity Verlet algorithm was applied in our MD simulations in order to follow the atomic scale motion. We chose $`\mathrm{\Delta }t`$ = 1 fs or 2 fs for the time step, depending on the temperature. These values are the typically used time steps for such simulations. To control the temperature of the systems we applied the velocity-rescaling method.

For calculating the inter-atomic forces in a-Se we used tight-binding (TB) Molina99  and self-consistent tight-binding (SCF-TB) Lomba00  models. The TB parametrization Molina99  has been introduced for disordered selenium following the techniques developed by Goodwin et al Goodwin89 . It was thoroughly tested by molecular dynamics calculations in liquid and amorphous phases and the results were compared to experiments and to ab initio calculations. The agreement with experiments and ab initio  calculations are rather good apart from the fact that the number of coordination defects in the solid and liquid phases are higher than the experimentally measured values. The authors improved their TB Hamiltonian by including the Hubbard correction Lomba00 . This implies that either the algorithm have to be made self-consistent or perturbation theory must be applied. Our choice was the first alternative. The convergence for this SCF-TB calculations was controlled by the charge differences on the atomic sites. Convergence criteria were considered to be satisfied if the deviation of atomic charges between the actual and the previous iterations was less than 0.01 electron/atom.

## III Sample preparation

### III.1 Cluster geometries

Crystalline forms of selenium consist of chains and eight-membered rings. Typical bond lengths are around 2.35 Å while most of the bond angle values are around $`103^{}`$. It is very likely that these local arrangements can be found in non-crystalline forms of selenium as well. The initial configuration of the eight-membered ring in our simulation had bond lengths of 2.38 Å and bond angles of 102. Dihedral angles were equal to 100 (See Fig. 1. top panel). For the eighteen-membered selenium chain (with 1-dimensional periodic boundary condition) these values were: 2.36 Å, 100, and 98, respectively. Every Se atom had two first-neighbors, i.e. there was no coordination defect.

### III.2 Glassy structure

Glassy selenium networks having 162 atoms (with three dimensional periodic boundary conditions) were prepared by “cook and quench” technique (from 3000 K to 20 K). The average quenching rate was $`410^{12}`$ K/s. In our simulation after quenching the periodic boundary constraint in the z direction was released and and the system was allowed to relax for another 100 ps. In this way the system could change its volume in the z direction. Four models have been prepared with different initial topologies and densities. The snapshot of one of our models can be seen in Fig.2.

## IV Results

### IV.1 Photo excitation of small Se clusters

Before illumination (photo excitation) the individual ring and the chain were relaxed for 4 ps at T=500 K. During this period the structures were stable. When a photon is absorbed an electron from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO) is transferred. This is a simple model of photo-excitation when an electron is shifted from the valence band to the conduction band (electron-hole pair creation). In our study electron-hole pair recombination was also simulated by taking back the electron from the conduction band to the valence band. The lifetime of the excitation was 200 fs.

After excitation one bond length in the ring started to increase and bond-breaking occurred. Two snapshots of this process can be seen in Figure 1. A similar result was published Shimojo98  by a Japanese group for S8. They performed MD simulation within the framework of density functional theory in the local density approximation. In our second MD simulation we investigated the linear chain structure. The same procedure was performed to model the excitation. An electron from HOMO to LUMO was transferred. Very similar result was obtained; a bond inside the chain was broken immediately after a HOMO electron was excited.

### IV.2 Photo excitation of glassy Se networks

The excitation procedure for the glassy networks was different from the method applied for small clusters. The hole procedure was carried out at the temperature of 20 K. Instead of a single electron excitation, five electrons were transferred from the valence band to the conduction band. The first excited electron was taken from the HOMO, the second excited electron was taken from the level below it and so on. A difficulty applying this procedure appeared during the simulation. The HOMO and the LUMO levels shifted into the middle of the gap after excitation and even level crossing was observed. After that point the result obtained by MD simulation was uncertain and useless. It could cause incorrect conclusions. To avoid this problem we applied the following calculation method: we separately calculated the influences of electron and hole creation and annihilation.

Our systems were open only along the z direction and periodic boundary conditions were applied in the other two directions. To measure the length (i.e. volume) of the sample we have proceeded in the following way. We have considered 40 atoms in z direction closest to the either ends of the sample. But we neglected those 10-10 atoms at the ends of these clusters, which had the largest (smallest) z coordinated to avoid the surface effects. The length of the sample was defined as the average z coordinates of 30 rest atoms at the left end of the sample minus the average z coordinates of 30 rest atoms at the right end.

The total length of four samples with average density of 4.5 g/cm3 is shown in Fig. 3. Arrows demonstrate the electron excitations and deexcitations. Only one excitation out of five caused volume decrease. The total influence is volume expansion. In the opposite process three out of five decrease, one increases the length. Final conclusion is that photo-induced volume expansion occurred during illumination and remained after finishing illumination. A typical atomic scale rearrangement after one photo-induced electron excitation is shown in Fig. 5.

More interesting is the number of coordination defects. Basically, selenium atoms have two nearest-neighbors. Average numbers of one-fold coordinated and threefold-coordinated atoms can be seen in Fig 4. (Practically, the number of twofold coordinated atoms remained the same in our simulations.) During illumination the number of one-coordinated defects was increasing while it was opposite for threefold coordinated selenium atoms. During the deexcitation the inverse process occurred. It is striking that there are more threefold defects and less one-fold coordinated atomic sites at the end of the simulation compared to initial configuration. All these observations suggest the atomistic explanation of volume changes in chalcogenide systems.

## V Summary

We have developed a molecular dynamics computer code to investigate the time development of ring and chain-like structures after photo-excitation. Furthermore, we investigated the photo-induced volume change in glassy networks. Our results suggest photo-induced bond breaking and the change of the number in coordination defects in these conditions. This might be an explanation of photo-induced volume changes in a-Se and other chalcogenide glasses.

## VI Acknowledgments

This work has been supported by the OTKA Fund (Grant No. T043231) and by Hungarian-British and Hungarian-Japanese intergovernmental S and T Programmes (No. GB-17/2003 and No. JAP-7/2000). The computer simulations were partly carried out at the Tokyo Polytechnics University (Japan). We are indebted to Prof. T. Aoki (TPU) for providing us this possibility.


format: text

text:
# Effects of two-site composite excitations in the Hubbard model

## Abstract

The electronic states of the Hubbard model are investigated by use of the Composite Operator Method. In addition to the Hubbard operators, two other operators related with two-site composite excitations are included in the basis. Within the present formulation, higher-order composite excitations are reduced to the chosen operatorial basis by means of a procedure preserving the particle-hole symmetry. The positive comparison with numerical simulations for the double occupancy indicates that such approximation improves over the two-pole approximation.

url\]http://www.scs.sa.infn.it

The analysis of highly correlated electron systems still reports many unsolved issues in spite of the large number of efforts which have been made for several decades. In general, difficulties come from the treatment of the collective excitations emerging in these systems, especially near the Mott-Hubbard transition. In this paper, we study the electronic states of the Hubbard model by use of the Composite Operator Method , which has shown to be capable to describe the physics of strongly correlated systems. Within projection techniques, going beyond the two-pole approximation, which implies the use of a local basis, is a very hard task. As first step along this direction, we here first used a four component operatorial basis with two non-local fields extending over two-sites.

The $`d`$-dimensional Hubbard Hamiltonian reads as follows,

$$H=\underset{ij\sigma }{}(t_{ij}\mu \delta _{ij})c_\sigma ^{}(i)c_\sigma (j)+U\underset{i}{}n_{}(i)n_{}(i),$$

where $`c_\sigma ^{}(i)`$ and $`c_\sigma (i)`$ are creation and annihilation operators of electrons with spin $`\sigma `$ at the site $`i`$, respectively. $`n_\sigma (i)=c_\sigma ^{}(i)c_\sigma (i)`$, $`\mu `$ is the chemical potential, $`t_{ij}=2dt\alpha _{ij}`$, $`\alpha [𝐤]=[\alpha _{ij}]=\frac{1}{d}_{i=1}^d\mathrm{cos}(k_ia)`$, $`a`$ is the lattice constant, $``$ is the Fourier transform, $`U`$ is the on-site Coulomb repulsion. We define the following operatorial basis,

$$\psi _\sigma (i)=\left(\begin{array}{c}\xi _\sigma (i)\hfill \\ \eta _\sigma (i)\hfill \\ \xi _{s\sigma }(i)\hfill \\ \eta _{s\sigma }(i)\hfill \end{array}\right),$$
(1)
where $`\xi _\sigma (i)=c_\sigma (i)\left(1n_\sigma (i)\right)`$ and $`\eta _\sigma (i)=c_\sigma (i)n_\sigma (i)`$ describe the transitions $`n(i)=0`$ $``$ 1 and 1 $``$ 2, respectively. $`\xi _{s\sigma }(i)`$ and $`\eta _{s\sigma }(i)`$ are

$$\begin{array}{c}\xi _{s\sigma }(i)=n_\sigma (i)\xi _\sigma ^\alpha (i)+c_\sigma ^{}(i)c_\sigma (i)\xi _\sigma ^\alpha (i)\hfill \\ +c_\sigma (i)\eta _\sigma ^\alpha (i)c_\sigma (i)\hfill \\ \eta _{s\sigma }(i)=n_\sigma (i)\eta _\sigma ^\alpha (i)+c_\sigma ^{}(i)c_\sigma (i)\eta _\sigma ^\alpha (i)\hfill \\ +c_\sigma (i)\xi _\sigma ^\alpha (i)c_\sigma (i)\hfill \end{array}$$
(2)
with $`c_\sigma ^\alpha (i)=_j\alpha _{ij}c_\sigma (j)`$. These operators describe two-site composite excitations which were not included in the previous work , and are eigenoperators of the interaction term of the Hamiltonian, similarly to $`\xi _\sigma (i)`$ and $`\eta _\sigma (i)`$.

The equations of motion of the basis read as

$$\begin{array}{c}i\frac{}{t}\xi _\sigma (i)=\mu \xi _\sigma (i)zt\left[c_\sigma ^\alpha (i)+\xi _{s\sigma }(i)+\eta _{s\sigma }(i)\right]\hfill \\ i\frac{}{t}\eta _\sigma (i)=(\mu +U)\eta _\sigma (i)+zt\left[\xi _{s\sigma }(i)+\eta _{s\sigma }(i)\right]\hfill \\ i\frac{}{t}\xi _{s\sigma }(i)=\mu \xi _{s\sigma }(i)\hfill \\ zt\left[\frac{2}{z}\eta _\sigma (i)\xi _{s\sigma }^\alpha (i)2\eta _{s\sigma }^\alpha (i)+DJ_\xi (i)\right]\hfill \\ i\frac{}{t}\eta _{s\sigma }(i)=(\mu +U)\eta _{s\sigma }(i)\hfill \\ zt\left[\frac{1}{z}\eta _\sigma (i)\xi _{s\sigma }^\alpha (i)+DJ_\eta (i)\right].\hfill \end{array}$$

where $`z=2d`$ is the coordination number. $`DJ_\xi `$ and $`DJ_\eta `$ are 3-site irreducible operators. By *irreducible* we mean that all local and two-site contributions have been subtracted. Hereafter, we will neglect $`DJ_\xi `$ and $`DJ_\eta `$. It is worth mentioning that the particle-hole symmetry of the model is preserved by this approximation. Then, the thermal retarded Green’s function $`G(𝐤,\omega )=[\psi (i)\psi ^{}(j)]`$ can be expressed as,

$$G(𝐤,\omega )=\underset{n=1}{\overset{4}{}}\frac{\sigma _n(𝐤)}{\omega E_n(𝐤)},$$
(3)
where

$$\begin{array}{c}E_1(𝐤)=\mu zt\alpha [𝐤]\hfill \\ E_2(𝐤)=\mu zt\alpha [𝐤]+U\hfill \\ E_3(𝐤)=\mu +zt\alpha [𝐤]zJ_U\hfill \\ E_4(𝐤)=\mu +zt\alpha [𝐤]+U+zJ_U,\hfill \end{array}$$
(4)
with $`J_U=1/z(\sqrt{4zt^2+(U/2)^2}U/2)`$.

The spectral densities $`\sigma _n(𝐤)`$ contain two correlation functions: $`\mathrm{\Delta }=\xi _{}^\alpha (i)\xi _{}^{}(i)\eta _{}^\alpha (i)\eta _{}^{}(i)`$ and $`p=\frac{1}{4}n_\nu (i)n_\nu ^\alpha (i)(c_{}(i)c_{}(i))^\alpha c_{}^{}(i)c_{}^{}(i)`$. $`\mathrm{\Delta }`$ can be directly computed in terms of the Green’s function. $`p`$ and $`\mu `$ are self-consistently evaluated through the constraint $`\xi (i)\eta ^{}(i)=0`$ and the equation defining the electron number density $`n`$.

To test the reliability of the present approximation, we calculated the double occupancy $`D=n_{}(i)n_{}(i)`$, and compared our results with the numerical data obtained by the Lanczos method for a 18-site system . The agreement is considerably good over the whole $`U/t`$ range and the results show a clear improvement over the ones obtained by the two-pole approximation . The details of the formulation and more extensive comparisons with numerical simulations will be presented elsewhere.


text:
# SB2 and eclipsing binaries with GAIA and RAVE

## 1. GAIA

GAIA is the ambitious flagship mission approved by ESA for a launch not later than 2012, and possibly already by Q2 of 2010. The mission is designed to obtain extremely precise astrometry (in the micro-arcsec regime), multi-band photometry and medium/high resolution spectroscopy for a large sample of stars. The goals call for astrometry and photometry (5 broad and 11 medium bands covering the optical range) to be complete over the whole sky to $`V=20`$, and spectroscopy complete to $`V=17.5`$. The spectroscopy is primarily meant to provide the radial velocities and thus the 6th component in the phase-space to complement the other five provided by astrometry. Each star will be visited around a hundred times during the five year mission life-time, in a fashion similar to the highly successful Hipparcos scanning mode. The spectra, photometry and astrometry acquired at each transit epoch will be individually accessible. GAIA has three telescopes/fields of view, roughly equally spaced over a great circle on the sky, so that every $``$2 hours along the 6 hours spin period, the same star will cross the field of view of Astro-1 (where astrometric position and magnitude in 5 broad bands are measured), of Astro-2 (identical to Astro-1), and of the Spectro-Photometric telescope (where spectra and photometry in $``$11 medium bands are collected).

Hipparcos performed photometry in three bands. $`H_P`$, a sort of white light covering the Johnson’s $`B`$ and $`V`$ bands, has been the most accurate, with individual measurements accurate to 0.01 mag for solar type stars brighter than $`V=7.5`$. The other two bands, $`V_T`$ and $`B_T`$, were significantly less accurate. GAIA will carry instead a fully featured photometric system, composed of a set of 5 broad and $``$11 medium bands, which are still under fine optimization. Figure 1 illustrates the expected errors on the magnitudes derived at each passage in the GAIA fields of view for three sample broad bands, three sample medium bands and for the photometric reading accompanying the white light astrometric measurement. Each epoch measurement in any band and for any star brighter than $`V=12`$ (irrespective of the spectral type) will be accurate to better than 0.01 mag. In some bands (like in 2B56 or 2F59) the same limit is reached at $`V=14`$. In white light, all individual measurements of stars brighter than $`V=18.5`$ will be more accurate than 0.01 mag.

Field stars have colors corresponding to $``$G0 at $`V=10`$, $``$K0 at $`V=15`$ and cooler spectral types at fainter magnitudes. Therefore, optimization of spectroscopic observations has to be performed for G-K stars. Low metallicity stars of the Halo that covers a large fraction of the sky should be properly observed as well. Wavelength region satisfying these constraints is that of the CaII triplet in the far red (cf. Figure 2). This spectral region is interesting even for hotter stars, because it includes the head of the Paschen series and strong multiplets #1 and 8 of NI (for a detailed discussion see Munari 1999, 2002). Figure 3 illustrates the expected errors of GAIA radial velocities for single field stars of G-K spectral types and for a single transit. All stars brighter than $`V=14`$ will have epoch radial velocities accurate to 10 km sec-1 or better, all those brighter than $`V=11.6`$ to 1 km sec-1 or better.

## 2. GAIA and the eclipsing binaries

Previous reviews of SB2 eclipsing binaries as seen by GAIA that nicely complement the present one have been given by Zwitter (2002, 2003), Zwitter & Munari (2004) and Milone (2003).

Minimal velocity separation for a firm recovery of the SB2 nature for solar type stars is of the order of 25-30 km sec-1. Scaling Figure 3 to binary star case we see that radial velocity of each star in an SB2 binary with $`V13`$ and $`L_2/L_1`$0.5 will be derived with an error not exceeding 10 km sec-1. With $``$100 spectra per star randomly distributed in phase it may be expected that 1/3 of them will contribute to precise definition of semi-amplitudes $`K_1`$, $`K_2`$. If their accuracy is better than 10 km sec-1, the resulting $`K_1`$ and $`K_2`$ will be accurate to 1 km sec-1 or better. Figure 1 shows that at $`V=13`$ the epoch photometry in any of the 5+11 bands is accurate to 0.01 mag and in white light to an amazing 0.0006 mag for any of the considered stellar types. Therefore all SB2 eclipsing binaries with $`V13`$ will be accurately observed by GAIA.

It is important to check if GAIA will collect enough photometric observations to properly constrain the shape of the light-curve, and of the eclipse phases in particular. The measurements in white light and broad bands are essentially simultaneous and altogether $``$100 observations scattered over the five year mission lifetime will be collected (the exact number depends on the ecliptic coordinates and is generally comprised between 60 and 200). GAIA will also collect $`170`$ measurements in the medium bands (again nearly simultaneous in all bands). Moreover, RVS sky mapper measurements obtained just before spectroscopic observations will provide additional $``$100 photometric points in a further photometric band covering the 8480-8740 Å spectroscopic wavelength range. Therefore, combining the 2B80 broad band, the 2F89 medium band and the RVS sky mappers, an average of about 370 independent photometric measurements will be obtained at wavelengths close to the Cousins $`I_\mathrm{C}`$ band, and 270 measurements for all other bands.

Such a number of photometric points should suffice for most of the cases. We have embarked in a long term project to assess the accuracy of GAIA investigations of SB2 eclipsing binaries (Munari et al. 2001, Zwitter et al. 2003, Marrese et al. 2004, Milone et al. 2004 to be submitted). Twelve systems have been studied so far, combining ground-based radial velocities acquired in the GAIA wavelength range with Hipparcos $`H_P`$ photometry (to mimic GAIA photometric mapping). Typical light (and radial velocity) curves are presented in Figure 5, where orbital solutions are over-plotted. Such $`H_P`$ mapping is generally sufficient to constrain the radii to 2-6%, depending of the total number of points and the width of the eclipse phases. With GAIA providing 3– to 4–times more points, the fraction of SB2 eclipsing binary systems for which the orbital solution will have formal errors in the 1-2% range should be remarkable. Similar conclusions were reached by Niarchos and Manimanis (2003) from their investigation in a GAIA-like mode of four over-contact systems.

Hipparcos derived a frequency of 0.8% of eclipsing binaries among the stars it surveyed. Applying this to star counts in the Galaxy, we may estimate that GAIA will discover 64 000 eclipsing binaries with $`V13`$. Applying Carquillat et al. (1982) conclusion that 1/4 of EBs are SB2, we end up with about 16 000 SB2 eclipsing binaries brighter than $`V=13`$ for which GAIA should provide orbital solutions formally accurate to 1-2%. This is a fantastic number compared to $``$100 systems studied at similar accuracy by ground-based observations so far. These 16 000 stars will have their parallax measured to 5 micro-arcsec accuracy which corresponds to an error $``$0.5% in their distance if closer than 1 Kpc. Such extremely accurate parallaxes will be sed to directly constrain the values of several parameters (e.g. absolute luminosity) and physics used in orbital modeling, so that the atmospheric structure and global parameters of eclipsing binaries could be refined to absolute errors well below 1%.

## 3. RAVE and the SB2 binaries

RAVE (RAdial Velocity Experiment) is an ambitious program to conduct an all-sky survey (complete to $`V`$=15.5 mag) of radial velocities and metallicities of 50 million stars using the 1.2 m UK Schmidt Telescope of the Anglo-Australian Observatory (AAO), together with a northern counterpart, over the period 2006 – 2010. A key feature of the RAVE project is a pilot program aiming to observe $`10^5`$ stars using the existing 6dF facility in unscheduled bright time over the period 2003–2005. This pilot program begun on April 11, 2003 and by the time this conference was held it logged the observations of 20 000 stars. The 6dF spectrograph brackets the CaII triplet region adopted for GAIA (c. Figure 2). Spectra are taken at a resolution of $`R=8000`$.

The target list of the pilot program includes a large fraction of the $`\mathrm{118\hspace{0.17em}000}`$ Hipparcos stars that are accessible from the southern hemisphere as well as some of the $`\mathrm{2\hspace{0.17em}539\hspace{0.17em}913}`$ stars of the Tycho-2 catalog. The survey focuses on stars in the color range $`0.4<BV<0.8`$. For these stars useful photometric parallaxes can be derived if meaningful trigonometric parallaxes are not available.

The RAVE main survey should utilize a new Echidna-style multi fiber spectrograph at the UK-Schmidt telescope and its northern counterpart. In the UK-Schmidt planned version it consists of a 2250-spine fiber array covering the full field of view (40 deg2). The fiber positions are reconfigurable in $`5`$ min. A detailed description of RAVE and in particular of its science goals, aiming to derive the kinematics and formation history of the Milky Way, can be found in Steinmetz (2003).

RAVE observes any given target star only once, even if with remarkable S/N. Therefore it is not possible to discover SB1 binaries or trace the motion of components in SB2 systems. However, the respectable resolution and high S/N of RAVE spectra allow it to detect all SB2 binaries with a suitable luminosity ratio and velocity separation at the time of the observation. Given the huge number of stars surveyed by RAVE, this would put useful statistical constraints (on velocity separation and luminosity ratio as a function of spectral type and metallicity) on the population of binary stars in the Milky Way. Figure 7 presents the results of simulations carried to out to infer the SB2 detectability by RAVE. They show that all F–M SB2 binaries with a luminosity ratio $``$0.5 and observed when velocity separation of their components is $``$35 km sec-1 will be recognized as such. And new SB2 binaries are already being discovered in the early data secured by RAVE in 2003!

## References

Carquillat J.M., Nadal R., Ginestet N., Pedoussaut A. 1982, A&A 115, 23

Marrese, P.M. et al. 2004, A&A 413, 635

Milone, E.F. 2003, in ‘GAIA Spectroscopy, Science and Technology’, U.Munari ed., ASP Conf. Ser. 298, pag. 303

Munari, U. 1999, in ‘GAIA Workshop’, V.Straižys ed., Baltic Astron. 8, 73

Munari, U. 2001, in GAIA: A European Space Project”, O.Bienaymé and C.Turon. eds., EDP Sciences, pag. 39

Munari, U. et al. 2001, A&A 378, 477

Munari, U. et al. 2003, in ‘GAIA Spectroscopy, Science and Technology’, U.Munari ed., ASP Conf. Ser. 298, pag. 275

Niarchos, P.G., Manimanis, V.N. 2003, A&A 405, 263

Steinmetz, M., 2003, in ‘GAIA Spectroscopy, Science and Technology’, U.Munari ed., ASP Conf. Ser. 298, pag. 381

Zwitter, T. 2002, in IAU Coll. 187, ASP Conf. Ser. 279, pag. 31

Zwitter, T. 2003, in ‘GAIA Spectroscopy, Science and Technology’, U.Munari ed., ASP Conf. Ser. 298, pag. 329

Zwitter, T. et al. 2003, A&A 404, 333

Zwitter, T., Munari, U. 2004, in ‘Environments and evolution of double and multiple stars’, IAU Coll. 191, RevMexAA, C. Scarfe ed., in press (astro-ph/0306019)


text:
# Freezing and folding behavior in simple off-lattice heteropolymers

## I Introduction

Many biological heteropolymers (in particular proteins) can adopt a single, sequence dependent conformation (the fold, or native state) under physiological conditions. This specificity allows groups of atoms within such a molecule to be placed accurately in space, allowing the vast range of precise biochemical functions within living organisms. Upon altering the environment around the protein (e.g. pH, solvent quality, temperature, etc.), the molecule can be made to lose this single conformation (i.e. unfold or denature), adopting a fluctuating, disordered state. This *folding transition* has been shown to be both reversible and first-order for many proteins Anfinsen (1973), and is still poorly understood.

The modern theoretical viewpoint on the thermodynamics of protein folding (recently reviewed by Pande et al. Pande et al. (2000)) uses concepts from the theory of spin glasses. This approach was introduced by Bryngelson and Wolynes Bryngelson and Wolynes (1987), who used the random-energy model (REM) Derrida (1980) to describe the low temperature behavior of heteropolymers. This simple model assumes the energies of the possible conformations of a heteropolymer to be independent and uncorrelated, and predicts a second order transition (no latent heat) to a glassy state dominated by a low energy conformation. The more complicated theory of Shakhnovitch et al. Shakhnovich and Gutin (1989) uses another trick from the spin glass repertoire, the replica approach, to arrive at the same conclusions.

The low energy conformation occupied in the glassy state is not necessarily the native state. Depending upon details of conditions, sequence and sample history, the heteropolymer may well become “misfolded”; that is, locked in to a non-native, low energy conformation. Indeed, most randomly selected sequences should have several closely-spaced conformations at the low end of their energy spectrum, and will not reliably or reproducibly fold into only one of these. Some annealing of the heteropolymer *sequence* (i.e. design work) should be needed to “drive down” the energy of a chosen (target) conformation well below the closely-spaced region of the energy spectrum, preventing sampling of competing conformations and allowing robust folding Pande et al. (2000). Simulations Bastolla et al. (2000); Pande et al. (2000) and numerical studies Shakhnovich and Gutin (1990) using lattice heteropolymer models support these predictions and the approximations behind the models.

However, the random energy approximation is essentially mean field in nature, and cannot accurately treat systems with persistent correlations Pande et al. (1996). Further, lattice models artificially constrain the configurations available to a polymer. While lattice homopolymer results have been instructive at high temperatures Grosberg and Khokhlov (1994); Flory (1969), they do not predict the low temperature freezing transition seen in the simulations of Zhou et al. Zhou et al. (1996, 1997), where an isolated homopolymer freezes into a single, apparently unique crystalline state via a first-order transition.

We have performed Monte Carlo simulations using ten random sequence realizations of a simple off-lattice heteropolymer model to test these predictions in continuum systems. We report that, rather than exhibiting glassy behavior at low temperature, all ten sequences studied pass through a first-order transition to an ordered low energy state, dominated by a single conformation, without need for sequence design work.

## II Model

Our heteropolymers are modeled as freely-jointed chains of $`N`$ hard-sphere monomers of diameter $`\sigma `$, each assigned a type; we denote the type of monomer $`i`$ by $`s_i`$. The model is a simple extension of that used by Zhou et al. Zhou et al. (1996, 1997) and is shown schematically in Fig. 1.

Bonded monomers interact via the potential

$$u_{i,i+1}(r)=\{\begin{array}{ccccccc}\mathrm{}& & & & r& <& (1\delta )\sigma \\ 0& & (1\delta )\sigma & <& r& <& (1+\delta )\sigma \\ \mathrm{}& & (1+\delta )\sigma & <& r& & \end{array}$$
(1)
The parameter $`\delta `$ is a bond-length flexibility parameter; this allows for conversion to molecular dynamics simulation and has been shown to make no phenomenological difference to homopolymer folding results Zhou et al. (1996, 1997). We use $`\delta =0.1`$, allowing a ten percent variation in bond length.

If monomers $`i`$ and $`j`$ are not directly bonded (that is, $`\left|ij\right|>1`$), then they interact via a square-well potential:

$$u_{ij}(r)=\{\begin{array}{ccccccc}\mathrm{}& & & & r& <& \sigma \\ ϵ_{s_is_j}& & \sigma & <& r& <& \lambda \sigma \\ 0& & \lambda \sigma & <& r& & \end{array}$$
(2)
where $`r`$ denotes the distance between monomers $`i`$ and $`j`$, $`\lambda \sigma `$ is the square well width, and $`ϵ_{ss^{}}`$ is an *interaction matrix* giving square-well interaction strengths (square-well depths) between the monomer types. In this work, we use a modified hydrophobic-polar, or *mHP* matrix 111An unmodified hydrophobic-polar matrix has polar-type monomers interacting as hard spheres. This form for the interaction leads to P-type rich sections of the heteropolymer sequence forming long loops off of a central spherical nucleus rich in hydrophobic-type monomers Chertovich et al. (2002); this behavior is not seen in real proteins.. This is given by:

$$ϵ^{mHP}=ϵ^{}\left(\begin{array}{cc}1& \frac{1}{2}\\ \frac{1}{2}& 0\end{array}\right)$$
(3)
where $`ϵ^{}`$ is the characteristic energy scale. With this model, monomers with $`s_i=0`$ are hydrophobic-like (H-type), preferring type-symmetric interactions, while type $`s_i=1`$ monomers are polar-like (P-type), preferring type-asymmetric interactions. Similar interaction matrices are widely used in on-lattice protein folding simulations and numerical work Pande et al. (1997); Li et al. (1996); Lau and Dill (1990) and are believed to be good approximations to real intraprotein interactions Li et al. (1997). In the dense, globular heteropolymer state, we expect this interaction matrix to lead to configurations with a core of hydrophobic monomers surrounded by polar monomers, as seen in globular proteins.

We use reduced units scaled by the hard-sphere radius $`\sigma `$ and the energy scale $`ϵ^{}`$, such that

$$\begin{array}{ccc}r^{}& =& r/\sigma \\ E^{}& =& E/ϵ^{}\\ T^{}& =& k_BT/ϵ^{}\end{array}$$
(4)
Polypeptides of length greater than approximately 30 amino acids can form stable folded domains (e.g. the 35 amino acid WW domain Ilsley et al. (2002)). We have chosen to study $`N=64`$ length heteropolymers; this length is short enough for simulation in reasonable time, while long enough that equivalent real proteins can fold. We have studied ten random sequences interacting via the *m*HP interaction matrix, each with a $`1:1`$ H:P composition ratio. If we consider P-type monomers to map onto charged and polar amino acid residues, and H-type monomers to map to all other amino acid residues, this composition ratio is similar to that found in real proteins Fukuchi and Nishikawa (2001).

## III Simulation Method

Monte Carlo simulations were performed for a single polymer in isolation. As such, we have not used periodic boundary conditions. We have implemented crankshaft, pivot and continuum configurational bias (CCB) polymer moves Sadus (2002). Bond length fluctuation is allowed both in the CCB move implementation, as well as through standard Monte Carlo particle displacement moves performed on the monomers. CCB moves regrow end segments of up to four monomers. Moves are performed with relative probabilities in the ratio $`N:(1/2):(1/4):(1/4)`$ for particle displacement, crankshaft, CCB and pivot moves respectively. Initial equilibration runs last $`3\times 10^9`$ Monte Carlo move attempts; data collection runs last $`1\times 10^9`$ move attempts.

The Monte Carlo simulations also use parallel tempering (also known as replica exchange) Smit and Frenkel (2002). Several Monte Carlo simulations of the same system at differing temperatures are run in parallel; these simulations attempt to exchange configurations at regular intervals. This helps to maintain the ergodicity of the simulations; an individual system which has become stuck at low temperature in a metastable “trap” may be promoted to higher temperatures, where it can explore a wider range of phase space, aiding its equilibration. Our parallel tempering setup uses ten individual simulations across the temperature range $`T^{}=0.15\mathrm{}1.5`$, distributed such that the ratio between consecutive temperatures is constant Kofke (2002). Configuration swap moves are attempted every $`1000`$ Monte Carlo move attempts.

The simulations sample the configurational energy $`E`$, the radius of gyration $`r_g`$, the two-particle density $`\rho ^{(2)}(r)`$, and the instantaneous contact map $`𝐜`$ for the system every $`10^5`$ Monte Carlo move attempts. The two-particle density is defined as

$$\rho ^{(2)}(r)=\frac{1}{N}\underset{i}{}\underset{ji}{}\delta \left(r\text{-}r_{ij}\right)$$
(5)
where $`r_{ij}`$ is the distance between monomers $`i`$ and $`j`$ Hansen and McDonald (1986).

The instantaneous contact map $`𝐜`$ reflects a heteropolymer conformation; it is an $`N\times N`$ matrix, with each element $`c_{ij}`$ equal to one if monomers $`i`$ and $`j`$ are within each others square well, and zero otherwise. The instantaneous contact map is therefore given by:

$$c_{ij}=\{\begin{array}{ccc}0& & r_{ij}<\lambda \sigma \\ 1& & r_{ij}>\lambda \sigma \end{array}$$
(6)
The average contact map $`c_{ij}`$ gives us information on the structure adopted by the system at given temperature. We can use it to calculate a *normalized thermal average overlap*:

$$𝒬_n=\frac{{\displaystyle \underset{ij}{}}c_{ij}^2}{{\displaystyle \underset{ij}{}}c_{ij}}$$
(7)
This normalized overlap parameter can take values between zero and one. Values close to one indicate that the system is dominated by only a very small number of very similar conformations; small values indicate that the system is exploring a variety of different conformations.

Energy and radius of gyration data from the individual simulations within a parallel tempering run are tied together using self consistent multiple histogram extrapolation Ferrenberg and Swendsen (1989). This allows the histograms of energy $`E`$ and other observables $`\{A\}`$ recorded from the individual simulations to be combined to estimate $`\mathrm{\Omega }(E,\{A\})`$, the density of states of the system. From this, we can interpolate the behavior of the system across, and to a certain extent beyond, the simulated temperature and parameter range. We use the extrapolated distribution at given temperature to calculate the constant volume specific heat $`c_v^{}`$, given by

$$c_v^{}=\frac{E^2E^{}^2}{T^2}$$
(8)
Histogram extrapolation methods rely on good sampling across the range of interest by the individual simulations. As such, simulations around first-order transitions with a significant free energy barrier between phases can be problematic. If simulations on either side of the transition temperature sample only the phase which is most stable at that temperature, but no simulation actually observes the system passing back and forth between phases, the extrapolation process can “blur out” the transition, leaving a signature only in an increase in the error estimate and the response functions (e.g. specific heat) across the region of the transition. Worse still, if the system exhibits strong metastability, this can introduce a large systematic error into any estimate for the position of the transition.

To avoid this problem, we have used multicanonical sampling Berg and Neuhaus (1992) to directly confirm or deny phase coexistence where we suspect the existence of a first-order transition. We have used an *iterative multicanonical scheme* Smith and Bruce (1995). Repeated biased Monte Carlo simulations refine our estimate for the density of states of the system, with the aim of sampling evenly in energy across the suspected range of coexistence. The bias can then be removed from the results of these simulations, and the refined estimate for the density of states can be spliced back into the wider simulation results. The individual simulations used in this procedure are performed using the data collection run parameters stated above.

## IV Simulation Results

Plots of $`E^{}`$, $`c_v^{}`$ and $`r_g^{}`$ for our sequences, interpolated across the studied temperature range, are shown in Fig. 2.

The phenomenology does not vary considerably with sequence. Above a temperature of $`T^{}0.9`$, the heteropolymer exists in a disordered coil state - a high energy structure that rapidly increases its radius of gyration and decreases its heat capacity (energy fluctuations) with temperature. Between $`T^{}0.3`$ and $`T^{}0.9`$, we see the expected low temperature collapsed globule state with a small radius of gyration (order of $`0.05\sigma `$) and with relatively high heat capacity (order of $`2.0k_B`$). Around $`T^{}0.26`$, however, we see sharp spikes in the specific heat, associated with sudden drops in both conformational energy and radius of gyration with decreasing temperature, suggestive of first-order transitions. We have used multicanonical sampling to study these apparent transitions; for all sequences, we have found temperatures at which double peaked probability distributions occur for conformational energy, indicating two-phase coexistence. These probability distributions are shown in Fig. 3(a).

The sequence dependence of the probability distributions is small.

In Fig. 3 we plot coexistence temperatures for the sequences studied against the energy difference $`\mathrm{\Delta }E^{}/T^{}`$ between the coexisting phases. The quantities are clustered around $`T^{}0.27`$ and $`\mathrm{\Delta }E/k_BT56`$, but otherwise appear uncorrelated. For comparison, the typical latent heat $`\mathrm{\Delta }H`$ of the folding transition for a globular protein is of order of $`100k_BT`$ Cooper et al. (2001).

### IV.1 Structural character of the low energy phase

We present two-particle densities $`\rho ^{(2)}`$ data, calculated using Eq. (5), for sequence 4 in Fig. 4;

these plots are typical of the results for all ten sequences studied. These data were collected before the multicanonical algorithm was used, and may exhibit some degree of metastability. We first consider $`\rho _{HH}^{(2)}\left(r^{}\right),`$the hydrophobic-hydrophobic two-particle densities, in Fig. 4(a). We see that below the transition temperature, the first three neighbor shells are clear and distinct, separated by regions with near-zero hydrophobic particle density. The next three neighbor shells are much more sparsely populated (limited by the diameter of the collapsed state), but remain distinct, separated by sharp minima. These data show a low energy phase structure with local translational order between hydrophobic monomers. For $`\rho _{HH}^{(2)}\left(r^{}\right)`$ above the transition, the first three neighbor shells are still clear up to $`T^{}=0.7`$, but not distinctly separated. The next three neighbor shells have blurred together into one single peak. Above the transition temperature, the hydrophobic monomers have no translational order, as expected in a disordered globular phase. We show typical low-energy phase and globular phase configurations in Fig. 5;

the orientation of the low-energy configuration has been chosen to clearly show the hydrophobic monomers arranged in hexagonal crystal planes, whereas no crystal structure exists in the globular configuration.

Figures 4(b) and (c) show hydrophobic-polar and polar-polar two-particle densities $`\rho _{HP}^{(2)}\left(r^{}\right)`$ and $`\rho _{PP}^{(2)}\left(r^{}\right)`$, respectively. Neither show any sign of translational order, either above or below the transition. The polar monomers remain in a disordered state across the temperature range studied.

We show the thermal average contact maps in Fig. 6.

Note that these data are collected from equilibrated parallel tempering runs before the multicanonical algorithm is used, and as such exhibit metastability. At low temperatures (below $`T^{}0.32`$), we see that the contact maps are dominated by very persistent contacts, suggesting that the low temperature structure is dominated by very few configurations. Contact maps collected at higher temperatures show large low persistence regions, as expected for a disordered globular state.

To quantify this, plots of the normalized overlap parameter, calculated from these contact maps, are shown in Fig. 7.

As with the contact maps, these results exhibit metastability around the transition. For all sequences, we see a large increase on cooling in $`𝒬_n`$ from values of $`𝒬_n0.3`$ in the disordered globule phase to values of $`𝒬_n0.85`$ in the low energy phase. This confirms that the low energy phase is dominated by small fluctuations around a single conformation.

## V Discussion and Conclusions

Our results indicate that, for all of the ten random heteropolymer sequences which we have studied, there exists a low temperature transition between the familiar collapsed liquid-like globule phase and a low energy phase. These transitions do not have the glassy character expected from replica theory, showing instead double-peaked energy distribution functions with substantial associated latent heat. In each case, the low energy phase is more compact than the disordered globule, having a lower $`r_g`$, exhibits local translational order and crystal planes, as shown by the two-particle density function $`\rho ^{(2)}`$ (similar to the low energy phases seen in off-lattice homopolymer simulations Zhou et al. (1996, 1997)), and is not glassy, being dominated by a single conformation of the polymer backbone, as shown by the normalized thermal average overlap parameter $`𝒬_n`$. In analogy to the more familiar phase behavior of bulk systems, we describe the transition as a freezing transition, and the low-energy state as a frozen state. We note the obvious parallels with protein folding transitions, where polypeptide heteropolymers undergo a first-order phase transition to a collapsed single conformation, exhibiting repeating structural units.

The large computational effort associated with simulating high-density, low temperature continuum polymer states has limited the number of sequences we have been able to study. Nevertheless, we can use the Newcombe-Wilson score method (method 3 in Ref. Newcombe (1998)) to estimate with 95% confidence that *at least* 72.25% of 64 length 1:1 HP composition ratio sequences will demonstrate an equivalent first-order freezing transition to a non-glassy single conformation state.

We suggest that the reason that we see different behavior from equivalent lattice models, which undergo low temperature glass transitions, is a matter of symmetries. As evidenced by the $`\rho ^{(2)}`$ data (see Fig. 4), the disordered globular state for our sequences has no translational order, with hydrophobic monomers able to occupy any position, within the bounds of bond lengths and hard-sphere overlap. Passing below the freezing temperature, the frozen states of our sequences gain crystalline-like order 222We say crystalline-like rather than crystalline since it is possible that, as in case of two dimensional crystals, this translational order is not long-ranged, and in the large $`N`$ limit would be destroyed by fluctuations. However, this will be unimportant for polymers which freeze into objects with short length scales, which we note includes real proteins as well as the simple model heteropolymers studied here. ; hydrophobic monomers must occupy lattice points, or points very near them. The spatial symmetry of the system has been broken in a manner which, in three dimensions, requires a first-order phase transition. In the lattice model, the approximation is that spatial symmetry is already broken for the system at all temperatures, and there is no possibility of such an order-disorder transition. Further, the symmetry of the lattice is fixed in advance; if the frozen state of the off-lattice system is on a different lattice, or if some distortion of the lattice is necessary for stability, the lattice model cannot accurately portray the ordered frozen state. With cooling, the lattice system will remain in a disordered state, eventually forming a glass.

The crystalline nature of the frozen states seen here do not resemble the ordered structures seen in biological heteropolymers (i.e., proteins); we see no signature of large scale helices or $`\beta `$-sheets in the contact maps. We believe this to be due to the complete flexibility of the chains in the model. It has been suggested Maritan et al. (2000) that helices and $`\beta `$-sheet like “saddle” structures are optimal shapes for a stiff stringlike object, in the same manner that the ubiquitous face centered cubic crystal structure is the optimal arrangement for sets of spheres. We are currently investigating this matter for both homo- and heteropolymeric systems.

It can be seen from the contact maps (see Fig. 6) that the frozen configurations are strongly sequence dependent. However, the thermodynamic behavior of the heteropolymers (Fig. 2) has very little (though still statistically significant) sequence dependence. We believe that the thermodynamic behavior, in the HP model at least, has stronger composition dependence than sequence dependence. Our 1:1 H:P composition ratio (with its freezing temperatures bunched around $`T^{}0.27`$) lies between the limits of the pure hydrophobic homopolymer (with a freezing temperature of $`T^{}=0.33`$ Zhou et al. (1996)), and the pure polar homopolymer with its effectively hard-sphere monomers, which should not freeze in isolation.

It is generally suggested that the fraction of heteropolymeric sequences capable of freezing to unique conformations should be small, around $`10^9`$ Alberts et al. (2002). Dill and coworkers estimate the chance of finding a protein sequence which will fold to a given approximate structure as being of order $`10^{10}`$ Lau and Dill (1990); Chan and Dill (1991); Dill (1999). Recent experimental work, using a large random-sequence protein library to find sequences which bind to ATP, estimates the fraction of protein sequences which are “functional” as around $`10^{11}`$ Keefe and Szostak (2001). For our system, we have estimated that over 70% of 64-length 1:1 composition ratio HP sequences freeze to unique conformations. Similar over-estimation is also seen for lattice models; the impressive direct enumeration performed by Li et al. Li et al. (1996) on an HP lattice model suggests that as many as 4.75% of all 27-length HP sequences have unique ground state configurations.

We believe such over-estimation to be due to a lack of sources of frustration in our model, as compared to real protein systems. Our model has only two monomer types, no size polydispersity or side chains, completely flexible bonds, and no orientational dependence in its interaction potential. Any of these factors could significantly lower the number of model sequences which exhibit folding-like behavior.

In summary, we have demonstrated that a simple continuum heteropolymer model can exhibit a first-order transition to an ordered low temperature state dominated by a single conformation, without glassy behavior. This behavior is similar to that seen in protein folding. It would be of interest to see whether extensions to this model yield further protein-like character, for example secondary structure, or instead introduce frustrations which prevent the observed behaviors.

###### Acknowledgements.

This work is supported by the BBSRC (grant reference B17005).


text:
# Incomplete Photonic Bandgap as Inferred from the Speckle Pattern of Scattered Light Waves

## Abstract

Motivated by recent experiments on intensity correlations of the waves transmitted through disordered media, we demonstrate that the speckle pattern from disordered photonic crystal with incomplete band-gap represents a sensitive tool for determination the stop-band width. We establish the quantitative relation between this width and the angualar anisotropy of the intensity correlation function.

Introduction. A wave propagating in a random medium undergoes multiple scattering and forms a complicated intensity pattern, commonly referred to as a specle pattern. It is described in statistical terms, with the help of probability distributions and correlation functions. The intensity correlation function, $`\delta I(𝒓)\delta I(𝒓^{})`$, where $`\delta I(𝒓)`$ is the deviation from the averaged intensity at point $`𝒓`$, contains a short range term, $`C_1(𝒓,𝒓^{})`$, which oscillates on a scale of the wave length and exponentially decays beyond the mean free path, $`l`$. It also contains small long-range terms, which become dominant for $`|𝒓𝒓^{}|l`$. Although the theory of speckle patterns was developed some time ago, it is only recently that the first experimental measurements of the spatial correlator $`C_1(𝒓,𝒓^{})`$ were reported both for microwaves  and optical waves. These experiments were carried out for isotropic disordered media; macroscopic isotropy was also assumed in the existing theories. In this isotropic situation comparison with the theory has enabled the authors to infer the effective refractive index, which is the only relevant parameter of the medium in the absence of disorder.

The main message of this paper is that, in a medium with underlying spatial structure, the pattern of intensity correlations exhibits a vastly richer behavior as compared to the isotropic case. Moreover, additional features in $`C_1(𝒓,𝒓^{})`$ carry quantitavive information about this structure. As an example, we consider a disordered incomplete-bandgap photonic crystal and demonstrate how the band-structure parameters can be extracted from the angular anisotropy of the correlator $`C_1`$.

Superficially, it may seem that a wave with frequency $`\omega `$, arriving from a distant source, after having been scattered by many random inhomogeneities, will lose all information about the crystal band-structure. Our point, though, is that $`C_1(𝒓,𝒓^{})`$ is essentially a local object: It is determined by scattering events in the vicinity of the points $`𝒓`$, $`𝒓^{}`$, and it is not sensitive to the “prehistory”, i.e., scattering events experienced by the wave before arrival to the point $`𝒓`$. In other words, $`C_1(𝒓,𝒓^{})`$ can serve as a ”microscop” for observation of the local interference picture (on a scale smaller than $`l`$), thus, revealing the band-structure of the inherent photonic crystal. Such a microscop is particularly well suited for determination of the band-structure of realistic photonic crystals, since it is not affected by the long-range disorder. On the other hand, long-range disorder is a generic feature of realisitic crystals, such as silica-based synthetic opals.

Disorder in synthetic opals. Silica-based synthetic opals, which are the fcc self-assemled arrangement of almost monodispersed silica spheres, play a distinguished role in fabrication of photonic band gap materials, the potential of which for light manipulation was first appreciated in Refs.. This is because the opals constitute a template subsequently infiltrated with high refraction index material, while the spheres are selectively removed (see e.g. ). Due to this distinguished role, the vast majority of experimental studies of light propagation in photonic crystals has been carried out on opals.

In opal photonic crystals, the contrast of the dielectric constant is weak, so that photonic bandgaps are not only incomplete, but the corresponding stop-bands are relatively narrow. For this reason, the photonic band structure of this materials is strongly obscured by the disorder. As a result, any reliable determination of the stop-band either from disorder-broadened reflection maxima or from transmission minima, which are also srongly affected by the disorder, is highly ambiguous. In fact, in early studies, actual extraction of the stop-band width from the data was based on the “thumb rules”.

Later studies have yielded a deeper insight into the microscopic origin of the disorder in opals. Namely, they indicated that one should distinguish three types of the disorder

(i) Short-range disorder due to point defects and spread in the sphere diameters. This type of the disorder causes the mean free path $`l15\mu m`$.

(ii) Stacking faults in $`[111]`$ direction. This disorder has a one-dimensional character, and thus might give rise to the 1D in-gap states within the incomplete stop-band. In fact, experimentally observed deep and sharp transmission minimum near the stop-band center has been accounted for these states. Subsequent interpretation of reflection and transmission as well as diffraction data  also relied on the prominent role played by the stacking faults.

(iii) macroscopic domains that are $`50100\mu m`$ in size. Such domains inavoidably emerge in course of self-assembly of thick enough photonic crystals. Due to their presence the stop-bands are strongly inhomogeneously broadened – a serious complication for determination of photonic band structure.

Summarizing, experimental studies have led to significant improvement of the understanding of the light propagation in strongly-disordered incomplete-bandgap photonic crystals. However, no quantitative information about the stop-band width could be inferred from these experiments. Below we demonstrate that this width naturally emerges in the angular dependence of the intensity correlator $`C_1`$.

Correlator $`C_1`$ in a disordered photonic crystal. The short-range correlator $`C_1(𝒓,𝒓^{})`$, is obtained by factorizing the product of four fields $`\psi (𝒓)`$ in the average $`I(𝒓)I(𝒓^{})`$. This leads to the relation $`C_1(𝒓,𝒓^{})=|C_\psi (𝒓,𝒓^{})|^2`$, where

$$C_\psi (𝒓,𝒓^{})=\frac{\psi (𝒓)\psi ^{}(𝒓^{})}{I(𝒓)^{1/2}I(𝒓^{})^{1/2}}$$
(1)
is the field correlation function. It follows from the Bethe-Salpeter equation that

$$C_\psi (𝒓,𝒓^{})=\frac{4\pi }{l}𝑑𝒓_1G(𝒓,𝒓_1)G^{}(𝒓^{},𝒓_1),$$
(2)
where $`G(𝒓,𝒓_1)`$ is the disorder-averaged Green’s function. The physical meaning of Eq. (2) is that correlation between points $`𝒓`$, $`𝒓^{}`$ is established via scattering on an intermediate impurity, at $`𝒓_1`$. The integral Eq. (2) can be reduced to

$$C_\psi (𝒓,𝒓^{})=\frac{4\pi c^3}{ϵ^{3/2}\omega }\text{Im}G(𝒓,𝒓^{}),$$
(3)
which has the following intuitive explanation: field correlation between a pair of points is due primarily to waves that are scattered in the vicinity of one point and arrive to the vicinity of the other point. The average amplitude to arrive from $`𝒓^{}`$ to $`𝒓`$ is $`G(𝒓,𝒓^{})`$, so that $`C_\psi (𝒓,𝒓^{})`$ should be a linear combination of $`G(𝒓,𝒓^{})`$ and $`G^{}(𝒓,𝒓^{})`$. Eq. (3) is the right linear combination since it reduces to unity for $`𝒓=𝒓^{}`$. Using Eq. (3) we present the field correlator as

$$C_\psi (𝒓,𝒓^{})=\frac{4\pi c^3}{ϵ^{3/2}\omega }\text{Im}\underset{\mu ,𝒌}{}\frac{\psi _{\mu ,𝒌}(𝒓)\psi _{\mu ,𝒌}^{}(𝒓^{})}{\omega ^2\omega _{\mu ,𝒌}^2+i\eta },$$
(4)
where $`\eta =c\omega /ϵ^{1/2}l`$, with $`ϵ`$ standing for the background dielectric constant, and $`\psi _{\mu ,𝒌}^{}(𝒓)`$, $`\omega _{\mu ,𝒌}`$ are the photonic crystal eigenmodes and eigenfrequencies, respectively. They are characterized by the wave vector $`𝒌`$ of the first Brillouin zone and the band index $`\mu =\pm 1`$.

Suppose that stop-band corresponds to the direction of propagation along the $`z`$-axis. Then, the anisotropy of $`C_\psi (𝒓,𝒓^{})`$ is expected within a narrow angular interval, $`\theta =\mathrm{arctan}(z^{}z)/(\rho ^{}\rho )\gamma ^{1/2}`$, where $`\gamma =\delta \omega /\omega 1`$ is the dimensionless frequency width of the stop-band. For such $`\theta `$, the main contribution to the sum in Eq. (4) comes from the small domain of $`𝒌`$-space around the Bragg condition $`k_zQ_B/2`$, where $`2\pi /Q_B`$ is the period along $`z`$ (see Fig. 1). Within this domain the dispersion law, $`\omega _{\mu ,𝒌}`$, can be simplified as

$$\frac{ϵ\omega _{\mu ,𝒌}^2}{c^2}=k_{}^2+\frac{Q_B^2}{4}+\mu Q_B\sqrt{\frac{1}{L_B^2}+\left(k_z\frac{Q_B}{2}\right)^2},$$
(5)
where we have introduced the Bragg length, $`L_B=4\omega /Q_B\delta \omega `$, of the decay along $`z`$ exactly at the Bragg condition. The corresponding eigenmodes, $`\psi _{\mu ,𝒌}`$ can be also simplified

$$\psi _{\mu ,𝒌}=𝒰_\mu (z)\mathrm{exp}\left[i(k_zQ_B/2)z+i𝒌_{}𝝆\right],$$
(6)
where $`𝒰_\mu `$ are the Bloch functions. It is seen from Eq. (5) that, for a given $`\omega _{\mu ,𝒌}=\omega `$ and the same $`k_z`$, the values of $`k_{}`$ are very different for upper $`\mu =1`$ and lower $`\mu =1`$ bands. As a result, the contribution to $`C_\psi `$ from the upper band is exponentially smaller. Then, keeping only $`\mu =1`$ terms, we can express the field-field correlation function as

$$C_\psi (𝒓,𝒓^{})=\frac{c^3𝒰_1(z)𝒰_1(z^{})}{2\pi ^3ϵ^{3/2}\omega }𝑑k_zk_{}𝑑k_{}𝑑\varphi \text{Im}\left\{\frac{\mathrm{exp}\left[iz(k_zQ_B/2)+i\rho k_{}\mathrm{cos}\varphi \right]}{\omega ^2\omega _{1,𝒌}^2+i\eta }\right\}.$$
(7)
Substituting Eq. (5) into the integral Eq. (7), and performing the integration over $`k_{}`$ and $`\varphi `$ we arrive at

$$C_\psi (𝒓,𝒓^{})=\left(\frac{\gamma }{4}\right)𝒰_1(z)𝒰_1(z^{})\text{Re}\left\{𝑑q\mathrm{exp}\left[iq\frac{\gamma }{4}\right]H_0\left(\frac{\theta \gamma ^{1/2}}{2}\sqrt{\mathrm{\Delta }+\sqrt{1+q^2}}\right)\right\},$$
(8)
where $`H_0`$ is a Hankel function of a zero order. The dimensionless length $``$ and “detuning” $`\mathrm{\Delta }`$ are defined as

$$=|𝒓𝒓^{}|Q_B,\mathrm{\Delta }=\frac{2}{\delta \omega }\left(\omega \frac{cQ_B}{2ϵ^{1/2}}\right)+i\frac{L_B}{2l}.$$
(9)
For $`\theta \gamma ^{1/2}1`$, i.e. $`|𝒓𝒓^{}|L_B`$, the asymptotic expansion of the Hankel function can be used, so that the integral Eq. (8) takes the form

$$C_\psi (𝒓,𝒓^{})=𝒰_1(z)𝒰_1(z^{})\left[\frac{\gamma ^{3/4}}{2(\pi )^{1/2}}\right]\text{Re}\left\{e^{i\pi /4}𝑑q\frac{\mathrm{exp}\left[iF(\theta ,q)\right]}{\left(\mathrm{\Delta }+\sqrt{1+q^2}\right)^{1/4}}\right\},$$
(10)
where $`F(\theta ,q)=\gamma q/4+\theta (\gamma ^{1/2}/2)\left[\mathrm{\Delta }+\sqrt{1+q^2}\right]^{1/2}`$. Function $`F(\theta ,q)`$ has a special point at which the first and the second derivatives with respect to $`Q`$, for some value of $`\theta =\theta _c`$, are zero. The position of this critical point, $`(\theta _c,q_c)`$, depends on the detuning $`\mathrm{\Delta }`$ and can be parametrized by a complex number $`w`$ as

$`\mathrm{\Delta }={\displaystyle \frac{4}{27}}w^3w,\theta _c=2\gamma ^{1/2}\left({\displaystyle \frac{w}{3}}\right)^{3/2},q_c={\displaystyle \frac{1}{3}}\left(4w^29\right)^{1/2}`$For $`\mathrm{\Delta }=0`$ we have $`w=3^{3/2}/2`$, $`\theta _c=(3^{3/4}/2^{1/2})\gamma ^{1/2}`$, and $`q_c=2^{1/2}`$. To avaluate the integral Eq. (10) for $`\theta `$ close to the critical point, $`\theta \theta _c`$, the function $`F(\theta ,q)`$ should be expanded up to the third order in $`(qq_c)`$

$$F(\theta ,q)=\varphi _0+\frac{\delta \theta }{\theta _t}+\delta q\frac{\delta \theta }{\theta _d}\frac{1}{3}(\delta q)^3,$$
(11)
where $`\delta \theta =\theta \theta _c`$, $`\delta q=(9/4w)(6/w)^{1/3}(qq_c)`$, $`F_0=\gamma \left(4w^29\right)^{3/2}/108`$, and characteristic angles $`\theta _t`$ and $`\theta _d`$ are defined as

$$\theta _t=\frac{6}{\gamma ^{1/2}}\left(\frac{3}{4w^39w}\right)^{1/2},\theta _d=6\left(\frac{2w}{3\gamma }\right)^{1/6}.$$
(12)
Substituting Eq. (11) into Eq. (10) and performing integration over $`\delta q`$, we obtain

$$C_\psi (𝒓,𝒓^{})=\left[\frac{2\gamma ^{1/5}}{3}\right]^{5/6}𝒰_1(z)𝒰_1(z^{})\text{Re}\left\{\frac{\pi ^{1/2}w^{1/3}}{(4w^29)^{1/4}}\mathrm{exp}\left(i\varphi _0+i\frac{\delta \theta }{\theta _t}\right)\text{Ai}\left(^{2/3}\frac{\delta \theta }{\theta _d}\right)\right\},$$
(13)
where Ai is the Airy function.

Discussion. It is easy to see that Eq. (13) has a form similar to the field distribution near the caustics, which is the envelope of rays reflected by a curved surface. This similarity is not accidental. Indeed, near the caustics the light rays cross over from geometrically allowed directions to the “shadow”. Eq. (13) also describes a crossover from the directions $`\theta \theta _c\gamma ^{1/2}`$ of free propagation to the angular domain $`\theta \gamma ^{1/2}`$ where propagation is forbidden due to the stop-band, which plays the role of “shadow” in the $`𝒌`$-space. Similarly to the case of caustics, the behavior of $`C_\psi `$ near the crossover $`\theta =\theta _c`$ is described by the following exponential decay at $`\theta <\theta _c`$

$`C_\psi (\theta )\mathrm{exp}\left[\left({\displaystyle \frac{\theta _c\theta }{\theta _d}}\right)^{3/2}\right]\mathrm{cos}\left[{\displaystyle \frac{\theta _c\theta }{\theta _t}}+\varphi _0\right],`$which, in the case of intensity correlation serves as an envelope of the fast oscillations. For $`\theta >\theta _c`$ this envelope oscillates itself as $`\mathrm{cos}\left[\left(\frac{\theta \theta _c}{\theta _d}\right)^{3/2}\pi /4\right]`$. Overall, there are three angular scales in the correlator $`C_\psi `$, namely, the critical angle, $`\theta _c`$, the period of the envelope, $`\theta _d/^{2/3}`$, and the period, $`\theta _t/`$, of the oscillations. At spatial distances $`|𝒓𝒓^{}|L_B`$ all three scales are of the same order $`\gamma ^{1/2}`$. For large $`|𝒓𝒓^{}|`$ we have $`\theta _c\theta _d/^{2/3}\theta _t/`$. Numerical examples, illustrating separation of the three scales is given in Fig. 2. Outside the stop-band, $`\theta \gamma ^{1/2}`$ the correlator $`C_1`$ saturates at the value $`\mathrm{exp}\left(|𝒓𝒓^{}|/l\right)`$, governed by the mean free path. Obviously, our main result Eq. (13) applies in the vicinity of each of the equivalent Bragg directions.

When large-scale disorder (domains) broadens significantly the reflectivity peak, $`𝑹(\theta )`$, the dependence $`C_1(\theta )`$ with all fine details remains unchanged. Even if disorder is purely short-range, the intensity correlations are much more sensitive to the periodic background than the reflectivity. To illustrate this, in Fig. 3 we present the dependence $`C_1(\theta )`$ together with $`𝑹(\theta )`$, which near the Bragg condition is given by

$$𝑹=\left|\frac{ϵ^{1/2}1+\left(ϵ^{1/2}+1\right)\left[g(\theta )\sqrt{1+g^2(\theta )}\right]}{ϵ^{1/2}+1+\left(ϵ^{1/2}1\right)\left[g(\theta )\sqrt{1+g^2(\theta )}\right]}\right|^2,$$
(14)
where $`g(\theta )=L_B/2l+4i\theta ^2/ϵ\gamma `$. The dependencies $`𝑹(\theta )`$, $`C_1(\theta )`$ were calculated for the effective refractive index $`ϵ^{1/2}=1.7`$ and $`L_B=l`$. Similar to $`C_1`$, the width of $`𝑹(\theta )`$ is $`\theta \gamma ^{1/2}`$ However, $`𝑹(\theta )`$ is strongly smeared when $`L_B=l`$, whereas $`C_1(\theta )`$ contains oscillations. Three oscillations are well pronounced for the distances $`|𝒓𝒓^{}|`$ equal to $`2l`$, and $`3l`$, as shown in Fig. 3. Since these values are comparable to $`L_B`$, the argument of the Hankel function in Eq. (8) is $`1`$. Therefore, $`C_1=|C_\psi |^2`$ was calculated directly from (8) without using the asymptotical form Eq. (10). Therefore, even when $`L_Bl`$, comparison of the saturation values of $`C_1(\theta )`$ at two $`|𝒓𝒓^{}|`$ allows to determine the value of the mean free path. Then $`L_B`$ and, correspondingly, the stop-band width, $`\delta \omega =2c/ϵ^{1/2}L_B`$, can be inferred from the oscillations in the following way. Knowledge of the oscillation period, $`\pi \theta _t(w)/`$ together with frequency and mean free path, allows to deduce the parameter $`w`$, then the detunig $`\mathrm{\Delta }`$, defined by Eq. (9), and, finally, $`L_B=2l\left\{\text{Im}\mathrm{\Delta }\right\}`$.

Conclusion. In the present paper we have demonstrated that the intensity correlations in a disordered photonic crystal with incomplete bandgap uncover the underlying band structure that is obscured by the disorder, so that conventional approaches applicable to the “clean” photonic crystals fail.

It is important to emphasize that the very phenomenon of intensity correlation is a disorder induced effect and that the function $`C_1(𝒓,𝒓^{})`$ is calculated under the condition of diffusive propagation of the wave inside the sample. Note that, conceptually, with $`C_1`$ playing the role of a microscope, there is an analogy between the correlation of the fluctuations and the experimental setup, in which the source of light is located inside the crystal and the directionality of the outgoing light is studied. On the contrary, the long-range contributions to the intensity correlator come from scattering trajectories that explore all directions inside the crystal, and, due to isotropization, do not carry information about the Bragg directions.

Note in conclusion, that sensitivity of a number of disorder-induced effects, such as coherent backscattering  and directionality of the exiting diffusing light, to the periodic background has already been discussed in the literature. This sensitivity has also been employed  for the estimation of the stop-band width. In this regard, the intensity correlations constitute a much more accurate tool. Indeed, in the shape of the coherent backscattering cone, the properties of periodic structure enter only through the convolution with cooperon, which describes the large-scale diffusive motion, whereas the correlator $`C_1`$ comes from short distances. Another advantage of using intensity correlations over traditional methods is that $`C_1`$ is insensitive to the long-range disorder due to domains that mask the photonic band structure. Earlier it was demonstrated  that, isolating a single domain and measuring reflection from this domain, results in a drastic narrowing of the reflection spectra. We emphasize, that correlation analysis of the speckle pattern allows to get rid of inhomogeneous broadening in a natural way.

Finally, an attractive avenue for the future work comes from proposed infiltrating opals with nematic liquid crystals to manupulate the photonic band gap externally. While in reflectivity measurements reorientation of the nematic director manifests itself in a simple shift of the optical Bragg reflection peak, the structure of the correlator $`C_1`$ might reveal the effect of birefringence on the photonic bandstructure.

Acknowledgements. One of the authors (B.S.) acknowledges the hospitality of the University of Utah. The work was supported by the Army Research Office under grant No. DAAD 19-03-1-0290 and by the Petroleum Research Fund under grant No. 37890-AC6.


text:
# Untitled Document

COMMISSIONS 27 AND 42 OF THE IAU

INFORMATION BULLETIN ON VARIABLE STARS

Number 5xxx

Konkoly Observatory Budapest 00 Month 200x HU ISSN 0374 – 0676

RAPID VARIATIONS IN V2275 CYG (NOVA CYG 2001#2)

GARNAVICH, P.M.1,2, MACDONALD, A.J.1, WU, B.3, EASTERDAY, S.M.1, LIBAL, A.J.1, PALUMBO, A.1, QUINN, M.A.1

1 Department of Physics, University of Notre Dame, Notre Dame, IN 46566 2 email: pgarnavi@nd.edu 3 Department of Electrical Engineering, University of Notre Dame, Notre Dame, IN 46566

Nova Cygni 2001#2 was discovered by A. Tago and K. Matayama on Aug. 18 (Nakamura 2001) at a magnitude of 6.6 . The brightness decay from maximum was one of the fastest ever recorded and it shows characteristics of recurrent novae (Kiss et al. 2002). Time resolved photometry in 2002 Oct. by Balman et al. (2003) revealed large amplitude variations with a period of 8 or 11 hours which might be associated with its orbital period.

We observed V2275 Cyg beginning on HJD2452930.61 (2003 Oct 18) with the 1.8-m Vatican Advanced Technology Telescope (VATT). This was 790 days after maximum light. Observations continued the next night starting at HJD2452931.59. The CCD was binned 2$`\times `$2 providing a scale of 0.4 arcsec/pixel. $`V`$-band exposures were 120 seconds with 30 second readout time covering about 10 hours over the two nights. The resulting images were bias subtracted and flat-fielded and instrumental magnitudes were measured using aperture photometry. Three comparison stars near the nova (Figure 1) were also measured and their positions listed in Table 1. Approximate $`V`$-band magnitudes for the stars were estimated from zero-point and airmass coefficients measured at the VATT on earlier runs with standard magnitude errors estimated to be $`\pm 0.05`$ mag.

The VATT light curves are shown in Figure 2 and were obtained by subtracting the instrumental magnitude of Star 1 from the nova instrumental magnitude. The light curve for Star 2 relative to Star 1 is also shown and demonstrates both stars were constant over the observing run. Variations in individual measurements show an RMS scatter of 0.01 mag per exposure for $`V18`$ mag stars.

V2275 Cyg clearly shows light variations with a full amplitude of 0.2 mag and a period of about 20 minutes as well as a slower brightening trend on both nights. At the time of the observations V2275 Cyg varied between $`18.0<V<18.5`$ mag. The long-period variation suggests a period $`>7`$ hours and may be the same phenomenon seen by Balman et al. (2003). The short period variation has not been previously observed and may a quasi-periodic oscillation (QPO) or a stable periodicity such as a spinning white dwarf. Power-spectrum analysis gives a period of 0.410$`\pm 0.005`$ hours for the variation. A plot of this period against the normalized light curve is shown in Figure 3 and indicates that this period was stable over the two nights of data. This supports the possibility that the short-term light variations come from reprocessing of light from a asynchronous spinning white dwarf, but more data is needed to confirm the stability of the period.

We acknowledge travel assistance from the University of Notre Dame Department of Physics. This work based on observations with the VATT: the Alice P. Lennon Telescope and the Thomas J. Bannan Astrophysics Facility.

Table 1. Comparison Stars

| Star | RA (2000) | Dec (2000) | V mag |
| --- | --- | --- | --- |
| 1 | 21:03:04.80 | +48:45:51 | 15.50 |
| 2 | 21:03:01.96 | +48:45:37 | 17.75 |
| 3 | 21:03:04.38 | +48:45:31 | 15.33 |

References:

Balman, S., et al. 2003, IAU Circ, 8074

Kiss, L.L., et al. 2002, A&A, 384, 982

Nakamura, A. 2001, IAU Circ, 7686


text:
# References

THE UNITED FORCES IN THE NATURE OF MATTER

Rasulkhozha S. Sharafiddinov

Institute of Nuclear Physics, Uzbekistan Academy of Sciences, Tashkent, 702132 Ulugbek, Uzbekistan

The united rest mass and charge of a particle correspond to two form of the same regularity of the unified nature of its ultimate structure. Each of them contains the electric, weak, strong and the gravitational contributions. As a consequence, the force of an attraction among two neutrinos and force of their repulsion must be defined from point of view of any of existing types of the actions. Therefore, to understand the nature of the micro world interaction at the fundamental level one must use the fact that each of the four well known forces includes in self both a kind of the Newton and a kind of the Coulomb components. The opinion has been speaked out that the existence of the gravitational parts of the united rest mass and charge implies the availability of the fifth force in the nature of matter.

In studying the fundamental structure of matter such characteristics as the mass and charge play a large role. At the same time their nature remains thus far not finally established. Usually it is accepted that there is no any connection between the mass and charge. However, according to the hypothesis of field mass based on the classical model of a hard electron , a particle all the mass is strictly electric.

Our analysis of the behavior of massive neutrinos in a nucleus field shows  clearly that between the mass of a Dirac neutrino and its electric charge there exists an intimate interconnection. Such a sharp dependence giving the possibility to investigate the compound structure of charge quantization law  and reflects the fact that each of existing types of charges testifies in favor of the availability of a kind of the inertial mass . Thereby this mass - charge duality of matter explains the coexistence of the united rest mass and charge for the same particle. At the account of earlier findings, nonweak  and undiscovered properties of the neutrino they have the form

$$m_\nu ^U=m_\nu ^E+m_\nu ^W+m_\nu ^S+\mathrm{},$$
(1)
$$e_\nu ^U=e_\nu ^E+e_\nu ^W+e_\nu ^S+\mathrm{}.$$
(2)
Here the indices $`E,`$ $`W`$ and $`S`$ imply that both mass and charge of the neutrino contain as well as the electric, weak and strong components.

Exactly the same one can as the development of these sights include in the discussion the gravitational mass and charge. This procedure, however, meets with many problems. One of them states that the united rest mass of a particle contains the part that corresponds to its gravitational charge.

In conformity with such contributions of the gravitational mass $`m_\nu ^G`$ and charge $`e_\nu ^G,`$ we not only recognize  that

$$m_\nu ^U=m_\nu ^E+m_\nu ^W+m_\nu ^S+m_\nu ^G,$$
(3)
$$e_\nu ^U=e_\nu ^E+e_\nu ^W+e_\nu ^S+e_\nu ^G$$
(4)
but also need elucidate what neutrino united rest mass and charge say about the unified force in the nature of matter. The answer to this question one can obtain by studying the ideas of each of existing types of charges and masses. All they therefore will be illuminated in the present work.

According to our presentations about the structural properties of rest mass and charge, the force of gravity of the Newton $`F_{N_{\nu \nu }}`$ among two neutrinos and force of the Coulomb $`F_{C_{\nu \nu }}`$ between themselves may also be expressed from point of view of any of all possible types of the actions. In other words, each of these forces becomes the function as well as of corresponding components of the united masses or charges of interacting objects.

One can define their structure in the limit of the electric masses and charges as follows :

$$F_{N_{\nu \nu }}^E=G\left(\frac{m_\nu ^E}{R}\right)^2,F_{C_{\nu \nu }}^E=\frac{1}{4\pi ϵ_0}\left(\frac{e_\nu ^E}{R}\right)^2,$$
(5)
where $`G`$ is the constant of the gravitational action.

These forces for the weak masses and charges may have the form

$$F_{N_{\nu \nu }}^W=G\left(\frac{m_\nu ^W}{R}\right)^2,F_{C_{\nu \nu }}^W=\frac{1}{4\pi ϵ_0}\left(\frac{e_\nu ^W}{R}\right)^2.$$
(6)
At the availability of strong masses and charges they must behave as

$$F_{N_{\nu \nu }}^S=G\left(\frac{m_\nu ^S}{R}\right)^2,F_{C_{\nu \nu }}^S=\frac{1}{4\pi ϵ_0}\left(\frac{e_\nu ^S}{R}\right)^2.$$
(7)
The contributions explained by the gravitational masses and charges of interacting particles are written in the form

$$F_{N_{\nu \nu }}^G=G\left(\frac{m_\nu ^G}{R}\right)^2,F_{C_{\nu \nu }}^G=\frac{1}{4\pi ϵ_0}\left(\frac{e_\nu ^G}{R}\right)^2.$$
(8)
In a similar way the studyied forces may also be defined for the united masses and charges:

$$F_{N_{\nu \nu }}^U=G\left(\frac{m_\nu ^U}{R}\right)^2,F_{C_{\nu \nu }}^U=\frac{1}{4\pi ϵ_0}\left(\frac{e_\nu ^U}{R}\right)^2.$$
(9)
Comparison of (9) with the corresponding size from (5) - (8) leads us to the consequence of correspondence principle that

$$F_{N_{\nu \nu }}^U=F_{N_{\nu \nu }}^E+F_{N_{\nu \nu }}^W+F_{N_{\nu \nu }}^S+F_{N_{\nu \nu }}^G,$$
(10)
$$F_{C_{\nu \nu }}^U=F_{C_{\nu \nu }}^E+F_{C_{\nu \nu }}^W+F_{C_{\nu \nu }}^S+F_{C_{\nu \nu }}^G.$$
(11)
Insertion of (3) and (4) in (9) it would seem says of that equations (10) and (11) do not correspond to the reality at all. It is easy to observe, however, that this is not quite so. The point is that any particle with the united mass and charge come forwards in the system as the unified and the whole. Nobody is in force to separate its by part in the mass or charge type dependence.

Furthermore, if it turns out that just the compound structures of the united rest mass and charge establish the intraneutrino harmony of forces of the different nature , this can also confirm the fact that $`m_\nu ^U`$ and $`e_\nu ^U`$ are the multicomponent vectors, squares of which become equal to

$$\left|\stackrel{}{m_\nu ^U}\right|^2=\left(m_\nu ^E\right)^2+\left(m_\nu ^W\right)^2+\left(m_\nu ^S\right)^2+\left(m_\nu ^G\right)^2,$$

$$\left|\stackrel{}{e_\nu ^U}\right|^2=\left(e_\nu ^E\right)^2+\left(e_\nu ^W\right)^2+\left(e_\nu ^S\right)^2+\left(e_\nu ^G\right)^2$$

and that, consequently, the existence of solutions (10) and (11) is by no means excluded naturally.

The absence of one of forces $`F_{N_{\nu \nu }}`$ or $`F_{C_{\nu \nu }}`$ would imply that both do not exist at all . This becomes possible owing to the mass - charge duality of matter. In other words, $`F_{N_{\nu \nu }}`$ and $`F_{C_{\nu \nu }}`$ correspond to the most diverse form of the same action at the different distances.

According to the recent presentations about the nature of strong matter, the nuclear forces at a small distances have the character of the repulsion. In a large distances dependence can essentially appear their property of an attraction. Taking into account these facts and all what neutrino masses and charges say about structures of fundamental forces, we are led to the implication that each of them includes in self those parts which correspond to the masses and charges of interacting objects. Formulating more concretely, one can write the electric $`F_{\nu \nu }^E,`$ weak $`F_{\nu \nu }^W,`$ strong $`F_{\nu \nu }^S`$ and the gravitational $`F_{\nu \nu }^G`$ forces of the interaction between the particles in general form:

$$F_{\nu \nu }^E=F_{N_{\nu \nu }}^E+F_{C_{\nu \nu }}^E,$$
(12)
$$F_{\nu \nu }^W=F_{N_{\nu \nu }}^W+F_{C_{\nu \nu }}^W,$$
(13)
$$F_{\nu \nu }^S=F_{N_{\nu \nu }}^S+F_{C_{\nu \nu }}^S,$$
(14)
$$F_{\nu \nu }^G=F_{N_{\nu \nu }}^G+F_{C_{\nu \nu }}^G.$$
(15)
So, it is seen that any of the four well known forces contains both a kind of the Newton and a kind of the Coulomb parts. In this the unified regularity is said of the nature of these forces. Therefore, to create at the fundamental level a truly picture of the micro world interaction one must establish the compound structures of their naturally united gauge potentials.

If now the idea of any of forces (9) is accepted, it should be added that

$$F_{\nu \nu }^U=F_{N_{\nu \nu }}^U+F_{C_{\nu \nu }}^U.$$
(16)
Using (10), (11) and taking (12) - (15), we find that

$$F_{\nu \nu }^U=F_{\nu \nu }^E+F_{\nu \nu }^W+F_{\nu \nu }^S+F_{\nu \nu }^G.$$
(17)
Thus, it follows that if the gravitational rest mass and charge of a particle are unequal to its all the mass and charge, this will indicate to the existence of the fifth force of the interaction which come forwards in the nature as the unified and the united.