Mykyta Chubynsky – Research

Past research (mostly pre-2007)

(back to recent research)

Rigidity and self-organization of elastic networks
Methods for accelerated computer simulation of materials
Other research

Rigidity and self-organization of elastic networks

(back to the top)

My main research topic both in graduate school and during my first two postdocs has been study of rigidity of elastic networks using methods of rigidity theory. Elastic networks can be thought of as networks of elastic springs connecting a set of sites. Mathematical rigidity theory deals with certain properties of these networks, such as the number of zero-frequency motions, sets of sites that do not deform in such motions (rigid clusters), stressed regions, etc., that depend only on the network geometry and often just on the topology, whereas the particular values of the spring constants and lengths are unimportant. This latter fact enabled the design of a very fast algorithm for analyzing network rigidity, the pebble game, developed by Don Jacobs, that essentially just does counting of elastic constraints on all length scales simultaneously. This algorithm was used, in particular, to study the phenomenon of rigidity percolation on diluted random elastic networks: as the constraint concentration increases, at some point a percolating rigid cluster emerges. My work in the area has been done in collaboration with M. F. Thorpe, N. Mousseau, J. C. Phillips, D. J. Jacobs, A. J. Rader and M.-A. Brière and has involved several subtopics:

Rigidity percolation on Bethe lattices Often, various statistical physics models can be solved analytically on Bethe lattices, which are networks with sites connected completely at random, regardless of the distance, and thus have no finite rings in the thermodynamic limit. This is the case for rigidity percolation as well. Following previous work in 2D (i.e., for networks with sites having two degrees of freedom), we have solved the rigidity percolation problem for 3D networks with both first- and second-neighbour constraints (bond-bending networks) which are interesting as models of covalent glasses. Unlike in both usual (connectivity) percolation on Bethe lattices and rigidity percolation on regular 3D bond-bending lattices, the rigidity transition on Bethe lattices is first order. We have obtained the phase diagram of the problem as a function of both the bond concentration and the "chemical order" parameter characterizing the propensity of sites with different number of neighbours (modeling different valences of atoms in a glass network) to bind together. See publication 3.

Self-organization, self-organized criticality and the intermediate phase in elastic networks The simplest percolation models involve studying randomly diluted networks where one starts from a regular lattice and removes bonds completely at random. But there is also interest in correlated percolation models. Many studies of rigidity percolation are inspired by covalent glass networks, and these networks are, of course, not totally random. Our self-organization models try to capture one aspect of this non-randomness, the tendency to avoid excessive mechanical stress. Being based on the pebble game constraint counting algorithm, these models are necessarily simplistic, but we believe they capture some important features of covalent glasses; indeed, there has been some experimental confirmation of our main qualitative results. The idea is to start with an empty lattice and then add bonds one by one as usual, but testing each bond for whether it would cause stress in the network and rejecting those that do, until introduction of stress becomes inevitable, after which bond insertion continues completely at random. It turns out that in this model, there are two separate percolation transitions, a rigidity transition and a stress transition, with a rigid but stress-free intermediate phase between them. More recently, recognizing that the original procedure is essentially an aggregation process, as bonds are only inserted and never removed, and thus is biased towards certain networks, we have introduced an additional equilibration step, consisting of repeated removal and reinsertion of bonds after every bond addition, maintaining the stress-free character of networks. If the equilibration is long enough, the result is an unbiased distribution of stress-free networks; this corresponds to equilibrium in the T → 0 limit for any statistical-mechanical model in which there is an energy cost associated with stress. In this variant of the model, the intermediate phase is still present, but it is unusual in that at any point within the intermediate phase, the probability that there is rigidity percolation is non-trivial (between 0 and 1). The distribution of finite cluster sizes is always power-law within the whole intermediate phase; also, percolation in the intermediate phase can often be created (destroyed) by a single bond addition (removal). This suggests that the system is in the critical state on the verge of percolation throughout the intermediate phase, an attribute of self-organized criticality, but in an equilibrium system. We have also calculated the bond-configurational entropy cost of self-organization using a Monte Carlo procedure that somewhat resembles thermodynamic integration; this cost is just a few percent of the total bond-configurational entropy, thus confirming that self-organization effects are expected to be observed in real systems close to the glass transition temperature. As mentioned, there is indeed some experimental evidence of these effects. See publications 4, 7, 10, 15, 18. More recently, I have also extended the model with equilibration above the upper boundary of the intermediate phase, by equilibrating while keeping the number of stress-causing constraints to a minimum, and also considered a finite-temperature analog of the model. These results are still unpublished.

Exactly mean-field conductivity in a correlated resistor network model As an offshoot of our work on self-organization in rigidity percolation, we have considered its direct analog in usual, connectivity percolation, the loopless percolation model. This model has been studied extensively before, as it is the q → 0 limit of the q-state Potts model; the novelty of our work is that by analogy with our rigidity model, we also consider what happens when bonds continue to be inserted at random after the network can no longer remain loopless (i.e., above the spanning tree limit). We have been able to first show numerically and then prove that when links in such a network are replaced by resistances, the conductivity changes exactly linearly as a function of bond concentration, just like in mean-field theory. The proof is based on the Kirchhoff theorem giving an expression for the conductance of a resistor network in terms of the number of trees that can be built on it. See publication 12.

Algorithms for rigidity analysis of general 3D networks While the theory underlying the pebble game algorithm is valid for 3D bond-bending networks, it is invalid in general for 3D elastic networks that are not bond-bending. But the question of practical interest is how frequent the errors would be in a pebble-game-type algorithm based on this theory. We have constructed such an algorithm, as well as a slower but exact algorithm based on network relaxation that we use to test the new pebble game. It turns out that at least in some cases of interest, errors in the new pebble game are extremely rare and thus the algorithm can be regarded as operationally exact in these cases. However, we have also developed a procedure that generates systematically the rare counterexamples and thus makes it possible to estimate their frequency. See publication 19.

Rigidity percolation transition in diluted 3D central-force networks One case when the pebble game for 3D networks is essentially exact is randomly diluted central-force regular lattices. Applying the pebble game to this case, we have found that the percolation transition is first order on the BCC and FCC bond-diluted networks and the BCC site-diluted networks, but second order on the FCC site-diluted network, contrary to what one might expect based on universality arguments. The results for first order transitions are especially spectacular: there are only very small rigid clusters (≲ 10 sites) below the transition, and then, upon a single bond addition, a huge percolating cluster taking up more than 90% of the network emerges. While the evidence we have is rather strong, the non-universality claim is certainly extraordinary enough that our results should be obtained by other means in the future. Exotic alternative explanations are possible. See publication 19.

Rigidity percolation in chemically ordered networks This is yet another case when the rigidity transition is first order. Take a bond-bending 3D network all sites of which have 3 neighbours (are 3-coordinated), but which otherwise is arbitrary. Then start breaking up bonds into two adding a site in the middle of the bond, thus introducing 2-coordinated sites. Never put two 2-coordinated sites next to each other, until that becomes necessary, after which never put three 2-coordinated sites next to each other, etc. In this way, maximum "chemical order" is achieved. We have proved that in this case, the whole network initially consists of a single rigid cluster, up until there are just enough 2-coordinated sites to always have a 2-coordinated site next to a 3-coordinated site and vice versa (at this point the fraction of 2-coordinated sites is 60%). This single cluster situation persists until 6 more sites are added, but when one more site is added, the huge rigid cluster breaks up into as many tiny clusters as there are sites in the network. This is thus the sharpest rigidity transition imaginable. See publication 11.

(back to the top)


Methods for accelerated computer simulation of materials

(back to the top)

In many condensed matter and molecular systems, motions of constituent atoms span a wide range of time scales. For instance, in protein molecules, the time scale of the folding process exceeds the fastest time scale of vibrational motions of atoms by at least 6 orders of magnitude and often more. Straightforward molecular dynamics (MD) simulations run into trouble, as the time step should correspond to the fastest time scale, but the total simulation time should be as long as the longest time. One possible approach is based on the fact that the configuration space of these systems can be separated into basins of attraction to potential energy minima, and the dynamics consists of high-frequency vibrations within these basins interrupted by much more rare jumps between them. One could then coarse-grain the dynamics by concentrating on the interbasin jumps (or activated events). This would still reproduce the correct long-time dynamics, provided that the relative frequencies of possible jumps are reproduced correctly. There have been many methods developed mostly by the computational chemistry community that allow one to find the pathways and then calculate the rates of the transitions. Most of these methods assume that the initial and the final states are known at least roughly. In many situations, this is not the case: for a given initial state, many different jumps are possible and the final states are not known a priori. Barkema and Mousseau have developed a method, the activation-relaxation technique (ART), for generating different possible jumps. In this method, the system starts at a minimum and gets displaced in a random direction until reaching a certain threshold (the "basin boundary") defined in terms of the lowest eigenvalue of the Hessian; once the threshold is reached, the system moves in the direction of the eigenvector corresponding to this eigenvalue relaxing in the perpendicular direction; this will often bring the system in the vicinity of a saddle point, thus completing the activation part, and then the relaxation brings the system to another minimum. Coupled with methods for calculating the rates, ART should allow, at least in principle, finding the rates of all jumps, executing them with appropriate relative frequencies and thus reproducing the dynamics. However, this approach requires the generation of a reasonably complete list of possible transitions for each visited minimum, which is computationally expensive. The problem arises from the fact that the probabilities of generating certain events in the ART method and other similar methods are not known in advance. My work has involved some attempts towards solving this problem and has been done in collaboration with N. Mousseau, G. T. Barkema and H. Vocks.

Combining intrabasin MD with ART: POP-ART The properly obeying probability ART (POP-ART) method modifies ART in such a way that the transitions between basins occur with relative probabilities obeying the detailed balance condition appropriate for the given temperature. This guarantees the correct distribution over different basins. The simulation is carried out at a constant energy (i.e., in the microcanonical ensemble), although coupling with a thermostat is possible. First of all, the random push within the initial basin is replaced by an MD run — this guarantees that points on the basin boundary (defined as in the original ART) are sampled with the correct probabilities. Once the boundary crossing occurs, the activation is carried out. This activation step is similar to that in ART, but it is continued until another basin boundary is reached and modified so as to ensure reversibility. If the activation step is viewed as a coordinate transformation, the detailed balance condition requires that the Jacobian of this transformation be equal to 1, or else an additional acceptance/rejection step is required. The Jacobian can be calculated on the fly during the activation. We have calculated the Jacobian for a model potential showing that when it is averaged over possible trajectories, it corresponds to the free energy change between the basins, having both energetic and entropic contributions. We have tested the method comparing it to MD for crystalline Si with an interstitial using the Stillinger-Weber potential showing that different interstitial configurations are reproduced with correct probabilities; simulations for a vacancy diffusion showed that the method can be much more efficient than MD at low enough temperatures. See publication 14.

Using memory of transitions in event-based simulations Even when the dynamics is coarse-grained by eliminating vibrations and concentrating on activated events, the sampling of the configuration space may still be inefficient. This is because in complex systems, there are often groups of basins, often called metabasins, such that transitions between basins within the metabasin are relatively frequent, but transitions to outside the metabasin are much more rare. The system will then visit the same states over and over again before moving elsewhere. For optimization problems, when the only purpose is finding the global energy minimum, there is a well-known approach, tabu search, in which a memory of recently visited states is introduced and returning to these states is prohibited for some time afterwards. In general, this is not a suitable approach when not just the lowest-energy state needs to be found, but reproducing the correct distribution over the states at a given temperature is required. However, using a toy model, we have found that if the memory of states is replaced by the memory of transitions, the obtained distribution is nearly exact. See publication 16.

(back to the top)


Other research

(back to the top)

Over the years, I have also been involved in research on other topics. While in some cases this has not resulted in anything publishable, this has been a useful experience, as it has given me a chance to familiarize myself with very different areas of physics.

Non-critical phase matchings in nonlinear optics This was one of my undergraduate research projects, under supervision of N. E. Kornienko from the Kiev Shevchenko University. Nonlinear processes in optics, such as the second harmonic generation or the sum and difference frequency generation, are only efficient when the waves involved phase-match so that constructive interference is achieved everywhere along the path. This would be difficult to achieve in an isotropic medium with dispersion, but is possible in an anisotropic crystal, when sending the light in a particular direction and using waves of different polarizations. In some cases, the matchings are non-critical, meaning that the matching condition is approximately satisfied over a wide range of frequencies and/or angles, which is desirable for some applications. My work involved calculating conditions of such matchings for a particular material. See publication 1.

Cluster expansion methods for electronic structure of substitutional alloys Electronic structure calculations for periodic solids utilize the periodicity by using Bloch's theorem and considering a single unit cell. For amorphous solids, this is, of course, impossible. In the intermediate case of a substitutional alloy having a periodic lattice but with a disorder in the occupancy of the lattice sites, there is still no perfect periodicity, but one can consider a periodic solid with parameters found in a self-consistent way as the zeroth approximation. The standard approach, the coherent potential approximation (CPA), uses the Green's function formalism and finds the self-consistent potential by requiring the vanishing averaged single-site scattering matrix. One can then improve the result by developing cluster expansions — one particular way is described in publication 2. As my undergraduate project, under supervision of S. P. Repetsky (the page is in Ukrainian) from the Kiev Shevchenko University, I have applied the CPA to a toy model with a continuous, rather than discrete, set of possible on-site parameters, and also have made an attempt to develop a self-consistent cluster expansion for this case.

Linear conversion of waves at a plasma boundary In a plasma in a magnetic field, several modes of electromagnetic waves are possible, with dispersion relations depending on both the angle between the direction of propagation and the magnetic field and the polarization. Just like in a birefringent crystal an ordinary wave can convert partially into an extraordinary wave upon reflection from the surface, conversion between different modes is possible at the plasma boundary. As yet another undergraduate project, under supervision of K. P. Shamrai from the Kiev Institute for Nuclear Research, I have considered analytically the conversion of helicon waves at the plasma boundary in the limits of both a very sharp and a very fuzzy boundary.

Stochastic models of population dynamics During my stay at Arizona State University, I had a chance to interact with T. J. Newman (now at Dundee), whose research is in the area of applications of statistical physics to biological systems and especially the effect of fluctuations in population dynamics, biochemical reactions, etc. I was involved, in particular, in the study of properties of the stochastic logistic model of population dynamics. Unfortunately, it turned out that most results obtained by me were obtained before, so this research has not resulted in any publications.

Kinetic Monte Carlo simulations of pattern formation due to a standing wave While this was not a research project that I can call my own, I participated in it extensively by aiding a graduate student in developing a kinetic Monte Carlo code for simulating pattern formation on a surface irradiated by a standing electromagnetic wave due to both inhomogeneous heating of the surface and steering of atoms being deposited on the surface by the electromagnetic field.

(back to the top)