Molecular dynamic simulations of surface morphology and pulsed laser deposition growth of lithium niobate thin films on silicon substrate*
Liu Yuea), Zhu Hao-Nanb), Pei Zi-Donga), Kong Yong-Faa),c)†, Xu Jing-Juna)‡
The MOE Key Laboratory of Weak-Light Nonlinear Photonics and TEDA Applied Physics School, Nankai University, Tianjin 300457, China
School of Mathematical Science, Peking University, Beijing 100871, China
R & D Center, Taishan Sports Industry Group, Leling 253600, China

Corresponding author. E-mail:

Corresponding author. E-mail:

*Project supported by the National Basic Research Program of China (Grant No. 2011CB922003), the International S&T Cooperation Program of China (Grant No. 2013DFG52660), the Taishan Scholar Construction Project Special Fund, China, and the Fundamental Research Funds for the Central Universities, China (Grant Nos. 65030091 and 65010961).


The molecular dynamic simulation of lithium niobate thin films deposited on silicon substrate is carried out by using the dissipative particle dynamics method. The simulation results show that the Si (111) surface is more suitable for the growth of smooth LiNbO3 thin films compared to the Si(100) surface, and the optimal deposition temperature is around 873 K, which is consistent with the atomic force microscope results. In addition, the calculation molecular number is increased to take the electron spins and other molecular details into account.

Keyword: 68.35.−p; 61.72.uf; 31.15.xv; lithium niobate; silicon; dissipative particle dynamics
1. Introduction

Lithium niobate (LiNbO3, or LN) has been intensively studied due to its interesting electro-optic, acousto-optic, elasto-optic, piezoelectric, pyroelectric, and nonlinear optical properties. Moreover, it can be easily doped with laser-active ions, which allows the simple fabrication of low-loss optical waveguides. In industry, LiNbO3 has been considered as the “ silicon of nonlinear optics” or “ silicon of photonics[1] due to its versatility. Because silicon is the backbone of integrated electronics, [2] growing a LiNbO3 thin film on silicon would be helpful to integrate optic and electric devices on one chip.

Several methods have been reported to fabricate LiNbO3 thin films on silicon, such as molecular beam epitaxy (MBE), pulsed laser deposition (PLD), and liquid phase epitaxy (LPE).[3, 4] Since the crystal structure of LiNbO3 is a distorted perovskite, and the lattice parameters are a = b = 5.15  Å , c = 13.863  Å , while silicon belongs to the diamond type, and the lattice parameter is a = 5.431  Å . The lattice mismatch, which is 5.2%, makes it difficult to grow high quality LiNbO3 thin films on the Si substrate. In experiment, the LiNbO3 film quality is sensitive to many parameters, such as temperature, substrate orientation, substrate morphology, target, and chamber atmosphere.[57] So far, the ideal growth conditions of high quality LiNbO3 crystal are still under debate.

Recently, several methods have been used to simulate thin film growth, [8] such as the molecule dynamics, Monte Carlo, and kinetic Monte Carlo methods.[912] These simulations have revealed the microscopic dynamics and the influence of the temperature on assembling particles. However, the calculation is slow and the time scale is short with limited particles considered.

The dissipative particle dynamics (DPD) method was developed by Koelman and Hoogerbrugge, [13, 14] which allows longer-time scale calculation and larger system simulation.[18] The DPD method has been used in studying a surfactant among water moleculars.[15, 16] Ryjkina et al.[15] and Kuo et al.[16] reported that the DPD method can also be used in liquid molecular aggregation.

In the PLD thin film growth process, the relative positions and velocities of the particles satisfy the proper conservation laws and symmetries, and the momentum transferred from the particles also satisfies the linear and angular conservations, which meet the demands of the DPD method. In this work, we apply the DPD method to the solid molecular field, which shows that the PLD process can be reproduced by the DPD technique. As a result, we find the optimal deposition temperature and substrate.

2. Simulation method and experiment

The lithium niobate molecule and the silicon atom are chosen to be the DPD particles to understand the lithium niobate morphology on the silicon substrate. Assume that the lithium niobate molecules are in stoichiometry (Li-Nb-O3) and electric neutrality under ideal target and oxygen atmosphere, and the total number is 1 × 106 in one simulation. The LiNbO3 molecular lattice parameters are a = b = 0.515  nm, c = 1.3863  nm, and the average parameter is ā = (a + b + c)/3 = 0.808  nm. In the two-dimensional simulation of a 100  nm × 100  nm film, the average thickness of the simulated film is h = 0.808 × 106/100/100 = 80.8  nm. Consider the ideal condition which does not include the defects and supposes that the particles at the boundaries are in balance. The thickness is at the same level with the atomic force microscope (AFM) results and the previous results (50  nm).[13] Approximately, the whole progress in PLD is designed to deposit on the silicon surface of 100  nm × 100  nm at different temperatures. The unreconstructed Si(111) and Si(100) surfaces are selected as the deposition surface.

The final position of a DPD particle is determined by the force exerted on it. The time evolution of the process follows Newton’ s law. A DPD particle is affected by the neighbor DPD particles via three non-bonded forces: a conservative force FC, a dissipative force FD, and a random force FR. Therefore, the effective force Fi acting on particle i is given by[14]

where i, j represent the DPD particles of the lithium niobate molecule and silicon atom, respectively. The FC is a repulsive central force with a maximum value aij as shown in Ref.  [15]

where aij, eij, rij, rc are the maximum repulsive parameter, the unit vector, the distance between particles i and j, and the cutoff radius, respectively.

To calculate FC, aij must be calculated. Here the interaction parameter is obtained by considering the mean pair contact interaction energy ⟨ Eij(T)⟩ between lithium niobate molecules and silicon atoms. First, we calculate the mixing energy between two DPD particles i and j according to[16]

where Z′ s are the coordination numbers of each DPD particle pair. The averaged pair interaction energy ⟨ E(T)⟩ is optimized. The electron state is set as the ground state, with default spin and B3LYP based on the density function theory (DFT) using Gaussian software. The basis set is 6-31G and no charge. Therefore, the electron spin is considered and the energy between particles i and j is optimized. These calculated averaged mixing energies are used to obtain the Flory– Huggins interaction parameter χ ij(T), [17]

The χ parameter is used to determine the maximum repulsive force aij between particles i and j, [17]

where the terms are affected by the interaction between the pure DPD particles i (aii = 15kBT).[18] The force and are responsible for the conservation of the total momentum in the system. To calculate and , the parameters are derived from Groot and Warren.[19] The whole programming is realized in Mathematica.

The LiNbO3 thin films are grown on Si(111) and Si(100) at 873  K for 2  h, and in-situ annealed for 1  h by the PLD. The target is LiNbO3 single crystal, the gas pressure is 30  Pa, the excimer laser wavelength is 248  nm, and the pulse width is 25  ns. The thin films’ morphology is investigated by atomic force microscopy (AFM) at room temperature.

3. Results and discussion

The deposition on Si(111) and Si(100) substrates at five different temperatures (723  K, 773  K, 823  K, 873  K, and 923  K) is simulated. Figure  1 shows the simulated and experimental results of lithium niobate deposited on different silicon surfaces at 873  K. Figures  1(a) and 1(c) are both on the Si(111) surface, and figures  1(b) and 1(d) are on Si (100). The simulations and experiments are consistent with each other on the island sizes and heights. As labeled by the red lines, the island edge is chosen to be the lowest line of the simulation picture and AFM results. Based on this, the growth conditions can be optimized to obtain a smooth surface. In the DPD simulation, the Si(111) surface leads to smaller island size of 55– 65  nm, which is similar with the AFM results (40– 90  nm) in Fig.  1(c). On the Si(100) surface, the edges of the simulated LiNbO3 islands cannot be observed, suggesting that the simulated LiNbO3 island is larger than 100  nm on average, as shown in Fig.  1(b), which is consistent with the AFM image in Fig.  1(d). The island heights are almost the same (∼ 50  nm) on Si(111) and Si(100), both experimentally and theoretically. The results show that LiNbO3, in terms of morphology, is more suitable to be deposited on the Si(111) surface compared with the Si(100) one.

The simulation results are meaningful, which indicate that the DPD method is feasible to simulate the thin film growth of complex oxides. The simulation results of LiNbO3 thin films grown on Si(111) substrate at different temperatures are shown in Fig.  2. The smallest island size is found to be optimized (55  nm) at 873  K, while most are between 60  nm and 64  nm at other temperatures. At 873  K, the film can becomes smoother because lower temperatures concentrate LiNbO3 molecules on silicon and higher temperatures reduce the repulsive force in Eq.  (4). Moreover, 873  K was reported as the proper temperature to grow smooth LiNbO3 thin films on Si(111).[7] For the morphology, the simulation result of 873  K is consistent with the reference.

Fig.  1. (a), (b) DPD simulation results for lithium niobate deposited on different silicon surfaces at 873  K; (c), (d) AFM images of LiNbO3 thin films deposited on silicon substrates by LMBE at 873  K: (a) and (c) on Si(111); (b) and (d) on Si(100).

Fig.  2. The DPD simulation results for lithium niobate deposited on Si(111) surface at different temperatures: (a) 673  K, (b) 723  K, (c) 773  K, (d) 873  K, (e) 923  K.

4. Summary

The DPD simulation suggests that the Si(111) surface is more suitable to produce smooth LiNbO3 thin films than the Si(100) surface. In addition, the optimal deposition temperature is found to be around 873  K, which is consistent with the experimental results. The results indicate that the DPD method is reliable for the simulation of thin film growth of complex oxides, which provides a new path for calculating the defects and the speed distribution of particles.

1 Kösters M, Sturman B, Werheit P, Haertle D and Buse K 2009 Nat. Photon. 3 510 DOI:10.1038/nphoton.2009.142 [Cited within:1] [JCR: 27.254]
2 Morales A M and Lieber C M 1998 Science 279 208 DOI:10.1126/science.279.5348.208 [Cited within:1]
3 Park Y B, Min B, Vahala K J and Atwater H A 2006 Adv. Mater. 18 12 DOI:10.1002/(ISSN)1521-4095 [Cited within:1] [JCR: 14.829]
4 Baumann R C, Rost T A and Rabson T A 1990 J. Appl. Phys. 68 2989 DOI:10.1063/1.346435 [Cited within:1] [JCR: 0.71]
5 Lee S H, Song T K, Noh T W and Lee J H 1995 Appl. Phys. Lett. 67 43 DOI:10.1063/1.115486 [Cited within:1] [JCR: 3.794]
6 Chaos J A, Gonzalo J, Afonso C N, Perričre J and Garćia-Gonźalez M T 2001 Appl. Phys. A 72 6 DOI:10.1007/s003390000683 [Cited within:1] [JCR: 1.545]
7 Wang X C, Ye Z Z and Zhao B H 2007 Solid State Commun. 142 694 DOI:10.1016/j.ssc.2007.04.036 [Cited within:2] [JCR: 1.534]
8 Rode S, Holscher R, Sanna S, Klassen S, Kobayashi K, Yamada H, Schmidt W G and Kuhnle A 2012 Phys. Rev. B 86 075468 DOI:10.1103/PhysRevB.86.075468 [Cited within:1]
9 Yamamoto R and Onuki A 1998 Phys. Rev. Lett. 81 4915 DOI:10.1103/PhysRevLett.81.4915 [Cited within:1] [JCR: 7.943]
10 Zheng S, Zhu W, Gao Y F, Stocks G M and Zhang Z 2010 Appl. Phys. Lett. 96 071913 DOI:10.1063/1.3310066 [Cited within:1] [JCR: 3.794]
11 Haberland H, Insepov Z and Moseler M 1995 Phys. Rev. B 51 11061 DOI:10.1103/PhysRevB.51.11061 [Cited within:1]
12 Zhang C and Li M 2010 Appl. Phys. Lett. 96 073112 DOI:10.1063/1.3327335 [Cited within:1] [JCR: 3.794]
13 Warrender J M and Aziz M J 2003 Mat. Res. Soc. Symp. Proc. 749w311 [Cited within:2]
14 Hoogerbrugge P J and Koelman J M 1992 Europhys. Lett. 19 155 DOI:10.1209/0295-5075/19/3/001 [Cited within:2] [JCR: 2.26]
15 Koelman J M and Hoogerbrugge P J 1993 Europhys. Lett. 21 363 DOI:10.1209/0295-5075/21/3/018 [Cited within:3] [JCR: 2.26]
16 Ryjkina E, Kuhn H, Rehage H, Felix M and Peggau J 2002 Angew. Chem. Int. Ed. 41 6 [Cited within:3]
17 Kuo M Y, Yang H C, Hua C Y, Chen C L, Mao S Z, Deng F, Wang H H and Du Y R 2004 ChemPhysChem 5 4 DOI:10.1002/cphc.200301036 [Cited within:2] [JCR: 3.349]
18 Espanol R and Warren P 1995 Europhys. Lett. 30 191 DOI:10.1209/0295-5075/30/4/001 [Cited within:2] [JCR: 2.26]
19 Groot R and Warren P 1997 J. Chem. Phys. 107 4423 DOI:10.1063/1.474784 [Cited within:1] [JCR: 1.578]