3
Melting of a two-dimensional system
We might look at the results of the last chapter from a slightly different perspective. By actively shearing the system, we put energy into it and, so to say, increased its temperature.
Eventually, this led to a phase transition where the system changed from an ordered – or crystalline – to an unordered – or liquid – state. Therefore, a more straightforward approach would be to directly put the system in contact with a heat bath. Naturally, two different questions come to mind. Obviously, in the limit of very low temperatures, the system should be found in the ordered state of a solid. At high temperatures on the other hand, we would expect to find a liquid-like behavior. Thus, the first question we might ask is at which temperature a phase transition of the system occurs. The second question is closely related to the first, as we might want to investigate the microscopic mechanism that accompanies the phase transition. Hence, besides a characterization of the state of the system, we will also focus on particle dynamics.
Melting of a two-dimensional system
−4 −2 0 2 4
r
10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 101 102
Φi(r,µ)
−3 −2 −1 0 1 2 3
r
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
Φ(r,µ)
µ= 3.0 µ= 6.0 µ= 9.0 µ= 12.0
a) b)
Figure 3.1: Shape of the Yukawa potential for different values of the screening parameter µ. (a) Single-particle potential. (b) Potential for a particle in a 1-dimensional chain with particles fixed at r =k and k ∈ Z0. With increasing µ the potential becomes steeper close to particle positions and flat in the intermediate regions, thus resembling a solid-core interaction.
between the particles in units of the lattice spacing. Since for now all particles are chosen to be identical,σij just gives an multiplicative offset and we setσij = 1. A largeµcorresponds to a short interaction range, mimicking hard cores, whereas the potential shifts towards a Coulomb potential for µ → 0. This dependence is also shown in figure 3.1, both for the pairwise potential and a 1-dimensional array of interacting particles. By convention, the potential energy is split between the two interacting particles so that the per-particle energy is defined as
Φi=X
j6=i 1
2Φij. (3.2)
Hence, the full potential energy is given by E= Φ =X
i
Φi= X
i,j<i
Φij. (3.3)
Physical units and scaling
The Yukawa potential describes a shielded electrostatic interaction, hence the parameters in eqn. (3.1) correspond to physical units as follows. The interaction strengthσij represents the Coulombic interaction between particles,
σij = QiQj
4πǫ0r0, (3.4)
whereQi are particle charges, ǫ0 is the vacuum permittivity and r0 represents the lattice spacing. The screening parameter, on the other hand, is given by
µ= r0
λ, (3.5)
38
3.1 Definition of the potential
where λis the Debye-length.
We measure energies, including thermal energy, in units of the mean interaction strength hσiji. The unit of masses is m, and the relevant length-scale is the lattice constant r0. Accordingly, the unit of time is τ0= mr20/hσiji1/2
. Forces
The two-particle forces are given by the negative of the gradient, Fij =−Fji=−∇iΦij = (µrij + 1)e−µrij
rij3 (xi−xj) (3.6) The Jacobian, on the other hand, is constructed by matrices
Jij =−(µrij + 1)e−µrij
rij3 I2+ µ2r2ij+ 3µrij + 3e−µrij
rij5 (xi−xj)⊗(xi−xj) (3.7a) for i 6= j, where ⊗ denotes the dyadic product. The diagonal part Jii is given by the negative sum over all matrices Jij in a row or column,
Jii=−X
j6=i
Jij =−X
j6=i
Jji. (3.7b)
As a consequence, the inter-particle Jacobian is symmetric, which is favorable when investigating the spectra. Within this definition, unstable modes of the spectrum have positive eigenvalues, whereas stable modes have negative ones.
When evaluating the potential numerically, we naturally have to cut off the potential at an appropriate radius rcut. As a condition we chose
Φij(µ, rcut)≤1×10−8.
The critical radius where this is fulfilled is shown in figure 3.2. It allows us to optimize stationary states up to a precision of kFk <10−8. Consequently, we considered a cut-off radius ofrcut= 8.0 to be sufficient for a screening ofµ >2. Belowµ= 1 the critical radius quickly diverges as the potential approaches its Coulomb interaction limit. Therefore, we restrict ourselves to the range 1≤µ≤12
System definitions
We consider a two-dimensional system. Unless stated otherwise, it is subject to periodic boundary conditions. Since the potential is purely repulsive, the system forms a Wigner crystal and its minimum energy configuration is given by a hexagonal lattice (Wigner, 1934; Meissner et al., 1976). Box-sizes are chosen such as to accommodate the hexagonal lattice with lattice constant r0, Lx = Nxr0 the width and Ly = √
3/2Nyr0 the height.
Correspondingly, N = NxNy is the total number of particles and n = 2N the system dimension. This results in a number density of ρ= 2/√
3≈1.15. In general, we consider systems withN = 2500 particles and choose a screening ofµ= 3.0, if not stated otherwise.
Melting of a two-dimensional system
2 4 6 8 10 12
µ 0
2 4 6 8 10 12 14
rcrit
Φ(µ, rcrit) = 1×10−8
Figure 3.2: Critical interaction range of the Yukawa potential as a function of the screening parameter µ. A cut-off ofrcut= 8.0 is considered sufficient forµ >2.
In some special cases we will use different boundary conditions, e.g. in the low-dimen-sional model system in the next chapter and when investigating the connection to elasticity theory in chapter 5. In either case, a fraction of the particles is kept fixed at positions of the underlying hexagonal lattice.
Natural oscillations
In thermal equilibrium and at low temperatures, particles in a crystalline configuration will oscillate about their lattice positions. The corresponding frequency can be estimated by the Einstein frequency ωE. It is defined as the frequency at which a single particle in the crystal vibrates, with all other particles being held fixed at their respective lattice sites. Thus, the other particles effectively constitute a potential well. The frequency can then be extracted by the harmonic approximation of the single-particle potential,
Φ1(x) =X
j6=1 1
2Φ1j ≈E0,1+ 12ω2E u2x+u2y
. (3.8)
Here,E0,1 is the particle’s potential energy at the lattice position, and ux and uy are the displacements in thex- andy-direction, respectively. Comparing eqn. (3.8) to the Taylor expansion of the potential,ω2E is found to be the negative eigenvalue of the single-particle Jacobian, eqn. (3.7b), which is degenerate due to rotational invariance of the single-particle potential in the crystalline configuration1. Moreover, since the submatrices Jii are given by the single-particle Jacobians, and the trace of a matrix equals the sum of its eigenvalues, the Einstein frequency is readily obtained by
ωE2 = 1
2Ntr (J). (3.9)
In the case of µ= 3, this yields an Einstein frequency of ωE2(µ= 3) = 2.098.
1We will return to this issue later on in chapter 5.5.
40
3.1 Definition of the potential
Scalings in other publications
Throughout publications addressing the Yukawa system, different scalings are used. First, instead of measuring thermal energy in units of the interaction energy, often the opposite is done. Thus, an interaction strength is defined as
Γ = σij
kBT. (3.10)
As a consequence, the temperature of the system is kept fixed and the interaction strength is varied in order to find e.g. phase transitions. The physical motivation is found in the experiments, where it is easier to adjust the interaction between particles instead of increasing the temperature of the heat bath. Moreover, changing the temperature in general also affects the Debye length, which is closely related to the screening parameterµ.
However, this is only of importance when actually simulating a specific physical system.
With our choice of σij = 1, the relation between Γ and the temperature in our system is readily given by
Γ = 1/kBT, (3.11)
i.e. Γ is the inverse temperature.
Another issue is the choice of the relevant length-scale. In place of the lattice spacing r0, frequently the Wigner-Seitz radius b is considered. In a two-dimensional system it is defined as the radius of the circle whose area equals the average area per particle in the system,
LxLy N = 1
ρ =πb2. (3.12)
It is thus an indicator of the density of the system. The relation (3.12) is essential when comparing to results in other publications, since it changes both the screening parameterµ as well as the interaction strength. For a box commensurate with the lattice, i.e.Lx=Nxr0 and Ly =√
3/2Nyr0, we find a direct relation betweenr0 and b:
r0 = s 2
ρ√ 3 =
s2π
√3b≈1.9b. (3.13)
Accordingly, our system is characterized by a Wigner-Seitz radius of b ≈0.53r0, and the parameters µand Γ have to be rescaled by
µ= 1.9µb and Γ = Γb
1.9, (3.14)
where the index bdenotes the parameters rescaled to the Wigner-Seitz radius.
A second length scale may be introduced when treating the system as a liquid instead of a crystal. The average distance between particles is then directly related to the particle density,
rp =p ρ−1 =
s√
3/2r20N
N =
s√ 3
2 r0. (3.15)
Melting of a two-dimensional system
Similarly to the Wigner-Seitz radius, quantities have to be rescaled according to µ= 1.075µp and Γ = Γp
1.075. (3.16)
Here, the indexp refers to the parameters rescaled to the average particle distance.