MONTE CARLO SIMULATIONS

          Simulations of the DHS were carried out for a fixed number of particles per unit volume (r* - the number density) and a fixed ratio, m*2,  between the dipolar energy and the thermal energy (m*2 = m2/(s3kBT)).  Technically, this means that the simulations were carried out in the canonical ensemble (or NVT). 

          Given a pair of values (r*, m*) one chooses the number of dipolar particles N, and places  them in a cube (square) of volume (area) equal to N/r*.  This initial configuration can be, e.g.,  random (random positions for the spheres and random orientations for the dipoles ) or ordered (e.g. spheres on a regular lattice and dipoles oriented in the same direction). The final results do not depend on this choice. After this, the Monte Carlo simulation is started and proceedes as follows:

        

 (1)  Calculate the energy E (dipolar energy) of the particles in the initial configuration. Set i  = 1.  

(2)  Choose  a particle i  and randomly generate a new position and orientation for this particle. Calculate the energy E’ of the new configuration.

(3) If the energy of the new configuration is lower than the energy  of the initial configuration (E’ < E), then change the initial configuration to the new one. This means that  new configurations that lower the energy are accepted. Set E=E’ and go to (5).

(4) If E’ > E, generate a random number x with a value between 0 and 1. Calculate the number y = exp(-(E’-E)/ kBT).  If y>x then change the initial configuration to the new one and set E=E’. Otherwise, reject the new configuration.  Since E/ kBT is proportional to m*2 , this means that the probability of accepting new configurations with higher  energy is smaller for systems where the dipolar energy dominates over the thermal energy . Go to (5)

(5) If i is equal to N then set i=1. Otherwise, set i = i + 1. Go to (2).

          After several cycles, the system reaches thermodynamic equilibrium: the mean value of its state variables that are not fixed a priori remains constant. Then the simulation continues and  starts exploring  different configurations of the system (microscopic states) that correspond to that equilibrium state (macroscopic state).  This allows the calculation of  thermal averages of relevant quantities. Technically, a simulation makes a direct estimate of the canonical partition function of the system, since it samples the most probable configurations for a given set of (NVT).

          A good quality simulation of a sufficiently large DHS system for a sufficiently long time (number of cycles) is a very demanding task, only possibe if care is taken in calculating the energy (inclusion of long range corrections) and various technical tricks (bias, cluster moves, use of tables to calculate functions, etc.) are used.

In the following figure, we represent typical  configurations (snapshots of the system in equilibrium) obtained for several sets of r* and m*  and for N=5776 in a 2-dimensions. 

Each point in the boxes represents a particle and one immediately sees that they have a strong tendency to form clusters. Most of the particles are aggregated linearly but some belong to branched clusters. Each dipole (not represented in the figures) shows, as expected, a strong tendency to align (“head-to-tail”) with its neighbors.  At the end of a simulation for a pair r*, m*  we keep the positions and orientations of the N dipoles of 100 to 300 configurations, i.e. “pictures” like the ones showed in the last figure. 

In the next section we explain how the structure of the DHS is analysed.