ABSTRACT

We report on molecular dynamics simulations performed using microcanonical ensemble to predict the melting of argon particles in nanometer size range 10 nm and to investigate the effect of the time step integration. We use the LennardJones potential functions to describe the inter-atomic interactions, and the results are evaluated by using caloric curves of the melting phenomenon. Thermodynamic properties, including the total energy, Lindemann parameter, kinetic and potential distribution’s functions, are used to characterize the melting process. The data shows bimodal behavior only in a certain interval of integration time step Δt, while the internal energy increases monotonically with the temperature. For the other time step values, the back bending disappears. We claim that negative specific heat is related to a possible decrease of entropy in an isolated system; this can be interpreted as a result of the internal interactions, especially attractive process and specific relaxation time.

1. Introduction

One of the great challenges in cluster physics in recent years has been the identification and the characterization of critical behaviour and phase transitions, including solidto liquid and liquid-to-gas phase transitions. Since clusters are particles of finite size, one is confronted with the general question of how to detect and/or characterize such a transition in a finite system, a question of interest for many microscopic or mesoscopic systems such as, for instance, melting and vaporization of metallic clusters, Bose condensation of quantum fluids, and nuclear liquidto-gas transition [1-5].

Haberland and co-workers [6] reported the first experimental determination of a caloric curve for the solid-toliquid like transition (melting) of a small cluster, i.e., a sodium cluster ion consisting of 139 atoms. A beam of cluster ions was generated with a canonical distribution of internal energy thus fixing the temperature. One cluster size was selected (thus switching to microcanonical system), irradiated by photons and the photofragmentation pattern (whose positions can be related to the energy) is measured as a function of cluster temperature. Similarly, Bachels et al. [7] reported a caloric curve for a free tin cluster distribution (without mass selection) impinging on a sensitive pyroelectric foil. The interpretation of their experiment was, however, later questioned [8].

Furthermore, Haberland and colleagues [9,10] were able by extending their method to map out the energy distribution of the fragments, demonstrating, by the bimodal shape observed, a possible existence of a negative heat capacity for melting a 147-atom sodium cluster. Haberland and colleagues had, in addition, constructed the caloric curve across the liquid-to-gas like transition [11] using known data of the atomic gas.

M. Farizon and colleagues [12-16] reported experimental measurements of caloric curves for solid or liquid to gas like transition of small charge selected hydrogen clusters. On the other hand, some theoretical and numerical simulations are reported concerning first order transition of argon clusters. In fact, K. L. Nierholz reported the structural and thermodynamic properties of Ar clusters, and derived many thermodynamic properties like caloric heat entropy and internal energy. The obtention of the caloric curves for the microcanonical and the canonical ensemble by molecular dynamics simulations and histogram methods exhibits bimodal behavior (back bending).

The analysis shows that this bimodal behavior is observed only in the case of microcanonical ensemble, while in the case of canonical description, the internal energy increases monotonically with the temperature.

All the reported papers related to the caloric curves present the back bending as a specific property of microcanonical description. For finite size systems, the thermodynamic ensembles are not equivalent.

For macroscopic objects, melting occurs at some well defined temperature; however, this is no longer true for small particles or clusters. At low temperatures, the atoms in a cluster or in a large piece of matter make only small amplitude vibrations around a fixed position. It takes a lot of energy to push an atom from its position in a solid. If the temperature increases, atoms in the cluster can affect neighboring sites and start a diffusive motion.

All properties change with the size of a cluster. The transition from the atom or molecule to the bulk is often quite smooth and the asymptotic behavior well understood. This is not the case for some thermal properties. Large and irregular fluctuations are observed, e.g. in the melting temperature, even for clusters containing several hundred of atoms.

2. Theoretical Approach

According to the Boltzmann’s definition, for an isolated system in equilibrium in microcanonical ensemble (NVE), the entropy is given by [17]

(1)

where k_B is the Boltzmann constant and Ω(NVE) is the number of discrete microstates in the configurational Г- space consistent with (NVE).

The entropy of any thermodynamic system in equilibrium is given by [18]

(2)

where p_i is the probability of the microstate i, and the summation is over all possible microstates compatible with the system.

In 1988 Tsallis introduced a new definition for entropy [19,20]

(3)

where 0 ≤ q ≤ 1 (q = 1 in the sense of q → 1).

The non extensivity of the entropy observed for small systems is not a consequence of the Tsallis definition. The entropies investigated [21] are all in the framework of the Boltzmann-Gibbs definition. Indeed, all of the non extensivities and non intensities observed in this paper are triggered by conversion of some of the external potential energy into the internal potential energy of the system, an effect that is significant only in small systems.

Entropy, like internal energy, becomes non extensive in small systems. However, there are two differences here:

• First, the non extensivity in the entropy is not as profound as in the internal energy.

• Second, while the internal energy decreases (becomes more negative) due to non extensivity, the entropy increases and, hence, is super extensive.

When two systems are combined, the entropy of the combination of the combined system, S_q(A + B) is given by

(4)

If two systems A and B are independent in the sense of the theory of probability, then q = 1 and the entropy is simply the sum of the entropies of its constituents; therefore, the entropy should be an extensive parameter. However, when the systems A and B are dependent, then q < 1 and the entropy as defined by Equation (4) is non extensive, regardless of the size of the system.

Yi-Fang Chang [22] reported that, according to the Boltzmann and Einstein fluctuation theories, all possible microscopic states of a system are equally probable in thermodynamic equilibrium and the entropy tends to a maximum value finally. It is known from statistical mechanics that fluctuations of the entropy may occur [23], while fluctuations always decrease the entropy [24].

When internal interactions exist among subsystems, the statistical independence and equal-probability are unavailable. If fluctuations are magnified [23,24] and the order parameter comes to a threshold value, phase transition will occur. In this case, the entropy may decrease in an isolated system, at least within a certain time. A self-organized structure whose entropy is smaller will be formed.

The Internal Energy of a system is [25] of a system is

(5)

where is the additive part of the particle energy in the state s, in most cases it and E are the kinetic energy; W′_ss and U′_ss are the absolute values of the attraction and repulsion energies of particles in the states s and s′, respectively.

When the probability changes with time, the entropy of a system composed of two subsystems changes also with time, and the entropy would be defined as [26]:

(6)

where. This shows that the entropy decreases with the internal interaction. Not only is this conclusion the same with the conditioned entropy on ρ₁ and ρ₂, but also it is consistent with the systems theory in which the total may not equal the sum of parts.

3. MD Simulation

The melting-like transition of Ar_N has been investigated through micro-canonical MD simulations using a LennardJones potential,

(7)

where r_ij is the distance between the particles, and ε and σ are the potential parameters of the LJ potential chosen appropriately to represent particular system of interest, which would be appropriate for solid argon.

The calculations are performed in dimensionless variables by scaling the energy in ε, length in σ, mass in M, temperature in ε/k_B, pressure in ε/σ³, time in (Mσ²/ε)^1/2.

Here we have chosen M = 0.66 × 10^{−22 g} is the atomic mass, and k_B is the Boltzman constant. ε/k_B = 120 K, σ = 3.84 Å, ε = 1.65324 × 10⁻¹⁴ erg, ε/σ³ = 287 bar and (Mσ²/ε)^1/2 = 2.44 ps.

Initially we assign a face-centred-cubic lattice to the positions of the atoms. To start the algorithm, velocities are all set to zero. Sometimes, instead of the velocities being scaled, they are in any case, one has to check the velocity distribution after the equilibration phase has been reached to make sure that it has the equilibrium Maxwell-Boltzmann form. At the end of this equilibration period, all memory of the initial configuration should have been lost.

In the present work, we adopt the latter approach to perform molecular dynamics simulations of melting-like transition of Ar_N, employing the velocity form of the Verlet algorithm [27].

(8)

The important parameter to choose in an MD simulation is the time increment Δt [28]. In a microcanonical ensemble simulation, the total energy of the system must be conserved. If Δt is too large, steps might become too large and the particle may enter the classically forbidden region where the potential energy is an increasing function of position. This can occur when two particles collide or when a particle hits the “wall” imposed by the external potential. Entering the classically forbidden region means that the new potential energy has become higher than the maximal value allowed. In this case, the total energy has increased, and this phenomenon keeps occurring for large step sizes until the total energy diverges. So, depending on the available total energy, Δt should be chosen small enough so that the total energy remains constant at all times, but not so small that it would require an extremely large number of steps to perform the simulation. The optimal value of Δt is usually found by trial and error. One femto (10⁻¹⁵) second is a good trial guess, but the optimal value really depends on the initial energy and the kind of potential considered.

For the micro-canonical ensemble (constant energy), the mean kinetic energy for different total energies is determined. Then the temperature T is defined by the equipartition theorem; then for practical purposes, it is common practice to define an instantaneous temperature T(t), proportional to the instantaneous kinetic energy EK(t) by a relation analogous to the one above.

(9)

With the kinetic energy EK and the number of atoms N, and K_B is the Boltzmann constant, and 3N − 6 represents the total number of internal degrees of freedom.

The total internal energy of a system can be written as a sum of kinetic and potential energy contributions, and the temperature of the system is proportional to its average kinetic energy.

The cluster’s configuration has a lower configurational energy EP which nearly compensates the effects of the potential truncation at 2.5σ.

Temperature changes are usually achieved by enabling a device in the code that brings the system to the desired temperature TD by rescaling the velocities. In the velocity verlet algorithm discussed at the end of this may be accomplished by replacing the equation

(10)

by

(11)

where T_D is the desired temperature, and T(t) the instantaneous temperature.

Switching to liquid-like state it is increasing its configurational energy; hence the kinetic energy has to be reduced to keep the total energy constant. Because of this, the same mean kinetic energy can be obtained at different total energies.

An observable commonly studied in the context of melting transitions is Lindemann’s parameters <δ> [29],

(12)

where N is the number of molecules in a cluster and r_ij is the distance between i^th and j^th Argon atoms.

This parameter measures the root mean square of the distance between two atoms averaged over all pairs. Even short isomer fluctuations with a subsequent return to the ground state can leave the cluster reordered; i.e., a previously nearest neighbour pair r_ij may become a secondor third-neighbour pair after the fluctuation. Such a reordering leads to a notable increase in, allowing for the clear identification of a melting transition. Note that the Lindemann criterion of melting, which measures the atomic fluctuations with respect to their equilibrium positions, is usually employed for bulk systems but is less well suited for cluster systems.

The second observable of interest is the specific heat, which can be obtained from the thermodynamic averages of the potential energy EP and its square.

The second observable of interest is the specific heat, which can be derived from the total energy expression.

4. Simulation Results

During the first 2500 MD steps, for N = 256 Argon atoms, positions are defined on a lattice assuming CFC crystal structure witch correspond to the most stable one at T = 0 with Leonnard Jones potential. Initial velocities have been taken to be zero, or are assigned taking them from a Maxwell distribution at a certain temperature T.

To ensure a reasonable numerical stability, the basic time increment is taken to be Δt = 0.008. The actual time will be fairly small since only a limited number of integration steps are possible.

Of course, such an initial state will not correspond to an equilibrium condition. However, the steady state is reached after a time of the other of 500 time steps Δt (Δt = 0.008), and one should wait for the system to reach equilibrium under clean constant energy conditions before collecting data. 1000 additional time steps are needed to determine the mean values of the calorific curve, heat capacity, root-mean-square (RMS) bond length fluctuations, and melting temperatures corresponding to the desired temperature T_D = 1.3476. The system is driving to this desired temperature by rescaling the velocities in the velocity Verlet algorithm in order to keep the system under control. All quantities are given in reduced units.

For Δt = 0.008 and at the indicated desired temperature TD, the kinetic energy distribution of the N Ar atoms is plotted in Figure 1. This curve exhibits Maxwell distri-

bution shape.

For a step of integration Δt = 0.008 and then by varying the temperature ΔT_D = 0.5, for each value of the temperature T_D are averaged physical quantities such as temperature, kinetic energy, potential energy, total energy, specific heat, and Lindemann factor. Then we draw for integration step Δt = 0.008 in Figure 2 the evolution of the average kinetic energy, potential, and total according to the desired temperature.

Starting from T = 0, the caloric curve increases roughly linearly, but near the melting temperature there is an acute change in slope characterized by a sudden jump analogous to bulk melting behavior. This jump is expected to be smoothed to a finite with small systems [27-30].

In our case, we observe that this jump occurs at the reduced temperature T_melt (T_melt =1.45) which corresponds to 174 K considered as the liquid-solid transition.

According to the study by Lutsko et al. [31], a factor should be introduced between the simulated structure of an isolated cluster and the actual (thermodynamic) melting points when periodic boundary conditions are enforced for bulk materials. The thermodynamic melting point is approximately 1.4 - 1.5 times its structural counterpart.

For a given same values of the integration time step, we observe, in addition, that the microcanonical internal energy exhibits an S-shaped curvature in the coexistence region where both solid-like and liquid-like states can occur. This can be understood by considering the constant energy ET = EP + EK, the cluster switches between solid and liquid along time scale. In solid-like state, the configuration of the cluster has a lower configurational energy. Switching to a liquid-like state, it is increasing its total energy; hence, the kinetic energy has to be reduced to keep the total energy constant. Because of this, the same mean kinetic energy can be obtained at different total energies.

With, we have three regions;

Region I

Region II

Region III

The potential curve is characterized by two distinct parts, where the solid part corresponds to Region I and the liquid part to Region III. Those two parts are separated by a gap in which the cluster fluctuates between the two phases; this region of coexistence of solid and liquid phases represents two behaviors where the potential energy distribution has a dent with an inverted curvature, and shows a bimodal structure.

Then the change in the potential energy can be interpreted as responsible for this back bending phenomenon.

Then we draw for integration step Δt = 0.008 in Figure 3 the evolution of the average specific heat according to the desired temperature. In the transition region, we observe a double peak structure that also can be interpreted by the coexistence of the solid and the liquid phase.

In our simulations of melting, we monitor the phase changes of clusters using atom-resolved root-mean-square (RMS) fluctuation of the interatomic distances.

For the solid state, in which atoms are fixed in lattice sites, the value of this parameter is typically smaller than 0.1, considerably smaller than that for the liquid state with its mobile atoms. This parameter corresponds to a time scale in which atomic transitions between sites for the solid cluster state are improbable, while these transitions for the liquid cluster state are effective.

Figure 4 demonstrates this for the 256-atom LennardJones cluster with the argon pair interaction parameters.