# Giant piezoelectricity in group IV monochalcogenides with ferroelectric AA
layer stacking

Seungjun Lee  Department of Electrical and Computer Engineering, University of
Minnesota, Minneapolis, Minnesota 55455, USA  Hyeong-Ryul Kim  Department of
Physics and Research Institute for Basic Sciences, Kyung Hee University, Seoul
02447, Korea  Wei Jiang  Department of Electrical and Computer Engineering,
University of Minnesota, Minneapolis, Minnesota 55455, USA  Young-Kyun Kwon
<EMAIL_ADDRESS>Department of Physics and Research Institute for Basic
Sciences, Kyung Hee University, Seoul 02447, Korea  Department of Information
Display, Kyung Hee University, Seoul 02447, Korea  Tony Low<EMAIL_ADDRESS>Department of Electrical and Computer Engineering, University of Minnesota,
Minneapolis, Minnesota 55455, USA

###### Abstract

The piezoelectricity of group IV monochalcogenides (MXs, with M = Ge, Sn and X
= S, Se) has attracted much attention due to their substantially higher
piezoelectric coefficients compared to other 2D materials. However, with
increasing layer number, their piezoelectricity rapidly disappears due to the
antiferroelectric stacking order, severely limiting their practical
applications. Using first-principles calculations, we investigated the
piezoelectricity of MXs with the ferroelectric AA stacking configuration,
which has recently been stabilized in experiments. We found that AA-stacked
MXs have a ferroelectric ground state with the smallest lattice constant among
other stacking configurations, resulting in a giant piezoelectric coefficient,
which is the first demonstration of a strategy where the piezoelectric
coefficients can increase with the number of layers. This can be attributed to
a strong negative correlation between the lattice constant along the armchair
direction and the piezoelectric coefficient, and spontaneous compressive
strain stabilized in ferroelectric AA stacking configuration.

##  I Introduction

An interesting connection between mechanical energy and electricity in
piezoelectric materials makes them highly valuable in various applications,
such as energy harvesting, sensors, and optoelectronic devices.
wang2006piezoelectric ; wang2007nanopiezotronics ; bao2021piezophototronic
Moreover, two-dimensional piezoelectric materials (2DPMs) are of practical
interest compared to their bulk counterparts because they are robust against
depolarization fields, enabling a reduction of device size.
hinchet2018piezoelectric ; zhang2021piezotronics ; jin2023emerging  Therefore,
many efforts contributed in the last decade have led to the discovery of
various types of monolayer or one-layer (1L) piezoelectric materials such as
hexagonal (h-) group III-IV materials, ares2020piezoelectricity  h-group II
oxides, blonsky2015ab  H-phase transition metal dichalcogenides (TMDs),
zhu2015observation ; alyoruk2015promising  In2Se3, zhou2017out ;
xue2018multidirection  CuInP2S6, liu2016room ; brehm2020tunable  and group IV
monochalcogenides (MXs, with M = Ge, Sn and X = S, Se). Fei2015 ;
gomes2015enhanced ; wang2017two ; higashitarumizu2020purely ; xu2022van  Among
the various candidates, MXs, in particular, have received much attention due
to their piezoelectric coefficients, which surpass those of other two-
dimensional materials by orders of magnitude.

Despite these promising advancements, the practical utilization of 2DPMs has
been limited because ground-state multilayer stacking configurations often
lose their piezoelectric response. li2013probing ; zhu2015observation  The
challenge arises from the interlayer dipole-dipole interactions, which
naturally favor an antiparallel order, resulting in the restoration of spatial
inversion symmetry. Nevertheless, recent advances in materials handling and
various growth techniques have paved the way for the experimental
stabilization of metastable stacking configurations, yoo2019atomic ;
higashitarumizu2020purely  resulting in novel ferroelectric (FE) order in
2DPMs. For example, interlayer twisting woods2021charge ; yasuda2021stacking
or sliding xu2022van ; meng2022sliding ; deb2022cumulative  give rise to
additional out-of-plane FE order, which is not possible in a 1L or minimum
energy stacking configuration. However, despite the fact that in-plane
ferroelectricity is also crucial for piezoelectric response, there is limited
understanding of the effect of interlayer interactions on the in-plane FE
order and the corresponding piezoelectric response.

Through first-principles calculation based on density functional theory (DFT),
we investigated the stacking-dependent piezoelectric responses of MXs. We
first carefully evaluated the piezoelectric coefficients of monolayer MX to
understand the large variations in their piezoelectric coefficients reported
previously. We found that an in-plane lattice constant (LC) along the armchair
direction primarily controls the piezoelectric coefficients of MX due to the
unique puckered structures. We further revealed that, the choice of van der
Waals (vdW) correction is particularly crucial for MXs, as it significantly
affects their in-plane lattice constants not only in multilayer cases but also
in the 1L limit, which is not expected for other 2D materials. Therefore, we
investigated the piezoelectric coefficients of AA and AC stacked MX with
carefully chosen vdW corrections. We found a significant enhancement of the
piezoelectric coefficients in AA stacked MX regardless of the choice of vdW
corrections. The physical origin of this enhancement is a spontaneous
compressive strain emerging in the AA stacking. Our results provide a deeper
understanding of the piezoelectric responses of MXs and pave the way for a
novel approach to optimize the piezoelectricity in 2DPMs through stacking
configuration.

##  II Computational details

The first-principles DFT calculation Kohn1965  was carried out using Vienna ab
initio simulation package (VASP). Kresse1996  We employed the plane wave basis
to expand the electronic wavefunctions with a kinetic energy cutoff of 500 eV.
The projector-augmented wave pseudopotentials Blochl1994 ; Kresse1999  were
used for the valence electrons, and the exchange-correlation (XC) functional
was treated within the generalized gradient approximation of
Perdew–Burke–Ernzerhof (PBE). Perdew1996  A sufficiently large vacuum region
($>20~{}\textrm{\AA}$) was included in the unitcell to avoid any spurious
interlayer interactions. The atomic basis was carefully relaxed until the
Helmann-Feynman force acting on every atom was smaller than 0.001 eV/Å, which
is very crucial to get a converged piezoelectric coefficient. The Monkhorst-
Pack 21$\times$21$\times$1 $k$-mesh was used to sample the Brillouin zone. To
comprehensively investigate the effects of interlayer interactions, we used
various types of vdW interactions including, simple Grimme methods, grimme-d2
; grimme-d3  as well as various non-local vdW-DF functionals. vdwdf ; vdwdf2 ;
optB88 ; revvdwdf2

The spontaneous polarization, $P$, of various MXs were calculated in the
framework of the modern theory of polarization based on the Berry phase.
spaldin2012beginner  Then, the planar elastic stiffness coefficients $C_{ij}$
and piezoelectric strain coefficients $e_{ijk}$ were evaluated as Fei2015

$C_{ij}=\frac{1}{A_{0}}\frac{{\partial}^{2}U}{\partial{\varepsilon}_{ii}\partial{\varepsilon}_{jj}},$
(1) $e_{ijk}=\frac{\partial P_{i}}{\partial{\varepsilon}_{jk}},\\\ $ (2)

where ${\varepsilon}$, $U$, and $A_{0}$ are the strain, the total energy and
the unitcell area. The effective thickness of the 1L MXs was treated as half
of the out-of-plane LC of their bulk counterparts, which was used to evaluate
an effective bulk polarization. Then, the piezoelectric stress coefficients
$d_{ij}$ can be calculated as

$d_{11}=\frac{e_{111}C_{22}-e_{122}C_{12}}{C_{11}C_{22}-C^{2}_{12}},$ (3)
$d_{12}=\frac{e_{122}C_{11}-e_{111}C_{12}}{C_{11}C_{22}-C^{2}_{12}}.$ (4)

##  III Results and discussion

Figure 1:  (a) Top and side views of a typical monolayer (1L) group IV monochalcogenides (MXs with M = Ge or Sn and X = S or Se) crystal structure. (b) Lattice constants, $a$ and $b$, of 1L SnSe obtained from different van der Waals (vdW). (c) Spontaneous polarization ${P}$ contour map of 1L SnSe, calculated by the PBE functional without vdW correction. Equilibrium lattice constants shown in (b) were overlaid for comparison. (d) 2D piezoelectric strain coefficient ($e_{ijk}$) contour map of 1L SnSe calculated by the derivative of (c) with strain along $a$. As a result, the upper left and lower right triangles represent $e_{111}$ and $e_{122}$, respectively. (e) ${P}$ and $e_{ijk}$, and (f) piezoelectric stress coefficients ($d_{ij}$) of 1L SnSe obtained from different vdW functionals.  Table 1: Lattice constants, $a$ and $b$, spontaneous polarization, $P$, piezoelectric strain coefficients, $e_{ijk}$, planar elastic stiffness coefficients, $C_{ij}$, and piezoelectric stress coefficients, $d_{ij}$ of 1L SnSe calculated by various functionals.  vdW |  $a$ (Å) |  $b$ (Å) |  $P$ (pC/m) |  $e_{111}$ (nC/m) |  $e_{122}$ (nC/m) |  $C_{11}$ (N/m) |  $C_{22}$ (N/m) |  $C_{12}$ (N/m) |  $d_{11}$ (pm/V) |  $d_{12}$ (pm/V)
---|---|---|---|---|---|---|---|---|---|---
none (PBE) | 4.383 | 4.292 | 195.7 | 3.17 | 0.48 | 21.47 | 44.73 | 17.6 | 205.01 | -69.98
none (PBE) Fei2015  | 4.35 | 4.24 |  | 3.49 | 1.08 | 19.88 | 44.49 | 185.7 | 250.57 | -80.31
none (PBE) guo2020small  | 4.41 | 4.29 |  | 2.81 | 0.52 | 23.06 | 42.82 | 18.89 | 175.32 | -65.11
none (PBE) wang2017two  |  |  | 181 |  |  |  |  |  |  |
D3 | 4.509 | 4.265 | 281.7 | 1.50 | -0.36 | 25.47 | 42.06 | 21.02 | 112.29 | -64.61
D2 | 4.305 | 4.248 | 150.6 | 5.31 | 1.46 | 17.85 | 41.64 | 17.24 | 439.27 | -146.81
optB88 | 4.384 | 4.284 | 192.5 | 2.74 | 0.27 | 25.02 | 45.00 | 20.26 | 164.7 | -68.15
optB88 fang2016layer  | 4.41 | 4.27 |  | 2.35 | 0.76 | 20 | 40.4 | 17 | 158.2 | -47.7
rev-vdW-DF2 | 4.331 | 4.271 | 153.8 | 3.91 | 0.59 | 26.61 | 47.6 | 20.19 | 202.8 | -73.62

Figure 1(a) shows the crystal structure of a typical 1L MX. We set $x$\- and
$y$-axes as in-plane directions with corresponding LCs of $a$ and $b$, and
$z$-axis as out-of-plane direction. It has puckered atomic structures with a
mirror plane of $M_{y}$ and belongs to $Pmn2_{1}$ space group.
gomes2015enhanced  More importantly, its broken $M_{x}$ symmetry leads to a
non-zero spontaneous polarization ${P}$ along $x$ or armchair direction.
Indeed, the spontaneous polarization of 1L MXs can be understood as a
rearrangement of the atomic basis along the armchair direction. It has also
been reported that the LC along the armchair direction plays a significant
role in determining its ${P}$. wang2017two  Therefore, we begin with our
discussion on the LCs of 1L MXs using SnSe as an exemplary material and will
carefully re-examine its piezoelectric coefficients.

Although the LCs of bulk MXs have been experimentally reported,
chattopadhyay1986neutron ; PhysRevB.16.832 ; Murgatroyd2020  those of 1L MXs
have not yet been done but are only available through DFT calculations. In
most previous studies, Fei2015 ; wang2017two ; xu2022van  the LCs of 1L MXs
have been obtained only with the PBE XC functional without vdW corrections,
which were considered only in the multilayer case. This is the most common
approach to determine the LCs of various 2D materials, since the effect of the
vdW correction may be negligible in the 1L case. (See Fig. S1 in Supplementary
Information (SI) for more details.) However, we found that the LC of 1L SnSe
strongly depends on the choice of vdW correction, as shown in Fig. 1(b). For
example, the LC ,$a$, obtained along the armchair direction ranges from 4.509
Å (PBE+D3) to 4.305 Å (PBE+D2), representing a deviation of almost 5%. Such
deviation is attributed to its puckered structure with a weak stiffness along
the armchair direction.

We found that a small change in the LCs of 1L MXs would cause a large
variation in their piezoelectric coefficients. To understand this result, we
calculated ${P}$ of 1L SnSe as a function of the two in-plane lattice
constants a and b using PBE XC, which is shown as a contour map in Fig. 1(c),
where other calculated LCs using other functions are superimposed. The contour
map is symmetric with respect to the $b=a$ line with $C_{4z}$ rotational
symmetry due to the degenerate ground states for $a>b$ and $a<b$. As clearly
visualized in Fig. 1(c), $P$ is positively correlated with $a$, and the aspect
ratio between $a$ and $b$ also has some secondary effects. As the LC
decreases, 1L SnSe loses its $P$ and becomes paraelectric. Therefore, the vdW
correction, which yields a smaller (larger) LC along the armchair direction,
may underestimate (overestimate) ${P}$. For example, the $P$ value of 281.7
pC/m estimated with the PBE+D3 correction is almost twice the value of 150.6
pC/m estimated with the PBE+D2 correction. Figure 1(d) shows the contour map
of the 2D piezoelectric strain coefficients, $e_{ijk}$, computed by taking the
first derivative of $P$ shown in Fig. 1(c) with respect to strain
$\varepsilon_{jk}$, as explained in Computational details and Fig S2 in SI. We
found significantly stronger $e_{111}$ near the boundary between the
paraelectric and FE phases due to the rapid emergence of non-zero ${P}$, and a
small enhancement in $e_{122}$, although much smaller than $e_{111}$. This
result clearly indicates that a compressive strain along the armchair
direction may play a major role in the enhancement of $e_{111}$ of SnSe.
Simultaneously, it also implies that the fluctuation of the LC according to
the choice of vdW functionals makes a large uncertainty in predicting its
piezoelectric coefficients.

For a more systematic study, we evaluated the piezoelectric coefficients of 1L
SnSe using various vdW functionals, as summarized in Fig. 1(e, f), and Table
1. Due to the strong correlation between ${P}$ and $a$, both the strain and
stress piezoelectric coefficients ($e_{ijk}$ and $d_{ij}$) also strongly
depend on the choice of functional. For example, $e_{111}$ and $d_{11}$ range
from from 5.31 nC/m to 1.50 nC/m, and from 439.27 pm/V to 112.29 pm/V,
respectively. As predicted above, the smaller $a$ leads to the higher
$e_{ijk}$ and $d_{ij}$. We further confirmed that our results are consistent
with previously reported studies, and the observed variation in both $e_{ijk}$
and $d_{ij}$ is simply understood as the result of variation in the relaxed
$a$. Therefore, we emphasize that the choice of computational options is
extremely important in the quantitative analysis and the vdW functional should
be used consistently in both 1L and multilayer cases.

Figure 2: Calculated lattice constants of bulk group IV monochalcogenides,
obtained using various van der Waals corrections (red dots). Experimental
values obtained from Ref. [40; 41; 42] are also shown as red solid lines for
comparison.

A careful investigation in the previous section indicates that there is a
significant uncertainty in determining the exact value of the piezoelectric
coefficients of 1L MXs without experimental input of LCs. Therefore, before
investigating the piezoelectric coefficients of multilayer MXs, we calculated
the LCs of bulk MXs, and then compared them with experimentally available
values to validate which vdW correction is most appropriate among various vdW
functionals. Figure 2 shows the LCs of four different bulk MXs calculated with
seven different functionals (red dots) together with their experimental values
(red solid lines) for comparison. Among various functionals, the rev-vdW-DF2
method shows better agreement with experimental values than the other
functionals. (See also Figure S3 and Table S1 in SI) Therefore, here, we
primarily use the rev-vdW-DF2 vdW correction in investigating piezoelectric
coefficients of 1L and multilayer MXs. For better understanding, we also
repeated all calculations using the Grimme-D3 vdW corrections, with all
results presented in SI. Note that, all physical results are qualitatively
consistent when using the same vdW functionals for both 1L and multilayer
cases, regardless of which vdW functionals were used. However, the
conventional approach (i.e., vdW functional is only considered in multilayers)
may lead to misinterpretation of the results, which will be discussed later.

Figure 3: (a) Side views of antiferroelectric (AFE) and ferroelectric (FE)
stacking configurations of bilayer MX and (b) top views of four possible in-
plane sliding configurations. In (a), blue arrows indicate the direction of
the spontaneous polarization of each layer. In (b), grey and blue dashed
rectangles visualize unitcells of each layer. (c) Contour maps of the
spontaneous polarization $P$ of the FE 2L SnSe in (c) AA and (d) AC stacking
configurations calculated by PBE with the rev-vdW-DF2 vdW correction.
Equilibrium lattice constants of 1L, 2L, 3L, and bulk AA and AC SnSe are
overlaid on (c) and (d), respectively.

Figure 3(a) and (b) show antiferroelectric (AFE) and FE stacking orders and
four possible sliding configurations. Thus, there are a total of eight
possible stacking configurations. Although it is known that the ground state
stacking configuration of MXs is the AFE AB stacking, xu2022van ;
higashitarumizu2020purely  it has also been predicted that the FE order can be
stabilized in both the AA and AC stacking configurations. xu2022van  In
addition, it has been experimentally confirmed that the AA-stacked FE SnS can
exist on a mica substrate below a certain critical thickness.
higashitarumizu2020purely  Therefore, here, we first focus on 2L AA and AC
stacked SnSe to understand the FE order in multilayer MXs.

Figure 3(c) and (d) show the $P$ contour maps of FE 2L SnSe in AA and AC
stacking configurations, respectively. Both contour maps are nothing but those
of their 1L counterparts shown in Fig. 1(c) with twice larger $P$ values,
implying that the ${P}$ of 2L SnSe is mainly due to the intralayer
contribution, regardless of the stacking configuration. Therefore, similar to
the 1L case, the equilibrium LCs primarily determine their ${P}$. The obtained
LCs of multilayer SnSe were overlaid in Fig. 3(c) and (d). The LCs of 2L AA
SnSe were evaluated to be $a=4.302$ Å and $b=4.263$ Å, whereas those of 2L AC
phase are $a=4.340$ Å and $b=4.251$ Å. Especially we focus on change in $a$
along the armchair direction when stacked from 1L to 2L. As shown in Fig. 3(c)
and (d), where the LCs of SnSe with different number of layers are overlaid,
AA stacking significantly reduces $a$, while AC stacking does not change $a$
considerably, compared to the 1L case. Thus, 2L AA SnSe has a smaller ${P}$
and may have a larger $e_{ijk}$ than 2L AC SnSe. It is worth noting that the
reduction of LCs of 2L AA SnSe can be understood as a result of dipole-dipole
interactions. Since parallel dipole-dipole interaction is inherently
unfavorable, relaxed structures can form in ways that favor either a decrease
in ${P}$ or an increase in the ${P}$-${P}$ distance. Notably, this result
remains robust across different vdW functionals and is thus applicable to
other MXs, emphasizing the generalizability of the observed phenomenon. (See
Table S2 in SI)

Figure 4: 2D piezoelectric strain coefficients (a) $e_{111}$ and (b)
$e_{122}$, and piezoelectric stress coefficients (c) $d_{11}$ and (d)
$d_{12}$) of 1L, 2L, and 3L SnSe with AA, AB, and AC stacking configurations,
calculated with the rev-vdW-DF2 vdW correction (black). For comparison, the
results of the Grimme-D3 corrections (grey) and PBE (blue) are also shown for
the AA SnSe and 1L SnSe, respectively.

Figure 4 shows the variation of four available piezoelectric coefficients
($e_{111}$, $e_{122}$, $d_{11}$, and $d_{12}$) of SnSe with the number of
layers up to 3L in AA, AB, and AC stacking configurations. (See Table S3 in SI
for more details) For the AB stacking configuration, due to its AFE order, the
piezoelectric coefficients of the even-numbered layers are calculated to be
zero, and those of the odd-numbered layers are also smaller than those of the
monolayer. On the other hand, the 2D piezoelectric strain coefficients
($e_{111}$ and $e_{122}$, in the unit of nC/m) increase with the number of
layers in both AA and AC stacking configurations. Moreover, AA stacking
exhibits higher piezoelectric coefficients than AC stacking, which is
consistent with our prediction based on the LCs. For example, the $e_{111}$ of
3L AA SnSe was calculated to be 13.5 nC/m which is larger than that of AC SnSe
(8.53 nC/m), and also six times larger than that of 3L AB SnSe (2.39 nC/m).
Note that the effective bulk piezoelectric coefficient (C/m2) of 3L AA SnSe is
also 15% larger than that of 1L SnSe thanks to the smaller LC.

The other important coefficients, the piezoelectric stress coefficients
($d_{11}$ and $d_{12}$) of multilayer AA SnSe are comparable or slightly
larger than those of 1L SnSe, depending on the choice of vdW correction. For
all other AA MXs, we calculated $e_{111}$, $e_{122}$, $d_{11}$, and $d_{12}$
up to 2L using both rev-vdW-DF2 and Grimme-D3 corrections, as summarized in
Tables S4-S6 in SI. We consistently found that as the LC decreases, $e_{111}$
and $e_{122}$ continue to enhance, while $d_{11}$ and $d_{12}$ remain
comparable to those of 1L MXs. It is worth noting that although there is no
meaningful improvement in $d_{11}$ and $d_{12}$ in multilayer AA MXs compared
to 1L, these values are much larger than those of AB MXs. For example, the
$d_{11}$ of 3L AA SnSe (155.96 pm/V) is almost three times larger than that of
3L AB SnSe (55.2 pm/V), suggesting that AA stacking MX is a promising
candidate in piezoelectric applications.

As shown in Fig. 3(c), the LC $a$ of AA SnSe continues to shrink with
stacking, but the bulk AA SnSe loses $P$ and becomes paraelectric, implying
that there exists a critical thickness below which non-zero $P$ is
accommodated. In fact, a previous experiment has successfully grown AA SnS on
mica substrates, exhibiting room temperature ferroelectricity up to 15 layers,
higashitarumizu2020purely  indicating that the critical thickness is not too
thin and could be further optimized by adjusting the experimental conditions.
This observation supports our prediction of optimizing the piezoelectric
coefficients of multilayer MXs by utilizing AA stacking configuration.

Figure 5: Summarized piezoelectric coefficients of (a) group IV
monochalcogenides (MXs) and (b) various 2D materials. In (a), purple, dark-
blue, yellow, green, skyblue, and orange circles represent AA (stacking) SnSe,
AB SnSe, AC SnSe, AA SnS, AA GeSe, and AA GeS, respectively. In (b), yellow,
orange, and purple circles denote materials of hexagonal Group III-V (h
III-V), hexagonal Group II oxides (h II Oxides), and H phase transition metal
dichalcogenides (2H-TMD), respectively, obtained from Ref. [8]. Skyblue and
red circles surrounded by a blue dashed rectangle denote our results of 1L MXs
and multilayer MXs.

Finally, we summarize the piezoelectric coefficients of multilayer MXs, as
shown in Fig. 5(a). SnSe shows the highest $e_{111}$ and $d_{11}$, followed by
SnS, GeSe, and GeS. In all MXs, the AA stacking boasts higher $e_{111}$ than
any other stacking, while $d_{11}$ is more or less than that of the
corresponding 1L case regardless of the stacking configuration. In addition to
$e_{111}$ and $d_{11}$, other physical properties such as stability and
switching barrier between the FE and AFE phases are also important in
piezoelectric applications. Therefore, all MXs can be practically useful, and
utilizing AA stacking configuration is an efficient way to optimize their
piezoelectric response.

For comparison, we also summarized the piezoelectric coefficients of several
other 2D materials, as shown in Fig. 5(b). Consistent with the previous
reports, we observed that 1L MXs have exceptionally large piezoelectric
coefficients compared to the other 2D materials, even larger coefficients in
AA stacked MXs. It is worthy of note that the giant piezoelectric coefficients
observed in AA MXs not only stand out among 2D materials, but are also
comparable to bulk piezoelectric materials. Recent high-throughput DFT
calculations revealed that among 941 bulk piezoelectric materials, only 5% of
materials satisfy $|e_{ij}|>3~{}$C/m2, including experimentally-confirmed
giant piezoelectric materials such as BaTiO3 (3.49 C/m2), SrHfO3 (8.73 C/m2),
RbTaO3 (8.93 C/m2), and BaNiO3 (27.46 C/m2). de2015database  Using the
effective thickness of 1.74 nm (1.74 nm=1.5$c$, $c$ is a out-of-plane lattice
constant of bulk SnSe), the effective bulk $e_{111}$ value of 3L AA SnSe was
calculated to be 7.76 C/m2, which is even larger than that of BaTiO3. Within
the Grimme-D2 vdW correction, which provides the smallest LCs among all vdW
functionals, the $e_{111}$ of 2L AA SnSe was estimated to be 12.66 C/m2 (14.68
nC/m). It is even larger than not only that of 3L AA SnSe estimated by rev-
vdW-DF2, but also that of SrHfO3 and RbTaO3. We now conclude that our
theoretical calculations have undoubtedly revealed that there is a strong
correlation between the reduction of LCs and giant piezoelectric coefficients
in AA-stacked MXs, regardless of the vdW functional used, making it a
promising candidate for low-dimensional piezoelectric applications.

##  IV Conclusion

In summary, we have used first-principles calculations to study the
piezoelectric coefficients of multilayer MXs with AA and AC stacking
configurations, which are known to stabilize ferroelectric order. Our re-
examination of 1L MXs revealed the pivotal role of the van der Waals
interaction in ensuring accurate and reliable predictions of piezoelectric
coefficients. Through systematic DFT calculations, we found that the AA-
stacked configurations exhibit remarkably larger piezoelectric coefficients
compared to all reported layered materials, including their monolayer
counterparts. The origin of this enhancement lies in the compressive strain
along the armchair direction, which is spontaneously introduced in the
ferroelectric AA stacking. These findings not only provide a comprehensive
understanding of the piezoelectric behavior of MXs, but also suggest a novel
strategy for optimizing piezoelectricity in low-dimensional materials through
stacking configuration.

###### Acknowledgements.

S.L., W.J., and T.L. were supported by the National Science Foundation (NSF)
through the DMREF program under Award No. DMR-1921629. S.L. is also supported
by Basic Science Research Program through the National Research Foundation of
Korea (NRF) funded by the Ministry of Education (NRF-2021R1A6A3A14038837).
H.-R.K. and Y.-K.K. acknowledge financial support from the Korean government
through the NRF (NRF-2022R1A2C1005505, NRF2022M3F3A2A01073562). The
computational work was partially done using the resources of the KISTI
Supercomputing Center (KSC-2023-CRE-0053).

## References

  * (1) Wang, Z. L. & Song, J.  Piezoelectric nanogenerators based on zinc oxide nanowire arrays.  _Science_ 312, 242–246 (2006).
  * (2) Wang, Z. L.  Nanopiezotronics.  _Adv. Mater._ 19, 889–892 (2007).
  * (3) Bao, R., Tao, J., Pan, C. & Wang, Z. L.  Piezophototronic effect in nanosensors.  _Small Sci._ 1, 2000060 (2021).
  * (4) Hinchet, R., Khan, U., Falconi, C. & Kim, S.-W.  Piezoelectric properties in two-dimensional materials: Simulations and experiments.  _Mater. Today_ 21, 611–630 (2018).
  * (5) Zhang, Q., Zuo, S., Chen, P. & Pan, C.  Piezotronics in two-dimensional materials.  _InfoMat_ 3, 987–1007 (2021).
  * (6) Jin, C.-C., Liu, D.-M. & Zhang, L.-X.  An emerging family of piezocatalysts: 2D piezoelectric materials.  _Small_ 19, 2303586 (2023).
  * (7) Ares, P. _et al._ Piezoelectricity in monolayer hexagonal boron nitride.  _Adv. Mater._ 32, 1905504 (2020).
  * (8) Blonsky, M. N., Zhuang, H. L., Singh, A. K. & Hennig, R. G.  Ab initio prediction of piezoelectricity in two-dimensional materials.  _ACS nano_ 9, 9885–9891 (2015).
  * (9) Zhu, H. _et al._ Observation of piezoelectricity in free-standing monolayer MoS2.  _Nat. Nanotechnol._ 10, 151–155 (2015).
  * (10) Alyoruk, M. M., Aierken, Y., Çakır, D., Peeters, F. M. & Sevik, C.  Promising piezoelectric performance of single layer transition-metal dichalcogenides and dioxides.  _J. Phys. Chem. C_ 119, 23231–23237 (2015).
  * (11) Zhou, Y. _et al._ Out-of-plane piezoelectricity and ferroelectricity in layered $\alpha$-In2Se3 nanoflakes.  _Nano Lett._ 17, 5508–5513 (2017).
  * (12) Xue, F. _et al._ Multidirection piezoelectricity in mono-and multilayered hexagonal $\alpha$-In2Se3.  _ACS nano_ 12, 4976–4983 (2018).
  * (13) Liu, F. _et al._ Room-temperature ferroelectricity in cuinp2s6 ultrathin flakes.  _Nat. Commun._ 7, 12357 (2016).
  * (14) Brehm, J. A. _et al._ Tunable quadruple-well ferroelectric van der waals crystals.  _Nat. Mater._ 19, 43–48 (2020).
  * (15) Fei, R., Li, W., Li, J. & Yang, L.  Giant piezoelectricity of monolayer group IV monochalcogenides: SnSe, SnS, GeSe, and GeS.  _Appl. Phys. Lett._ 107, 173104 (2015).
  * (16) Gomes, L. C., Carvalho, A. & Neto, A. C.  Enhanced piezoelectricity and modified dielectric screening of two-dimensional group-IV monochalcogenides.  _Phys. Rev. B_ 92, 214103 (2015).
  * (17) Wang, H. & Qian, X.  Two-dimensional multiferroics in monolayer group IV monochalcogenides.  _2D Mater._ 4, 015042 (2017).
  * (18) Higashitarumizu, N. _et al._ Purely in-plane ferroelectricity in monolayer SnS at room temperature.  _Nat. Commun._ 11, 2428 (2020).
  * (19) Xu, B., Deng, J., Ding, X., Sun, J. & Liu, J. Z.  Van der Waals force-induced intralayer ferroelectric-to-antiferroelectric transition via interlayer sliding in bilayer group-IV monochalcogenides.  _Npj Comput. Mater._ 8, 47 (2022).
  * (20) Li, Y. _et al._ Probing symmetry properties of few-layer MoS2 and h-BN by optical second-harmonic generation.  _Nano Lett._ 13, 3329–3333 (2013).
  * (21) Yoo, H. _et al._ Atomic and electronic reconstruction at the van der waals interface in twisted bilayer graphene.  _Nat. Mater._ 18, 448–453 (2019).
  * (22) Woods, C. _et al._ Charge-polarized interfacial superlattices in marginally twisted hexagonal boron nitride.  _Nat. Commun._ 12, 347 (2021).
  * (23) Yasuda, K., Wang, X., Watanabe, K., Taniguchi, T. & Jarillo-Herrero, P.  Stacking-engineered ferroelectricity in bilayer boron nitride.  _Science_ 372, 1458–1462 (2021).
  * (24) Meng, P. _et al._ Sliding induced multiple polarization states in two-dimensional ferroelectrics.  _Nat. Commun._ 13, 7696 (2022).
  * (25) Deb, S. _et al._ Cumulative polarization in conductive interfacial ferroelectrics.  _Nature_ 612, 465–469 (2022).
  * (26) Kohn, W. & Sham, L. J.  Self-Consistent Equations Including Exchange and Correlation Effects.  _Phys. Rev._ 140, A1133–A1138 (1965).
  * (27) Kresse, G. & Furthmüller, J.  Efficient Iterative Schemes for $ab$ $initio$ Total-Energy Calculations Using a Plane-Wave Basis Set.  _Phys. Rev. B_ 54, 11169–11186 (1996).
  * (28) Blöchl, P. E.  Projector Augmented-Wave Method.  _Phys. Rev. B_ 50, 17953–17979 (1994).
  * (29) Kresse, G. & Joubert, D.  From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method.  _Phys. Rev. B_ 59, 1758–1775 (1999).
  * (30) Perdew, J. P., Burke, K. & Ernzerhof, M.  Generalized Gradient Approximation Made Simple.  _Phys. Rev. Lett._ 77, 3865–3868 (1996).
  * (31) Grimme, S.  Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction.  _J. Comput. Chem._ 27, 1787–1799 (2006).
  * (32) Grimme, S., Antony, J., Ehrlich, S. & Krieg, H.  A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu.  _J. chem. phys._ 132 (2010).
  * (33) Dion, M., Rydberg, H., Schröder, E., Langreth, D. C. & Lundqvist, B. I.  Van der Waals density functional for general geometries.  _Phys. Rev. Lett._ 92, 246401 (2004).
  * (34) Lee, K., Murray, É. D., Kong, L., Lundqvist, B. I. & Langreth, D. C.  Higher-accuracy van der Waals density functional.  _Phys. Rev. B_ 82, 081101 (2010).
  * (35) Klimeš, J., Bowler, D. R. & Michaelides, A.  Chemical accuracy for the van der Waals density functional.  _J. Condens. Matter Phys._ 22, 022201 (2009).
  * (36) Hamada, I.  van der waals density functional made accurate.  _Phys. Rev. B_ 89, 121103 (2014).
  * (37) Spaldin, N. A.  A beginner’s guide to the modern theory of polarization.  _J. Solid State Chem._ 195, 2–10 (2012).
  * (38) Guo, S.-D., Guo, X.-S., Zhang, Y.-Y. & Luo, K.  Small strain induced large piezoelectric coefficient in $\alpha$-AsP monolayer.  _J. Alloys Compd._ 822, 153577 (2020).
  * (39) Fang, W. _et al._ Layer dependence of geometric, electronic and piezoelectric properties of snse.  _arXiv preprint arXiv:1603.01791_ (2016).
  * (40) Chattopadhyay, T., Pannetier, J. & Von Schnering, H. G.  Neutron diffraction study of the structural phase transition in SnS and SnSe.  _J. Phys. Chem. Solids_ 47, 879–885 (1986).
  * (41) Grandke, T. & Ley, L.  Angular-resolved uv photoemission and the band structure of GeS.  _Phys. Rev. B_ 16, 832–842 (1977).
  * (42) Murgatroyd, P. A. E. _et al._ GeSe: Optical Spectroscopy and Theoretical Study of a van der Waals Solar Absorber.  _Chem. Mater._ 32, 3245–3253 (2020).
  * (43) De Jong, M., Chen, W., Geerlings, H., Asta, M. & Persson, K. A.  A database to enable discovery and design of piezoelectric materials.  _Sci. Data_ 2, 150053 (2015).


$Y$, whose filtration comes from the cellular filtration of $Y$. To set up
notation, we have the $E_{1}$-page of this spectral sequence

$E_{1}=\oplus_{n=0}^{3}\pi_{*}tmf\Longrightarrow\pi_{*+n}tmf\wedge Y.$

For a homotopy class $\beta$ in $\pi_{*}tmf\wedge Y$, we denote by $\alpha[n]$
the element that detects it in the $E_{1}$-page of the $tmf$-based
Atiyah–Hirzebruch spectral sequence, where $n$ is the Atiyah–Hirzebruch
filtration of $\beta$, and $\alpha$ is a class in $\pi_{*}tmf$. The stem of
$\beta$ is then the stem of $\alpha$ plus $n$.

###### Proof of 6.38..

In our Atiyah–Hirzebruch notation, we can rewrite the two $\nu$-extensions of
6.38 as

  1. (1)

$\nu\cdot\bar{\kappa}^{2}\kappa[3]=\Delta\eta\kappa\bar{\kappa}[1]$,

  2. (2)

$\nu\cdot\bar{\kappa}^{3}[2]=\Delta^{2}\nu\kappa[0]$.

We first prove (2), namely, that
$\nu\cdot\bar{\kappa}^{3}[2]=\Delta^{2}\nu\kappa[0]$. In $\pi_{*}tmf\wedge
C_{\eta}$, we have

$\nu\cdot\bar{\kappa}^{3}[2]=\langle\nu,\bar{\kappa}^{3},\eta\rangle[0]$

by Lemma 5.3 of [WX18]. By Moss’s Theorem and the differential
$d_{11}(\Delta^{2}\kappa)=\eta\bar{\kappa}^{3}$ in the elliptic spectral
sequence of $tmf$, we have

$\langle\nu,\bar{\kappa}^{3},\eta\rangle=\Delta^{2}\nu\kappa.$

Mapping this relation along the inclusion $C_{\eta}\rightarrow Y$ gives us
(2).

For (1), note that in $\pi_{*}tmf\wedge\Sigma C_{\eta}$, we have

$\nu\cdot\bar{\kappa}^{2}\kappa[3]=\langle\nu,\bar{\kappa}^{2}\kappa,\eta\rangle[1]$

by Lemma 5.3 of [WX18]. Since $\bar{\kappa}^{2}\kappa$ is $\nu$-divisible in
$\pi_{*}tmf$, we may shuffle

$\langle\nu,\bar{\kappa}^{2}\kappa,\eta\rangle=\langle\bar{\kappa}^{2}\kappa,\nu,\eta\rangle.$

By Moss’s theorem and the differential
$d_{5}(\Delta\kappa\bar{\kappa})=\nu\bar{\kappa}^{2}\kappa$ in $tmf$, we have

$\langle\bar{\kappa}^{2}\kappa,\nu,\eta\rangle=\Delta\eta\kappa\bar{\kappa}.$

Pulling back this relation along the quotient map $Y\rightarrow\Sigma
C_{\eta}$ gives (1).

Extensions (3) and (4) follow by duality. The fact that there are no other
exotic $\nu$-extensions is discussed below. ∎

Most possibilities for other exotic $\nu$-extensions are ruled out in a
straightforward way by analyzing $i_{*}$ and $p_{*}$, duality, the fact that
$\bar{\kappa}\nu=0$. However, the following two extensions require us to
analyze the classical Adams Spectral Sequence. The following proof depends on
checking algebraic extensions in

$\mathop{\mathrm{Ext}}\nolimits_{\mathcal{A}}((H\mathbb{F}_{2})^{*}(tmf\wedge
Y),(H\mathbb{F}_{2})^{*})$

using Bruner’s $\mathop{\mathrm{Ext}}\nolimits$-program [Bru]. See Figure 18
for classical Adams $E_{2}$-charts for $tmf\wedge V(0)$ and $tmf\wedge Y$, and
see [DFHH14, Chapter 13] for $tmf$.

######  Lemma 6.39.

In $\pi_{*}tmf\wedge Y$,

  1. (1)

$\nu a[31,1]=0$

  2. (2)

$\nu a[36,2]=0$

Dually, we have

  3. (3)

$\nu a[134,2]=0$

  4. (4)

$\nu a[139,3]=0$

###### Proof.

To show this, we need to prove that

  1. (1)

$\nu a[31,1]\neq a[34,6]$,

  2. (2)

$\nu a[36,2]\neq a[39,7]$.

In our Atiyah–Hirzebruch notation, we can rewrite these extensions as

  1. (1)

$\nu\cdot\kappa^{2}[3]\neq\kappa\bar{\kappa}[0]$,

  2. (2)

$\nu\cdot\Delta\nu^{3}[3]\neq\Delta\eta\kappa[0]$.

We give a proof for (1) that $\nu\cdot\kappa^{2}[3]\neq\kappa\bar{\kappa}[0]$
using the classical Adams Spectral Sequence. We consider the Adams Spectral
Sequence for $tmf\wedge Y$ and its subquotients. We will show that the Adams
filtration of $\kappa^{2}[3]$ is 7 and the Adams filtration of
$\kappa\bar{\kappa}[0]$ is 8. The fact that there is no such $\nu$-extension
follows from the algebraic fact that on the Adams $E_{2}$-page, the
$h_{2}$-multiple of the first element is not the second element, which is
checked by a computer program.

For the class $\kappa\bar{\kappa}[0]$, it is clear that the Adams filtration
of $\kappa\bar{\kappa}$ in $\pi_{34}tmf$ is 8, (it is detected by the element
$d_{0}g$,) and it maps nontrivially on the Adams $E_{2}$-pages along the map
$tmf\rightarrow tmf\wedge Y$. The image under this map, which we denoted by
$d_{0}g[0]$, is a permanent cycle. It cannot be killed due to filtration
reasons. Therefore, the class $\kappa\bar{\kappa}[0]$ is detected by
$d_{0}g[0]$ and, in particular, it has Adams filtration 8.

For the class $\kappa^{2}[3]$, we first consider the class $\kappa^{2}[1]$ in
$\pi_{29}tmf\wedge V(0)$. Since $\pi_{29}tmf=0,\pi_{30}tmf=0$, we have
$\pi_{30}tmf\wedge V(0)=0$. This forces three nonzero Adams differentials
eliminating the three elements in the Adams $E_{2}$-page for $tmf\wedge V(0)$.
In particular, we learn that $\kappa^{2}[1]$ in $\pi_{29}tmf\wedge V(0)$ is
detected by the only remaining element $j[0]$ in Adams filtration 7, and that
there is a nonzero $d_{3}$-differential from $(t-s,s)$-bidegrees $(31,6)$ to
$(30,9)$.

Considering the quotient map $tmf\wedge Y\rightarrow tmf\wedge\Sigma^{2}V(0)$,
we learn that $\kappa^{2}[3]$ is detected in Adams filtration at most 7.
Considering the induced map on the Adams $E_{2}$-pages, we also learn that it
is an isomorphism on the $(t-s,s)$-bidegrees $(31,6)$ and $(30,9)$. So, in
particular, the element in $(t-s,s)$-bidegree $(31,6)$ does not survive.
Therefore, $\kappa^{2}[3]$ is detected in Adams filtration exactly 7.

For (2), that $\nu\cdot\Delta\nu^{3}[3]\neq\Delta\eta\kappa[0]$, we use the
Adams spectral sequence again in a very similar way. We will show that the
Adams filtration of $\Delta\nu^{3}[3]$ is 8 and the Adams filtration of
$\Delta\eta\kappa[0]$ is 9. The fact that there is no such extensions then
follows as in (1).

For the class $\Delta\eta\kappa[0]$, it is clear that the Adams filtration of
$\Delta\eta\kappa$ in $\pi_{39}tmf$ is 9, (it is detected by the element $u$,)
and it maps nontrivially on the Adams $E_{2}$-pages along the map
$tmf\rightarrow tmf\wedge Y$. The image under this map, which we denoted by
$d_{0}g[0]$, is a permanent cycle. It cannot be killed due to filtration
reasons. Therefore, the class $\Delta\eta\kappa[0]$ is detected by $u[0]$, and
in particular it has Adams filtration 9.

For the class $\Delta\nu^{3}[3]$, we first consider the class $\Delta\nu^{3}$
in $\pi_{33}tmf$. The class $\Delta\nu^{3}$ in $\pi_{33}tmf$ is detected in
the Adams filtration 8. Considering the quotient map $tmf\wedge
Y\rightarrow\Sigma^{3}tmf$, we learn that $\Delta\nu^{3}[3]$ is detected in
Adams filtration at most 8. To show that it is detected in Adams filtration 8,
we will show that the only other element in lower filtration, the class in
$(t-s,s)$-bidegree $(36,7)$, supports a nonzero $d_{2}$-differential.

The maps in the zigzag

$\textstyle{tmf\wedge S^{1}}$$\textstyle{tmf\wedge
V(0)\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{tmf\wedge
Y}$

are isomorphisms in $(t-s,s)$-bidegrees $(36,7)$ and $(35,9)$ on Adams
$E_{2}$-pages. So the claimed nonzero $d_{2}$-differential follows from the
one in the Adams spectral sequence of $tmf$, from $(t-s,s)$-bidegrees $(35,7)$
and $(34,9)$. ∎

We now turn to the study of the $v_{1}$-extensions. First, recall the
discussion on $v_{1}$-self maps and $A_{1}$ from Section 2.6. The homotopy
groups of $tmf\wedge A_{1}$ are studied by the third author in [Pha18].
Furthermore, the knowledge of the homotopy groups of $tmf\wedge A_{1}$ are
sufficient to allow us to deduce much of the action of $v_{1}$ on the homotopy
groups of $tmf\wedge Y$, via the long exact sequence on homotopy of the
cofiber sequence

$tmf\wedge\Sigma^{2}Y\xrightarrow{v_{1}}tmf\wedge Y\rightarrow tmf\wedge
A_{1}.$

Since the outcome depends on the choice of the $v_{1}$-self-map, we call a
$v_{1}$-self-map of type $I$ if its cofiber is $A_{1}[01]$ or $A_{1}[10]$ and
of type $II$, otherwise. Again, see Section 2.6 for the definition of
$A_{1}[ij]$.

######  Lemma 6.40.

(a) For all $v_{1}$-self maps of $Y$, there exotic $v_{1}$-extensions and
those induced by $\bar{\kappa}$-linearity:

  1. (1)

$v_{1}a[9,1]=a[11,3]$

  2. (2)

$v_{1}a[15,1]=a[17,3]$

  3. (3)

$v_{1}a[30,2]=\bar{\kappa}a[12,2]$

  4. (4)

$v_{1}a[33,1]=a[35,3]$

  5. (5)

$v_{1}a[38,2]=\bar{\kappa}a[20,2]$

  6. (6)

$v_{1}\Delta^{2}a[9,1]=\Delta^{2}a[11,3]$

  7. (7)

$v_{1}a[99,1]=\bar{\kappa}^{3}a[21,3]$

  8. (8)

$v_{1}a[104,2]=\bar{\kappa}^{4}a[26,0]$

  9. (9)

$v_{1}a[105,1]=a[107,3]$

  10. (10)

$v_{1}(v_{1}a[108,2])=\bar{\kappa}^{3}a[52,0]$

  11. (11)

$v_{1}a[114,2]=\bar{\kappa}^{4}a[36,2]$

  12. (12)

either $v_{1}a[116,2]=\bar{\kappa}^{2}a[78,0]$ or
$v_{1}a[116,2]=\bar{\kappa}a[98,0]$

  13. (13)

$v_{1}\bar{\kappa}a[105,1]=\bar{\kappa}^{3}a[67,3]$

  14. (14)

$v_{1}a[129,1]=a[131,3]$

  15. (15)

either $v_{1}a[131,3]=\bar{\kappa}^{2}a[93,3]$ or
$v_{1}a[131,3]=\bar{\kappa}a[113,3]$.

  16. (16)

$v_{1}a[134,2]=\bar{\kappa}a[116,2]$

  17. (17)

$v_{1}\bar{\kappa}a[115,3]=\bar{\kappa}a[117,13]$

  18. (18)

$v_{1}(v_{1}a[139,3])=\bar{\kappa}^{3}a[83,3]$

  19. (19)

$v_{1}\bar{\kappa}a[120,3]=\bar{\kappa}a[122,14]$

  20. (20)

$v_{1}(v_{1}\bar{\kappa}a[124,0])=\bar{\kappa}^{4}a[68,2]$

  21. (21)

$v_{1}a[147,1]=\bar{\kappa}a[129,1]$

  22. (22)

$v_{1}a[152,2]=\bar{\kappa}a[134,2]$

  23. (23)

$v_{1}a[156,10]=a[158,16]$

  24. (24)

$v_{1}a[162,2]=\bar{\kappa}^{2}a[124,0]$

(b) For $v_{1}$-self-maps of type $I$, there are also the following
$v_{1}$-extensions, and those induced from these by $\bar{\kappa}$-linearity:

  1. (1)

$v_{1}a[68,2]=\bar{\kappa}^{2}a[30,2]$

  2. (2)

$v_{1}a[83,3]=\bar{\kappa}^{4}a[15,1]$

###### Proof.

For all parts, except for $(\ref{intrus1})$, $(\ref{intrus2})$,
$(\ref{intrus3})$, we see, by inspecting the relevant parts of the homotopy
groups of appropriate $tmf\wedge A_{1}[ij]$, that the targets of the stated
$v_{1}$ extensions are sent to zero via the natural map

$\pi_{*}(tmf\wedge Y)\rightarrow\pi_{*}(tmf\wedge A_{1}[ij]).$

Therefore, they are in the image of a $v_{1}$-multiplication and the stated
$v_{1}$-extensions are the only possibilities.

For part (9), consider

$\mathrm{sk}_{4}(tmf\wedge Y)=(tmf\wedge Y)/(tmf\wedge Y)_{5},$

where $(tmf\wedge Y)_{5}$ is the $5$th term in the $X(4)$-Adams tower of
$tmf\wedge Y$. It is a module over $\mathrm{sk}_{4}(tmf)$. Since
$\Delta^{4}\in\pi_{96}(\mathrm{sk}_{4}(tmf))$, this element acts on
$\pi_{*}\mathrm{sk}_{4}(tmf\wedge Y)$. We see that the induced map
$\pi_{*}(tmf\wedge Y)\rightarrow\pi_{*}\mathrm{sk}_{4}(tmf\wedge Y)$ sends
$a[9,1]$ and $a[11,3]$ to non-trivial elements, which we denote by the same
names. Furthermore, it sends $a[105,1]$ and $a[107,3]$ to elements detected by
the products $\Delta^{4}a[9,1]$ and $\Delta^{4}a[11,3]$. Since
$v_{1}a[9,1]=a[11,3]$ by part (1),

$v_{1}\Delta^{4}a[9,1]=\Delta^{4}v_{1}a[9,1]=\Delta^{4}a[11,1]$

in $\pi_{*}\mathrm{sk}_{4}(tmf\wedge Y)$. It follows that $v_{1}a[105,1]$ must
be detected by $a[107,3]$ in the $E_{\infty}$-term of the elliptic spectral
sequence of $tmf\wedge Y$. ∎

######  Remark 6.41.

We are left with two undecided $v_{1}$-extensions, namely (12) and (15) in
6.40. We expect that some of these unsettled $v_{1}$-extensions can be
resolved using comparison with the classical Adams Spectral Sequence for
$tmf\wedge V(0)$, $tmf\wedge Y$ and $tmf\wedge V(0)/v_{1}^{4}$. These will
soon appear in upcoming work of [BR].

Figure 18. Classical Adams Spectral Sequence $E_{2}$-pages for $tmf\wedge
V(0)$ (top) and $tmf\wedge Y$ (bottom) computed with Bruner’s
$\mathop{\mathrm{Ext}}\nolimits$-program [Bru].

Figure 19. $d_{11}$ to $d_{23}$ differentials in stems 46 to 120

Figure 20. $d_{11}$ to $d_{23}$ differentials in stems $120$ to $194$

Figure 21. Exotic extensions in the elliptic spectral sequence of $tmf\wedge
Y$. This records
$tmf_{*}Y\cong\widetilde{tmf}_{*+3}(\mathbb{R}P^{2}\wedge\mathbb{C}P^{2})$.
The zigzags denote exotic $v_{1}$-extensions that occur only for certain
choices of $v_{1}$ self-maps.

Figure 22. Exotic extensions in the elliptic spectral sequence of $tmf\wedge
Y$. This records
$tmf_{*}Y\cong\widetilde{tmf}_{*+3}(\mathbb{R}P^{2}\wedge\mathbb{C}P^{2})$.
The zigzags denote exotic $v_{1}$-extensions that occur only for certain
choices of $v_{1}$ self-maps.

## References

  * [Bau08] T. Bauer.  Computation of the homotopy of the spectrum tmf.  In Groups, homotopy and configuration spaces, volume 13 of Geom. Topol. Monogr., pages 11–40. Geom. Topol. Publ., Coventry, 2008.
  * [BBB+19] A. Beaudry, M. Behrens, P. Bhattacharya, D. Culver, and Z. Xu.  On the tmf-resolution of Z.  arXiv e-prints, page arXiv:1909.13379, September 2019.
  * [BE20a] P. Bhattacharya and P. Egger.  A class of 2-local finite spectra which admit a $v_{2}^{1}$-self-map.  Adv. Math., 360:106895, 40, 2020.
  * [BE20b] P. Bhattacharya and P. Egger.  Towards the $K(2)$-local homotopy groups of $Z$.  Algebr. Geom. Topol., 20(3):1235–1277, 2020.
  * [Bea15] A. Beaudry.  The algebraic duality resolution at p = 2.  Algebr. Geom. Topol., 15(6):3653–3705, 2015.
  * [Bea17] A. Beaudry.  Towards the homotopy of the $K(2)$-local Moore spectrum at $p=2$.  Adv. Math., 306:722–788, 2017.
  * [Beh12] M. Behrens.  The Goodwillie tower and the EHP sequence.  Mem. Amer. Math. Soc., 218(1026):xii+90, 2012.
  * [BG18] I. Bobkova and P. G. Goerss.  Topological resolutions in $K(2)$-local homotopy theory at the prime 2.  Journal of Topology, 11(4):917–956, 2018.
  * [BGH17] A. Beaudry, P. G. Goerss, and H.-W. Henn.  Chromatic splitting for the $K(2)$-local sphere at $p=2$.  arXiv e-prints, page arXiv:1712.08182, December 2017.
  * [BL01] A. Baker and A. Lazarev.  On the Adams spectral sequence for $R$-modules.  Algebr. Geom. Topol., 1:173–199, 2001.
  * [Bou79] A. K. Bousfield.  The localization of spectra with respect to homology.  Topology, 18(4):257–281, 1979.
  * [BR] R. R. Bruner and J. Rognes.  The Adams Spectral Sequence for Topological Modular Forms.  Mathematical Surveys and Monographs. American Mathematical Society, Providence, RI.  To appear.
  * [Bro10] K. S. Brown.  Lectures on the cohomology of groups.  In Cohomology of groups and algebraic $K$-theory, volume 12 of Adv. Lect. Math. (ALM), pages 131–166. Int. Press, Somerville, MA, 2010\.
  * [Bru] R. R. Bruner.  Cohomology charts software ext.1.9.3. available at http://www.rrb.wayne.edu/papers/. 2018.
  * [Bru78] R. Bruner.  Algebraic and geometric connecting homomorphisms in the Adams spectral sequence.  In Geometric applications of homotopy theory (Proc. Conf., Evanston, Ill., 1977), II, volume 658 of Lecture Notes in Math., pages 131–133. Springer, Berlin, 1978.
  * [DFHH14] C. L. Douglas, J. Francis, A. G. Henriques, and M. A. Hill, editors.  Topological modular forms, volume 201 of Mathematical Surveys and Monographs.  American Mathematical Society, Providence, RI, 2014.
  * [DM81] D. M. Davis and M. Mahowald.  $v_{1}$\- and $v_{2}$-periodicity in stable homotopy theory.  Amer. J. Math., 103(4):615–659, 1981.
  * [GH16] P. G. Goerss and H.-W. Henn.  The Brown-Comenetz dual of the $K(2)$-local sphere at the prime 3\.  Adv. Math., 288:648–678, 2016.
  * [GHM04] P. G. Goerss, H.-W. Henn, and M. Mahowald.  The homotopy of $L_{2}V(1)$ for the prime 3.  In Categorical decomposition techniques in algebraic topology (Isle of Skye, 2001), volume 215 of Progr. Math., pages 125–151. Birkhäuser, Basel, 2004.
  * [GHM14] P. G. Goerss, H.-W. Henn, and M. Mahowald.  The rational homotopy of the $K(2)$-local sphere and the chromatic splitting conjecture for the prime 3 and level 2.  Doc. Math., 19:1271–1290, 2014.
  * [GHMR05] P. G. Goerss, H.-W. Henn, M. Mahowald, and C. Rezk.  A resolution of the $K(2)$-local sphere at the prime 3.  Ann. of Math. (2), 162(2):777–822, 2005.
  * [GHMR15] P. G. Goerss, H.-W. Henn, M. Mahowald, and C. Rezk.  On Hopkins’ Picard groups for the prime 3 and chromatic level 2.  J. Topol., 8(1):267–294, 2015.
  * [Gre16] J. P. C. Greenlees.  Ausoni-Bökstedt duality for topological Hochschild homology.  J. Pure Appl. Algebra, 220(4):1382–1402, 2016.
  * [GS18] J. P. C. Greenlees and V. Stojanoska.  Anderson and Gorenstein duality.  In Geometric and topological aspects of the representation theory of finite groups, volume 242 of Springer Proc. Math. Stat., pages 105–130. Springer, Cham, 2018.
  * [Hen07] H.-W. Henn.  On finite resolutions of $K(n)$-local spheres.  In Elliptic cohomology, volume 342 of London Math. Soc. Lecture Note Ser., pages 122–169. Cambridge Univ. Press, Cambridge, 2007.
  * [HKM13] H.-W. Henn, N. Karamanov, and M. Mahowald.  The homotopy of the $K(2)$-local Moore spectrum at the prime 3 revisited.  Math. Z., 275(3-4):953–1004, 2013.
  * [IWX20] D. C. Isaksen, G. Wang, and Z. Xu.  More stable stems.  arXiv e-prints, page arXiv:2001.04511, January 2020.
  * [Mos70] R. M. F. Moss.  Secondary compositions and the Adams spectral sequence.  Math. Z., 115:283–310, 1970.
  * [MR99] M. Mahowald and C. Rezk.  Brown-Comenetz duality and the Adams spectral sequence.  Amer. J. Math., 121(6):1153–1177, 1999.
  * [Pha18] V.-C. Pham.  Homotopy groups of $E_{C}^{hG_{24}}\wedge A_{1}$.  arXiv e-prints, page arXiv:1811.04484, November 2018.
  * [Rav78] D. C. Ravenel.  A novice’s guide to the Adams-Novikov spectral sequence.  In Geometric applications of homotopy theory (Proc. Conf., Evanston, Ill., 1977), II, volume 658 of Lecture Notes in Math., pages 404–475. Springer, Berlin, 1978.
  * [Rav86] D. C. Ravenel.  Complex Cobordism and Stable Homotopy Groups of Spheres, volume 121 of Pure and Applied Mathematics.  Academic Press Inc., Orlando, FL, 1986.
  * [Rav87] D. C. Ravenel.  Localization and periodicity in homotopy theory.  In Homotopy theory (Durham, 1985), volume 117 of London Math. Soc. Lecture Note Ser., pages 175–194. Cambridge Univ. Press, Cambridge, 1987.
  * [Rez] C. Rezk.  Supplementary notes for MATH 512.  Online notes.
  * [Sil86] J. H. Silverman.  The Arithmetic of Elliptic Curves, volume 106 of Graduate Texts in Mathematics.  Springer-Verlag, New York, 1986.
  * [Sto12] V. Stojanoska.  Duality for topological modular forms.  Doc. Math., 17:271–311, 2012.
  * [Sto14] V. Stojanoska.  Calculating descent for 2-primary topological modular forms.  In An alpine expedition through algebraic topology, volume 617 of Contemp. Math., pages 241–258. Amer. Math. Soc., Providence, RI, 2014\.
  * [WX18] G. Wang and Z. Xu.  Some extensions in the Adams spectral sequence and the 51-stem.  Algebr. Geom. Topol., 18(7):3887–3906, 2018.


# Vertical-supercooling-controlled interfacial instability for a spreading
liquid film

Li Chen  Feng Wang  Yingrui Wang  Peng Huo  Yuqi Li  Xi Gu  Man Hu  Daosheng
Deng<EMAIL_ADDRESS>Department of Aeronautics and Astronautics, Fudan
University,Shanghai, 200433, China

###### Abstract

Thermal effect is essential to regulate the interfacial instabilities for
diverse technology applications. Here we report the fingering instability at
the propagation front for a spreading liquid film subjected to the
supercooling at the vertical direction. We find the onset timescale of
hydrodynamic instability is strongly correlated with that of the vertical
solidification process. This correlation is further validated in a non-uniform
geometry, demonstrating the capability of controlling fingering instability by
structure design. Moreover, based on the experimental observations, we propose
a physical mechanism by considering thermal Marangoni effect at the spreading
front, and the predicted wavelength from the linear stability analysis agrees
with experiments excellently. This work offers another valuable dimension by
gating the vertical temperature to exploit the interfacial stabilities and
steer liquid flow, consequently shedding light on the microfluidic cooling for
electronics, and the advanced functional fibers and fabrics.

††preprint:

Interfacial instabilities and the correlated intriguing patterns are not only
ubiquitous in nature, such as the snowflakes and patterns of bacterial
colonies langer1980instabilities ; ben2000cooperative , but also essential for
technology applications, such as the coating processing, microfluidics, and
Lab-on-a-Chip craster2009dynamics ; stone2004engineering ;
darhuber2005principles . Particularly, thermal effect plays an indispensable
role for the interfacial instabilities davis1987thermocapillary ;
schatz2001experiments . For example,during the unidirectional freezing,
Mullins-Sekerka instability occurs at the moving planar liquid-solid interface
mullins1964stability ; and under a unidirectional temperature gradient, the
thermocapillary instability appears boos1997thermocapillary .At the
microscopic scale, the heat transfer becomes even more pronounced for
interfacial stabilities, significantly influencing the fluid dynamics and the
subsequently produced patterns cazabat1990fingering ; kataoka1999patterning .
In spite of these extensive studies, heat transfer and temperature gradient
are mainly confined within the liquid film.

Recently, by surveying an impacted droplet on a cold substrate
Soninphysfluid1997 , diverse patterns have been clearly revealed, such as the
crack pattern changing from a 2D fragmentation to a hierarchical fracture
ghabache2016frozen , the self-peeling of impacting droplets RuiterNatphy2018 ,
and morphology evolution from conical tips to toroidal shapes hu2020frozen .
Fingering growth emerges in a binary droplet, which freezes from the outside
prior to the impacting on the cold surface kant2020pattern . At a microscopic
level, the nucleation process and the associated solidification or
crystalization during the droplet spreading and cooling have been revealed to
be responsible for the eventual morphology and the arrest of the contact line
Pallav2020fastfreezing ; Koldeweij2021Initial ; Grivet2022Contact ;
LambleyNatphys2023 . However, the freezing of impacted droplets and the
correlated spreading process are extremely complicated, because of the
inherent complex interplay of the rapid impacting hydrodynamics, the non-
uniform liquid thickness, the transient heat transfer, and the intricate
solidification.

In this Letter, different from the aforementioned frozen droplets, we explore
the spreading of liquid film with uniform thickness by injecting a viscous jet
at an elevated temperature into a Hele-Shaw cell while being subjected to the
supercooling along the walls. This setup allows us to thoroughly investigate
the effect of vertical heat transfer, as another degree of freedom to
precisely control interfacial instability, for the first time to be best of
our knowledge. We observe fingering instability at the spreading front, and
identify the onset of the instability is primarily determined by the
solidification process at the vertical direction. We propose a physical
mechanism of this instability, and the linear stability analysis agrees with
experiments excellently. Moreover, we demonstrate this vertical-supercooling-
induced fingering instability can be further engineered simply through
structure designing of a non-uniform geometry.

Figure 1: Spatiotemporal evolution of fingering instability for an injected
liquid at an elevated temperature into a Hele-Shaw cell subjected to the side-
wall supercooling. (a) Schematic diagram of the experimental setup; (b, c)
Snapshot for a stable interface at a higher $T_{w}\approx T_{m}$ ($T_{w}=67$
∘C) and fingering instability at a lower $T_{w}<T_{m}$ ($T_{w}=25\ ^{\circ}$C)
for a paraffin wax film in a glass cell ($Q=4$ ml/s, $h=$ 400 $\mu$m); (d)
Temporal evolution of the azimuthal undulations along the radius $r(\theta)$
for the interface in row (c).

_Fingering instability_ As shown in Fig. 1a, paraffin wax is melted into a
viscous state at an elevated temperature ($T_{l}=95\ ^{\circ}$C) above its
melting point ($T_{m}=69\ ^{\circ}$C). Subsequently, this viscous liquid is
injected into a Hele-Shaw cell consisting of two glass plates with a gap
thickness $h$ = 400 $\mu$m through a central inlet ($2r_{in}=2.5$ mm). The
flow rate of injection is varied by a syringe pump ($Q=0.25-4$ ml/s), and the
cooling process is controlled or gated by the wall temperature of Hele-Shaw
cell ($T_{w}$). At a higher $T_{w}=67\ ^{\circ}$C ($T_{w}\approx T_{m}$), when
the injected liquid spreads within the Hele-Shaw cell, the moving front is
stable, and the circular symmetry of the interface is preserved (Fig. 1b).

Figure 2: Scaling analysis for the onset of instability. (a) The temporal
evolution of $r$ dependent on flow rates $Q$; (b) Prior to the instability
onset, the spreading radius has a power law of $1/2$ with time,
$r/r_{c}\sim(t/t_{c})^{1/2}$; (c) At the onset of the instability, the
experimental data of $t_{c}$ agreeing with theory of Eq. (2) and (3)
excellently $t_{c}=t_{s}\sim h^{2}$; (d) Critical radius
$r_{c}\sim(Re)^{1/2}$.

However, at a lower $T_{w}=$ 25 ∘C ($T_{w}<T_{m}$) (Fig. 1c), the injected
viscous liquid accordingly undergoes the supercooling process, and an
interfacial front instability is identified with the pronounced fingering
patterns along the radial direction (Supplementary Movie 1). As shown by the
spatiotemporal snapshots, initially the liquid film spreads outward uniformly
($t=305,430$ ms), and the propagation front forms a nearly perfect circular
interface. Once its radius is beyond a critical value ($r>r_{c}$), some little
notches appear at the moving front and signify the onset of the fingering
instability ($t=555$ ms). As the fluid is continuously injected, the front
progressively propagates radially, and fingers further extend longer outward
($t=680,805$ ms). This fingering instability is quantitatively characterized
by temporal evolution of the azimuthal undulations along the interface
[$r(\theta)$] as shown by Fig. 1d. Note that the notches at the interface are
pinning solidification points during the subsequent growing process, and the
wavelength $\lambda$ between the nearest pinning solidification points is
associated with the circumferential length and the number of the observed
fingers ($N_{f}$), $\lambda_{exp}\sim 2\pi r_{c}/N_{f}$.

_Onset of the fingering instability_ A quantitative understanding for the
onset of the instability is obtained by checking the critical radius ($r_{c}$)
with various physical parameters (such as $Q,h$) (Supplementary Movie 2,
Supplementary Note 3). The temporal evolution of radius ($r$) for various flow
rates ($Q$) (Fig. 2a) indicates the turning point for the onset of instability
(a critical time $t_{c}$ or radius $r_{c}$), after which the fingering
instability further develops (the radius is denoted for the length of typical
fingers $r_{f}$). Prior to the instability, based on the mass conservation,
$Q=2\pi rhdr/dt$, the radius dependent on time is obtained,

$r\sim(Qt/\pi h)^{1/2}\sim t^{1/2}.$ (1)

As shown in Fig. 2b, the experimental data collapse into a master curve with a
$1/2$ scaling law, $r/r_{c}\sim(t/t_{c})^{1/2}$.

Because of the geometric confinement in the Hele-Shaw cell, the injected
viscous liquid at a high temperature $T_{l}$ is subjected to conductive
cooling due to the bottom and top plate at a low temperature ($T_{w}$) at the
vertical direction, likely resulting in the aforementioned instability. By
balancing the diffusive heat flux with the solidification rate
ghabache2016frozen ; hu2020frozen , the timescale of solidification process at
the vertical direction is obtained,

$t_{s}={\rho\mathcal{L}(h/2)^{2}}/{2\kappa\Delta T}\sim h^{2},$ (2)

where $\Delta T=T_{m}-T_{w}$, $\mathcal{L}$, $\rho$, $\kappa$ for the
supercooling, latent heat, density and thermal conductivity (Table I for the
materials parameters). Indeed, the experimental data of $t_{c}$ for the onset
of the instability agrees excellently with the $t_{s}$ for the solidification
timescale at the vertical direction in theory (Fig. 2c),

$t_{c}=t_{s}.$ (3)

Additionally, the scaling law holds for various gap thicknesses as well
($t_{c}=t_{s}\sim h^{2}$).

Further, from Eq. (1)-(3), the critical radius $r_{c}$ is attained,

$r_{c}\sim(Qh)^{1/2}(\rho\mathcal{L}/8\pi\kappa\Delta T)^{1/2}\sim Re^{1/2},$
(4)

where $Re=\rho v_{in}h/\mu$ ($v_{in}=Q/\pi r_{in}^{2}$, $\mu$ for viscosity)
for Reynolds number. This $r_{c}\sim Re^{1/2}$ scaling is consistent with the
experimental observation for paraffin wax (Fig. 2d). For other materials with
various physical and thermal properties (Table I in Supplementary and
Supplementary Note 4), the $1/2$ scaling law still holds (Fig. 1d) in eicosane
($\rm C_{20}H_{42}$) with the similar thermal properties of the paraffin wax.
But for the different thermal properties, such as liquid metal of Bi-Sn alloy
with a larger thermal conductivity and smaller latent heat, $r_{c}$ becomes
smaller to be comparable with $r_{in}$, hence the instability appears
immediately right after the injection of this liquid metal.

Figure 3: Geometry effect on the fingering instability. (a) Schematic diagram
of the non-uniform gap thickness with a small constant gradient ($\varphi\ll
1$); (b) Interface susceptible to instability in the narrower side, while
maintaining stability in the wider side; (c) $r$ at the interface as a
function of azimuth at different times; (d) Emergence of fingering
instability, once $\tilde{t}>2$ (as shaded in the green), consistent with the
experimental observation in (b). Orange dots in (b-d) marking the region of
interfacial instability.

_Geometry with a nonuniform gap_ Since the timescale of heat transfer along
the vertical direction determines the onset of instability [$t_{c}\sim h^{2}$
in Eq. (2)], we further directly demonstrate the effect of thickness in a
nonuniform geometry on the resultant instability. In this wedge-shaped Hele-
Shaw cell (Fig. 3a) al2012control , the constant depth gradient is extremely
small on the order of $O(10^{-3})$ expnonuniformdata .

When the molten paraffin wax is injected into this cell at $T_{w}=$ 25 ∘C and
$Q=$ 4 ml/s, the high-speed images show the pattern evolution (Fig. 3b and
Supplementary Movie 3), and the interfacial profile $r(t,\theta)$ depends on
time and the azimuthal angle (Fig. 3c). Initially the spreading profile
remains a nearly circular shape ($t=$ 160 ms). Afterwards, in the left-side
region with a narrower gap, the interface or the propagation front is more
susceptible to the instability because of rapid cooling, and fingering
instability becomes pronounced beyond a critical time ($t_{c}\sim 344$ ms).
Alternately, in the right-side region with a wider gap, the interface still
remains stable and smooth ($t=$ 560 ms). In this way, the asymmetry pattern is
formed, _i.e._ , the stability is persisted at the wider gap while the
instability is triggered at the narrower gap ($t=$ 720 ms).

This geometry effect on the pattern stability is further analyzed by a simple
scaling analysis. By combining equation Eq. (2) and the local height at the
interface which is described by
$h(r,\theta)=h_{in}+r(t,\theta)\rm{cos}\theta\rm{sin}\varphi$
[$h_{in}=(h_{n}+h_{w})/2$], the local solidification timescale at the
interface is obtained,

$t_{s}(r,\theta)=\rho\mathcal{L}h^{2}/8\kappa\Delta T\sim h(r,\theta)^{2}.$
(5)

We establish the stability diagram for the interface at propagation front by
defining dimensionless time $\tilde{t}(\theta)=t/t_{s}$ (Fig. 3d). The
interface is stable for $\tilde{t}<2$, while fingering instability sets up for
$\tilde{t}>2$. The experimental data agrees reasonably well with this scaling
analysis, which is built on the heat transfer at the vertical direction by
addressing the solidification timescale during the supercooling process.

_Temperature visualization_ In order to unravel the role of the supercooling
process for this fingering instability, the spatial-temporal evolution of
temperature for the viscous liquid is directly recorded by a thermal camera
(Fig. 4a, Supplementary Movie 4, Supplementary Note 2). Fig. 4b presents
temperature distribution along the radial direction of the marked finger (the
black arrow in the snapshot at 296 ms). Adjacent to the inlet ($r\approx
r_{in}$), the injected liquid is immediately cooled down by contacting with
the sidewall; but in the spreading region ($r_{in}<r<r_{f}$), temperature
slightly decreases with the radius. In the vicinity of the propagation front
($r\approx r_{f}$), temperature drops dramatically, generating a pronounced
temperature gradient (d$T/$d$r$) around tens of K$/$mm (Fig. 4c).

Arising from this temperature gradient, the resultant surface-tension gradient
might play an essential role to be responsible for the emergence of the front
instability. Consequently, this thermal Marangoni force is balanced by the
viscous force
$\mathrm{d}\gamma/\mathrm{d}r=(\mathrm{d}\gamma/\mathrm{d}T)(\mathrm{d}T/\mathrm{d}r)\sim\mu
v_{\mathrm{M}}/h$, resulting in a characteristic Marangoni velocity of
$v_{\mathrm{M}}\approx(\mathrm{d}\gamma/\mathrm{d}r)h/\mu$. From the physical
parameters of paraffin wax Koller2017Liquid ; Klein2019Liquid ;
Yang2021Surface , $|\mathrm{d}\gamma/\mathrm{d}T|\approx 0.01\rm{\ mN/(m\cdot
K)}$ with a linear approximation, and $\mu=$ 1.8 mPa $\cdot$s,
$v_{\mathrm{M}}\approx 20\rm{\ mm/s}$ is estimated theoretically.

From the mass conservation, the velocity of propagation front or the interface
is $v_{f}=Q/2\pi rh\sim r^{-1}$, and experimental data agrees with the $-1$
power law as expected (Fig. 4d). However, once $v_{f}$ reaches down to
$v_{\textup{M}}$, indeed, the instability is triggered, which is reproducible
at various $Q$ and $h$ in experiments.

Figure 4: Spatiotemporal evolution of temperature. (a) Thermal images of the
viscous liquid of the paraffin wax at $T_{w}=24\ ^{\circ}$C, $Q=0.25$ ml/s,
$h=$ 270 $\mu$m; (b) Temperature profile along the radial direction of the
black arrow in the last snapshot of (a); (c) The pronounced temperature
gradient in the vicinity of the propagation front; (d) For various
experiments, once the front velocity $v_{f}$ reaching down to
$v_{\textup{M}}$, the onset of fingering instability occurs.

_Physical mechanism_ Based on this viewpoint, we propose a physical mechanism
of the observed radial fingering instability by taking into account the
thermal Marangoni force. As shown in the sketch of Fig. 5a, during the stage I
for the onset of the instability ($t<t_{c}$), as $v_{f}$ gradually slows down,
the thermal Marangoni force overcomes the viscous force, and the instability
at the liquid-air interface is triggered. Concurrently, the protrusion
associated with the instability leads to the formation of the pinning
solidification points as subjected to the supercooling locally, and the
spacing between the pinning solidification points is correlated with the
instability wavelength ($\lambda_{exp}$) to be explored hereafter.
Subsequently, during the stage II for the fingering growth ($t>t_{c}$), as the
flow is continuously injected into the cell, the liquid further spreads
outward passing through the pinning solidification points. By this fashion,
the fingers are further extending radially.

In order to validate the proposed mechanism, numerical simulation of the
incompressible Navier-Stokes equations is performed (Supplementary Note 5 and
Supplementary Movie 5). For the stage I (Fig. 5b), the liquid-gas interface
gradually undergoes the circumferential instability as driven by the gradient
of surface tension. Indeed, the Marangoni-induced flow during the onset of the
instability is presented by the magnified view of velocity fields. For the
stage II (Fig. 5c), the fingers grow continuously and extend radially between
the adjacent pinning solidification points, which are formed during the onset
of instability.

_Linear stability analysis_ In order to quantitatively characterize the
underlying mechanism, we perform a linear stability analysis to highlight the
prominent thermal Marangoni effect based on a depth-averaged 2D governing
equations (Supplementary Note 7) al2012control ; de2021marangoni . The
Marangoni flow induced by the surface tension gradient is required to be
explicitly expressed in the dimensionless boundary condition at the liquid-air
interface, _i.e._ , the tangential surface stresses should be balanced by
viscous stresses associated with fluid motion de2021marangoni .

Then based on the linear stability analysis, the dispersion relation is
obtained for the instability growth rate ($\sigma$) dependent on the
dimensionless wavenumber ($k=2\pi r_{in}/\lambda$, $r_{in}$ for the
characteristic length scale),

$\frac{k}{\left(k+\sqrt{k^{2}+\chi^{2}}\right)^{2}}\cdot\frac{1}{\sqrt{k^{2}+\sigma
Pe}}=-\frac{2}{Ma\tilde{T}^{gr}}-\frac{2}{f_{p}RePe}.$ (6)

Here $Ma=-({d\gamma}/{dT})\cdot[{r_{in}(T_{l}-T_{w})}/{\alpha\mu}]$ is for
Marangoni number, $Pe=r_{in}U_{0}/\alpha$ for Peclet number ($\alpha$ for the
thermal diffusivity), and $\tilde{T}^{gr}$ for temperature gradient across the
interface, $\chi$ and $f_{p}$ for a constant coefficient (Supplementary Note
7). Here the radial temperature gradient is obtained from the experimental
observation in Fig. 4c and is applied for Marangoni number and linear
stability analysis (Supplementary Note 6).

Qualitatively, since the left-hand side of equation (6) is positive and the
last term $-{2}/{f_{p}RePe}$ at the right-hand side is negative, the
dispersion relationship holds only if $\tilde{T}^{gr}<0$ for temperature
decreasing with the front propagation at the interface, otherwise the
interface is stable. Correspondingly, a critical magnitude of temperature
gradient ($\tilde{T}_{cr}^{gr}$) is necessary for the possible instability by
setting left side of equation (6) to be
zero,$|\tilde{T}_{cr}^{gr}|=f_{p}RePe/Ma$.

Quantitatively, by numerically solving equation (6), a phase diagram for the
interfacial stability is obtained, where the growth rate
($\sigma(k,\tilde{T}^{gr})$) relies on the wavenumber and temperature gradient
at the interface (Fig. 5d). A blue line for marginal stability ($\sigma=0$)
divides the phase into two regimes: the top-left regime in grey for the stable
interface ($\sigma<0$); the bottom-right regime for the unstable interface
(colorbar for the value of $\sigma>0$ ). The fastest growth rate associated
with the maximum instability $\sigma_{fast}=max\ [\sigma(k,\tilde{T}^{gr})]$
is indicated by a red line. Specially, a dispersion relationship $\sigma(k)$
(Fig. 5e) is presented for a given typical temperature gradient, such as
$|\tilde{T}^{gr}|=32\mathrm{\ K/mm}$ in our experiments, and the wavenumber of
$k_{fast}$ and $k_{cr}$ corresponds to the growth rate of $\sigma_{fast}$ and
$\sigma=0$, respectively.

Additionally, based on the Stefan model for the interfacial temperature
GielenJFM2020 ; hu2020frozen ; lyu2021JCP  (Supplementary Note 8) and scaling
law of $r_{c}$ of Eq. (4), the scaling of temperature gradient at the
interface for the onset of instability is attained,

$\tilde{T}^{gr}|_{r=r_{c}}=\left.\frac{\mathrm{d}T}{\mathrm{~{}d}r}\right|_{r=r_{c}}\sim
r_{c}^{-1}\sim(Qh)^{-1/2}.$ (7)

By combining the above Eq. (7) with the measured temperature gradient at a
given $Q=0.25$ ml/s (Fig. 4c), the relationship between temperature gradient
and flow rate ($Q$) is obtained, $\tilde{T}^{gr}|_{r=r_{c}}(Q)$. Then from the
dispersion relationship of Eq. (6) together with the calculated
$\tilde{T}^{gr}|_{r=r_{c}}(Q)$, the linear stability analysis further predicts
wavelength $\lambda_{cr}$ and $\lambda_{fast}$, corresponding to $\sigma=0$
and $\sigma_{fast}$ (Fig. 5f).

Indeed, as shown in Fig. 5f, the observed wavelength $\lambda_{exp}$ is
greater than $\lambda_{cr}$, as expected. Furthermore, $\lambda_{exp}$ is
quantitatively comparable with the fastest instability growth mode, falling
into the regime ($\lambda_{fast},3\lambda_{fast}$). Considering the simplified
model for such complicated coupling processes of hydrodynamics, heat transfer
and phase change, this remarkable agreement between theory and experiment
clearly validates the proposed mechanism based on the thermal Marangoni
instability which is determined by the ambient cooling at the vertical
direction.

After the onset of instability in stage I, the subsequent growth of fingers
during stage II can be so complicated (due to solidification, heat transfer,
thermal expansion and other factors) that the interface might be deformed into
the irregular morphology, which is beyond the scope of this work.

Figure 5: Fingering instability mechanism and linear stability analysis.(a)
Schematic of the onset of instability (stage I) and subsequently the fingers
growing between the pinning solidification points (stage II); (b,c) Numerical
simulation showing the snapshots for the stage I (b) and stage II (c); (d)
Phase diagram of the growth rate $\sigma$ as a function of $k$ and d$T$/d$r$,
the blue and red curves for $\sigma=0$ and $\sigma_{fast}$; (e) A typical
dispersion relationship, and $k_{fast}$ and $k_{cr}$ for $\sigma_{fast}$ and
$\sigma=0$; (f) Wavelength $\lambda$ dependent on injection flow rates $Q$,
lines for the model, squares for experiments and error bars for five
reproducible experiments.

_Outlook_ This work presents that the strategy of the vertical supercooling
and the structure design to control and steer the liquid flow and the
subsequently correlated instabilities for the potential technological
implications, such as microfluidic cooling for semiconductor devices
microfluidicscoolingnature . Particularly, this supercooling process of the
molten jet and liquid film is inherent within the state-of-the-art thermal
drawing process, which involves heating and the subsequent cooling stage, and
consequently sheds light on the diverse fluid instabilities to fabricate
versatile functional structures for the advanced functional fibres and fabrics
AbouraddyNatureMat2007 ; KaufmanNature2012 ; DengAFM2020 ; DengNatureCom2022 .

###### Acknowledgements.

D. D. is grateful to Prof. Howard Stone, and Prof. Martin Bazant for the
insightful and supportive comment on this work. We also acknowledge Prof. K.
L. Chong for providing the reference de2021marangoni  for the linear stability
analysis. This work is supported by the funding from the National Program in
China and startup in Fudan University.

## References

  * (1) J. S. Langer, Instabilities and pattern formation in crystal growth.  Rev. Mod. Phys. 52, 1 (1980).
  * (2) E. Ben-Jacob, I. Cohen, H. Levine, Cooperative self-organization of microorganisms.  Adv. Phys. 49, 395 (2000).
  * (3) R. V. Craster, O. K. Matar, Dynamics and stability of thin liquid films.  Rev. Mod. Phys. 81, 1131 (2009).
  * (4) H. A. Stone, A. D. Stroock, A. Ajdari, Engineering flows in small devices: microfluidics toward a lab-on-a-chip.  Annu. Rev. Fluid Mech. 36, 381 (2004).
  * (5) A. A. Darhuber, S. M. Troian, Principles of microfluidic actuation by modulation of surface stresses.  Annu. Rev. Fluid Mech. 37, 425 (2005).
  * (6) S. H. Davis, Thermocapillary instabilities.  Annu. Rev. Fluid Mech. 19, 403 (1987).
  * (7) M. F. Schatz, G. P. Neitzel, Experiments on thermocapillary instabilities.  Annu. Rev. Fluid Mech. 33, 93 (2001).
  * (8) W. W. Mullins, R. Sekerka, Stability of a planar interface during solidification of a dilute binary alloy.  J. Appl. Phys. 35, 444 (1964).
  * (9) W. Boos, A. Thess, Thermocapillary flow in a hele-shaw cell.  J. Fluid Mech. 352, 305 (1997).
  * (10) A. Cazabat, F. Heslot, S. Troian, P. Carles, Fingering instability of thin spreading films driven by temperature gradients.  Nature 346, 824 (1990).
  * (11) D. E. Kataoka, S. M. Troian, Patterning liquid flow on the microscopic scale.  Nature 402, 794 (1999).
  * (12) S. Schiaffino and A. A. Sonin, Molten droplet deposition and solidification at low Weber numbers, _Phys. Fluids_ 9, 3172 (1997).
  * (13) E. Ghabache, C. Josserand, T. Séon, Frozen impacted drop: from fragmentation to hierarchical crack patterns.  Phys. Rev. Lett. 117, 074501 (2016).
  * (14) J. de Ruiter, D. Soto, K. K. Varanasi, Self-peeling of impacting droplets.  Nat. Phys. 14, 35 (2018).
  * (15) M. Hu, F. Wang, Q. Tao, L. Chen, S. M. Rubinstein, D. Deng, Frozen patterns of impacted droplets: from conical tips to toroidal shapes.  Phys. Rev. Fluids 5, 081601(R) (2020).
  * (16) P. Kant, H. Müller-Groeling, D. Lohse, Pattern formation during the impact of a partially frozen binary droplet on a cold surface.  Phys. Rev. Lett. 125, 184501 (2020).
  * (17) P. Kant, R. B. J. Koldeweij, K. Harth, M. A. J. van Limbeek, D. Lohse, Fast-freezing kinetics inside a droplet impacting on a cold surface.  Proc. Natl. Acad. Sci. U.S.A. 117, 2788-2794 (2020).
  * (18) R. B. J. Koldeweij, P. Kant, K. Harth, R. de Ruiter, H. Gelderblom, J. H. Snoeijer, D. Lohse, M. A. J. van Limbeek, Initial solidification dynamics of spreading droplets.  Phys. Rev. Fluids 6, L121601 (2021).
  * (19) R. Grivet, A. Monier, A. Huerre, C. Josserand, T. Séon, Contact line catch up by growing ice crystals.  Phys. Rev. Lett. 128, 254501 (2022).
  * (20) H. Lambley, G. Graeber, R. Vogt, _et al._ , Freezing-induced wetting transitions on superhydrophobic surfaces. _Nat. Phys._ 19, 649 (2023).
  * (21) T. T. Al-Housseiny, P. A. Tsai, H. A. Stone, Control of interfacial instabilities using flow geometry.  Nat. Phys. 8, 747 (2012).
  * (22) Geometry parameters of the non-uniform gap as following: $h_{n}=$ 0.1 mm, $h_{w}=$ 0.7 mm, $l=$ 15 cm, $\varphi\sim\Delta h/\Delta l=(h_{w}-h_{n})/l=4\times 10^{-3}\ll 1$.
  * (23) T. M. Koller, T. Klein, C. Giraudet, J. Chen, A. Kalantar, G. P. van der Laan, M. H. Rausch, A. P. Fröba, Liquid viscosity and surface tension of n-dodecane, n-octacosane, their mixtures, and a wax between 323 and 573 K by surface light scattering.  J. Chem. Eng. Data 62, 3319 (2017).
  * (24) T. Klein, S. Yan, J. Cui, J. W. Magee, K. Kroenlein, M. H. Rausch, T. M. Koller, A. P. Fröba, Liquid viscosity and surface tension of n-hexane, n-octane, n-decane, and n-hexadecane up to 573 K by surface light scattering.  J. Chem. Eng. Data 64, 4116 (2019).
  * (25) J. Yang, J. Wu, S. Bi, Surface tension measurements by pendant drop method of 10 pure long-chain alkanes and alcohols for temperatures up to 573.15 K.  J. Chem. Eng. Data 66, 2615 (2021).
  * (26) R. A. Lopez, de la Cruz, C. Diddens, X. Zhang, D. Lohse, Marangoni instability triggered by selective evaporation of a binary liquid inside a Hele-Shaw cell.  J. Fluid Mech. 923, A16 (2021).
  * (27) M. V. Gielen, R. de Ruiter, R. B. Koldeweij, D. Lohse, J. H. Snoeijer, H. Gelderblom, Solidification of liquid metal drops during impact.  J. Fluid Mech. 883, A32 (2020).
  * (28) S. Lyu, K. Wang, Z. Zhang, A. Pedrono, C. Sun, and D. Legendre, A hybrid VOF-IBM method for the simulation of freezing liquid films and freezing drops.  J. Comput. Phys. 432, 110160 (2021).
  * (29) R. Van Erp, R. Soleimanzadeh, L. Nela, G. Kampitsis, E. Matioli, Co-designing electronics with microfluidics for more sustainable cooling.  Nature 585, 211 (2020).
  * (30) A. Abouraddy, M. Bayindir, G. Benoit, S. Hart, K. Kuriki, N. Orf, O. Shapira, F. Sorin, B. Temelkuran, Y. Fink, Towards multimaterial multifunctional fibres that see, hear, sense and communicate.  Nat. Mater. 6, 336 (2007).
  * (31) J. J. Kaufman, G. Tao, S. Shabahang, E.-H. Banaei, D. S. Deng, X. Liang, S. G. Johnson, Y. Fink, A. F. Abouraddy, Structured spheres generated by an in-fibre fluid instability.  Nature 487, 463 (2012).
  * (32) B. R. Xu, S. Q. Ma, Y. Z. Xiang, J. Zhang, M. F. Zhu, L. Wei, G. M. Tao, D. S. Deng, In-fiber structured particles and filament arrays from the perspective of fluid instabilities, _Adv. Fiber. Mater._ 2, 1 (2020).
  * (33) W. Esposito, L. Martin-Monier, P. Piveteau, B. R. Xu, D. S. Deng, F. Sorin, Controlled filamentation instability as a scalable fabrication approach to flexible metamaterials, _Nat. Commun._ 13, 6154 (2022).