# Spin polarization in Lateral two-dimensional Heterostructures

S. Hannan Mouasvi1, and H. Simchi2 1 Electrical Engineering Department,
Islamic Azad University, Central Tehran Branch, Tehran, Iran  2 Department of
Physics, Iran University of Science and Technology, Narmak, Tehran 16844, Iran
<EMAIL_ADDRESS>
###### Abstract

In this work, we study the spin polarization in the $MoS(Se)_{2}-WS(Se)_{2}$
Transition metal dichalcogenide heterostructures by using the non-equilibrium
Green’s function (NEGF) method and a three-band tight-binding model near the
edges of the first Brillouin zone. Although it has been shown that the
structures have no significant spin polarization in a specific range of energy
of electrons, by applying a transverse electric field in the sheet of the
metal atoms, shedding light on the sample, and under a small bias voltage, a
significant spin polarization in the structure could be created. Besides, by
applying a suitable bias voltage between leads and applying the electric field
a noticeable spin polarization can be found even without shedding the light on
the heterostructures.

  * •

November 2020

##  1 Introduction

By decreasing the minimum feature size of electronic devices, the charge-based
ones encounter some limitations including power dissipation and the
manufacturing cost. There are different strategies for overcoming these kinds
of limitations. Beyond the Moore strategy, scientists have tried to substitute
not only novel materials instead of silicon but also the spin degree of
freedom instead of the charge of an electron. After discovering the single
atomic layer of graphite, called graphene [1], due to some graphene
boundaries, scientists tried to discover and use the graphene-like two
dimensional materials as an alternative of the conventional three dimensional
semiconductors for manufacturing electronic devices. Transition metal
dichalcogenides (TMDs) with general formula as $MX_{2}$ (M=Mo, W and X=S, Se)
are one group of them. TMDs have fascinating and unique properties such as
strong spin-orbit coupling (SOC) in the valence band edge, the number of
layers dependent bandgap, higher absorption coefficient per unit thickness as
compared to conventional semiconductors [2], long spin lifetime [3, 4, 5, 6],
and strong excitonic effects [7]. Regarding the strong spin-orbit coupling,
significant conduction and valence band spin splitting, and very long electron
spin lifetimes, TMDs are one of the main candidate for manufacturing the spin-
based devices [8, 9, 10] and many scientific groups are trying to improve the
specifications of TMDs. For instance, Guo et al., have used doping techniques
to enhance the spin splitting in monolayer of $WS_{2}$ [3]. It has been shown
that the spin relaxation time is enhanced by a small in-plane magnetic field
in monolayer of $WSe_{2}$ [11]. Komsa et al., tuned the bandgap and in
consequence the spectral region of $MoS_{2}/MoSe_{2}/MoTe_{2}$ compounds by
alloying [12]. Based on the assertion of another group, the size and direct-
indirect bandgap in the $WX_{2}$ (X=S, Se, and Te) monolayer can be altered by
strain engineering technique [13]. Stability enhancement of $MoS_{2}$ by
changing its phase from 2H to 1T through Li intercalation is the subject of
Ref.[14]. Forming stacked $MX_{2}$ (M=Mo, W and X=S) heterostructures as it
has been asserted in Ref.[15] is a way to expedite electron-hole separation
for photonic applications. Besides, some references have reported a
correlation between the formation of heterostructures and altering the
photocatalytic and optical properties of TMDs [7, 16].
In the process of fabricating electronic devices, a homojunction is created by
doping an impurity in the semiconductor materials whilst a heterojunction, as
the basis of heterostructures, is made by using two semiconductors with
different energy gaps. Generally, TMD heterostructures are divided into two
lateral and vertical groups. At least two layers of different TMD monolayers
are needed to be stacked and form a vertical (planar) heterostructure.
Similarly, If the layers stitch beside each other, a lateral (in-plane)
heterostructure can be formed. In this regard, various groups have been
reported practical achievement of sharp and defect-free lateral
heterostructures via different methods such as lateral epitaxial and thermal
CVD processes [6, 14, 17]. This is while improving the interfacial properties
of the heterostructures by means of various methods is still in progress [6].
A wide range of TMD based applications has been reported so far, e.g. in-plane
transistors, diodes, p-n photodiodes, and complementary metal-oxide-
semiconductor (CMOS) inverters [6]. Furthermore, optoelectronic applications
of van der Waals (vdW) and lateral heterostructures are increasing due to the
rapid interlayer charge transfer [18] and type engineering ability [6]. Heo et
al., have reported a rotation-based interlayer photoexcitation in
$MoS_{2}/WS_{2}$ monolayer stacks due to the in-direct to direct transitions
[19]. Likewise, the enhancement of photoexcited charge carriers separation and
collection through photovoltaic structure under lateral and vdW
heterostructures was carried out by Atwater et al [2].
A fine review on spintronic in two dimensional materials and their
heterostructures can be found in Ref.[20]. It has been shown that the
heterostructures can significantly improve the SOC [20]. Some examples are
graphene with a WS2 substrate, TIs/graphene vdW heterostructure, and MoS2 on
the graphene channel [20]. Zhang et al., have reported high-spin polarization
at room temperature in 2D layers by reducing its carrier lifetime via the
construction of vdW heterostructures. A near unity degree of spin polarization
has been observed in PbI2 layers with the formation of type-I and type-II band
aligned vdW heterostructures with monolayer of $WS_{2}$ and $WSe_{2}$ [21].
Zhou et al., have studied the vdW graphene/hafnene heterostructure with
different stacking configurations and shown that the distinct electronic
distribution and spin-polarized characteristic of the vdW heterostructures
come from the intensity of orbital overlap induced by diverse stacking
configurations [22]. Deyan et al., have reviewed and discussed some of the
distinctive effects observed in ferromagnetic junctions with prominent 2D
crystals such as graphene, hexagonal boron nitride, and transition metal
dichalcogenides and showed how spin interface phenomena at such junctions
affects the observed magnetoresistance in devices [23]. By referring to
another available references like Refs.[24, 25, 26], one can study more about
different applications of TMD heterostructures as well as prefect spin and
valley polarization in some TMD-based monolayers.
Using the NEGF method, spin polarization in lateral $MoX_{2}-WX_{2}$ (X=S, Se)
heterostructures has been studied by our team in this paper. Besides, the
domination of W (Mo)-atom orbitals on the negative (positive) energy range of
electrons for both types of spin degree of freedom, and the main role of
$d_{{z}^{2}}$ orbital of both W- and Mo-atoms in the quantum conductance of
the device under the equilibrium condition (no bias voltage between leads)
have been revealed by our studies. Under equilibrium condition and without
applying a transverse electric field and shedding light on the structure, no
significant spin polarization was seen when the energy of an electron was less
than -0.1 eV. On the other hand, under the non-equilibrium condition plus
applying the field and shedding light on the sample, in a specific range of
electron energy, due to the breaking of the time-reversal symmetry, a
significant spin polarization was seen in both $MoSe_{2}-WSe_{2}$ and
$MoS_{2}-WS_{2}$ lateral heterostructures. In addition, it has been shown, the
$d_{{z}^{2}}$ orbital of edge atoms plays the main role in the quantum
conductance of the device under the cited condition. In consequence, we show
that a significant spin polarization can be seen in some TMD heterostructures
under special conditions and give a general positive answer to the probable
future application of spin degree of freedom in the TMD heterostructures [24,
25, 26], as it has been demonstrated that a heterostructure can improve the
SOC [27]. It is worth noting, up to our best knowledge, the spin polarization
in lateral TMD heterostructures haven’t been studied theoretically and
experimentally. In consequence, we can not compare our results with others.
The structure of this article is as follows. The material parameters and
calculation methods are provided in sections 2 and 3, respectively. Section 4
includes the results and discussion, and in section 5 one can read the
summary.

##  2 Material parameters

The honeycomb configuration of the $MoS(Se)_{2}-WS(Se)_{2}$ TMD
heterostructures which are illustrated in Fig.1 are composed of three sections
including the left lead, right lead, and channel. In our calculations, the
channel area which is divided into two $MoS(Se)_{2}$ and $WS(Se)_{2}$ semi-
infinite nanoribbons at the left and right side appropriately, consists of 72
metal atoms. Note that the mismatch between the two semi-infinite nanoribbons
is negligible [28]. Fig.1 (c) presents the next nearest neighborhood atoms
respect to a central metal atom. The magnitude of each displacement vector,
$\overrightarrow{R}_{i}$(i=1,…,6), in that figure equals to 3.19
$\mathring{A}$ (3.34 $\mathring{A}$) for $MoS(Se)_{2}-WS(Se)_{2}$
heterostructure. It must be emphasized, in the numerical calculations, we only
consider the metal atoms with a triangular lattice structure, and in
consequence, electrons only hop between the metal atoms in six different
directions.

Figure 1: (Color online) Schematic of $MoS(Se)_{2}-WS(Se)_{2}$ TMD
heterostructure $(a)$ top view, $(b)$ left lead side view, $(c)$ next nearest
neighborhood atoms respect to each central metal atom.
$\overrightarrow{R}_{i}$(i=1,…,6) are the displacement vectors. The lattice
constant and the length of each $\overrightarrow{R}_{i}$ vector are equal to
3.19 $\mathring{A}$ (3.34 $\mathring{A}$) in $MoS(Se)_{2}-WS(Se)_{2}$
heterostructure.

##  3 Calculation method

Exploiting a three-band Hamiltonian model [29, 30, 31, 32, 33, 34] and the
NEGF method, besides the spin polarization studies in lateral two- dimensional
$MoS(Se)_{2}-WS(Se)_{2}$ TMD heterostructures, it has been shown, at the two K
and K’ corners of first Brillouin zone (BZ), the Bloch states in monolayer
$MX_{2}$ (M=Mo, W and X=S, Se) mostly consist of $d_{{z}^{2}}$ orbitals for
the conduction band minimum (CBM), and $d_{xy}$ and $d_{{x}^{2}-{y}^{2}}$ for
the valence band maximum (VBM) [5, 35, 36, 37, 38, 39]. Additional information
including the hopping matrix between metal atoms in each displacement
direction and on-site matrix of each metal atom is provided in appendix A. In
order to calculate the spin polarization and spin current which are defined as
Eq.1 and Eq.2 respectively, one should evaluate the transmission probability,
T +(-), of spin up (down) electrons through the NEGF method [40, 41, 42, 43].

$\displaystyle{SP=\frac{T^{+}-T^{-}}{T^{+}+T^{-}}}$ (1)
$\displaystyle{I}^{+(-)}(E,V)=\frac{q}{\pi{}\hslash{}}\int{T}^{+(-)}(E,V)[\mathit{f}_{1}(E,V)-\mathit{f}_{2}(E,V)]dE$
(2)

Figure 2: (Color online) The self-consistent method for calculating the effect
of shedding light on the heterostructure.

In Eq.2, $\hslash{}$, q, E, V, and $\mathit{f}_{1,2}$ refer to the reduced
Planck constant, elementary charge of electron, energy, potential, and Fermi
distribution function of leads, respectively. Near the equilibrium condition,
by shedding the light on the channel, some electrons are excited to higher
energy levels and their energy increases by $\hslash{}\omega{}$ where
$\omega{}$ is the frequency of the incident light. After a while, they come
back to lower energy levels due to losing the same energy. In consequence, we
are calculating the effect of the incident light on the quantum conductance of
a two levels system [40] in which the energy of its electrons is equal to
$E\pm\hslash{}\omega{}$. Using Fermi’s golden rule [44], one can calculate the
electron-photon coupling matrix and consequently the self-energy of incident
light to implement the influence of it on the system. This could be done by
adding the mentioned self-energy to the Green’s function of the electrons
similar to the effect of leads and their self energies (Appendix B). Althuogh
the detailed information about incident light self-energy calculations exists
in Appendix B, a brief explanation is provided in the continuation of this
paragraph. In the first step, the electron and hole correlation functions and
the initial value of carrier’s density are calculated by using the Green’s
function of the channel and self-energies of leads without shedding the light
(zero self-energy of the photon). By means of the coupling matrix and
correlation functions, one can calculate the initial value of photon self-
energy which is used for evaluating the new amount of the Green’s function of
the channel and the self-energy of the photon. Using Poisson’s equation, the
loop is repeated until the convergence is achieved i.e., it is a self-
consistent method (Appendix B) in accordance with Fig.2.

Since the transverse electric field is applied in the plane of metal atoms,
its effect is added to the diagonal elements of the Hamiltonian matrix as a
potential term. The mentioned term is calculated by multiplying the strength
of the field to the y-coordinate of each atom while the lower (upper) edge of
the nanoribbon is considered as the reference line. In addition, the bias
voltage between leads changes the Fermi distribution function at left (right)
lead equal to $\mathit{f}(E+V_{DS}/2)(\mathit{f}(E-V_{DS}/2))$ where $V_{DS}$
is bias voltage. The difference between these two distributions i.e.,
$[\mathit{f}(E+V_{DS}/2)-\mathit{f}(E-V_{DS}/2)]$ should be considered in
calculating the transmission probability [40]. It should be noted, both the
bias voltage and shedding light, place the system in the non-equilibrium
condition. Thus, the self-consistent method should be used for calculating the
transmission probability and spin polarization [40]. It is important to note
that we used the optimal values for energy of photons, intensity of light, and
the intensity of $E_{T}.y$ which are 0.12 eV, $20nW/\mu m^{2}$, and 0.2 eV,
respectively in our mentioned numerical calculations. In other words, we found
that, by altering the mentioned parameters, the spin polarization could be
changed. Therefore, first, we calculated the optimum values for them aimed at
getting the maximum spin polarization, and then the optimal values were
exploited in our calculations.

##  4 Results and discussion

Fig.3 (a) and (b) show the spin up localized density of states (LDOS) for
$MoS_{2}-WS_{2}$ and $MoSe_{2}-WSe_{2}$ TMD heterostructures under equilibrium
condition, respectively. As these figures show, the LDOS of Mo (W)-atom
dominates on the positive (negative) range of energy of electrons. Although by
using density functional theory (DFT) Kang et al., have shown that the highest
occupied molecular orbital (HOMO) belongs to W-atom and the lowest unoccupied
molecular orbital (LUMO) belongs to Mo-atom [45], because of considering the
molecular orbitals of metal atoms, our results differ from their results.
Density of state (DOS) which is the summation of LDOSs, is shown in the dashed
line in Fig.3 (a) and (b) for $MoS_{2}-WS_{2}$ and $MoSe_{2}-WSe_{2}$
heterostructures accordingly. Based on the figures, one can say that the DOS
is zero for -0.1 to 0.2 eV electron energy range in $MoS_{2}-WS_{2}$ and -0.02
to 0.02 eV in $MoSe_{2}-WSe_{2}$ heterostructures. In consequence, it is
expected that a transmission gap will be seen in the curve of the quantum
conductance.

Figure 3: (Color online) Spin up localized density of states $(a)$
$MoS_{2}-WS_{2}$ and $(b)$ $MoSe_{2}-WSe_{2}$, under equilibrium condition,
i.e., without light, electric field, and bias voltage.

Figure 4: (Color online) Quantum conductance of $MoS(Se)_{2}-WS(Se)_{2}$ for
$(a)$ spin up and $(b)$ spin down electrons, uner equilibrium condition.
Figure 5: (Color online) Spin polarization for $MoS(Se)_{2}-WS(Se)_{2}$
without light, electric field, and bias voltage.

To calculate the quantum conductance of $MoS_{2}-WS_{2}$ and
$MoSe_{2}-WSe_{2}$ TMD heterostructures for both spins in accordance with
Fig.4 (a) and (b) respectively, we used the DOS of the system. As the
mentioned figures show, when the energy of an electron,
$E_{e}\geqslant-0.25eV$, a small difference is seen between the conductance of
spin-up electrons and spin-down ones. Correspondingly, a transmission gap
equals to 0.3 (0.04) eV is seen in $MoS_{2}-WS_{2}$ $(MoSe_{2}-WSe_{2})$
heterostructure. This implied that both nanoribbons are narrow bandgap
semiconductors. Based on the curve of spin polarization under equilibrium
condition as Fig.5, for $E_{e}\geqslant-0.25$eV, spin polarization is not
zero, generally. For better understanding the role of atomic orbitals on the
quantum conductance, we show the magnitude of LDOSs for spin-up electrons as
circles on the lattice position of each atom (Fig.6). As the figure shows, the
$d_{{z}^{2}}$ orbital of all atoms has the main role in the quantum
conductance not only in W-side but also in Mo-side. According to Fig.6 (a) and
(b), the appeared difference in the quantum conductance curves can be
attributed to the different roles and the radius of blue circles on the
position of W and Mo atoms which are placed at the upper edge of a nanoribbon.
Fig.7 shows the LDOS after shedding the light on the sample. Comparing Fig.7
with Fig.6, one can conclude that the LDOS decreases due to the occupation of
states by excited electrons, probably. It should be noted, in order to have a
good view of LDOS changes in the presence of light as well as make the
mentioned comparison possible, the circles’ diameters are illustrated about
four times larger than the calculated amounts.

Figure 6: (Color online) The magnitude of LDOSs for spin up electrons $(a)$
$MoS_{2}-WS_{2}$ and $(b)$ $MoSe_{2}-WSe_{2}$ under equilibrium condition. The
$d_{{z}^{2}}$ orbital is shown in blue color, $d_{xy}$ orbital in red color
and $d_{{x}^{2}-{y}^{2}}$ in green color.

Figure 7: (Color online) The magnitude of LDOSs after shedding the light for
spin up electrons $(a)$ $MoS_{2}-WS_{2}$ and $(b)$ $MoSe_{2}-WSe_{2}$ under
equilibrium condition. The $d_{{z}^{2}}$ orbital is shown in blue color,
$d_{xy}$ orbital in red color and $d_{{x}^{2}-{y}^{2}}$ in green color.

Figure 8: (Color online) Spin current versus bias voltage $(a)$
$MoS_{2}-WS_{2}$ and $(b)$ $MoSe_{2}-WSe_{2}$ nanoribbon. Here, the
$E_{T}=0.2eV$.

Creating spin polarization in a structure could be done by implementing an
external factor such as electric field [46], magnetic field [25, 47, 48],
exchange field [49], or bias voltage [50] and consequently breaking time-
reversal symmetry (TRS). With that mean, first, we apply a bias voltage at a
range from -4 to 4 volts and a transverse electric field equal to 0.2 eV and
calculate the spin-dependent current. As the spin-dependent current versus
bias voltage curve in Fig.8 shows, for some values of bias voltage, the spin-
dependent current is not zero and therefore, the spin polarization has a non-
zero value. A negligible value of the spin polarization for both
$MoS_{2}-WS_{2}$ and $MoSe_{2}-WSe_{2}$ heterostructures in energies from
0.41eV to 0.45eV has been achieved (see Fig.5). Moreover, as presented in
Fig.8, for $-0.1V\leqslant V_{DS}\leqslant 0.1V$, the spin-up current is equal
to spin down current approximately. Now, a question can be asked. Is it
possible under these conditions i.e., $0.41eV\leqslant E_{e}\leqslant 0.45eV$
and $-0.1V\leqslant V_{DS}\leqslant 0.1V$, in which the non-zero spin
polarization changes to a significant spin polarization value under shedding
the light on the channel and applying a transverse electric field ($E_{T}$)?

As illustrated in Fig.9, which is the curve of spin polarization versus the
energy of electrons, when $V_{DS}=0.1V$, $E_{T}.y=0.2eV$, and the channel is
under shedding the light with $20nW/\mu m^{2}$ intensity and 0.12 eV energy, a
significant and almost constant spin polarization is seen for $0.42eV\leqslant
E_{e}\leqslant 0.44eV$. It means that the answer to the above question is
positive when the $MoSe_{2}-WSe_{2}$ heterostructure is considered. Whilst, in
accordance with Fig.10, with shedding light on the channel of $MoS_{2}-WS_{2}$
heterostructure, a non-zero spin polarization can be seen.

To clarify the subject with regard to the influence of spin-orbit coupling in
VBM and CBM degeneracy, spin-valley dependency, and optical selection rule in
$W(M)X_{2}$ monolayers [51], one can say that the light excites the electrons
from non-degenerate VBM to non-degenerate CBM based on the spin-valley
dependency and optical selection rule. In the excitation process by the light
the difference between the density of spin up and spin down electrons changes
and in consequence, it is expected that a non-zero spin polarization is seen.
As Fig.9 shows, when the energy of electron, $E_{e}$, is approximately in a
range from 0.42 to 0.44 eV, the spin polarization in a $MoSe_{2}-WSe_{2}$
heterostructure has a nearly constant positive value (0.6 to 0.8) due to the
domination of spin-up carriers while $E_{T}.y=0.2eV$, $V_{DS}=0.1V$, and light
with the previous specifications exists on the channel. In accordance with the
previous one, the $MoS_{2}-WS_{2}$ TMD heterostructure has no spin
polarization under the same condition while the channel is dark. In other
words, by applying $E_{T}$ and $V_{DS}$ simultaneously and without shedding
the light on the channel, the difference between the densities of spin up and
spin down electrons is negligible, whereas with and without $V_{DS}$ the spin
polarization is created while the light and $E_{T}$ are present. Similarly,
the light excites the electrons based on the spin-valley dependency and
optical selection rule such that the difference between the density of spin-up
electrons and spin-down electrons changes. It should be noted that the perfect
spin polarization is seen with $E_{T}$, shedding the light, and without bias
for $E_{e}=0.444eV$ (Fig.10).

Figure 9: (Color online) Spin polarization in $MoSe_{2}-WSe_{2}$ under light
shedding, 0.2 eV electric field, and various bias voltage. The intensity of
light is $20nW/\mu m^{2}$ and its energy is 0.12 eV. d (l) means dark (light).
Figure 10: (Color online) Spin polarization for $MoS_{2}-WS_{2}$ under light,
0.2eV electric field, and various bias voltage. The intensity of light is
$20nW/\mu m^{2}$ and its energy is 0.12 eV. d (l) means dark (light).

##  5 Summary

Exploiting the NEGF method and three-band Hamiltonian model, we studied the
spin polarization in $MoS_{2}-WS_{2}$ and $MoSe_{2}-WSe_{2}$ lateral
heterostructures. Under the equilibrium conditions, It was shown that non-
significant spin polarization happens when the energy of an electron ($E_{e}$)
is from 0.42 to 0.44 eV, in both heterostructures, whereas by applying a
transverse electric field, bias voltage, and shedding light on the
$MoSe_{2}-WSe_{2}$ structure, a significant spin polarization was seen for
$0.42eV\leqslant E_{e}\leqslant 0.44eV$. Conversely, it was shown that a
significant spin polarization happens only by shedding the light on the bias-
less heterostructure while a transverse electric field exists.

##  Appendix A Tight-binding details

The Bloch states in monolayer $MX_{2}$ (M=Mo, W and X=S, Se) mostly consist of
$d_{{z}^{2}}$ orbitals for conduction band minimum (CBM), and $d_{xy}$ and
$d_{{x}^{2}-{y}^{2}}$ for the valence band maximum (VBM) [5, 35, 36, 37, 38,
39]. Therefore, by considering the base functions $|\Phi_{1}>=d_{{z}^{2}}$,
$|\Phi_{2}>=d_{xy}$, $|\Phi_{3}>=d_{{x}^{2}-{y}^{2}}$ the next nearest tight-
binding Hamiltonian [35] can be written as equation 6.

$\displaystyle H^{NN}(K)=\left(\begin{array}[]{ccc}h_{0}&h_{1}&h_{2}\\\
h_{1}^{*}&h_{11}&h_{12}\\\ h_{2}^{*}&h_{12}^{*}&h_{22}\\\ \end{array}\right)$
(6)

Each element of $H^{NN}$ is a $3\times 3$ matrix that represents the hopping
integrals [52]. These hopping integrals change in the different displacement
directions which are shown in Fig.1 (c). The hopping matrices are as 10 to 18
relations [52].

$\displaystyle h(R_{1})=\left(\begin{array}[]{ccc}t_{0}&t_{1}&t_{2}\\\
-t_{1}&t_{11}&t_{12}\\\ t_{2}&-t_{12}^{*}&t_{22}\\\
\end{array}\right)=h(R_{4})^{\dagger}$ (10) $\displaystyle
h(R_{6})=\left(\begin{array}[]{ccc}t_{0}&\frac{1}{2}t_{1}-\frac{\sqrt{3}}{2}t_{2}&-\frac{\sqrt{3}}{2}t_{1}-\frac{1}{2}t_{2}\\\
-\frac{1}{2}t_{1}-\frac{\sqrt{3}}{2}t_{2}&\frac{1}{4}t_{11}+\frac{3}{4}t_{22}&\frac{\sqrt{3}}{4}t_{11}-t_{12}-\frac{\sqrt{3}}{4}t_{22}\\\
\frac{\sqrt{3}}{2}t_{1}-\frac{1}{2}t_{2}&-\frac{\sqrt{3}}{4}t_{11}+t_{12}+\frac{\sqrt{3}}{4}t_{22}&\frac{3}{4}t_{11}+\frac{1}{4}t_{22}\\\
\end{array}\right)=h(R_{3})^{\dagger}$ (14) $\displaystyle
h(R_{2})=\left(\begin{array}[]{ccc}t_{0}&\frac{1}{2}t_{1}+\frac{\sqrt{3}}{2}t_{2}&\frac{\sqrt{3}}{2}t_{1}-\frac{1}{2}t_{2}\\\
-\frac{1}{2}t_{1}+\frac{\sqrt{3}}{2}t_{2}&\frac{1}{4}t_{11}+\frac{3}{4}t_{22}&-\frac{\sqrt{3}}{4}t_{11}-t_{12}+\frac{\sqrt{3}}{4}t_{22}\\\
-\frac{\sqrt{3}}{2}t_{1}-\frac{1}{2}t_{2}&\frac{\sqrt{3}}{4}t_{11}+t_{12}-\frac{\sqrt{3}}{4}t_{22}&\frac{3}{4}t_{11}+\frac{1}{4}t_{22}\\\
\end{array}\right)=h(R_{5})^{\dagger}$ (18)

Therefore, the three-band TB Hamiltonian will be rewritten as relation 22:

$\displaystyle H=\left(\begin{array}[]{ccc}h_{D}&\cdots&h(R_{n})\\\
\vdots&\ddots&\vdots\\\ h(R_{n})^{\dagger}&\cdots&h_{D}\\\ \end{array}\right)$
(22)

where,

$\displaystyle h_{D}(up)=\left(\begin{array}[]{ccc}\epsilon_{1}&0&0\\\
0&\epsilon_{2}&i\lambda\\\ 0&-i\lambda&\epsilon_{2}\\\ \end{array}\right)$
(26) $\displaystyle h_{D}(down)=\left(\begin{array}[]{ccc}\epsilon_{1}&0&0\\\
0&\epsilon_{2}&-i\lambda\\\ 0&i\lambda&\epsilon_{2}\\\ \end{array}\right)$
(30)

The constants of the above equations are defined in Table 1 [34].

Table 1: GGA based three-band TB model parameters of monolayer MX2 in the unit of eV [34] | $\epsilon_{1}$ | $\epsilon_{2}$ | $t_{0}$ | $t_{1}$ | $t_{2}$ | $t_{11}$ | $t_{12}$ | $t_{22}$ | $\lambda$
---|---|---|---|---|---|---|---|---|---
$MoS_{2}$ | 1.046 | 2.104 | -0.184 | 0.401 | 0.507 | 0.218 | 0.338 | 0.057 | 0.073
$MoSe_{2}$ | 0.919 | 2.065 | -0.188 | 0.317 | 0.456 | 0.211 | 0.290 | 0.130 | 0.091
$WS_{2}$ | 1.130 | 2.275 | -0.206 | 0.567 | 0.536 | 0.286 | 0.384 | -0.061 | 0.211
$WSe_{2}$ | 0.943 | 2.179 | -0.207 | 0.457 | 0.486 | 0.263 | 0.329 | 0.034 | 0.228

##  Appendix B Light interaction

The retarded Green’s function $G^{R}$ is defined as:

$\displaystyle
G^{R}=[(E+i\eta)\mathbf{I}-H_{c}-\Sigma_{l1}-\Sigma_{l2}-\Sigma_{photon}]^{-1}$
(31)

Where $\eta$ equals $1.5\times 10^{-3}$, $\mathbf{I}$ is the unit matrix,
$\Sigma_{l1,2}$ are self-energy of leads, and $\Sigma_{photon}$ is photon
self-energy. The effect of leads on channel appears as shifting the energy
levels and broadening the energy level. It is assumed, due to the energy of
photons, $\hslash{}\omega{}$, the energy of electrons changes as
E+($\hslash{}\omega{}$) and E-($\hslash{}\omega{}$) under shedding the light
on the channel, and in consequence, similar to the leads the effect of photons
can be added to the Hamiltonian of the system as a self-energy. In order to
evaluate the mentioned self-energy, in the first instance, we assumed that the
self-energy of a photon in that equation is zero. This allows us to calculate
the initial value of $\rho$ by:

$\displaystyle\rho=dE[diag(\mathit{f}_{L}*G^{R}*\Gamma_{l1}*G^{A})+diag(\mathit{f}_{R}*G^{R}*\Gamma_{l2}*G^{A})]$
(32)

Where, $\mathit{f}_{L,R}$ and $\Gamma_{l1,2}$ are Fermi distribution function
and broadening matrix of electrodes, respectively. Using equation 33 one can
calculate the preliminary value of photon self-energy.

$\displaystyle\Sigma_{photon}=\frac{i[\Sigma_{photon}^{in}+\Sigma_{photon}^{out}]}{-2}$
(33)

Where $\Sigma_{photon}^{in}$ and $\Sigma_{photon}^{out}$ are defined [42] as:

$\displaystyle\Sigma_{photon}^{in}=\sum_{pq}\quad
M_{lp}*M_{qm}*[NG_{pq}^{n}(E-\hslash\omega)+(N+1)G_{pq}^{n}(E+\hslash\omega)]$
(34) $\displaystyle\Sigma_{photon}^{out}=\sum_{pq}\quad
M_{lp}*M_{qm}*[NG_{pq}^{p}(E+\hslash\omega)+(N+1)G_{pq}^{p}(E-\hslash\omega)]$
(35)

Where $M_{qm}$, N, $G^{n(p)}$, and $\hslash\omega$ are electron-photon
coupling matrix, number of photons, electron (hole) correlation function, and
energy of photons, respectively. With the aim of calculating the proper value
of photon self-energy, we have to calculate the potential of the system
through Poisson’s equation as:

$\displaystyle U=(-\rho)/P$ (36)

Where P is the discretized Poisson operator (appendix C). Finally, we compare
the new potential value with the initial one. If their difference is bigger
than a specific arbitrary value, we guess the new self-energy of photon and
substitute it in equation 31. This self- consistent loop continues until the
mentioned potential difference satisfies the mentioned requirements.

##  Appendix C Poisson’s equation

Poisson’s equation in Cartesian geometry has a form like equation 37.

$\displaystyle\frac{\partial^{2}u}{\partial
x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}+\mathit{f}=0$ (37)

With the aim of changing the Cartesian geometry to Triangular one, three new
vectors called u, v, and w vectors as shown in Fig.11 are defined. It should
be noted, the lattice of metal atoms is triangular and in consequence, the
angles $\alpha$ and $\beta$ are equal to 60 and 120 degrees, respectively.
Therefore, the Eq.37 can be written as [53]:

$\displaystyle\frac{\partial^{2}u}{\partial
x^{2}}+\frac{\partial^{2}u}{\partial
y^{2}}=\frac{2}{3}(\frac{\partial^{2}u}{\partial
u^{2}}+\frac{\partial^{2}u}{\partial v^{2}}+\frac{\partial^{2}u}{\partial
w^{2}})$ (38)

Figure 11: (Color online) The triangular lattice and definition of new vectors
u, v, and w.

Using the finite difference method, the Eq.38 can be discretized and equation
39 will be found [54]:

$\displaystyle-6\Phi_{i,j}+\Phi_{i,j-2}+\Phi_{i,j+2}+\Phi_{i+1,j+1}+\Phi_{i-1,j+1}+\Phi_{i+1,j-1}+\Phi_{i-1,j-1}=Q_{i,j}+O(h^{6})$
(39)

## References

  * [1] Novoselov K S, Geim A K, Morozov S V, Jiang D, Zhang Y, Dubonos S V, Grigorieva I V and Firsov A A 2004 science 306 666–669
  * [2] Jariwala D, Davoyan A R, Wong J and Atwater H A 2017 Acs Photonics 4 2962–2970
  * [3] Guo S, Zheng H, Wang Y and Zhang J 2019 AIP Advances 9 075304
  * [4] Liu W, Wong P K J and Xu Y 2019 Progress in Materials Science 99 27–105
  * [5] Zhu Z Y, Cheng Y C and Schwingenschlögl U 2011 Physical Review B 84 153402
  * [6] Ávalos-Ovando O, Mastrogiuseppe D and Ulloa S E 2019 Journal of Physics: Condensed Matter 31 213001
  * [7] Jing Y, Liu B, Zhu X, Ouyang F, Sun J and Zhou Y 2020 Nanophotonics 1 1675–1694
  * [8] Kośmider K, González J W and Fernández-Rossier J 2013 Physical Review B 88 245436
  * [9] Yang L, Sinitsyn N A, Chen W, Yuan J, Zhang J, Lou J and Crooker S A 2015 Nature Physics 11 830–834
  * [10] Wang Z, Zhao L, Mak K F and Shan J 2017 Nano letters 17 740–746
  * [11] Guimarães M H D and Koopmans B 2018 Phys. Rev. Lett. 120 266801
  * [12] Komsa H P and Krasheninnikov A V 2012 The journal of physical chemistry letters 3 3652–3656
  * [13] Amin B, Kaloni T P and Schwingenschlögl U 2014 Rsc Advances 4 34561–34565
  * [14] Chen P, Zhang Z, Duan X and Duan X 2018 Chemical Society Reviews 47 3129–3151
  * [15] Hong X, Kim J, Shi S F, Zhang Y, Jin C, Sun Y, Tongay S, Wu J, Zhang Y and Wang F 2014 Nature nanotechnology 9 682–686
  * [16] Amin B, Singh N and Schwingenschlögl U 2015 Phys. Rev. B 92 075439
  * [17] Wang J, Li Z, Chen H, Deng G and Niu X 2019 Nano-Micro Letters 11 48
  * [18] Jin C, Ma E Y, Karni O, Regan E C, Wang F and Heinz T F 2018 Nature nanotechnology 13 994–1003
  * [19] Heo H, Sung J H, Cha S, Jang B G, Kim J Y, Jin G, Lee D, Ahn J H, Lee M J, Shim J H et al. 2015 Nature communications 6 1–7
  * [20] Liu Y, Zeng C, Zhong J, Ding J, Wang Z M and Liu Z 2020 Nano-Micro Letters 1–26
  * [21] Zhang D, Liu Y, He M, Zhang A, Chen S, Tong Q, Huang L, Zhou Z, Zheng W, Chen M et al. 2020 Nature communications 11 1–8
  * [22] Zhou B, Wang X, Dong S, Zhang K and Mi W 2017 Carbon 120 121–127
  * [23] Dayen J F, Ray S J, Karis O, Vera-Marun I J and Kamalakar M V 2020 Applied Physics Reviews 7 011303
  * [24] Heshmati-Moulai A, Simchi H, Esmaeilzadeh M and Peeters F 2016 Physical Review B 94 235424
  * [25] Heshmati-Moulai A, Simchi H and Esmaeilzadeh M 2017 The European Physical Journal B 90 128
  * [26] Li M Y, Chen C H, Shi Y and Li L J 2016 Materials Today 19 322–335
  * [27] Liu Y, Zeng C, Zhong J, Ding J, Wang Z M and Liu Z 2020 Nano-Micro Letters 12 1–26
  * [28] Á valos Ovando O, Mastrogiuseppe D and Ulla S E 2019 Physical Review B 99 035107
  * [29] Cappelluti E, Roldán R, Silva-Guillén J, Ordejón P and Guinea F 2013 Physical Review B 88 075409
  * [30] Roldán R, López-Sancho M P, Guinea F, Cappelluti E, Silva-Guillén J Á and Ordejón P 2014 2D Materials 1 034003
  * [31] Ridolfi E, Le D, Rahman T, Mucciolo E and Lewenkopf C 2015 Journal of Physics: Condensed Matter 27 365501
  * [32] Fang S, Carr S, Cazalilla M A and Kaxiras E 2018 Physical Review B 98 075106
  * [33] Choukroun J, Pala M, Fang S, Kaxiras E and Dollfus P 2018 Nanotechnology 30 025201
  * [34] Liu G B, Shan W Y, Yao Y, Yao W and Xiao D 2013 Physical Review B 88 085433
  * [35] Bromley R, Murray R and Yoffe A 1972 Journal of Physics C: Solid State Physics 5 759
  * [36] Mattheiss L 1973 Physical Review B 8 3719
  * [37] Kadantsev E S and Hawrylak P 2012 Solid State Communications 152 909–913
  * [38] Lebegue S and Eriksson O 2009 Physical Review B 79 115409
  * [39] Ataca C, Sahin H and Ciraci S 2012 The Journal of Physical Chemistry C 116 8983–8999
  * [40] Datta S 2005 Quantum transport: atom to transistor (Cambridge university press)
  * [41] Tan S G and Jalil M B 2012 5 - spintronics and spin hall effects in nanoelectronics Introduction to the Physics of Nanoelectronics (Woodhead Publishing) pp 141–197
  * [42] Alam M and Voss P L 2015 Beilstein journal of nanotechnology 6 726–735
  * [43] Paulsson M 2002 Non equilibrium green’s functions for dummies: Introduction to the one particle negf equations (Preprint cond-mat/0210519)
  * [44] Dresselhaus M, Dresselhaus G, Cronin S and Filho A 2018 Solid State Properties: From Bulk to Nano (Springer Berlin Heidelberg)
  * [45] Kang J, Sahin H and Peeters F M 2015 The Journal of Physical Chemistry C 119 9580–9586
  * [46] Wu S, Ross J S, Liu G B, Aivazian G, Jones A, Fei Z, Zhu W, Xiao D, Yao W, Cobden D et al. 2013 Nature Physics 9 149–153
  * [47] Neumann A, Lindlau J, Colombier L, Nutz M, Najmaei S, Lou J, Mohite A D, Yamaguchi H and Högele A 2017 Nature nanotechnology 12 329–334
  * [48] Jelver L, Stradi D, Stokbro K, Olsen T and Jacobsen K W 2019 Physical Review B 99 155420
  * [49] Yang Y, Xu Z, Sheng L, Wang B, Xing D and Sheng D 2011 Physical review letters 107 066602
  * [50] Lee J, Wang Z, Xie H, Mak K F and Shan J 2017 Nature materials 16 887–891
  * [51] Kaviraj B and Sahoo D 2019 RSC Adv. 9 25439–25461
  * [52] Pavlović S and Peeters F 2015 Physical Review B 91 155410
  * [53] Evans D J 1991 International journal of computer mathematics 39 81–98
  * [54] Pickering W and Harley P 1988 Journal of computational physics (Print) 75


# Molecular dynamics simulation of thermal transport
across solid/liquid interface created by meniscus

L. Klochko  V. Mandrolko  Université de Lorraine, CNRS, LEMTA, 54000 Nancy,
France  Faculty of Physics, Taras Shevchenko National University of Kyiv,
64/13, Volodymyrs’ka St., Kyiv 01601,Ukraine  G. Castanet  Université de
Lorraine, CNRS, LEMTA, 54000 Nancy, France  G. Pernot  Université de Lorraine,
CNRS, LEMTA, 54000 Nancy, France  F. Lemoine  Université de Lorraine, CNRS,
LEMTA, 54000 Nancy, France  K. Termentzidis  Université de Lyon, CNRS, INSA-
Lyon, CETHIL UMR5008, F-69621, Villeurbanne, France  D. Lacroix  Université de
Lorraine, CNRS, LEMTA, 54000 Nancy, France  and M. Isaiev  Université de
Lorraine, CNRS, LEMTA, 54000 Nancy, France

###### Abstract

Understandings heat transfer across a solid/liquid interface is important to
develop new pathways to improve thermal management in various energy
applications. One of the important questions that arises in this context is
the impact of three-phase contact line between solid, liquid and gas on the
perturbations of the heat fluxes at the nanoscale. Therefore, this paper is
devoted to the investigations of features of thermal transport across
nanosized meniscus constrained between two solid walls. Different wetting
states of the meniscus were considered with molecular dynamics approach by the
variation of the interactional potential between atoms of the substrate and
the liquid. The effect of the size of the meniscus on the exchange of energy
between two solid walls was also invetigated. It was shown that the presence
of a three phase contact line leads to a decrease of the interfacial boundary
resistance between solid and liquid. Further, investigations with the finite
element method were used to link atomistic simulations with the continuum
mechanics. We demonstrate that the wetting angle and the interfacial boundary
resistance are the required key-parameters to perform multiscale simulations
of such engineering problems with an accurate microscale parametrization.

keywords: thermal transport, molecular dynamics, finite element method,
interfacial thermal resistance, meniscus

##  1 Introduction

Current miniaturization of devices is often related to different thermal
management issues such as overheating, hotspots and thermal interface effects.
The latters constitute bottlenecks for improving their stability, reliability,
and lifetime. Thus, tuning heat transfer at small scales, have a significant
importance for a wide range of applied fields. For instance, thermal transport
across solid/liquid interface is crucial for: energy applications, materials
elaboration, solidification/melting processes, cooling applications, etc. In
addition, such interfaces play an increasingly important role at the nanoscale
when the surface-to-volume ratio becomes significant. As some examples of such
role, the outstanding enhancement of energy transfer and thermal conductivity
in nanofluids [1, 2] or in porous system filled by liquids [3, 4] can be
mentioned. Such remarkable properties of the systems cause significant applied
interest for theirs application in the field of energy harvesting,
dissipation, and storage. However, the morphology features of the contact
between two different species separating a heat source and heat sinks result
in the presence of interfacial thermal resistance (ITR), which lead to the
degradation of the heat dissipation performance [5].

Despite the significant applied and fundamental interest of this issue,
mechanisms that rule the thermal resistance close to a solid/liquid interface
remain unclear. This is particularly the case while considering thermal
transport in nanofluids with significant amount of interfaces. In the latter
case, physic of transport is still under debates. Specifically, numerous
models have been already developed for the description of thermal conductivity
enhancement in nanofluids, however each of the model works well only for
specific nanofluids and concentration level of nanoparticles [6].

Going back to the basics, heat transfer without mass flow, at the macroscopic
level, can be described in the frame of Fourier’s law. In the case of
multiphase systems, a temperature jump ($\Delta T$) appears at the interface
between the two media. The latter is defined in terms of an ITR, $R=\Delta
T/J$, where $J$ is the normal heat flux across the interface.

From a microscopic point of view, the heat transfer from the solid to the
liquid depends on the characteristics of the vibrational modes presented in
the solid phase. For example, C. Ulises Gonzalez-Valle and B. Ramos-Alvarado
[7] showed that out-of-plane modes (perpendicular to the surface) of the
substrate contribute more significantly than in-plane-modes (parallel to the
surface) to the thermal transport across the interface. However, the
dependence of the ITR with the interaction parameters (wetting angle,
penetration length adhesion work) is the genuine difference with the
solid/solid ITR (Kapitza resistance) which mostly depends on bulk properties
of materials [8, 9].

The ITR can significantly affect the thermal transport in systems where there
is an important surface-to-volume fraction of liquid. Its value depends on the
interactions between atoms of solid and fluid in the vicinity of the interface
[10, 11]. In solid/liquid systems, interfacial properties also affect: the
wetting angle and the adhesion strength, while both parameters influence the
heat transfer. For instance, a well-pronounced dependence of the ITR with the
wetting angle was found [12, 13]. Based on the spectral analysis of the heat
flow in silicon-based systems with solid/liquid interface, it was shown that
the scaling law of ITR dependence on wetting angle is not unique [14, 15]
while liquid density depletion layer thickness [7, 16] complies with a more
generic behaviour. However, this parameter is not universal for all materials
[17], and further investigation is required to go deeper in the understanding
of interfacial thermal transport.

The three-phase contact line is the line separating three different phases for
instance solid, liquid and gas. The line often appears in the physical
processes like nucleate boiling [18], meniscus characteristics [18],
crystallization occurance [19] etc. The energy transport close to the three-
phase contact line is important for understanding: the heat transfer in phase
change materials [20, 21], the contact photothermally-induced bubbles grow
[22], or the various physical-based theranostic modalities [23]. It is clear
that the efficiency of energy transfer in this case depends on the physical
properties of all species as well as interactions between them. However, the
features of this interaction is even more poorly understood then the
interfacial thermal transport. Eventually, three-phase contact line can also
significantly contribute to the thermal transport at the nanoscale, especially
where phase change occurs. As an example, it was previously reported [24, 25]
that the contribution of the contact line on the thermal transport in the
nucleation process is crucial.

Thus, in order to describe this complex behaviour more accurately, one needs
to find scalable approaches, as the characteristic size of a nucleation site
of the water is in micron range. Thus, physical insight regarding features of
thermal transport close to the three-phase contact line is essential for
further improvement of the approaches in the energy field.

In this framework, the main goal of this work is to provide a detailed
analysis of the thermal transport close to the triple line. For this, we
considered heat transfer across a meniscus connecting two solid walls as a
model system. Firstly, molecular dynamics (MD) approach was adopted to mine
insight regarding features of energy propagation across the meniscus. The
dependence of the wetting angle as well as on the meniscus size was analyzed.
In a second step, results of MD simulations were compared with ones performed
by finite elements method (FEM) to find appropriate parametrization required
to build links between nanoscale and micro/macroscale approaches.

##  2 Model and Methods

###  2.1 Simulations details

To investigate heat transfer in solid/liquid systems, we performed non-
equilibrium molecular dynamics (NEMD) simulations with LAMMPS [26, 27].
Distinct systems were considered. Specifically, the following systems were
used: i) _“water layer confined between two Si layers”_ (or confined water),
where a slab of water totally fills the gap between two solid walls, and ii)
_a “water meniscus between two Si layers_ ”, where the water domain fills only
the middle part of the gap between the two Si layers, (see Fig. 1). For the
latter we performed simulations for different water drop sizes:
$6\text{a}_{0}$, $10\text{a}_{0}$, and $14\text{a}_{0}$, where
$\text{a}_{0}=0.543$ nm is the lattice parameter for Si. The initial systems
were equilibrated at $T=300$ K, then a temperature gradient set between
$T=270$ K to $T=330$ K was applied in Si regions. More details about the
models and simulations are provided in (SI).

Figure 1:  The initial configuration for: a) confined water (case i), b)
meniscus $6\text{a}_{0}$ (case ii), c) meniscus $10\text{a}_{0}$ (case ii),
and d) meniscus $14\text{a}_{0}$ (case ii). Periodic boundary conditions were
applied in all three dimensions.

The interaction between silicon atoms of the solid subsrate and oxygen atoms
of the water was simulated with a Lenard-Jones potential, given as:

$V(r)=4\varepsilon\left(\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right),$
(1)

where $r$ is the distance between atoms, $\varepsilon$ is the depth of the
potential, $\sigma$ is the distance at which potential is equal to zero. In
all our simulations, we take the distance from the Matthiessen’s mixing rule:

$\sigma=\sigma_{Si-O}=\frac{\sigma_{Si-Si}+\sigma_{O-O}}{2}.$ (2)

The value of the $\varepsilon$ parameter was varied to obtain the different
wetting state. The variation range of $\varepsilon$ was chosen between 10 to
21 meV, with step of 1 meV, that corresponds to a continuous variation of the
nanoscale wetting angle between 126 to 49∘ [28]. The snapshots of water atoms
in the meniscus for all the studied values of potential depth after procedure
of thermalisation are represented in Fig. 2.

Figure 2:  The thermalized configuration of water menisci ($6\text{a}_{0}$,
$10\text{a}_{0}$, and $14\text{a}_{0}$ from top to bottom) for different
parametrizations: a), d), g) $\varepsilon=10$ meV; b), e), h) $\varepsilon=16$
meV; c), f), i) $\varepsilon=21$ meV. Periodic boundary conditions were
applied in all three dimensions.

This range of $\varepsilon$ values allows us to mimic water/substrate
behaviour going from hydrophobic to hydrophilic cases respectively. For each
considered value of the $\varepsilon$ parameter, we performed five independent
NEMD simulations.

###  2.2 Heat flux evaluation

In order to determine the heat flux flowing through the water junction, two
independent methods were used. The first method, named energy balance method
(EBM) [29] is based on the added/extracted energy, $\left<E\right>$, in
thermostats (source/sink), and the heat flux was defined as the amount of
energy passing through a unit area, $A=L_{x}\times L_{y}$, per unit time $t$:

$J=\frac{\left<E\right>}{t\times A},$ (3)

where $L_{x}$, $L_{y}$ are the system sizes of the simulated box.

The second method called all atom contributed method (ACM) for the evaluation
of the heat flux makes possible to decompose different contributions of heat
transport [30, 31] and to consider more precisely the spatial distributions of
heat fluxes:

$\boldsymbol{J}=\frac{1}{V}\left[\sum_{i}e_{i}\boldsymbol{v}_{i}-\sum_{i}\hat{S}_{i}\boldsymbol{v}_{i}\right],$
(4)

where $e_{i}$, $v_{i}$, and $\hat{S}_{i}$ are the per-atom energy (sum of
potential and kinetic), velocity and stress tensor respectively, and $V$ is
the volume of a bin where the calculations were performed. The first term of
Eq. 4 is called convective term and relies to heat transfer due to matter
displacement. The second term, called virial, corresponds to heat transfer due
to stress mechanisms. Knowing the heat flux as well as temperature variations
is mandatory to evaluate ITR as it will shown hereafter.

##  3 Results and discussions

###  3.1 Meniscus morphology

Mass density maps of water menisci with different sizes (case ii) are
presented in Fig. 3. The distribution of the density inside the menisci is
also shown for three different wettability. For small $\varepsilon$
value,10meV, (i.e. large wetting angle, see Fig. 3a, d, g) the expected
hydrophobic case is recovered. With the increase of the $\varepsilon$
parameter the hydrophilicity increases as expected. In addition, we can
observe that the density of water at the edge of the meniscus (see Fig. 3a, d,
g)) is lower for small $\varepsilon$ values (light blue area). Only the
central part of the domain is filled with water with a density close to the
bulk one. Intuitively, this situation will not be favorable to heat transfer.
As $\varepsilon$ increases, the wetting angle decreases as well as the amount
of “low density water”, present on the different edges, for all sizes of
meniscus (Fig. 3). In the latter cases, energy coupling through the water
meniscus is more efficient as it will be demonstrated hereafter.

The shape of the water menisci was deduced from the isodensity line
corresponding to $\rho=0.5\rho_{bulk}$ ($\rho_{bulk}$ is the density of
confined water). The isodensity lines are depicted in the Fig. 3 by black
solid lines. Such lines were compared with the analytical predictions shown by
red lines. For the analytic description of the menisci shapes, we use the
assumptions of: “circular shape” of the water profile, a constrained volume of
the meniscus, and assumes that the wetting angle can be taken the same as for
a droplet with the same $\varepsilon$ value [28] (see details in (SI)). Such
theoretical profiles are also presented in Fig. 3 with the red lines. As one
can see, the excellent agreement between the results of simulations and
analytic modeling can be stated. The latter model will be used to define
system boundaries in FEM simulations. Density maps for all considered values
of wetting angle together with density profiles can be founded in (SI).

Figure 3:  MD density profiles of water menisci ($6\text{a}_{0}$,
$10\text{a}_{0}$, and $14\text{a}_{0}$ from top to bottom) for different
parametrizations: a), d), g) $\varepsilon=10$ meV; b), e), h) $\varepsilon=16$
meV; c), f), i) $\varepsilon=21$ meV. The black solid line connects the points
where the density of meniscus is equal to the half of the density of confined
water $\rho=(0.47-0.53)$ g/cm3. Red line corresponds to the analytic
evaluation of the meniscus shape based on the assumption of its circular
geometry (see (SI)).

###  3.2 Interfacial thermal boundary resistance between the liquid meniscus
and the solid substrate

Knowing the thermal energy going through across the interface is important to
estimate efficacy of heat transport performance. The related power can be
calculated as the energy required for maintaining the temperature difference
in thermostats per unit of time in EBM, $P=\left<dE\right>/dt$.

In the frames of ACM approach the power was calculated from the heat fluxes as
follows:

$P=\int_{\Sigma}\boldsymbol{J}\cdot\boldsymbol{ds},$ (5)

where $\boldsymbol{J}$ is calculated with the Eq. 4, $\boldsymbol{ds}$ is the
unit area taken as $\boldsymbol{ds}=dxdy\cdot\boldsymbol{e_{z}}$ in our
calculations, and $\Sigma$ is the cross-sectional area of the meniscus normal
to $z$ . In order to decrease numerical errors caused by atoms fluctuations,
we averaged the power over different slices in $z$ direction.

The Fig. 4a presents dependence of the exchanged power with the interaction
potential $\varepsilon$ for different initial meniscus sizes as well as for
confined water layer. One can see that the exchanged power shows a linear
trend with the increase of $\varepsilon$. In addition, the rate of the change
is increasing with the meniscus size. Eventually, as it can be seen in Fig. 4a
for the confined water, both EBM and ACM methods gives very similar results
(relative error lower than 3%). Beyond the fact that both approaches can be
used, it also confirms the correctness of the computed heat flux values.

Temperature jump as a function of $\varepsilon$ at the interface is shown in
Fig. 4b. The decrease of the temperature difference can be stated with
increasing of the strength of interaction. This is coherent with a better
water adsorption at Si interface and the occurrence of a higher water density
layer that improves the thermal contact. It is interesting to note, that for
the smaller $\varepsilon$, the temperature difference is almost the same for
all meniscus sizes as well as for confined water layer. In the same time,
there is a visible difference in the value of $\Delta T$ for higher
interactions. This may caused by the capillary nature of the edges of the
meniscus, which allow transverse modes propagate [32], that may significantly
improve transverse phonon interactions with the liquid.

Figure 4:  (a) Power balance and (b) temperature jump as a function of
$\varepsilon$ for case i) (blue line) and case ii) obtained by MD, FEM, and
ACM approaches.

The evaluated exchanged power and temperature jump at the interface give us
the possibility to calculate the absolute thermal resistance ($R_{th}$) at the
interface by:

$R_{th}=\frac{\Delta T}{P}.$ (6)

As can be seen from Fig. 5a, the value of absolute thermal resistance is
decreasing with increasing $\varepsilon$, showing rather similar behavior for
all meniscus sizes. Nevertheless, the difference of the $R_{th}$ is
significant for the considered sizes in the case of the weak interaction, and
this difference gets smaller for the high values of $\varepsilon$.
Additionally, it can be mentioned that while increasing the $R_{\text{th}}$
meniscus size the curve cowerage forward the curve corresponding to the
confined water layer.

The silicon/water contacted area ($A$) as a function of $\varepsilon$ shows
linear variation with increasing interaction parameter $\varepsilon$, and
obvious slope dependence with the size of meniscus. Those results are
represented in Fig. 5b. This information allows us to evaluate the interfacial
thermal resistance (ITR) with the use of the following relation:

$R=R_{th}\cdot A.$ (7)

Comparisons of the ITR values of $R$ vs. $\varepsilon$, both respectively
determined by EBM and ACM methods, are shown in Fig. 6. The computed ITR
values for confined water layer are similar to the literature results which
are typically in the range of $0.4\times 10^{-8}<R<2.3\times 10^{-8}$ m2 K/W,
see Barisik et al [33].

For both situations of confined water layer and menisci, $R$ decreases with
the increase of $\varepsilon$. This is naturally due to the growth of
interaction strength, as the adhesion behavior started to be more and more
hydrophylic. The latter ($\varepsilon=21$meV) creates a better coupling
between the water and the silicon molecules. ITR values for the meniscus (case
ii) are smaller than for confined water (case i). Mechanisms that play an
important role are thus the weak bonding of water molecules at the interface
and edge effects where water with “low density” surrounds the meniscus.
Occurrence of both lead to a degradation of the energy transfer. The
difference between the corresponding $R$ values for the meniscus and confined
water slightly decreases with increasing $\varepsilon$. The wetting angle was
determined based on already estimated values [28]. It should be noted that the
value of $R$ is almost the same for all meniscus sizes. Thus, the decreasing
of the ITR between cases ii) and i) can be connected with heat transfer
channels related to the presence of liquid/vapor interface.

Figure 5: Left graph (a) shows the dependence of the thermal resistance
normalized by the silicon/water contact area on different values of
$\varepsilon$ for: system with confined water layer (case i) and water menisci
with different sizes, $6\text{a}_{0}$, $10\text{a}_{0}$, and $14\text{a}_{0}$
(case ii). The contact area is shown in the right graph (b). Dependence of
silicon/water contact area with $\varepsilon$ for water meniscus with
different sizes, $6\text{a}_{0}$, $10\text{a}_{0}$, and $14\text{a}_{0}$ (case
ii) is shown in the left graph. Figure 6:  Dependence of interfacial thermal
resistance (ITR) with $\varepsilon$ for: system with confined water layer
(case i) and water meniscus with different sizes, $6\text{a}_{0}$,
$10\text{a}_{0}$, and $14\text{a}_{0}$ (case ii).

###  3.3 MD insights for upper scale approaches

MD simulations give interesting insights about ITR variations with the contact
angle at the microscale. However, engineering applications cannot be currently
described by MD tools. To tackle this issue, we propose to use parameters
resulting from the atomic scale modelling as inputs for FEM modelling. To
perform FEM calculations, we used the COMSOL® software (see subsection 1.2 in
SI).

For such simulations, we approximate the shape of meniscus edges with
isodensity lines (see (SI)) as it was mentioned above, which were presented in
Fig. 3 (red lines) and which match the isodensity lines computed in MD (black
lines).

In order to calculate temperature profile with the FEM, we assume that there
is no exchange on the lateral edge of the water domain. The wetting angle [28]
for chosen $\varepsilon$ at the contact points with solid was set as a
boundary condition. To find the position of the edges, we also constrained the
volume of the meniscus. More details about the procedure and resulted figures
are provided in (SI). In COMSOL®, silicon slabs were modelled by two
rectangles with lengths and spacing coinciding with the corresponding
dimensions used in MD. The temperatures of the lower and upper boundaries were
set at 330 and 270 K respectively. Periodic boundary conditions were applied.
As input data, we used ITR values obtained from MD. The thermal conductivity
of silicon and water were chosen equal to 220 and 0.86 W/(m K) [4].

To take into account the non-constant cross-sectional area of the drop, it was
decided to characterize the energy transfer as the total heat flux transferred
from hot to cold sides. The results are shown in Fig. 7. One can see from Fig.
7 that, at small $\varepsilon$ values, the heat flux distributions obtained
with FEM and MD are a bit different, although at $\varepsilon=$21 meV they
almost converge. This trend could be related to ITR results given in Fig. 6.
While $\varepsilon$ is increases, $R$ gets smaller as well as channels for
heat transport are better defined (dense water slab).

Figure 7:  Temperature profiles of the menisci ($6\text{a}_{0}$,
$10\text{a}_{0}$, and $14\text{a}_{0}$ from top to bottom) with heat flux
distributions (black arrows) for different hydrophobicities: a), c), e)
$\varepsilon=$10 meV, b), d), f) $\varepsilon=$21 meV. In each central plot is
shown the one-dimensional distribution of heat fluxes (left side) and
temperatures (right side). The regions where the thermostats with $T=270$ K
and $T=330$ K were applied are shown with blue and red areas respectively.
Brown lines the 2D distributions of isotemperature are shown.

The right part of Fig. 7 shows the temperature profile in the water meniscus
at $\varepsilon=10,16,21$ meV with heat flux distribution represented by black
arrows, and the colorbar displays the temperature range. Those arrows are
given along $x$ and $z$ axes of the system. Results obtained with FEM and MD
show similar trends (see Figs. 14–Figs. 16 in SI). While increasing
$\varepsilon$, exchanged heat flux rises and becomes more structured. Fluxes
obtained with MD exhibit more fluctuations than the one obtained with FEM due
to the inherent variations of MD method. Averaging over more than one
simulation reduces those fluctuations, improving the contact between hot and
cold baths by reducing the ITR, and also gives more clear patterns. As can be
seen from the 1D distributions, the temperature inside the meniscus varies
linearly, and the temperature difference increases with $\varepsilon$. For all
cases, the iso-temperature lines (red lines) are mostly parallel to the
interface.

The temperature jump at the interface as well as total exchanged power as
function of $\varepsilon$ calculated using FEM are presented in Fig. 4 for the
meniscus with the initial size equals to $10a_{0}$. As one can see, there is
an excellent agreement between the results obtained with two different
approaches.

##  4 Conclusions

To conclude, inspired by the impact of the three-phase contact line between
solid, liquid and gas on the thermal transport at the nanoscale, we
investigated features of thermal transport across the liquid meniscus
constrained between two solid walls. Firstly, we considered morphology of the
capillary interface for different interaction strengths between atoms of solid
and liquid and for different meniscus sizes. An analytic model allows us to
describe the shape of the meniscus with respect of wetting angle and volume of
the liquid.

Further, linear trend of the exchanegd power across the meniscus to maintain
temperature difference between the walls as function of interactional
potential depth ($\varepsilon$) was observed for a considered range of wetting
angles varrying between hydrophobic and hydrophilic states. We found that the
temperature jump at the solid/liquid interface is decreasing with increasing
of $\varepsilon$.

Additionally, it was stated that a significant difference exists between
absolute thermal resistance of the smallest meniscus and confined water layer
for the hydrophobic case. The latter becomes small while the system is
reaching the hydrophilic state. With increasing of the meniscus size the
system behavior goes to the same as for the confined water layer.

Regarding the interfacial thermal resistance (ITR) between solid/liquid
interface, it was found that for the meniscus, the value of ITR shows the same
behavior as for the confined water layer – it decreases with increasing of
$\varepsilon$. However, the ITR for the meniscus is systematically lower then
one for the confined water layer. Furthermore, the size dependence of ITR for
all meniscus sizes is almost negligible. It may appear because of arising of
additional possible channel for heat transfer due to capillary nature of
liquid/vapour interface.

Finally, the ITR and wetting angle values evaluated with MD approach were used
to combine atomistic simulation with FEM. It was shown that with such
parametrization the results obtained with both approaches are consistent. This
gives us the possibility to link simulations at different scales.

Acknowledgements

This paper contains the results obtained in the frames of the project
“Hotline” ANR-19-CE09-0003. This work was performed using HPC resources from
GENCI-TGCC and GENCI-IDRIS (Grant 2020-A0080907186), in addition HPC resources
were partially provided by the EXPLOR centre hosted by the Universitté de
Lorraine. The authors thanks the ERASMUS+ for providing mobility grant for V.
Mandrolko.

## References

  * [1] H. M. R. Gonçalves, R. F. P. Pereira, E. Lepleux, L. Pacheco, A. J. M. Valente, A. J. Duarte, V. de Zea Bermudez,  _Small_ 2020, _16_ , 28 1907661.
  * [2] O. Mahian, E. Bellos, C. N. Markides, R. A. Taylor, A. Alagumalai, L. Yang, C. Qin, B. J. Lee, G. Ahmadi, M. R. Safaei, S. Wongwises,  _Nano Energy_ 2021, _86_ 106069\.
  * [3] P. Lishchuk, D. Andrusenko, M. Isaiev, V. Lysenko, R. Burbelo,  _International Journal of Thermophysics_ 2015, _36_ , 9 2428.
  * [4] M. Isaiev, X. Wang, K. Termentzidis, D. Lacroix,  _Applied Physics Letters_ 2020, _117_ , 3 033701.
  * [5] T. Kim, S. Kim, E. Kim, T. Kim, J. Cho, C. Song, S. Baik,  _Small_ 2021, _17_ , 38 2102128.
  * [6] N. Ali, J. A. Teixeira, A. Addali,  _Journal of Nanomaterials_ 2018, _2018_ 1\.
  * [7] C. U. Gonzalez-Valle, B. Ramos-Alvarado,  _International Journal of Heat and Mass Transfer_ 2019, _131_ 645\.
  * [8] Y. Zhang, D. Ma, Y. Zang, X. Wang, N. Yang,  _Frontiers in Energy Research_ 2018, _6_.
  * [9] P. L. Kapitza,  _Phys. Rev._ 1941, _60_ 354\.
  * [10] S. Merabia, S. Shenogin, L. Joly, P. Keblinski, J.-L. Barrat,  _Proceedings of the National Academy of Sciences_ 2009, _106_ , 36 15113.
  * [11] M. Frank, D. Drikakis,  _Entropy_ 2018, _20_ , 5 362.
  * [12] H. Acharya, N. J. Mozdzierz, P. Keblinski, S. Garde,  _Industrial & Engineering Chemistry Research_ 2011, _51_ , 4 1767.
  * [13] N. Shenogina, R. Godawat, P. Keblinski, S. Garde,  _Physical Review Letters_ 2009, _102_ , 15.
  * [14] S. Alosious, S. K. Kannam, S. P. Sathian, B. D. Todd,  _The Journal of Chemical Physics_ 2019, _151_ , 19 194502\.
  * [15] B. Ramos-Alvarado, S. Kumar,  _The Journal of Physical Chemistry C_ 2017, _121_ , 21 11380.
  * [16] C. U. Gonzalez-Valle, L. E. Paniagua-Guerra, B. Ramos-Alvarado,  _The Journal of Physical Chemistry C_ 2019, _123_ , 36 22311.
  * [17] C. U. Gonzalez-Valle, B. Ramos-Alvarado,  _ACS Applied Nano Materials_ 2021, _4_ , 4 3821.
  * [18] D. M. Basavaraja, R. Kumar, A. K. Das,  _Langmuir_ 2021, _37_ , 33 10006.
  * [19] A. Azimi Yancheshme, G. Momen, R. Jafari Aminabadi,  _Advances in Colloid and Interface Science_ 2020, _279_ 102155\.
  * [20] J. Qiu, D. Huo, Y. Xia,  _Advanced Materials_ 2020, _32_ , 25 1.
  * [21] R. Jilte, A. Afzal, S. Panchal,  _Energy_ 2021, _219_ 119564\.
  * [22] X. Yan, J. Xu, Z. Meng, J. Xie, H. Wang,  _Small_ 2020, _16_ , 29 1.
  * [23] D. Li, Y. Yang, D. Li, J. Pan, C. Chu, G. Liu,  _Small_ 2021, _2101976_ 1\.
  * [24] C. Kunkelmann, K. Ibrahem, N. Schweizer, S. Herbert, P. Stephan, T. Gambaryan-Roisman,  _International Journal of Heat and Mass Transfer_ 2012, _55_ , 7 1896.
  * [25] P. A. Raghupathi,  Ph.D. thesis, Kate Gleason College of Engineering, Rochester Institute of Technology,Rochester, New York, 2018.
  * [26] S. Plimpton,  _Journal of Computational Physics_ 1995, _117_ , 1 1\.
  * [27] Large-scale atomic/molecular massively parallel simulator,  https://www.lammps.org.
  * [28] M. Isaiev, S. Burian, L. Bulavin, M. Gradeck, F. Lemoine, K. Termentzidis,  _Molecular Simulation_ 2016, _42_ , 11 910.
  * [29] T. Ikeshoji, B. Hafskjold,  _Molecular Physics_ 1994, _81_ , 2 251.
  * [30] D. Surblys, H. Matsubara, G. Kikugawa, T. Ohara,  _Physical Review E_ 2019, _99_ , 5.
  * [31] P. Boone, H. Babaei, C. E. Wilmer,  _Journal of Chemical Theory and Computation_ 2019, _15_ , 10 5579.
  * [32] L. D. Landau, E. M. Lifshitz,  _Fluid Mechanics, Second Edition: Volume 6 (Course of Theoretical Physics)_ ,  Course of theoretical physics / by L. D. Landau and E. M. Lifshitz, Vol. 6. Butterworth-Heinemann, 2 edition, 1987.
  * [33] M. Barisik, A. Beskok,  _International Journal of Thermal Sciences_ 2014, _77_ 47\.
  * [34] F. H. Stillinger, T. A. Weber,  _Physical Review B_ 1985, _31_ , 8 5262.
  * [35] M. Orsi,  _Molecular Physics_ 2013, _112_ , 11 1566.
  * [36] D. Andrusenko, M. Isaiev, A. Tytarenko, V. Lysenko, R. Burbelo,  _Microporous and Mesoporous Materials_ 2014, _194_ 79\.
  * [37] R. J. Stoner, H. J. Maris,  _Physical Review B_ 1993, _48_ , 22 16373.
  * [38] R. J. Stevens, A. N. Smith, P. M. Norris,  _Journal of Heat Transfer_ 2005, _127_ , 3 315.

Supporting Information

##  1 The model and simulation approaches

###  1.1 Molecular dynamics calculation

Simulation domain has size $30\text{a}_{0}\times 10\text{a}_{0}\times
42\text{a}_{0}$ in $x$,$y$, and $z$ directions respectively, where
$\text{a}_{0}=0.543$ nm is the lattice parameter for Si. Both systems consist
of two silicon slabs with a slab of water. For the case of confined water
(case i), the slab of water covers all silicon surface, and for the meniscus
(case ii) it covers only the central part. All the systems consist of $67800$
atoms of Si in total, but with variety of water molecules: $15028$ for
confined water, $3179$ for meniscus with size $6\text{a}_{0}$, $5201$ for
meniscus with size $10\text{a}_{0}$, and $7225$ for meniscus with size
$14\text{a}_{0}$. The distance between two silicon slabs is set as $10a_{0}$.

Silicon was simulated as diamond lattice with lattice parameters of $0.543$ nm
and its surface is in contact with water with $(0,0,1)$ crystal plane. The
interaction between silicon atoms were modeled via the Stillinger-Weber
potential [34]. Atoms in the bottom and top of silicon slabs ($1.086$ $<$
$z_{1}<$ $2.172$ nm, $20.634$ $<z_{2}<21.72$ nm) were fixed with string forces
to maintain their position close to the initial. The simple sketch that shows
how input parameters are arranged compared to the system is showm in Fig. 8.

Figure 8:  Representation of the initial systems. Here $L_{\text{water}}$ is
the size of the initial systems cases i) $30a_{0}$, and ii) $6a_{0}$,
$10a_{0}$,$14a_{0}$, $d_{\text{water}}=10a_{0}$, $d_{\text{empty}}=2a_{0}$,
$d_{\text{fixed}}=2a_{0}$, $d_{\text{substrate}}=14a_{0}$,
$d_{\text{therm}}=3a_{0}$.

Water was simulated in frame of extended simple point charge model, $SPC/E$
[35]. The interaction potential was determined by a combination of the
Lennard-Jones, and the Coulombic potentials with cutoff distances $r_{c}=1.0$
nm for both of them. For long-range Coulomb forces, each atom is assigned a
charge such as $q(H)=+0.4238e$, and $q(O)=-0.8476e$. Usage of selected
parameters give us a neutral charge for water molecules. To keep the O-H bond
length, $l=0.096$ nm, and H-O-H angles, $109.46^{\circ}$, sustainable the
SHAKE algorithm was applied. For O-O interaction the values of
$\varepsilon=6.736$ meV, and $\sigma$ = 0.3166 nm were chosen. In case of O-Si
interaction, epsilon varied within 10 – 21 meV, and $\sigma$ = 0.26305 nm.

Firstly, the initial system was tempered at temperature $T=1$ K by using
Gaussian distribution for velocity field. Then the system evolved in the $NVE$
ensemble with a gradual temperature rescaling up to 300 K for 150 ps.
(Temperature rescaling to 10 K was applied for 50 ps, to 50 K for 50 ps and
finally to 300 K for another 50 ps). Further the systems were equilibrated in
$NVT$ ensemble for 150 ps.

After equilibration, temperature rescales to 270 K and 330 K were applied to
the regions 2.172 < $z$ < 3.801 nm and 19.005 < $z$ < 20.634 nm respectively.
Then the simulation was performed for 1 ns. The values of temperature,
density, potential and kinetic energy of particles and stress tensors were
obtained. The obtained values are averaged from all values calculated every
0.01 ps. The binning of the system was 0.2 nm. In the region with water,
calculations were performed with binning 0.1 nm for more precise calculations
of water density.

###  1.2 Usage of FEM procedure

In order to analyze system by using finite element method (FEM) we performed
simulation in COMSOL® software in steady-state mode with convergence criteria
= 0.001 (or relative tolerance). All the boundaries that confined the droplet
and silicon slabs were thermally insulated (except those to which a
temperature gradient from 270 to 330 K was applied). Within simulations a
triangular grid with maximum resolution (or “extremely fine” resolution) was
constructed. The estimated time that we spend for one seed is about $t=2\sim
5$ sec.

##  2 The droplet shape derivation

Curves corresponding to the boundaries of the droplets were obtained based on
the circle approximation (see Fig. 9):

$(x-x_{0})^{2}+(z-z_{0})^{2}=R^{2}_{\text{circle}},$

where $x_{0}$ and $z_{0}$ – coordinates of the center of the circle. If $d$ is
the distance between silicon slabs, and $\phi$ is the wetting angle:

$\displaystyle d/2=R_{\text{circle}}\cos\left(\phi\right),$
$\displaystyle(x-x_{0})^{2}+(z-z_{0})^{2}=R^{2}_{\text{circle}},$
$\displaystyle(x-x_{0})^{2}+(z-z_{0})^{2}=d^{2}/(4\cos^{2}\left(\phi\right)),$
$\displaystyle
x=x_{0}\pm(d^{2}/(4cos^{2}\left(\phi\right)-(z-z_{0})^{2})^{1/2},$

where the roots with “$\pm$” and sign corresponds to a hydrophobic (+) and
hydrophilic (-) cases respectively. $x_{0}$ was calculated based on the fact
that the volume $v$ of the droplet for all cases remains constant.

Figure 9:  A simple sketch that shows connection between curves corresponding
to the boundaries of the droplets and the system, that we used to obtain
curves in Figs. 3, 17\- 19. The input parameters that are shown here are well
described in equation above.

We calculated the volume of the droplet based on the number of water molecules
and its density.

$\frac{v}{2b}=\int_{0}^{x_{0}\pm(d^{2}/(4\cos^{2}\left(\phi\right))-(z-z_{0})^{2})^{1/2}}\int_{-d/2}^{d/2}dxdz,$

$b$ is the thickness of droplet which is not used for curve fitting because we
are considering a 2D system. After calculating the integral, we obtain:

$\displaystyle\frac{v}{2b}=x_{0}d$
$\displaystyle\pm\frac{1}{8}\Big{[}\frac{d^{2}}{\cos^{2}\left(\phi\right)}\cdot\arcsin\left(\frac{\left(2z_{0}+d\right)\cos\left(\phi\right)}{d}\right)-$
$\displaystyle-\frac{d^{2}}{\cos^{2}\left(\phi\right)}\cdot\arcsin\left(\frac{\left(2z_{0}-d\right)\cos\left(\phi\right)}{d}\right)+$
$\displaystyle+\left(d-2z_{0}\right)\sqrt{4z_{0}\left(d-z_{0}\right)+d^{2}\tan^{2}\left(\phi\right)}+$
$\displaystyle+\left(d+2z_{0}\right)\sqrt{-4z_{0}\left(d+z_{0}\right)+d^{2}\tan^{2}\left(\phi\right)}\Big{]}$

Thus, we obtained an analytical equation of the curve bounding the drop based
on the distance between the silicon slabs and the wetting angle.

##  3 Temperature profiles of the meniscus with heat fluxes distributions

The same temperature difference between two silicon slabs (60 K) was applied.
Fig. 11 shows the temperature profiles with heat fluxes obtained using FEM
(left) and MD (right) methods. The central plot shows one-dimensional
temperature distributions in the middle of the system obtained from both of
these approaches (solid green line represents FEM values and purple dots
represents MD values). As it can be seen, the fluxes on the graph for MD are
less structured than for FEM. This is expected due to fluctuations of physical
quantities that are present in the MD simulation. Distribution for MD was
obtained by averaging five simulations, each of them separately gives much
more disordered fluxes. Averaging over more simulations would give the more
directional fluxes with an increased computational costs.

The increase of the density peak near the surface with increasing
$\varepsilon$ is represented in Fig. 10 as well. Moreover, from Fig. 10, it
can be seen how water penetrates deeper into silicon for high $\varepsilon$.
Increasing the $\varepsilon$ creates a stronger coupling between water and
silicon, which increases the depth of water penetration into the silicon
lattice, and the density of the surface layer of water, which in turn reduces
the mismatch in the density of vibration states and provides better heat
transfer.

Figure 10: Density profile of confined water obtained at different values of
$\varepsilon$.

(a)  (b)

(c)  (d)

(e)  (f)

(g)  (h)

Figure 11:  Temperature profiles of the meniscus with size $6\text{a}_{0}$,
($\varepsilon$ = 10 - 21 meV) with heat fluxes distributions obtained from MD
simulations.

(a)  (b)

(c)  (d)

(e)  (f)

(g)  (h)

Figure 12:  Temperature profiles of the meniscus with size $10\text{a}_{0}$,
($\varepsilon$ = 10 - 21 meV) with heat fluxes distributions obtained from MD
simulations.

(a)  (b)

(c)  (d)

(e)  (f)

(g)  (h)

Figure 13:  Temperature profiles of the meniscus with size $14\text{a}_{0}$,
($\varepsilon$ = 10 - 21 meV) with heat fluxes distributions obtained from MD
simulations.

(a)  (b)

(c)  (d)

(e)  (f)

(g)  (h)

Figure 14:  Temperature profiles of the meniscus with size $6\text{a}_{0}$,
($\varepsilon$ = 10 - 21 meV) with heat fluxes distributions obtained from FEM
simulations.

(a)  (b)

(c)  (d)

(e)  (f)

(g)  (h)

Figure 15:  Temperature profiles of the meniscus with size $10\text{a}_{0}$,
($\varepsilon$ = 10 - 21 meV) with heat fluxes distributions obtained from FEM
simulations.

(a)  (b)

(c)  (d)

(e)  (f)

(g)  (h)

Figure 16:  Temperature profiles of the meniscus with size $14\text{a}_{0}$,
($\varepsilon$ = 10 - 21 meV) with heat fluxes distributions obtained from FEM
simulations.

(a)  (b)

(c)  (d)

(e)  (f)

(g)  (h)

Figure 17:  The density profiles of meniscus with size $6\text{a}_{0}$ for
different parametrization, ($\varepsilon$ = 10 - 21 meV). Values of density
are shown in colorbar. The black solid line connects the points where the
density of meniscus is equal to the half of the density of confined water. Red
lines are the analytically obtained borders of droplet.

(a)  (b)  (c)

(d)  (e)  (f)

(g)  (h)  (i)

(j)  (k)  (l)

Figure 18:  The density profiles of meniscus with size $10\text{a}_{0}$ for
different parametrization, ($\varepsilon$ = 10 - 21 meV). Values of density
are shown in colorbar. The black solid line connects the points where the
density of meniscus is equal to the half of the density of confined water. Red
lines are the analytically obtained borders of droplet.

(a)  (b)

(c)  (d)

(e)  (f)

(g)  (h)

Figure 19:  The density profiles of meniscus with size $14\text{a}_{0}$ for
different parametrization, ($\varepsilon$ = 10 - 21 meV). Values of density
are shown in colorbar. The black solid line connects the points where the
density of meniscus is equal to the half of the density of confined water. Red
lines are the analytically obtained borders of droplet.