Stepwise Homogeneous Melting of Benzene Phase I at High Pressure

: We investigate, using molecular dynamics simulations, the spontaneous homogeneous melting of benzene phase I under a high pressure of 1.0 GPa. We ﬁnd an apparent stepwise transition via a metastable crystal phase, unlike the direct melting observed at ambient pressure. The transition to the metastable phase is achieved by rotational motions, without the di ﬀ usion of the center of mass of benzene. The metastable crystal completely occupies the whole space and maintains its structure for at least several picoseconds, so that the phase seems to have a local free energy minimum. The unit cell is found to be unique—no such crystalline structure has been reported so far. Furthermore, we discuss the inﬂuence of pressure control on the melting behavior.


Introduction
Benzene is a renowned small and simple molecule. Regardless, its polymorphism [1][2][3] and local liquid structure [4][5][6] have been extensively explored until now. Five crystalline phases of pure benzene have been reported so far experimentally [7][8][9][10][11][12][13][14], and computational studies predicted many other potential crystalline structures [1,[15][16][17][18]. Although the phase diagram shows the most stable phase at a given condition, there can be rich intermediate phases on the transition pathway from the initial to the finally prevailed phase [19][20][21][22][23]. Ostwald argued that a phase transition can proceed via an intermediate metastable phase due to the reduction in the surface energy of nucleation (Ostwald's step rule) [24] in contrast to the direct nucleation described in the classical nucleation theory [25]. Indeed, a metastable phase or structure has been shown to play a key role in the melting of colloids [26], copper and aluminum [27], and ice [28]. Further, under high pressures, the melting pathway can be much more complicated, because of the competition between potential energy, entropy, and packing effect (negligible at ambient pressure).
The existence of intermediate metastable phase in the melting of benzene phase I is still controversial, although the phase transition dynamics of benzene has been investigated [29][30][31]. Tohji et al. claimed the existence of a premelting stage at 10 K below melting point using powder X-ray diffraction [30]. Furthermore, they predicted that a plastic crystal transiently appears in the melting. This premelting stage was later disproven by the investigation over the complete temperature range of 4 to 280 K [31]. Very recently, we performed the molecular dynamics (MD) simulations for the homogeneous melting of benzene phase I crystal near the limit of superheating and statistically demonstrated that there is no intermediate transient state between the crystal and liquid phases [32]. While these studies have been conducted under ambient pressure, the melting dynamics under high pressure have not been explored.
In this study, we perform MD simulations of the homogeneous melting of benzene phase I at a high pressure of 1.0 GPa and observe the stepwise transition via a metastable crystal phase. We show that the formation of the metastable phase is achieved by rotational motions, and the unit cell structure differs from any other crystalline phases reported so far.

Force Field
Benzene was modelled with full atomistic detail using CHARM22 [33,34]. This force field reproduces the crystal structure of benzene phase I [29,35] and the physical properties of liquid benzene, such as density and self-diffusion coefficient [36]. The intermolecular nonbonded interactions were described by Lennard-Jones plus Coulomb potential. The intermolecular interactions were truncated at 1.20 nm. The Lennard-Jones parameters for cross-interactions were obtained using Lorentz-Berthelot mixing rules and the long-range Coulombic interactions were evaluated using particle-mesh Ewald algorithm [37].

MD Simulations
MD simulations were performed using GROMACS 2019 package [38], in which the equations of motion are integrated using leapfrog algorithm with a time step of 2 fs. The temperature T for equilibration was controlled using a Berendsen thermostat [39] with a damping constant of 1.0 ps, while a Nose-Hoover thermostat [40,41] was used for the production runs. The pressure p was isotropically controlled by a Berendsen barostat [39] for both equilibration and production runs, with a damping constant of 2.0 ps. The pressure was set to 1.0 GPa for all the simulations. The periodic boundary condition was applied in all three directions.
To evaluate the influence of pressure control on the melting behavior, we performed some additional simulations using anisotropic pressure control, where the box dimensions and the box angles can change.

Crystals
The structure of benzene phase I was obtained from Cambridge Structural Database (refcode: BENZEN) [42]. Phase I is the most stable crystalline phase up to 1.2 GPa [9,14]. The unit cell consists of four benzene molecules with space group Pbca. The lattice parameters are a = 0.7460 nm, b = 0.9666 nm, and c = 0.7034 nm at 270 K [11]. Lattice axes a, b, and c corresponded to x, y, and z axes, respectively, in this study.
For the spontaneous homogeneous melting, we replicated the unit cell to manufacture an approximately cubic system with a box size of 5.21 nm × 4.78 nm × 4.84 nm. The number of molecules, N, was 980. To equilibrate the crystal structure, energy minimization was first carried out using the steepest descent algorithm. Subsequently, an isochoric isothermal (NVT) MD simulation was performed at 270 K for 500 ps. Thereafter, an isobaric isothermal (NpT) MD simulation for 5 ns at 270 K and 0.1 MPa was performed, followed by a 5 ns NpT MD simulation at 270 K and 1.0 GPa. Furthermore, an NpT MD simulation at 500 K and 1.0 GPa was performed for 5 ns. The final configuration obtained was used as the initial structure for the following heating simulations.
To determine the equilibrium melting temperature we used the two-phase method [43][44][45][46]. First, we developed two identical boxes of benzene phase I crystal containing 448 molecules with a box size of 2.68 nm × 6.45 nm × 2.86 nm. We randomly erased 48 molecules from one of these crystal structures, and we melted it by performing a short NVT MD simulation at 800 K. Then, we joined these crystal and liquid boxes. The box dimension was 2.68 nm × 13.5 nm × 2.86 nm. To equilibrate the solid-liquid coexistence structure, energy minimization was first carried out using the steepest descent algorithm. Subsequently, an NVT MD simulation was performed at 200 K for 20 ps. Thereafter, the system was gradually heated up by performing each 20 ps NpT MD simulations at two temperatures of 200 and

Results and Discussion
We initially evaluated the equilibrium melting temperature Tm of benzene phase I at 1.0 GPa in our model using the two-phase method [43][44][45][46]. The initial configuration for production runs was the coexistence between liquid and phase I crystal, as shown in Figure 1A. We performed each 10 ns NpT MD simulation at 430K, 435K, and 440 K and at 1.0 GPa. Figure 1B depicts the time evolution of the volume for these runs. At Tm, the phase I crystal and liquid exhibit the same stability in free energy, and the total volume does not change. On the other hand, the volume increases at T > Tm due to melting, while it decreases at T < Tm. The results show that Tm is 435 ± 5 K at 1.0 GPa for our computational model, which is approximately 20 K lower than the experimental estimation of 457 K [9]. We also computed Tm at 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, and 1.2 GPa in the same manner. Figure 1C shows that the melting curve for our computational model qualitatively reproduces the experimental result [9,47].  [9]. The panel C also presents the homogeneous melting temperature (red sharp), and the stable (filled diamonds) and unstable (open diamonds) conditions for the metastable phase, estimated in this study.
Benzene phase I, equilibrated at 500 K, was heated up in the increments of 10 K at 1.0 GPa, and an NpT MD simulation was performed for 10 ns at each T. We first observed the melting of crystal at 590 K. This result indicates that the limit of superheating achieved by our heating rate (1 K/ns) is 590 K, which is considerably higher than the equilibrium melting temperature of 435 K for our computational model. The large superheating for homogeneous melting has also been reported for other materials both experimentally [26,48,49] and computationally [28,[50][51][52][53][54][55], because of the absence of an apparent trigger, such as surfaces and impurities, to initiate the disordering. To observe the melting under moderate conditions, we focused on the melting at 584 K, heated up from 580 K. At 584 K, 6 of 95 independent trajectories melted within 10 ns when we gave different momenta to the same configuration. Figures 2A and 2B present the time evolution of potential energy (PE) and density for a typical melting trajectory at 584 K. We find a plateau between 490 and 510 ps (shown by green shade), where both the PE and density remain constant. The preservation of these properties over a time  [9]. The panel C also presents the homogeneous melting temperature (red sharp), and the stable (filled diamonds) and unstable (open diamonds) conditions for the metastable phase, estimated in this study.
Benzene phase I, equilibrated at 500 K, was heated up in the increments of 10 K at 1.0 GPa, and an NpT MD simulation was performed for 10 ns at each T. We first observed the melting of crystal at 590 K. This result indicates that the limit of superheating achieved by our heating rate (1 K/ns) is 590 K, which is considerably higher than the equilibrium melting temperature of 435 K for our computational model. The large superheating for homogeneous melting has also been reported for other materials both experimentally [26,48,49] and computationally [28,[50][51][52][53][54][55], because of the absence of an apparent trigger, such as surfaces and impurities, to initiate the disordering. To observe the melting under moderate conditions, we focused on the melting at 584 K, heated up from 580 K. At 584 K, 6 of 95 independent trajectories melted within 10 ns when we gave different momenta to the same configuration. Figure 2A,B present the time evolution of potential energy (PE) and density for a typical melting trajectory at 584 K. We find a plateau between 490 and 510 ps (shown by green shade), where both the Crystals 2019, 9, 279 4 of 10 PE and density remain constant. The preservation of these properties over a time range implies the existence of a metastable phase. Figure 3A-C depicts the molecular structures at each stage of 100, 490, and 575 ps, respectively. In contrast to the disordered liquid structure in Figure 3C, layered structures along the y axis are depicted in Figure 3A,B. However, their molecular orientations are clearly different. The same structure as that in Figure 3B was also observed in all the other melting trajectories at 584 K, although the structure sometimes did not fill the whole space of the simulation box and its lifetime was in the order of several 10 picoseconds.
Crystals 2019, 9,x FOR PEER 4 range implies the existence of a metastable phase. Figure 3A-C depicts the molecular structures at each stage of 100, 490, and 575 ps, respectively. In contrast to the disordered liquid structure in Figure 3C, layered structures along the y axis are depicted in Figure 3A,B. However, their molecular orientations are clearly different. The same structure as that in Figure 3B was also observed in all the other melting trajectories at 584 K, although the structure sometimes did not fill the whole space of the simulation box and its lifetime was in the order of several 10 picoseconds.  Although the metastable structure exists for a very short time in the melting pathway at 584 K, the same structure might be (meta) stable at different pT conditions. To explore such a pT region, we performed NpT MD simulations for each 10 ns at pressures ranging from 0.8 to 1.2 GPa and at temperatures ranging from 450 to 550 K. At a given pT condition, we gave 10 different momenta to range implies the existence of a metastable phase. Figure 3A-C depicts the molecular structures at each stage of 100, 490, and 575 ps, respectively. In contrast to the disordered liquid structure in Figure 3C, layered structures along the y axis are depicted in Figure 3A,B. However, their molecular orientations are clearly different. The same structure as that in Figure 3B was also observed in all the other melting trajectories at 584 K, although the structure sometimes did not fill the whole space of the simulation box and its lifetime was in the order of several 10 picoseconds.  Although the metastable structure exists for a very short time in the melting pathway at 584 K, the same structure might be (meta) stable at different pT conditions. To explore such a pT region, we performed NpT MD simulations for each 10 ns at pressures ranging from 0.8 to 1.2 GPa and at temperatures ranging from 450 to 550 K. At a given pT condition, we gave 10 different momenta to Although the metastable structure exists for a very short time in the melting pathway at 584 K, the same structure might be (meta) stable at different pT conditions. To explore such a pT region, we performed NpT MD simulations for each 10 ns at pressures ranging from 0.8 to 1.2 GPa and at temperatures ranging from 450 to 550 K. At a given pT condition, we gave 10 different momenta to the structure obtained at 490 ps in the melting trajectory ( Figure 3B). The metastable structure was found to be preserved for 10 ns in the all trajectories at several thermodynamic conditions: 450 K and 0.8, 0.9, 1.0, 1.1, and 1.2 GPa; 500 K and 0.9, 1.0, 1.1, and 1.2 GPa; 550 K and 1.1 and 1.2 GPa. At the other conditions, the metastable structure melted into liquid. These stable and unstable T-p conditions for the metastable structure are plotted in Figure 1C. These results imply that the metastable structure can be a "phase" rather than a transiently appearing fragile structure. Its free energy is likely to be between phase I and liquid benzene, in accordance with Ostwald's step rule [23].
In contrast to the short lifetime of the metastable phase, the waiting times for the phase transition from phase I to the metastable phase range from hundreds of picoseconds to several nanoseconds at 584 K. The significant difference in their lifetimes indicates that the free energy barrier from the metastable phase to liquid is much lower than that from phase I to the metastable phase at 584 K.
To characterize the self-diffusion over melting, we computed the mean square displacement (MSD) of center of mass (COM) of benzene from the initial configuration at 584 K (t = 0 ps), using the expression of | → r (t) − → r (0)| 2 ( Figure 2C). The time-dependent MSD shows a small hump at around 490 ps, arising from the re-organization of molecules to form the metastable phase. The difference in the MSD of COM is negligibly small (0.04 nm 2 ). The MSD does not change between 490 and 510 ps, where the PE and density are also preserved. After 510 ps, the MSD gradually increases, during which a liquid nucleus forms and erodes the metastable phase. Finally, the MSD increases rapidly after 530 ps, at which point the liquid structure covers approximately half the space of the simulation box. Figure 4 shows the radial distribution functions (RDFs) between COMs of benzene molecules in the three phases. The results show that the metastable phase maintains a long-range order structure in comparison with the liquid phase. Although the metastable phase exhibits broader peaks than phase I, the peaks in the two phases are located at a similar radial distance. This indicates that there is only a minor change in the position of COM between phase I and the metastable phases. Hence, the analyses of MSD and RDF demonstrate that the transition from phase I to the metastable phase arises from the rotational motion, and COM in the metastable phase exhibits an ordered structure. Crystals 2019, 9,x FOR PEER 5 the structure obtained at 490 ps in the melting trajectory ( Figure 3B). The metastable structure was found to be preserved for 10 ns in the all trajectories at several thermodynamic conditions: 450 K and 0.8, 0.9, 1.0, 1.1, and 1.2 GPa; 500 K and 0.9, 1.0, 1.1, and 1.2 GPa; 550 K and 1.1 and 1.2 GPa. At the other conditions, the metastable structure melted into liquid. These stable and unstable T-p conditions for the metastable structure are plotted in Figure 1C. These results imply that the metastable structure can be a "phase" rather than a transiently appearing fragile structure. Its free energy is likely to be between phase I and liquid benzene, in accordance with Ostwald's step rule [23].
In contrast to the short lifetime of the metastable phase, the waiting times for the phase transition from phase I to the metastable phase range from hundreds of picoseconds to several nanoseconds at 584 K. The significant difference in their lifetimes indicates that the free energy barrier from the metastable phase to liquid is much lower than that from phase I to the metastable phase at 584 K.
To characterize the self-diffusion over melting, we computed the mean square displacement (MSD) of center of mass (COM) of benzene from the initial configuration at 584 K (t = 0 ps), using the expression of ⟨| ⃗ − ⃗ 0 |⟩ ( Figure 2C). The time-dependent MSD shows a small hump at around 490 ps, arising from the re-organization of molecules to form the metastable phase. The difference in the MSD of COM is negligibly small (0.04 nm 2 ). The MSD does not change between 490 and 510 ps, where the PE and density are also preserved. After 510 ps, the MSD gradually increases, during which a liquid nucleus forms and erodes the metastable phase. Finally, the MSD increases rapidly after 530 ps, at which point the liquid structure covers approximately half the space of the simulation box. Figure 4 shows the radial distribution functions (RDFs) between COMs of benzene molecules in the three phases. The results show that the metastable phase maintains a long-range order structure in comparison with the liquid phase. Although the metastable phase exhibits broader peaks than phase I, the peaks in the two phases are located at a similar radial distance. This indicates that there is only a minor change in the position of COM between phase I and the metastable phases. Hence, the analyses of MSD and RDF demonstrate that the transition from phase I to the metastable phase arises from the rotational motion, and COM in the metastable phase exhibits an ordered structure. We next characterize the molecular orientation in the metastable phase. Figure 3B demonstrates that benzene molecules in every other layer of the metastable phase (we name them M1 and M2 layers) tend to have the same alignment. We project the probability density of the normal vector on a sphere using HEALPix algorithm [56][57][58], in which the order parameter ( , ) is used to define the vector orientation ( Figure 5A). In phase I, benzene molecules are most likely to orient at (90°, 40°) and (90°, 140°), as shown in Figure 5B. It is to be noted that the molecular orientations of (90°, 40°) and (90°, 140°) are identical to (90°, −140°) and (90°, −40°), respectively. The molecular orientation in We next characterize the molecular orientation in the metastable phase. Figure 3B demonstrates that benzene molecules in every other layer of the metastable phase (we name them M1 and M2 layers) tend to have the same alignment. We project the probability density of the normal vector on a sphere using HEALPix algorithm [56][57][58], in which the order parameter (θ, ϕ) is used to define the vector orientation ( Figure 5A). In phase I, benzene molecules are most likely to orient at (90 • , 40 • ) and (90 • , 140 • ), as shown in Figure 5B. It is to be noted that the molecular orientations of (90 • , 40 • ) and (90 • , 140 • ) are identical to (90 • , −140 • ) and (90 • , −40 • ), respectively. The molecular orientation in phase I is slightly different from that obtained experimentally [11] (see the green circles in Figure 5B); the difference probably arises from the limitations of the force field or the effect of high pressure. In the metastable phase, there are distinct orientations in both M1 and M2 layers ( Figure 5C,D). More specifically, the molecules in the M1 layer are most likely to orient at (90 • , 0 • ) and (90 • , 180 • ) (or (90 • , −180 • )), while those in the M2 layer are most likely to orient at (90 • , 90 • ) and (90 • , −90 • ). These preferential orientations are obviously different from those for phase I. These results indicate that all molecules change the orientations in the transition from phase I to the metastable phase. In the previous study, we have shown that at ambient pressure the flipping motion played an important role in the formation of the critical nucleus [32]. Hence, the rotational motion triggers the melting transition of phase I at both the pressures. Further, it is primarily the high pressure that induces the formation of the metastable phase in the melting dynamics, because there is no intermediate metastable phase at ambient pressure. phase I is slightly different from that obtained experimentally [11] (see the green circles in Figure  5B); the difference probably arises from the limitations of the force field or the effect of high pressure.
In the metastable phase, there are distinct orientations in both M1 and M2 layers (Figures 5C and  5D). More specifically, the molecules in the M1 layer are most likely to orient at (90°, 0°) and (90°, 180°) (or (90°, −180°)), while those in the M2 layer are most likely to orient at (90°, 90°) and (90°, −90°). These preferential orientations are obviously different from those for phase I. These results indicate that all molecules change the orientations in the transition from phase I to the metastable phase. In the previous study, we have shown that at ambient pressure the flipping motion played an important role in the formation of the critical nucleus [32]. Hence, the rotational motion triggers the melting transition of phase I at both the pressures. Further, it is primarily the high pressure that induces the formation of the metastable phase in the melting dynamics, because there is no intermediate metastable phase at ambient pressure. Furthermore, we find that the probability densities in the metastable phase are more broadly distributed than phase I. The loose orientation results in the increase of conformational entropy, which may compensate for the increase in PE (Figure 2A) and stabilize the metastable phase in term of free energy. Thus, based on the preferential orientation in each layer and the negligible diffusivity of COM, we conclude that the metastable phase is a crystal.
Using the same method, we evaluated the probability density for the two-fold axis of the benzene molecule, which is in parallel to the molecular plane and determined its most probable orientation. The average box size was calculated between 490 and 510 ps. Based on the repeating structural pattern, we determined that the unit cell is tetragonal and contains two molecules with a dimension of 0.51 nm x 0.51 nm x 0.96 nm ( Figure 6). The space group is P42nm. Interestingly, as far as we know, the obtained unit cell structure is different from any other crystalline structures of benzene reported so far [1,[7][8][9][10][11][12][13][14][15][16][17][18] It is worth investigating the influence of pressure control on the melting behavior. In the above simulations, we used the isotropic pressure control, where the dimension of box was isotropically changed. For comparison, we performed some NpT MD simulations using anisotropic pressure Furthermore, we find that the probability densities in the metastable phase are more broadly distributed than phase I. The loose orientation results in the increase of conformational entropy, which may compensate for the increase in PE (Figure 2A) and stabilize the metastable phase in term of free energy. Thus, based on the preferential orientation in each layer and the negligible diffusivity of COM, we conclude that the metastable phase is a crystal.
Using the same method, we evaluated the probability density for the two-fold axis of the benzene molecule, which is in parallel to the molecular plane and determined its most probable orientation. The average box size was calculated between 490 and 510 ps. Based on the repeating structural pattern, we determined that the unit cell is tetragonal and contains two molecules with a dimension of 0.51 nm × 0.51 nm × 0.96 nm ( Figure 6). The space group is P4 2 nm. Interestingly, as far as we know, Crystals 2019, 9, 279 7 of 10 the obtained unit cell structure is different from any other crystalline structures of benzene reported so far [1,[7][8][9][10][11][12][13][14][15][16][17][18] formed phase III structure was preserved for 10 ns. Note that on the pT phase diagram, the phase III locates at very high pressures over 4 GPa and does not have the phase boundary with phase I [9,14]. Our results imply that the phase III can be a metastable phase between phase I and liquid. Second, we heated up the phase I crystal under the anisotropic pressure control until it melted at 580K. However, we did not observe a stepwise homogeneous melting via phase III at the limit of superheating. The phase III may appear on the melting pathway at lower temperatures. These comparisons indicate that the metastable structure that we found in this study is likely to appear under the isotropic pressure control, which can be realized by the recently developed experimental techniques [59][60][61].

Summary
We carried out the MD simulations to investigate the homogeneous melting of benzene phase I under a high pressure of 1.0 GPa. We found a stepwise change in the PE, density, and MSD in the melting trajectories. The metastable phase preserves the layered structure as phase I but the molecular orientations are obviously different from those of phase I. The formation of the metastable phase is attributable to the flipping motion of benzene rather than the diffusion of COM. Although the flipping motion has been previously shown to lead the formation of the critical nucleus, no intermediate metastable phase was observed at ambient pressure. Furthermore, the probability density of the normal vector demonstrates the preferential orientation in the metastable phase. Their orientations are more relaxed than those in phase I, indicating the gain of conformational entropy compensates the increase of potential energy. We further determined the unit cell and found that such a crystalline structure has not been reported so far.
It is worth investigating the influence of pressure control on the melting behavior. In the above simulations, we used the isotropic pressure control, where the dimension of box was isotropically changed. For comparison, we performed some NpT MD simulations using anisotropic pressure control, which allows the deformation of the simulation box. First, we performed 10 independent NpT MD simulations at 500 K for each 10 ns under the anisotropic pressure control, initiated from the metastable structure obtained under the isotropic pressure control. In all the trajectories, the metastable structure transformed into a different crystal form-phase III with line defect [16]. The formed phase III structure was preserved for 10 ns. Note that on the pT phase diagram, the phase III locates at very high pressures over 4 GPa and does not have the phase boundary with phase I [9,14]. Our results imply that the phase III can be a metastable phase between phase I and liquid. Second, we heated up the phase I crystal under the anisotropic pressure control until it melted at 580K. However, we did not observe a stepwise homogeneous melting via phase III at the limit of superheating. The phase III may appear on the melting pathway at lower temperatures. These comparisons indicate that the metastable structure that we found in this study is likely to appear under the isotropic pressure control, which can be realized by the recently developed experimental techniques [59][60][61].

Summary
We carried out the MD simulations to investigate the homogeneous melting of benzene phase I under a high pressure of 1.0 GPa. We found a stepwise change in the PE, density, and MSD in the melting trajectories. The metastable phase preserves the layered structure as phase I but the molecular orientations are obviously different from those of phase I. The formation of the metastable phase is attributable to the flipping motion of benzene rather than the diffusion of COM. Although the flipping motion has been previously shown to lead the formation of the critical nucleus, no intermediate metastable phase was observed at ambient pressure. Furthermore, the probability density of the normal vector demonstrates the preferential orientation in the metastable phase. Their orientations are more relaxed than those in phase I, indicating the gain of conformational entropy compensates the increase of potential energy. We further determined the unit cell and found that such a crystalline structure has not been reported so far.

Funding:
We acknowledge the support from the Japan Society for the Promotion of Science (KAKENHI 18K19060) and Inoue Foundation for Science (Inoue Science Research Award).