Thermal Boundary Characteristics of Homo-/Heterogeneous Interfaces

TL;DR

This study uses ReaxFF MD simulations and phenomenological theory to analyze thermal boundary resistance at solid interfaces, revealing how chemical reactivity influences TBR and enabling design of multi-layered structures.

Contribution

It combines ReaxFF MD and phenomenological theory to characterize thermal boundary resistance at reactive and non-reactive interfaces, providing new insights into interface design.

Findings

Reactive interfaces have twice the TBR of non-reactive ones.

A stable amorphous layer forms at reactive interfaces, increasing TBR.

Temperature profiles differ between homogeneous and heterogeneous interfaces.

Abstract

The interface of two solids in contact introduces a thermal boundary resistance (TBR), which is challenging to measure from experiments. Besides, if the interface is reactive, it can form an intermediate recrystallized or amorphous region, and extra influencing phenomena are introduced. Reactive force field Molecular Dynamics (ReaxFF MD) is used to study these interfacial phenomena at the (non-)reactive interface. The non-reactive interfaces are compared using a phenomenological theory (PT), predicting the temperature discontinuity at the interface. By connecting ReaxFF MD and PT we confirm a continuous temperature profile for the homogeneous non-reactive interface and a temperature jump in case of the heterogeneous non-reactive interface. ReaxFF MD is further used to understand the effect of chemical activity of two solids in contact. The selected Si/SiO$_2$ materials showed that the TBR of the reacted interface is two times larger than the non-reactive, going from $1.65\times 10^{-9}$ to $3.38\times 10^{-9}$ m$^2$K/W. This is linked to the formation of an intermediate amorphous layer induced by heating, which remains stable when the system is cooled again. This provides the possibility to design multi-layered structures with a desired TBR.

Figures

Tables

Table 1. Table 1: Comparison of volumetric expansion coefficients ( α v subscript 𝛼 𝑣 \alpha_{v} ) obtained from the present ReaxFF MD simulations, and reported from literature (in parentheses).

Element	${𝜶}_{𝒗}\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟓}}$ (m3/m3K)
Pt	2.0 (2.7 Bolz (1973))
Ni	3.2 (3.9 Kollie (1977))

Table 2. Table 2: Thermal conductivities for different system sizes of Pt and Ni.

System	Conductivity	Conductivity
	Pt (W/mK)	Ni (W/mK)
3$\times$3$\times$24	$4.0\pm 0.4$	$7.7\pm 0.5$
3$\times$3$\times$36	$5.9\pm 0.3$	—
3$\times$3$\times$60	—	$17.2\pm 0.6$
3$\times$3$\times$84	$10.4\pm 0.3$	$21.7\pm 0.6$
3$\times$3$\times$108	$15.0\pm 0.5$	$24.4\pm 0.5$
3$\times$3$\times$132	$16.3\pm 0.3$	—
3$\times$3$\times$156	$19.2\pm 0.5$	—
Extrapolated	$49.8\pm 10.5$	$74.4\pm 9.2$
Literature	71.6 Lide (2003)	90.7 Lide (2003)

Table 3. Table 3: Thermal Boundary Resistance (TBR) of Si/SiO 2 systems, with  T C subscript 𝑇 𝐶 T_{C} and T H subscript 𝑇 𝐻 T_{H} (K) as thermostatted temperature for the heat sink and source respectively.

Case	${𝑻}_{𝑪}$ (K)	${𝑻}_{𝑯}$ (K)	Temperature Jump, $𝚫T$ ($K$)	TBR $\left(\frac{{\text{m}}^{\mathrm{𝟐}}\text{K}}{\text{W}}\right)$
1	100	150	12.8	$1.65\times {10}^{-9}$
2	300	350	9.7	$2.40\times {10}^{-9}$
3	800	850	7.8	$2.34\times {10}^{-9}$
4	1700 $\to$ 100	1700 $\to$ 150	16.6	$3.38\times {10}^{-9}$

Table 4. Table 4: Results of BM-eos fitting to different ReaxFF.

	Reference	${𝑩}_{\mathrm{𝟎}}$	${𝑽}_{\mathrm{𝟎}}$	${𝑬}_{\mathrm{𝟎}}$	${𝒓}^{\mathrm{𝟐}}$	Figure
	Force Field	(GPa)	$($Å${}_{}{}^{\mathrm{𝟑}})$	(kcal/mol)
Pt	Mueller et al. (2010)	240	61.52	$-$532.4	1.0	Figure 2a L1
	Sanz-Navarro et al. (2008)	179	62.91	$-$534.3	1.0	Figure 2a L2
	Nielson et al. (2005)	166	64.52	$-$560.6	1.0	Figure 2a L3
Ni	Mueller et al. (2010)	155	46.96	$-$414.9	1.0	Figure 2b L2
	Nielson et al. (2005)	167	48.21	$-$369.8	1.0	Figure 2b L1
Si	Kulkarni et al. (2012)	144	272.1	$-$1617	0.92	Figure 2c
	Fogarty et al. (2010)	165	273.6	$-$1611	0.88
	Narayanan et al. (2011)	292	252.9	$-$1729	0.95
	Newsome et al. (2012)	216	291.2	$-$1675	0.53
	Nielson et al. (2005)	295	265.7	$-$2244	0.84
	Rahnamoun and van Duin (2014)	235	268.3	$-$2206	0.95
	Zou and Van Duin (2012)	334	255.2	$-$1728	0.95
	Pitman and Van Duin (2012)	289	252.3	$-$1730	0.95
	Castro-Marcano and van Duin (2013)	334	255	$-$1728	0.95
SiO2	Kulkarni et al. (2012)	35	242.1	$-$2128	0.83	Figure 2d
	Fogarty et al. (2010)	500	244.6	$-$1847	0.60
	Narayanan et al. (2011)	33	260.0	$-$1818.6	0.86
	Newsome et al. (2012)	34	238.1	$-$1793	0.98
	Nielson et al. (2005)	47	244.7	$-$1860	0.83
	Rahnamoun and van Duin (2014)	234	247.2	$-$1828	0.78
	Zou and Van Duin (2012)	331	263.4	$-$1841	0.19
	Pitman and Van Duin (2012)	273	258.1	$-$1837	0.16
	Castro-Marcano and van Duin (2013)	331	263.4	$-$1840	0.19

Table 5. Table 5: Thermal conductivity of Pt-systems in vacuum in z-direction, and different sizes in x- and y-direction.

System	Thermal Conductivity (W/mK)
$3\times 3\times 32$	9.7
$5\times 5\times 32$	8.0
$10\times 10\times 32$	10.7

Table 6. Table 6: Lattice constant for Pt and Ni.

System	Lattice [a; b; c] (Å)	Deviation from Literature (%)
Pt - literature Wyckoff (1963)	3.9231; 3.9231; 3.9231	-
Pt - EM ReaxFF	3.9473; 3.9473; 3.9473	+0.6
Pt - EM ReaxFF + vacuum in z-direction	3.9412 ; 3.9412 ; —	+0.5
Pt - interface	3.9192; 3.9192; —	$-$0.1
Ni - literature Swanson et al. (1953)	3.5238; 3.5238; 3.5238	—
Ni - EM ReaxFF	3.6122; 3.6122; 3.6122	+2.5
Ni - EM ReaxFF + vacuum in z-direction	3.6048 ; 3.6048 ; —	+2.3
Ni - interface	3.5273; 3.5273; —	+0.1

Table 7. Table 7: Lattice constant for Si (bc8), and cristobalite SiO 2 .

System	Lattice [a; b; c] (Å)	Deviation From Literature (%)
Si (bc8) - literature Kasper and Richards (1964)	6.636; 6.636; 6.636	—
Si (bc8) - EM ReaxFF	6.4393; 6.4393; 6.4393	$-$2.9
Si (bc8) - EM ReaxFF + vacuum in z-direction	6.4383; 6.4046; —	$-$2.9; $-$43.5
Si (bc8) - interface	6.6273; 6.6273; —	$-$0.1
SiO2 - literature Nieuwenkamp (1935)	4.964; 4.964; 6.920	—
SiO2 - EM ReaxFF	5.0443; 5.0443; 7.0063	+1.0
SiO2 - EM ReaxFF + vacuum in z-direction	5.0443; 5.0443; —	+1.0
SiO2 - interface	4.9705; 4.9705; —	+0.1

Table 8. Table 8: Thermal conductivity of a 3 × \times 3 × \times 132 Platinum structure under different mechanically induced stresses.

Simulation	No Stress	1% Compression	1% Stretched
	k (W/mK)	k (W/mK)	k (W/mK)
1	15.2	14.6	13.4
2	17.8	13.9	14.6
3	16.7	15.3	17.4
4	15.6	16.5	16.5
Average	16.3 $\pm$ 1.2	15.1 $\pm$ 1.1	15.5 $\pm$ 1.8

Equations

$$TBR=R=\Delta T/Q$$

$${\text{n}}_1\cdot {\tau }_{U,1}(P,t)+{\text{n}}_2\cdot {\tau }_{U,2}(P,t)=0$$

$$\alpha_{v}=\dfrac{1}{V}\left(\dfrac{\partial V}{\partial T}\right)_{p}=\Bigg{[}\dfrac{\partial ln(V)}{\partial T}\Bigg{]}_{p}$$

$$k=\frac{1}{3}{C}_{v}v{\lambda }_L$$

$$\frac{1}{{\lambda }_L}=\frac{1}{{\lambda }_{\infty }}+\frac{4}{L}$$

$$\frac{1}{k_L}=\frac{12}{C_{v}v}\frac{1}{L}+\frac{1}{k_{\infty }}$$

$$t_{\text{non-dim}}=t/t_{eq}$$

Full text

\pubvolume

9 \issuenum5 \articlenumber663

\historyReceived: 15 March 2019; Accepted: 19 April 2019; Published: 26 April 2019 \TitleThermal Boundary Characteristics of Homo-/Heterogeneous Interfaces \AuthorKoen Heijmans 1\orcidA, Amar Deep Pathak 1\orcidB, Pablo Solano-López 2, Domenico Giordano 3, Silvia Nedea 1,* and David Smeulders 1 \AuthorNamesKoen Heijmans, Amar Deep Pathak, Silvia Nedea, Pablo Solano-LÓpez, Domenico Giordano, and David Smeulders

\corresCorrespondence: s.v.nedea@tue.nl

1 Introduction

Molecular characteristics of solids in contact play a key role in various fundamental studies related to heat transfer Samvedi and Tomar (2009), mechanical behavior Pilania et al. (2014), micro/nano-fluidics Kim et al. (2015), and catalysis Mueller et al. (2010). Besides, it plays an important role in applications related to semiconductors Majumdar (1993), microelectronics Schutze et al. (2001); Nunzi et al. (2008) and heat-shielding in re-entry vehicles for aerospace applications Giordano et al. (2017). In the latter case, insight into the thermal resistance at an interface of a multilayer structure is crucial for prediction and control of overheating of the thermal protection system used for the re-entry shuttle. The non-equilibrium effects of a hypersonic flow impinging on a solid interface requires detailed investigation of the boundary processes. Mass, momentum, energy transfer and chemical reactions on the interface are critical under these extreme conditions Dias et al. (2016). These processes can change the interfacial properties significantly compared to initial bulk properties. Furthermore, in case of chemical reactions at the interface (as shown in Figure 1), insight into the heat transfer in a small layer of material (a few molecular layers) is required.

The influence of the solid interfaces on thermal properties can be analyzed by the local temperature profiles on a molecular level. Experimental measurement of a temperature profile at the molecular level is extremely challenging. Therefore, computational models can be useful accurate tools to provide thermal insight and to establish the interfacial thermal correlations. In the context of building up a macroscopic theory of gas–surface interactions targeting the hypersonic re-entry flows, Giordano et al. Giordano et al. (2017) have proposed a Phenomenological-Theory (PT) to study heat transfer between two solids in contact. Another tool is Molecular Dynamics (MD) simulations, which has been used before to investigate thermal transport across solid interfaces Merabia and Termentzidis (2012); Mahajan et al. (2011); Deng et al. (2014); Lampin et al. (2012); Chen et al. (2012); Zushi et al. (2015); Schelling et al. (2002), and significant influences of the solid interface on the thermal conduction are reported. A schematic representation of these two methods is shown in Figure 1, with a continuum view for PT, and a molecular view for ReaxFF.

Though temperature profiles at the solid interface are investigated before, these studies and methods focus mainly on non-reactive interfaces. However, chemical activity at the interfaces can influence the thermal behavior of solids in contact. Therefore, reactive force field Molecular Dynamics (ReaxFF van Duin et al. (2001); Senftle et al. (2016)) is used, which is able to capture the chemical reaction and its influence on the surface transformation and temperature profile on molecular level.

In this study, we first consider the characterization of material properties like thermal expansion, thermal conductivity, and elastic properties using ReaxFF, to validate the force field. Thereafter, we build a generic non-reactive system, in which an interface is created based on the same material, for this, we considered two Platinum slabs (homogeneous Pt/Pt interface). Further, we create an interface between two different materials in contact, a non-reactive heterogeneous Platinum-Nickel interface (Pt/Ni). Accordingly, the temperature profiles from the ReaxFF MD simulations are compared with the macroscopic level PT based model of Giordano et al. Giordano et al. (2017). For the computation of the temperature profile with PT, relevant material properties are computed with ReaxFF MD and upscaled, to be used as input in the PT model. Platinum and Nickel were selected because they are well studied, non-reactive, monatomic, and have similar lattice size.

After we analyzed this generic model, we created the reactive heterogeneous interface of Silicon and Silicon-oxide (Si/SO2). This Si/SiO2 interface is of relevance for many applications in the semiconductor industry, as well as aerospace engineering Kulkarni et al. (2012). Because of its relevance, the Si/SiO2 interface is studied numerous times before Van Duin et al. (2003), including the thermal boundary resistance (TBR) based on experiments Hurley et al. (2011), or numerical methods like MD Deng et al. (2014); Lampin et al. (2012); Chen et al. (2012), acoustic and diffusive mismatch models (AMM, DMM) Hu et al. (2001), and phonon wave-package method Deng et al. (2014). With ranging values between 0.4–3.5 $\times 10^{-9}$ m2K/W, depending on the method and composition of the materials. Furthermore, Chen et al. showed a strong correlation between the coupling between the materials and the TBR. These coupling can be directly related to reactions happening at the interface. However, to our knowledge no systematic studies has been done on the influence of a reactive interface on the TBR. Therefore, we used ReaxFF to study the influence of the reactive interface. The Si/SiO2 system is kept at various temperatures, within the ReaxFF simulations, to increase/decrease the chemical activity. Accordingly, the TBR is computed. The TBR is defined as the temperature discontinuity at the interface ($\Delta T$) divided by the heat flux ($Q$) that crosses the interface, see Equation (1). The TBR is often referred as the Kapitza resistance Kapitza (1941), however, we kept the analogy of Peterson et al. Peterson and Anderson (1973)

[TABLE]

2 Methodology & Material Properties Estimation

2.1 Phenomenological-Theory (PT)

To study the heat transfer between two solids in contact, Giordano et al. Giordano et al. (2017) have proposed a Phenomenological-Theory (PT). With the aim of building up a macroscopic theory of gas–surface interactions targeting the hypersonic re-entry flows. They have remarked the lack of a physical principle justifying the standard temperature-continuity boundary conditions as a replacement of temperature-continuity, and have introduced tension continuity:

[TABLE]

where, n1 and n2 are normal unit vectors at a point of contact $P$ and time $t$. This macroscopic theory is founded on momentum conservation and represents a more physically motivated boundary condition. For the mathematical formulation of the phenomenological-theory and further details, readers can refer to the original paper Giordano et al. (2017). This method is used to compare the temperature profile non-reactive interfaces studied with MD. For the input of required material properties, the MD calculated values are used.

2.2 Reactive Force Field Molecular Dynamics

Molecular Dynamics (MD) is a computational method to obtain macroscopic and microscopic properties from approximated trajectories of individual particles. These approximated trajectories, obtained from Newton’s equations of motion, form an ensemble from which macroscopic properties of materials can be obtained Frenkel and Smit (2001). To capture the chemical change during a reaction, Reactive force field (ReaxFF) van Duin et al. (2001) is used. ReaxFF is computationally more expensive than the non-reactive force field, however, it allows bond formation and bond breaking during the simulations, which makes simulations of chemical reactions possible. According to ReaxFF the bond order between a pair of atoms can be obtained directly from the inter-atomic distance, which relation is used to mimic chemical change. The feature of bond formation and breaking allows the user not to give predefined reactions pathways, these should present themselves given the right temperatures and chemical environment. However, the accuracy of this relies directly on the training set and the weights that are used to parameterize the reactive force field. Therefore, we tested several available ReaxFF on their ability to predict relevant material characteristics for our study. ReaxFF is widely used in studying chemical activities at a molecular level Senftle et al. (2016); Pathak et al. (2016), including many Si/SiO2 systems Kim et al. (2014); Suek et al. (2018); Nielson et al. (2005); Kulkarni et al. (2012); Fogarty et al. (2010); Narayanan et al. (2011); Castro-Marcano and van Duin (2013); Pitman and Van Duin (2012); Zou and Van Duin (2012); Newsome et al. (2012); Van Duin et al. (2003).

For the in-silico characterization using ReaxFF MD methodology, we compute the material properties of relatively simple systems of Platinum (Pt), Nickel (Ni), Silicon (Si) and Silicon dioxide (SiO2). Furthermore, the elastic properties, thermal expansion coefficient, and thermal conductivity of the materials are important parameters for the phenomenological-theory, and thus the computed values are used as input parameters to upscale the molecular results up to the macroscopic level.

2.2.1 Force Field Selection

Selection of an appropriate force field is very important for in-silico characterization. The calculations of elastic properties, thermal expansion coefficients and the radial distribution function (RDF), guide as a selection criterion for the appropriate force fields. Supercells of $(5\times 5\times 5)$ Pt, $(5\times 5\times 5)$ Ni, $(3\times 3\times 3)$ Si and $(4\times 4\times 4)$ SiO2 are created, containing approximately 400–800 atoms, with initial volumes of 7547, 5469, 7890, and 10,913 Å3, respectively. Periodic Boundary Conditions (PBC) are applied in all directions. The unit cells of Pt Wyckoff (1963), Ni Swanson et al. (1953), Si Kasper and Richards (1964) and SiO2 Nieuwenkamp (1935) are taken from experimental crystallographic information files.

We have chosen three reactive force fields Mueller et al. (2010); Sanz-Navarro et al. (2008); Nielson et al. (2005) available for Pt and Ni and nine force fields for Si and SiO2 Nielson et al. (2005); Kulkarni et al. (2012); Fogarty et al. (2010); Narayanan et al. (2011); Castro-Marcano and van Duin (2013); Pitman and Van Duin (2012); Zou and Van Duin (2012); Newsome et al. (2012). These force fields are tested by deforming the crystals in the range of 0.86 to 1.16 times their initial volume. The resulting increase in potential energy is shown in Figure 2.

For clarity only the best performing force field is shown for Si and SiO2. The relation between volume and energy is found by integration of the pressure in the third order Birch-Murnaghan equation of state (BM-eos) Murnaghan (1944); Birch (1978); Colonna et al. (2011), this relation is fitted to the results of the deformed crystal. The resulting parameters ($B_{0}$ and $V_{0}$) are given in Appendix A and compared with the literature.

The force field developed by J.E. Mueller et al. Mueller et al. (2010) showed the best results, and was therefore chosen for further use including Pt and Ni. This force field was parameterized for studying hydrocarbon reactions on nickel surfaces. They included the equation of state (EOS) for different Ni bulk structures in the training. Both our and their calculations predicted the EOS, together with the lattice parameters, in close agreement with quantum mechanical calculations. We predicted the equilibrium volume ($V_{0}$) of Pt Wyckoff (1963) and Ni Swanson et al. (1953) unit cells within 1.9% and 7.0% deviations from their experimentally observed crystal structures. The deviations for the bulk modulus ($B_{0}$) of Pt and Ni are 9.8% Holmes et al. (1989) and 16.0% Chen et al. (2000) respectively. Furthermore, Mueller et al. computed cohesive energies in close agreement with experimental values. The force field developed by Kulkarni et al. Kulkarni et al. (2012) was chosen for further use including Si and SiO2. This force field is an extension to include gas–surface reactions between oxygen and silica into an existed force field developed by van Duin et al. Van Duin et al. (2003). This original force field was parametrized to include the chemistry of silicon and silicon oxides, and the interface between these materials. Previous work of Tian et al. Tian et al. (2017) also indicated that this original force field predicts the thermal conductivity of vitreous SiO2 in close agreement with experimental values. The force field of Kulkarni et al. is able to predict the equilibrium volume ($V_{0}$) within 7.5% Kasper and Richards (1964) and 30% Nieuwenkamp (1935) deviations for Si and SiO2, respectively. The bulk modulus ($B_{0}$) is within 22.9% Bernstein et al. (2000) and 5.7% Liu (1993) deviations for Si and SiO2, respectively. Therefore, this force field is selected for the study that includes Si and SiO2. These force fields are chosen for further investigation.

To validate further the applicability of the chosen force fields, we have obtained the Radial Distribution Functions (RDFs) of Pt–Pt, Ni–Ni, Si–Si (Si and SiO2) and Si–O (SiO2) pairs from ReaxFF MD simulations in periodic solid supercells as shown in Figure 3.

The sharp peaks in RDF elucidate the extent of ordering in the supercell, thus representing the solid phase. The locations of the peaks coincide with the position of neighboring atoms (represented by ‘$\ast$’) in the experimentally observed solid crystal Wyckoff (1963); Swanson et al. (1953); Kasper and Richards (1964); Nieuwenkamp (1935). Concluding that the selected force fields Mueller et al. (2010); Kulkarni et al. (2012) can capture the crystalline phase of Pt, Ni, Si, and SiO2.

The volumetric thermal expansion coefficient ($\alpha_{v}$) can be obtained from the slope of the natural logarithm of the volume ($\ln V$) versus imposed temperature ($T$) Mashreghi (2012):

[TABLE]

where $\alpha_{v}$ is the volumetric thermal expansion coefficient at constant pressure. We have varied the temperature over 250–500 K at atmospheric pressure in an NPT ensemble. The thermal expansion coefficients for Pt and Ni computed from ReaxFF MD and existing literature values are given in Table 1. The results on the molecular scale are in reasonable agreement with the bulk experimental value Bolz (1973); Kollie (1977).

2.2.2 Thermal Conductivity with Steady State NEMD Method

The thermal conductivity of a solid can be computed from Steady State Non-Equilibrium Molecular Dynamics (ss-NEMD). This method has been previously used Chantrenne and Barrat (2004); Schelling et al. (2002); Stackhouse and Stixrude (2010); Zhou et al. (2009); Heino and Ristolainen (2003); Tian et al. (2017), and it is based on imposing a temperature gradient over a system to estimate the thermal conductivity. A schematic view of this method is shown in Figure 4.

Two strongly coupled regions (using a Berendsen thermostat with damping constant $\tau$ = 100 fs) are created, one hot zone (red zone, $T_{H}$ = 330 K) and one cold zone (blue zone, $T_{C}$ = 300 K), which act as the heat source and sink, respectively. In between these two zones, there are weakly coupled regions (gray zone, $\tau$ = 105 fs), this damping constant proved to have a negligible low influence on the system. This results in a steady state temperature gradient ($dT/dx$) and an energy flux ($q$), in the weakly coupled regions between heat source and sink.

From the energy flux and the temperature gradient, the thermal conductivity ($k$) can be computed, following Fourier’s law. This intuitive principle makes NEMD well suited to study thermal conductivity of different matter, and investigate the influence of structural defects and solid interfaces Chantrenne and Barrat (2004). According to the kinetic theory, the thermal conductivity ($k$) is related to the mean free path ($\lambda$) of energy carriers (Equation (4)). If the characteristic length of the system is larger than the mean free path of carriers, thermal energy is transferred by multiple collisions. In this diffusive regime, the Fourier law is still valid. In cases when the characteristic length of the system is in the order of the mean free path, the energy carriers may travel ballistically between source and sink. This scattering in the heat source and sink introduces an extra limiting effect on the mean free path, and thus a reducing effect on the conductivity (Equation (4)). Thus, the conductivity equation must be corrected for the enhanced scattering effect Chen (2005):

[TABLE]

where, $C_{v}$ is the heat capacity, $v$ the energy carrier velocity, and $\lambda_{L}$ the corrected mean free path for a system of size $L$. This can be estimated from Matthiessen’s rule Dames and Chen (2004), which states that the corrected resistivity is the sum of the intrinsic scattering and the scattering due to impurity. Thus, the corrected mean free path can be expressed as a combined effect of the mean free path of bulk ($\lambda_{\infty}$) and length of the system ($L$) as:

[TABLE]

In a system with periodic boundary conditions, the average distance for an energy carrier to scatter with the heat source or sink is $L/4$ Schelling et al. (2002). Combining Equations (4) and (5), the thermal conductivity ($k_{L}$) of system size $L$ can be expressed as:

[TABLE]

In this equation, $\frac{1}{k_{\infty}}$ is the thermal conductivity of the bulk material. Bulk thermal conductivity ($k_{\infty}$) can be estimated by extrapolating the effective thermal conductivity obtained for small system sizes ($k_{L}$).

The thermal conductivities of Pt and Ni are computed using the ss-NEMD method for different system sizes. The total length of the systems varied from $3\times 3\times X$ (with $X$ from 24 to 156), including $3\times 3\times 6$ zones as the heat source and sink. The energy flux is taken from the average of the heat flux added and extracted by the two strongly coupled regions. The temperature gradient is computed over the weakly coupled region, and the system is equilibrated up to 1 ns with time steps ($\Delta t$) of 0.25 fs. The by ss-NEMD computed thermal conductivity values are presented in Table 2, and increases with system size for both Pt and Ni systems.

The thermal conductivities are extrapolated for long length by fitting the linear expression between $1/L$ and $1/k$ (Equation (6)) as shown in Figure 5. The fitted lines intersect the Y axis at $1/k_{\text{Pt}}$ = 0.021 mK/W and $1/k_{\text{Ni}}$ = 0.013 mK/W, which results in a bulk thermal conductivity of $k_{\text{Pt}}$ = 49.8 $\pm$ 10.5 W/mK and $k_{\text{Ni}}$ = 74.4 $\pm$ 9.2 W/mK. The small deviations in an individual system result in large deviations in the bulk thermal conductivity, due to the extrapolation Berendsen (2011). The computed thermal conductivities are approximately 35 % and 18 % lower than literature Lide (2003) values for Pt and Ni, respectively. This was expected because the ReaxFF formalism we used does not describe free electrons. The difference with experimental values can be explained by limitations of the ReaxFF method we used. ReaxFF MD is not able to model free electrons, thereby we underestimate the thermal conductivity. However, our aim is to compare models and confirm the temperature jump, not to compute exactly the thermal conductivities of Pt and Ni. There is a ReaxFF expansion including free electrons (e-ReaxFF) developed by Islam et al. Islam et al. (2016). At the moment of writing this e-ReaxFF concept does not include the studied materials, nonetheless, this might be interesting for future research.

The gradients of the extrapolated curve (Figure 5) obtained from ReaxFF MD simulations are 2.1 $\times$ 10*-9* m2K/W and 1.0 $\times$ 10*-9* m2K/W for Pt and Ni, respectively. The gradient should be equal to the theoretical gradient $12/C_{v}v$ (Equation (6)) obtained from the kinetic theory. By assuming a specific heat of $C_{v,\text{Pt}}$ = 29 $\times$ 105 J/m3K, and $C_{v,\text{Ni}}$ = 40 $\times$ 105 J/m3K Lide (2003), the computed velocities of thermal transport carriers are $v_{\text{Pt}}\approx 2~{}\times 10^{3}$ m/s, $v_{\text{Ni}}\approx 3~{}\times~{}10^{3}$ m/s. These values are found to be in agreement with the speed of sound in the lateral direction through Pt and Ni from literature Lide (2003). We also studied the final size effects perpendicular to the heat flow for platinum systems, see Appendix B. However, no finite size effects were observed for perpendicular directions, corresponding to the findings of Zhou et al. Zhou et al. (2009).

2.2.3 Building the Interfacial Molecular System

Schematic diagrams of the studied Pt/Pt, Pt/Ni and Si/SiO2 systems are given in Figure 6. The crystal structures of Pt Wyckoff (1963), Ni Swanson et al. (1953), Si, and SiO2 Nieuwenkamp (1935) are used to build the interfaces. For Si the bc8 form given by Kasper et al. Kasper and Richards (1964) is used, and we used the cristobalite SiO2 of Nieuwenkamp et al., these two where specifically chosen to closely match each-others crystal lattice. The top and bottom sections ($9\times 9\times 6$ Pt, $10\times 10\times 6$ Ni, $3\times 3\times 2$ Si and $4\times 4\times 2$ SiO2), are attached to strongly coupled thermostats ($\tau$ = 100 fs) and act as a heat source and heat sink. The intermediate regions ($9\times 9\times 9$ Pt, $10\times 10\times 10$ Ni, $3\times 3\times 25$ Si and $4\times 4\times 24$ SiO2) are weakly coupled ($\tau$ = 105 fs). The supercells are initially placed at a small distance and approach each other during an energy minimization to form the interface. From these energy minimized (merged) systems, the simulations are started. In the non-reactive systems, the top sections are kept at $T_{H}$ = 330 K and bottom sections are kept at $T_{C}$ = 300 K. In the reactive systems temperature values are varied to trigger a chemical reaction at the interface. To calculate the TBR, all the simulations are done over 1 ns, from which the last 0.75 ns are considered to obtain the data.

When an interface between two different materials is created artificial mechanical stresses are introduced by fitting the different lattice constants in one single periodic box. To restrict this to a minimum we carefully selected the materials, supercells, and orientation to create the interface. Thereby, the deformation of the crystals is limited to 0.1% compared to their literature value. Furthermore, we studied the influence of 1% deformation of Platinum on the thermal conductivity. The thermal conductivity for a compressed, as well as, a stretched crystal was lower, however for both cases within the standard deviation of the original system. The deformation of the crystals, and the study towards the thermal conductivity can be found in Appendix C.

3 Results

ReaxFF MD simulations are carried out to understand the temperature discontinuity across the solid interfaces of homogeneous (Pt/Pt), heterogeneous (Pt/Ni) and heterogeneous reactive materials (Si/SiO2). The computed material properties from the previous section are used to upscale the results from the molecular level to macroscopic phenomenological-theory level. For the thermal conductivity, the extrapolated value is used. A comparison is made between both methods for the non-reactive interfaces.

3.1 Non-Reactive Interfaces

In realistic experimental conditions, the thermostat takes time to set the desired temperature, thus the temperature of the heat source evolves with time. To study this effect on the final temperature profile in ReaxFF MD, we compared two different settings to increase the temperature of the heat source. One with a gradual temperature increase to $T_{H}$, and one with an instantaneously high temperature at $T_{H}$. See Appendix D, and Figure 10, for the result of the comparison between the two cases. We observe that the final temperature profile is almost the same for both cases. Thus in the following cases, we have initialized the temperature of the heat source instantaneously at high temperature (instant \Delta$$T).

For the non-reactive ReaxFF MD interface investigation, the systems given in Figure 6a,b are studied. The temperature profile evolution across the solid interfaces of Pt/Pt and Pt/Ni systems with ReaxFF MD is plotted in gray-scale after every 200 ps, which can be seen in Figure 7a,b. The light-gray to black lines represents respectively the earlier and later time periods.

The equilibrated temperature is represented by the solid brown and green lines, respectively Platinum and Nickel. The temperature profiles developed over time obtained from ReaxFF MD are compared with the temperature profile from phenomenological-theory as shown in Figure 7. The molecular level ReaxFF MD simulations and macroscopic level phenomenological-theory have a different time scale, thus to compare them, a non-dimensional time ($t_{\text{non-dim}}$) is defined as:

[TABLE]

where $t$ is the actual time and $t_{eq}$ is the time assumed the system is in steady state. An steady state time of 0.5 ns is assumed for the MD simulations. The red, blue and gray numbers in the figures represent these different transient $t_{\text{non-dim}}$ states. The temperature profile obtained from the phenomenological-theory evolves slowly with time when compared with the ReaxFF MD simulations, where the transient states are quickly converging and fluctuation around the equilibrium state. In the ReaxFF MD results, a continuous temperature profile is observed at the Pt/Pt interface while a temperature jump (discontinuity) is observed at the Pt/Ni interface as shown in Figure 7a,b, respectively. Similar temperature profiles are also observed from the phenomenological-theory. The ReaxFF MD method shows a temperature jump of approximately 39%, where the phenomenological-theory results in a 55% jump of the imposed temperature difference of 30 K. Since one method is based on molecular level and another one is based on macroscopic level theory, a slight discrepancy in the magnitude of the temperature jump can be expected. These results confirm that the temperature jump is observed at the solid interface between different materials for both molecular level and macroscopic level modeling.

3.2 Reactive Interfaces

At the surface of re-entry vehicles chemical reactions can occur, these reactions contribute to the heating of such vehicles Giordano et al. (2017); Kulkarni et al. (2012). Furthermore, these reactions form a small layer, and influence the heat and mass transport at the surface. To gain more fundamental knowledge of such a surface, we studied a reactive solid Si/SiO2 interface (see Figure 6c,d). The building of the physical system is similar to the two previous systems (Pt/Pt and Pt/Ni). The length between the heat source and sink is approximately 327 Å, and temperatures of the heat source and sink are varied to increase/decrease the chemical reaction at the solid interface Tromp et al. (1985). Four different cases are studied, for the first case (Case 1) the set temperatures for heat source ($T_{H}$) and sink ($T_{C}$) are 150 and 100 K, respectively. For the second case (Case 2), the temperatures are $T_{H}=$ 350 K and $T_{C}$ = 300 K, and for the third case (Case 3) the temperatures are $T_{H}$ = 850 K and $T_{C}$ = 800 K. These cases have an interface temperatures of approximately 125, 325, and 825 K. In the fourth case (Case 4) the complete system was heated to a high temperature (1700 K) for 3.5 ns to create a reactive region, and thereafter, cooled back to $T_{C}$ = 100 K, and $T_{H}$ = 150 K to stop the chemical activity completely again. After the cooling, a new steady state simulation was done at $T_{C}$ = 100 K, and $T_{H}$ = 150 K. From Figure 11 in Appendix E, one can see that the major part of the interface formation takes place within the first nanosecond. In terms the thickness of the interface, as well as, the depth of the oxygen penetration into the silicon surface only little changes are observed after the first nanosecond. Therefore, it was not needed to go for longer simulations to create the interface. This fourth case was chosen to get a distinct comparison between the non-reacted and reacted interface, at the same temperature (Case 1 and 4).

The resulted temperature profiles for Case 1–3 are plotted in Figure 8a. The temperature is made non-dimensional by taking 100–150 K as reference and divide by 100. The thick solid lines are fitted to the data, and extrapolated to the interface, to get the temperature jump. The initial interface is positioned at $z$ = 0 Å in the figures. A temperature jump is observed between Si and SiO2 for all the cases. It reiterates that there is a temperature discontinuity at the reactive heterogeneous solid interface as well. Case 1 (100–150 K) shows a clear jump, with a sharp temperature profile. When the temperature is increased to 300–350 K (case 2), the jump remains, however, it is less sharp. This is caused by some small deformation at the interface induced by the temperature. Moving to even higher temperatures 800–850 K (case 3), not only a temperature jump but also a drop of the temperature profile over the entire system is observed. This suggests a heat sink at the interface, due to energy consumption by the deformation of the crystals at the interface. This deformation has an impact on the heat transfer and results in an intermediate region between the Silicon and Silica crystals of a few ångströms.

To better study the effect and size of the intermediate region, the entire Si/SiO2 system was heated up to 1700 K to increase the reactivity and advance the formation of the intermediate region. Higher temperatures were also tested, however, these resulted in melting of the Silicon crystal, and separation of the two slabs. Lower temperatures would require more time to create a similar intermediate region. The heating process results in a larger intermediate amorphous region, where oxygen diffused up to 5 Å into the Silicon crystal, and deformation of both materials is visible up to 10 Å from the interface. After the heating, the system was cooled back to 100–150 K, this temperature was chosen to stop the chemical activity as far as possible. The final equilibrated temperature profile is compared with the temperature profile of case 1, where no reactive region is had been present at the interface. This comparison shown in Figure 8b, and a closer profile around the interface is shown in Figure 9a–c. From Figure 9c, the thicker interfacial region for the reacted interface can be clearly observed, compared to the non-reacted clean interface (Figure 9b). For the heated interface (case 4), the temperature jump ($\Delta T$) is larger, however, is has become less sharp than case 1 and smoothed over the formed intermediate region.

The thermal boundary resistances (TBR) are given in Table 3. The calculated value for the low temperature (case 1) is in good agreement with Deng et al. Deng et al. (2014), who found a value of $1.48~{}(\pm 0.46)\times 10^{-9}$ m2K/W using NEMD, and $1.37~{}(\pm 0.42)\times 10^{-9}$ m2K/W using phonon wave-package dynamics approach. The higher temperatures (case 2,3) approach the experimental results of Hurley et al. Hurley et al. (2011), who measured a resistance of $2.3\times 10^{-9}$ m2K/W. The reacted interface, which includes an amorphous SiO2 region, is in good agreement with the DMM results of Hu et al.Hu et al. (2001), who found a resistance of $3.5\times 10^{-9}$ m2K/W, for amorphous SiO2 with Si. The temperature jump at the interface for the 300–350 K and 800–850 K temperatures (case 2,3) is smaller, however, the calculated TBR is higher, caused by the deformation of the interface which acts as an extra heat sink. The TBR for the reacted interface (case 4) is more than twice the TBR of the non-reacted clean interface (case 1) at the same temperature, caused by the amorphous Si–SiO2 intermediate region.

4 Conclusions

ReaxFF MD is known to capture the physical and chemical phenomena under various conditions Senftle et al. (2016); Pathak et al. (2016). We have chosen various force fields for Pt/Ni and Si/SiO2 systems, which can mimic their material properties. The selected force fields predict the equilibrium volumes Wyckoff (1963); Swanson et al. (1953); Kasper and Richards (1964); Nieuwenkamp (1935) and bulk modulus Holmes et al. (1989); Chen et al. (2000); Bernstein et al. (2000); Liu (1993) of respectively Pt, Ni, Si and SiO2 in close agreement with experiments/theory. To validate further, thermal expansion and thermal conductivity coefficients of Pt and Ni are estimated using ReaxFF MD. The thermal expansion coefficients are found to be in reasonable agreement with experiments. The thermal conductivity of a solid material is size-dependent on the molecular level. Thus, we have obtained the thermal conductivity of Pt and Ni for various system sizes and extrapolated to a very long length to determine the bulk thermal conductivity.

The elastic and thermal properties, obtained from ReaxFF MD, served as input parameters for the macroscopic level phenomenological-theory (PT) Giordano et al. (2017). Temperature profiles of non-reactive interfaces, obtained from both methods, are compared. In this comparison, we have reported the temperature profiles across a homogeneous (Pt/Pt), and heterogeneous (Pt/Ni) interface. Temperature continuity is observed at the solid homogeneous interface of Pt/Pt. The temperature profile of the molecular level simulation is faster at equilibrium than the phenomenological-theory. The temperature profiles between Pt/Ni has a discontinuity at the interface observed in both molecular and macroscopic level. The temperature jump obtained from the molecular level calculation is 18% lower than the one obtained from PT calculations. The discrepancy between the two models in the temperature jump for Pt/Ni is minimal, and can be explained by the fact that the length- and time-scale for both calculations are different, and/or the length dependence of the thermal conductivity. We can conclude that both models, the molecular level ReaxFF MD simulations and the PT, predict a temperature discontinuity across the solid boundary if the materials are not the same.

The ReaxFF MD methodology can capture chemical reactions, therefore, interesting insights could be obtained for solid pairs which can form a reactive interface. For this purpose, the Si/SiO2 pair was chosen and the heat source and heat sink were varied to increase/decrease the chemical reaction at the interface. Three distinct solids (Si, amorphous reacted Si–SiO2 interface, and SiO2) have been observed. The thermal boundary resistance (TBR) is computed at the Si/SiO2 interface for the different systems, providing us with information of the TBR over interfaces with different chemical activity. It can be concluded that the reacted amorphous region at the interface introduces extra resistivity, compared to the non-reactive clean interface. Showing the opportunity to control the thermal resistivity of a multi-layered system by controlling the interfacial reactive regions.

\authorcontributions

conceptualization, S.N. and D.G.; methodology, S.N. and D.G.; validation, K.H., A.P. and P.S.; investigation, K.H., A.P. and P.S.; resources, D.S.; writing–original draft preparation, K.H., A.P. and S.N.; writing–review and editing, K.H., A.P., S.N., P.S., D.G. and D.S.; supervision, S.N., D.G. and D.S.; funding acquisition, D.S.

\funding

This research received no external funding.

Acknowledgements.

The authors would like to thank the Netherlands Organization for Scientific Research (NWO) for access to the national high performance computing facilities (Dossiernr: 17092 6026). \conflictsofinterestThe authors declare no conflict of interest. \appendixtitlesyes

Appendix A Birch-Murnaghan Equation of State Fitting Si & SiO2

To select the correct force field, a BM-EOS study is done. The results of the fitting are presented in Table 4, with $r^{2}$ as correlation coefficient of the fitting. The results are compared with literature, which gives $V_{0,\text{Pt}}=60.38~{}\text{\AA}^{3}$ Wyckoff (1963), $V_{0,\text{Ni}}=43.76~{}\text{\AA}^{3}$ Swanson et al. (1953), $V_{0,\text{Si}}=292~{}\text{\AA}^{3}$ Kasper and Richards (1964), $V_{0,\text{SiO}_{2}}=171~{}\text{\AA}^{3}$ Nieuwenkamp (1935), $B_{0,\text{Pt}}=266$ GPa Holmes et al. (1989), $B_{0,\text{Ni}}=185$ GPa Chen et al. (2000), $B_{0,\text{Si}}=98$ GPa Hopcroft et al. (2010), and  $B_{0,\text{SiO}_{2}}=36$ GPa Liu (1993). The ReaxFF of Mueller et al.  Mueller et al. (2010) and Kulkarni et al.  Kulkarni et al. (2012) proves to be the best applicable for respectively the Pt/Ni system and the Si/SiO2 system.

Appendix B Finite Size Effects, Perpendicular to the Heat Flow

Platinum systems ($3\times 3\times 32$, $5\times 5\times 32$, $10\times 10\times 32$) were used to study the finite size effects in perpendicular direction to the heat flow. The system are placed in vacuum in z-direction and have periodic boundary conditions in x- and y-direction, the given crystal sizes are including heat source and sink. No final size effects are observed in perpendicular direction to the heat flow.

Appendix C Influence of Mechanical Deformation of Slabs

When a heterogeneous interface with two different materials, and thus different lattice parameters, is created the materials are compressed or stretched to fit both materials within the same periodic box. This introduces extra mechanical stresses in the crystals. To restrict this to a minimum we have chosen the materials and supercells in such a way that these artificial deformations are kept to a minimum. The lattice parameters are given in Tables 6 and 7 for the literature value, after an energy minimization in ReaxFF, after an energy minimization in ReaxFF with a vacuum in z direction, and the size used in this work. The lattice parameters are compared with the literature value and the error is given in the last column, this shows that the values are all within 3%, and there is only 0.1% deformation in this work compared to the lattice from literature.

Influence of Stress on Thermal Conductivity

To form an interface with different materials, the materials are slightly stressed to match each other lattice constants. One of the two materials was slightly compressed, and the other one slightly stretched to form the interface. In the Tables 6 and 7 amount of deformation is shown for the materials used in this work, which are within $\pm 0.1\%$. To gain more knowledge on the effect of these stresses on the heat transport across the material, we computed the thermal conductivity for 3 $\times$ 3 $\times$ 132 Platinum structures with lattices corresponding to the literature value, 1% compressed structures, a 1% stretched structures. The computed thermal conductivities for the different deformation are given in Table 8, and are within each other’s standard deviation. Thereby, we conclude that we can neglect the effect of stress in this study.

Appendix D Comparison of Gradual and Instant Induced Temperature

In realistic experimental conditions, the thermostat takes time to set the desired temperature, thus the temperature of the heat source evolves with time. To investigate the effect of heating the systems on the final temperature distribution in ReaxFF-MD, we have compared a gradual temperature increase to $T_{H}$ and an instantaneously temperature at $T_{H}$. For the gradual temperature setting, we have increased the temperature of the heat source ($T_{H}$ = 330 K) in the steps of 5 K per 0.1 million iterations (25 ps). After 0.6 million iterations (150 ps), the temperature profile of the gradual temperature rise system is compared with the system in which temperature of hot zone was instantaneously set at $T_{H}$ = 330 K (instant \Delta$$T) as shown in Figure 10, the comparison was done over a total range of 1 million iterations. We observe that the equilibrated temperature profile is almost the same for both cases. Thus in the following cases, we have initialized the temperature of the heat source instantaneously at high temperature (instant \Delta$$T).

Appendix E Development of the Amorphous Si/SiO2 Interface

To create a reactive formed interface between Si and SiO2, the system was heated for to high temperatures (1700 K) for 3.5 ns. Thereafter, it was cooled back to $T_{C}$ = 100 K, and $T_{H}$ = 150 K to stop the chemical activity completely again. The development of the interface over the first 3 ns can be observed in the snapshots in Figure 11, also the number of oxygen atoms are counted for these snapshots.

\reftitle

References