Simulation of relaxed bubble collision
                                                                 Raz Carmi

The project is written for MATLAB 5.2 .
In order to run the simulation write:  start

The files are:  start.m  start.mat  simulate.m  picabst.m  picabst.mat  picsett.m  picsett.mat  picresu.m  picresu.mat  picalgo.m  picalgo.mat  texabst.m  texsett.m  texresu.m  texalgo.m

Spontaneous formation of topological defects in non-equilibrium phase transitions was extensively studied in the past few years. This issue is relevant to the early universe theories, superfluids 4He and 3He, liquid crystals, superconductors and Bose-Einstein condensation. A theory of bubble collision of domains having different values of the phase belonging to the complex order parameter can help to predict the net spontaneous flux which should be measured in a rapid superconducting phase transition. The simulation is based on a 2D hexagonal lattice where each point represents a domain of order parameter with random phase and constant amplitude. The matrix element (a point in the lattice) has a real value (the phase). Between each three points a quantized vortex, antivortex or homogeneous site may be created according to the three random phases and a condition called ''the relaxed geodesic rule'' (described by an appropriate function with independent random variable). The geodesic rule is the assumption of minimal gradient and the ''relaxing'' is the probability for deviation from this rule. The program should find the random array of vortices and antivortices, and the net flux (which is a measurable quantity in a real experiment). The relaxing strength can be changed in different simulations and the influence on the net flux and on the total vortex density can be investigated.

Parameter setting:
There are four setting options:
1.  Matrix dimension: determines the size of the square matrix for the simulation.
     A real superconductor sample is typically analog to a matrix of the order of
     magnitude of 1000 to 10,000 square but it also depends on the rapidness of the
     quench. A simulation with such a matrix will execute for a long while.
     Here, we are limited to smaller systems in order to show the behavior of the results
     in a reasonable period of time.

2. No. of quenches: determines the number of independent simulation runs which are
     included in the statistics. Each quench is analog to a cooldoun of the system from
     the normal to the superconducting state and then counting the spontaneous
     vortices. Where this number is higher the accuracy of the statistics is better, but it
     takes longer to complete the simulation.

3.  Max. relaxing and No. of relaxing points: determine the maximal relaxing strength
     and the number of the relaxing points which will appear in the comparison
     between the net and the total flux. The scale for the relaxing strength is chosen to
     be equal to 1 for "very high" relaxing. The first result to appear in the graph
     (relaxing=0) is always the reference point: net flux = total flux =1 (as a relative

The main result of the simulations is that the increasing of the net flux (vortices minus antivortices) is much sensitive to the relaxed geodesic rule than the increasing of the total flux. The maximal theoretical statistical value of the net flux is the square root of the total flux (similar to the standard deviation in a random walk). Here we can see that even when the relaxing strength is not too high the net flux is close to about 1/10 of the maximal value or even higher (providing the matrix dimension is not too small). This result has an experimental implication for superconductors.
The graphs show:
1.  A sample of the vortex array (40x40). When the relaxing strength start growing the density is higher and
     more double vortices/antivortices appear. Too many of these objects is not physically reasonable since
     they are unstable in a real system.
2. Distribution Histograms of the net and the total flux, containing the number of quenches defined in the
    setting. The average of the net flux is always around zero. The distribution of the total flux is very close to
    its average.
3. A comparison between the net flux and the total flux, both relative to their values without relaxing. Here,
    the net flux is actually the standard deviation of the distribution (red). The total flux is the average of the
    distribution (blue).

The important aspect of the algorithm is how to define properly the relaxed geodesic rule, while maintaining the correct physical meaning. The main steps are:
1. Each matrix element of the square array has a random value between 0 to 2pi, represents the phase of the
    complex order parameter.
2. The interactions between the matrix elements are based on a structure of hexagonal lattice, because a
    vortex may be created between 3 different phase values. That is, if the sum of the phase differences
    equals +-2pi.
3. In the simple geodesic rule the phase difference is always minimized. If the geodesic rule is relaxed, there
    is some probability, for example, that a phase difference of pi/2 will be considered as -3pi/2.
4. Here, the program calculates the sum of the 3 phase difference values around each contact point in the
    lattice and checks if it is equal to 0 (homogeneous site), +-1 (vortex or antivortex), +-2 or more (double
    vortex/antivortex). Before the summation each phase difference is "relaxed" with a probability determined
    by the relaxing strength parameter.
5. For a physical meaning a non-symmetric probability is declared (in some mathematical steps). For
    example, if there is originally a vortex between three phase regions, the probability for homogeneous site
    is higher than for a double vortex (in an appropriate way), since in a real system a double vortex is more
    costly in energy.