# The perturbation method for the skew-symmetric strongly elliptic systems of
PDEs

. . Bagapsh1), 2)

For a Jordan domain with sufficiently smooth boundaries, the solution to the
Dirichlet problem for second order skew-symmetric strongly elliptic system
with constant coefficients and regular enough boundary data is constructed in
the form of a power series of a small parameter describing the perturbation of
the given system from the Laplace one. The coefficients of this series are the
functions that are determined sequentially as solutions to special Dirichlet
problems for the usual Laplace and Poisson equations. The obtained series
converges uniformly in the closure of the domain under consideration.

Keywords: elliptic system; Dirichlet problem; perturbation method.

11footnotetext: Federal Research Center ”‘Computer Science and Control”’ of
the Russian Academy of Sciences, Moscow, Russia; St.Petersburg State
University, St.Petersburg, Russia; Moscow Center for Fundamental and Applied
Mathematics, Moscow, Russia
e-mail<EMAIL_ADDRESS>This research was supported by the
Grants Council of the President of the Russian Federation, grant
MK-1204.2020.1, by Ministry of Science and Higher Education of the Russian
Federation, project 0705-2020-0047, and by the Theoretical Physics and
Mathematics Advancement Foundation BASIS.

##  1\. Introduction and the main result

In the present paper we are dealing with the Dirichlet problem for second
orderhomogeneous elliptic equations with constant _complex_ coefficients in
Jordan domains in the complex plane $\mathbb{C}$. Under certain regularity
assumptions on the boundaries of domains under consideration and on the
boundary functions, we present a solution representation in the form of a
special power series of a small parameter describing the perturbation of a
given equation from the Laplace one, and prove that the corresponding series
converges uniformly in the closure of the domain under consideration.

Consider the following partial differential equation

$af_{xx}+2bf_{xy}+cf_{yy}=0$ (1)

on the complex–valued function $f$ of the complex variable $z=x+iy$ with
constant coefficients $a$, $b$, $c\in\mathbb{C}$. Denote $a=a_{1}+ia_{2}$,
$b=b_{1}+ib_{2}$, $c=c_{1}+ic_{2}$. The equation (1) can be written as a
system of equations

$\left(A\frac{\partial^{2}}{\partial x^{2}}+2B\frac{\partial^{2}}{\partial
x\partial y}+C\frac{\partial^{2}}{\partial
y^{2}}\right)\left(\begin{matrix}u\\\
v\end{matrix}\right)=\left(\begin{matrix}0\\\ 0\end{matrix}\right)$ (2)

on the real and imaginary parts $u=\text{Re}f$, $v=\text{Im}f$, where

$A=\left(\begin{matrix}a_{1}&-a_{2}\\\ a_{2}&a_{1}\end{matrix}\right),\qquad
B=\left(\begin{matrix}b_{1}&-b_{2}\\\ b_{2}&b_{1}\end{matrix}\right),\qquad
C=\left(\begin{matrix}c_{1}&-c_{2}\\\ c_{2}&c_{1}\end{matrix}\right)$ (3)

are the constant real skew-symmetric matrices. Through all this paper the
equation (1) is assumed to be elliptic, which means that

$\det(A\xi^{2}+2B\xi\eta+C\eta^{2})\neq
0\quad\text{when}\quad\xi^{2}+\eta^{2}>0.$

Moreover, we will also assume that the equation (1) is strongly elliptic. The
latter property means that

$\det(A+2\alpha B+\beta C)\neq 0\quad\text{for
each}\quad\alpha\leqslant\beta^{2},$

see, for instance [1]–[4]. The conditions of ellipticity and strong
ellipticity given above are stated in terms of the system (2) that is
equivalent to the equation (1). It is worth to show how these conditions look
line in terms of the initial equation (1). In such terms the ellipticity
property means that both roots of the corresponding characteristic equation
$a\lambda^{2}+2b\lambda+c=0$ are not real, while the strong ellipticity
property means that these roots belong to the different half-planes of
$\mathbb{C}$ with respect to the real axis.

Using some suitable affine transformations of the independent variables ($x$
and $y$) and of the dependent variables ($u$ and $v$), and some suitable
linear combination of equations in system (2), one can transform the initial
system to the system corresponding to the equation $\mathcal{L}f=0$ on new
complex valued function $f$ of the new complex variable $z$ with the operator

$\mathcal{L}=\partial\overline{\partial}+\tau\partial^{2},\qquad\tau\in[0,1),$
(4)

see [4] and [5], where

$\overline{\partial}=\frac{1}{2}\left(\frac{\partial}{\partial
x}+i\frac{\partial}{\partial
y}\right),\qquad\partial=\frac{1}{2}\left(\frac{\partial}{\partial
x}-i\frac{\partial}{\partial y}\right)$

are the standard Cauchy–Riemann differential operator and its complex
conjugate. We will also use the notation $\overline{\partial}_{z}$ and
a$\partial_{z}$ for them. Note that when $\tau=0$ the equation
$\mathcal{L}f=0$ turns out to be the Laplace equation $\Delta f=0$, because
$\Delta=4\partial\overline{\partial}$.

We are interested in the question on solutions construction to the Dirichlet
problem for the equation $\mathcal{L}f=0$ in its classic setting. Let us
recall the statement of this problem.

######  Problem 1.

Let $\varOmega$ be a bounded simply connected domain, and let
$\varGamma=\partial\varOmega$ be its boundary. Given a function $h\in
C(\varGamma)$, find a function $f\in C(\overline{\varOmega})$ such that
$\mathcal{L}f=0$ in $\varOmega$ and $f=h$ on $\varGamma$.

The most important and intriguing question related with Problem 1 is to
describe such domain $\varOmega$ where the corresponding Dirichlet problem for
the equation $\mathcal{L}f=0$ is solvable for every continuous boundary
function $h$. This question remains unsolved in the general setting, even it
is not known, whether Problem 1 is solvable for an arbitrary boundary boundary
function $h$ data in a Jordan domain of general type. The most general known
result (at the best of our knowledge) is that any Dirichlet problem for all
$\mathcal{L}$ under consideration is always solvable in $C^{1}$–polygons, see
[6].

In this paper we are dealing with the special explicit method of solutions
constructing to Problem 1, which is based on the solution representation of
the form

$f=\sum\limits_{n=0}^{\infty}f_{n}\tau^{n}.$ (5)

Using this method we answer the question on solvability of Problem 1 with
Hölder boundary data in Jordan domains with sufficiently smooth boundaries,
and present new explicit formulae for such solutions.

Our main result, Theorem 1, states that if $\varOmega$ is a Jordan domain and
if the function $h\in C(\varGamma)$ is such that $h\circ\omega$ belongs to the
Hölder class $C^{\alpha}(\mathbb{T})$ for some $1/2<\alpha<1$ and some
conformal mapping $\omega$ of the unit disk $\mathbb{D}$ to $\varOmega$, then
for each $\tau\in[0,1)$ there exists a solution of the Dirichlet problem that
is represented in the form (5), converging uniformly on the closure
$\overline{\varOmega}$ of $\varOmega$.

Observe that the condition $h\circ\omega\in C^{\alpha}(\mathbb{T})$ is
satisfied, for instance, when $h\in C^{\alpha}(\varGamma)$ and $\varGamma$ is
a Diny smooth curve, where $\varGamma$ is the boundary of $\varOmega$, as
above. Indeed, by Lindelöf’s theorem $\omega^{\prime}\in
C(\overline{\mathbb{D}})$, so for each $z_{1}$, $z_{2}\in\mathbb{T}$ it holds

$|h\circ\omega(z_{1})-h\circ\omega(z_{2})|\leqslant[h]_{\alpha}\,|\omega(z_{1})-\omega(z_{2})|^{\alpha}\leqslant[h]_{\alpha}\,\|\omega^{\prime}\|_{C(\overline{\mathbb{D}})}^{\alpha}\,|z_{1}-z_{2}|^{\alpha},$

where
$[h]_{\alpha}=\sup_{\zeta_{1}\neq\zeta_{2}}|h(\zeta_{1})-h(\zeta_{2})|/|\zeta_{1}-\zeta_{2}|^{\alpha}$.

The paper has the following structure. In Section 2 we explain the
aforementioned method of solutions constructing in the form of representations
(5). Further, in Section 3 we prove the convergence of the series (5) under
assumptions to $h$ and $\varOmega$ mentioned above.

##  2\. Perturbation method

Let $\Omega$ be a Jordan domain with the boundary $\varGamma$. Take a function
$h\in C(\varGamma)$. We are going to obtain the function $f$ of the form (5)
such that $\mathcal{L}f=0$ and $f|_{\varGamma}=h$. To do this we substitute
the expansion (5) into the equation $\mathcal{L}f=0$ and equate to zero the
factors $f_{n}$ for the same powers of the parameter $\tau$. Moreover, we put
$f_{0}|_{\varGamma}=h$ and $f_{n}|_{\varGamma}=0$ for $n\geq 1$. Therefore,
the conditions on $f_{n}$, $n\geq 0$, are reduced to the following sequence of
boundary value problems:

$\partial\overline{\partial}f_{0}=0\quad\text{in}\quad\varOmega,\qquad\qquad
f_{0}|_{\varGamma}=h$ (6)

and

$\partial\overline{\partial}f_{n}=-\partial^{2}f_{n-1}\quad\text{in}\quad\varOmega,\qquad\qquad
f_{n}|_{\varGamma}=0$ (7)

for $n\geqslant 1$.

Let $\omega\colon\mathbb{D}\to\varOmega$ be some conformal mapping from the
unit disk $\mathbb{D}=\\{z\in\mathbb{C}\colon|z|<1\\}$ onto $\varOmega$. If
$\varOmega$ is Jordan then, according to the Carathéodory extension theorem
[9, Theorem 2.6] that $\omega$ is extended to the eponymous homeomorphism from
$\overline{\mathbb{D}}$ onto $\overline{\varOmega}$. Let us move the problems
(6) and (7), $n\geq 1$, to the disk $\mathbb{D}$. Define the functions

$F:=f\circ\omega,\qquad H:=h\circ\omega\qquad F_{n}:=F_{n}\circ\omega.$

Then we have

$\mathcal{L}f=\frac{1}{|\omega^{\prime}|^{2}}\left[\partial\overline{\partial}F+\tau\partial\left(\frac{\overline{\omega^{\prime}}}{\omega^{\prime}}\partial
F\right)\right]=:\mathcal{M}F.$ (8)

Equations (5), (6) and (7) yield

$F=\sum\limits_{n=0}^{\infty}F_{n}\tau^{n},$ (9)

where

$\partial\overline{\partial}F_{0}=0\quad\text{in}\quad\mathbb{D},\qquad\qquad
F_{0}|_{\mathbb{T}}=H$ (10)

and

$\partial\overline{\partial}F_{n}=-\partial\left(\frac{\overline{\omega^{\prime}}}{\omega^{\prime}}\partial
F_{n-1}\right)\quad\text{in}\quad\varOmega,\qquad\qquad F_{n}|_{\mathbb{T}}=0$
(11)

for $n\geqslant 1$. Assuming that the functions $h$ and $F_{n-1}$ (for some
$n>1$) are sufficiently regular, we can write out the solutions for (10) and
(11):

$F_{0}(z)=\frac{1}{2i}\int_{\mathbb{T}}\partial_{\zeta}G(\zeta,z)h(\zeta)d\zeta,\qquad
F_{n}(z)=-\int_{\mathbb{D}}G(\zeta,z)\partial\left(\frac{\overline{\omega^{\prime}(\zeta)}}{\omega^{\prime}(\zeta)}\partial
F_{n-1}(\zeta)\right)d\mu,$ (12)

$n\geqslant 1$. Here

$G(\zeta,z)=\frac{2}{\pi}\log\left|\frac{\zeta-z}{1-\zeta\overline{z}}\right|$
(13)

is the Green’s function for the operator
$\partial\overline{\partial}=\Delta/4$ in the disk $\mathbb{D}$ and $\mu$ is
the planar Lebesgue measure.

Let us define the operator $\mathcal{P}$ acting on $C(\mathbb{T})$ and the
operators $\mathcal{K}$, $\mathcal{K}_{z}$, $\mathcal{K}_{\overline{z}}$
acting on $L_{p}(\mathbb{D})$ as follows:

$\mathcal{P}[\varphi(z)]:=\frac{1}{2i}\int_{\mathbb{T}}\partial_{\zeta}G(\zeta,z)\varphi(\zeta)d\zeta,\qquad\mathcal{K}[\varphi(z)]:=\int_{\mathbb{D}}\partial_{\zeta}G(\zeta,z)\varphi(\zeta)d\mu$

and

$\mathcal{K}_{z}[\varphi(z)]:=\text{v.p.}\int_{\mathbb{D}}\partial_{z}\partial_{\zeta}G(\zeta,z)\varphi(\zeta)d\mu,\qquad\mathcal{K}_{\overline{z}}[\varphi(z)]:=\text{v.p.}\int_{\mathbb{D}}\partial_{\overline{z}}\partial_{\zeta}G(\zeta,z)\varphi(\zeta)d\mu.$

The symbol $\text{v.p.}\int$ means that the corresponding integrals are
understood in the sense of Cauchy principal values. For the sake of brevity we
omit this symbol in what follows everywhere, where it will not cause any
misunderstanding. From (13) we obtain

$\displaystyle\mathcal{P}[\varphi(z)]=\frac{1}{2\pi}\int_{\mathbb{T}}\left(\frac{1}{\zeta-z}+\frac{\overline{z}}{1-\zeta\overline{z}}\right)\varphi(\zeta)d\zeta,$
(14)
$\displaystyle\mathcal{K}[\varphi(z)]=\frac{1}{\pi}\int_{\mathbb{D}}\left(\frac{1}{\zeta-z}+\frac{\overline{z}}{1-\zeta\overline{z}}\right)\varphi(\zeta)d\mu$

and

$\mathcal{K}_{z}[\varphi(z)]=\frac{1}{\pi}\int_{\mathbb{D}}\frac{\varphi(\zeta)d\mu}{(\zeta-z)^{2}},\qquad\mathcal{K}_{\overline{z}}[\varphi(z)]=\frac{1}{\pi}\int_{\mathbb{D}}\frac{\varphi(\zeta)d\mu}{(1-\zeta\overline{z})^{2}}.$
(15)

In terms of the operators defined above, the formulae (12) for constructing
the functions $F_{n}$ can be written in the short form

$F_{0}=\mathcal{P}[h],\qquad
F_{n}=\mathcal{K}[(\overline{\omega^{\prime}}/\omega^{\prime})\partial
F_{n-1}],\quad n\geqslant 1.$ (16)

We also need the partial sums of series (5) and (9) respectively:

$s_{m}=\sum_{n=0}^{m}f_{n}\tau^{n},\qquad S_{m}=\sum_{n=0}^{m}F_{n}\tau^{n}.$
(17)

Now we are ready to state our main result, which is the following theorem.

######  Theorem 1.

Let $\varOmega$ be a Jordan domain with the boundary $\varGamma$, let $\omega$
be some conformal mapping from $\mathbb{D}$ onto $\varOmega$, and let $h\in
C(\varGamma)$. Suppose $\varOmega$ and $h$ are such that $h\circ\omega\in
C^{\alpha}(\mathbb{T})$ for some $1/2<\alpha<1$. Then for every $\tau\in[0,1)$
the series (5), where $f_{n}=F_{n}\circ\omega^{-1}$, $n\geq 0$, and where
$F_{n}$ are defined by (16), converges uniformly on $\overline{\varOmega}$ to
the function $f\in C(\overline{\varOmega})$ that satisfies the equation
$\mathcal{L}f=0$ in $\varOmega$ and coincides with $h$ on $\varGamma$.

##  3\. Auxiliary lemmas and proof of Theorem 1

Let us recall, that $L_{p}(\mathbb{D})$ is the usual Lebesgue space of
functions considered with respect to the planar Lebesgue measure on
$\mathbb{D}$, and $W^{1}_{p}(\mathbb{D})$ is the standard Sobolev space of
functions in $\mathbb{D}$. We start with the following technical lemma.

######  Lemma 1.

If $\varphi\in C^{\alpha}(\mathbb{T})$, where $0<\alpha<1$, then
$\mathcal{P}\varphi\in W^{1}_{p}(\mathbb{D})$ with any exponent
$0<p<(1-\alpha)^{-1}$.

###### Proof.

Take $\psi=\mathcal{P}\varphi$. Since $\varphi\in C(\mathbb{T})$, then by the
properties of the Poisson integral, we have $\psi\in
C(\overline{\mathbb{D}})$, so obviously $\psi\in L_{p}(\mathbb{D})$. We are
going to prove the $L_{p}$–integrability of the first partial derivatives of
$\psi$. We represent the function $\psi$ as the sum
$\psi(z)=\psi_{1}(z)+\overline{\psi_{2}(z)}$ of holomorphic components

$\psi_{1}(z)=\frac{1}{2\pi
i}\int_{\mathbb{T}}\frac{\varphi(\zeta)d\zeta}{\zeta-z},\qquad\psi_{2}(z)=-\frac{1}{2\pi
i}\int_{\mathbb{T}}\frac{\overline{\varphi(\zeta)}d\zeta}{\zeta-z}-\frac{1}{2\pi}\int_{\mathbb{T}}\overline{\varphi(\zeta)}|d\zeta|.$

Since $\varphi\in C^{\alpha}(\mathbb{T})$ with $0<\alpha<1$, then $\psi_{1}\in
C^{\alpha}(\overline{\mathbb{D}})$, because of the Privalov theorem [7] for
Cauchy type integral. Put
$T(z,r):=\\{\zeta\in\mathbb{C}\colon|\zeta-z|=r\\}\subset\mathbb{D}$. Using
the Cauchy formula

$\psi_{1}(z)=\frac{1}{2\pi
i}\int_{T(z,r)}\frac{\psi_{1}(\zeta)d\zeta}{\zeta-z}$

we obtain the equality

$\partial\psi(z)=\psi_{1}^{\prime}(z)=\frac{1}{2\pi
i}\int_{T(z,r)}\frac{\psi_{1}(\zeta)d\zeta}{(\zeta-z)^{2}}=\frac{1}{2\pi
i}\int_{T(z,r)}\frac{\psi_{1}(\zeta)-\psi_{1}(z)}{(\zeta-z)^{2}}d\zeta,$

and, furthermore, the estimate

$|\partial\psi(z)|\leqslant\frac{1}{2\pi
i}\int_{T(z,r)}\frac{|\psi_{1}(\zeta)-\psi_{1}(z)|}{|\zeta-z|^{2}}|d\zeta|\leqslant\frac{1}{2\pi
i}\int_{T(z,r)}\frac{c_{\alpha}|\zeta-z|^{\alpha}}{|\zeta-z|^{2}}|d\zeta|=\frac{c_{\alpha}}{r^{1-\alpha}},$

where $c_{\alpha}=\sup_{\zeta\neq
z}|\psi_{1}(\zeta)-\psi_{1}(z)|/|\zeta-z|^{\alpha}$. Taking the limit
$r\to(1-|z|)$ we obtain $|\partial\psi(z)|\leqslant
c_{\alpha}(1-|z|)^{\alpha-1}$, see also [8, page 74] or [9, page 50]. This
means that $\partial\psi\in L_{p}(\mathbb{D})$ if $(1-\alpha)p<1$. Similarly,
the $L_{p}$–integrability of the derivative
$\overline{\partial}\psi=\overline{\psi_{2}^{\prime}}$ is established under
the same condition on $p$. The lemma is proved. ∎

Consider the operator

$K\varphi(z):=\frac{1}{\pi}\int_{\mathbb{C}}\frac{\varphi(\zeta)d\mu}{(\zeta-z)^{2}},$
(18)

which is a bounded one from $L_{p}(\mathbb{C})$ to itself for $p\in(1,\infty)$
by the celebrated theorem of Calderon and Zygmund [10]. We denote by
$\|K\|_{p}$ the norm of the operator $K$ in $L_{p}(\mathbb{C})$ and similarly
we denote the norms of operators acting in $L_{p}(U)$ with arbitrary domain
$U$. It is important for us that $\|K\|_{p}\to 1$ when $p\to 2$, see [11, page
89], [12, pages 5–6].

######  Lemma 2.

The operators (14), (15) have the following properties:

(i) $\mathcal{K}\colon L_{p}(\mathbb{D})\to L_{p}(\mathbb{D})$ is bounded for
$p>1$;

(ii) $\mathcal{K}_{z}\colon L_{p}(\mathbb{D})\to L_{p}(\mathbb{D})$ is bounded
for $p>1$, and $\|\mathcal{K}_{z}\|_{p}=\|K\|_{p}$;

(iii) $\mathcal{K}_{\overline{z}}\colon L_{p}(\mathbb{D})\to
L_{p}(\mathbb{D})$ is bounded for $p\geqslant 2$; moreover,
$\|\mathcal{K}_{\overline{z}}\|_{p}\leqslant\|K\|_{p}$ and, in particular,
$\|\mathcal{K}_{\overline{z}}\|_{2}=\|K\|_{2}$.

###### Proof.

The statement (i) follows from the fact that the kernel of the integral
$\mathcal{K}$ is the sum of two kernels with a weak singularity, that is they
have the growth of the first order.

(ii) Take $\varphi\in L_{p}(\mathbb{D})$, $p>1$. We have
$\mathcal{K}_{z}\varphi=K\varphi_{1}$, where the function $\varphi_{1}$
coinsides with $\varphi$ in the unit disk $\mathbb{D}$ and equals to zero
outside $\mathbb{D}$, thus
$\|\varphi_{1}\|_{L_{p}(\mathbb{C})}=\|\varphi\|_{L_{p}(\mathbb{D})}$.
Therefore, $\|\mathcal{K}_{z}\|_{p}=\|K\|_{p}$.

(iii) Making in (15) the change of variables $\zeta\to\xi=1/\overline{\zeta}$,
we obtain

$\mathcal{K}_{\overline{z}}[\varphi(z)]=\frac{1}{\pi}\int_{\mathbb{C}\setminus\overline{\mathbb{D}}}\frac{\varphi(1/\overline{\xi})}{\xi^{2}}\cdot\frac{d\mu}{(\overline{\xi}-\overline{z})^{2}}=\overline{K[\varphi_{2}(z)]},$
(19)

where $\varphi_{2}(z)=\overline{\varphi(1/\overline{z})}/\overline{z}^{2}$ for
$z\in\mathbb{C}\setminus\overline{\mathbb{D}}$ and $\varphi_{2}(z)=0$ for
$z\in\mathbb{D}$. Then, by means of reverse change
$\xi\to\zeta=1/\overline{\xi}$ we find

$\|\varphi_{2}\|_{L_{p}(\mathbb{D})}=\left(\int_{\mathbb{C}\setminus\overline{\mathbb{D}}}|\varphi_{2}(\xi)|^{p}d\mu\right)^{\frac{1}{p}}=\left(\int_{\mathbb{D}}|\xi|^{2p-4}\cdot|\varphi(\zeta)|^{p}d\mu\right)^{\frac{1}{p}}\leqslant\|\varphi\|_{L_{p}(\mathbb{D})}$

for $p\geqslant 2$ with the equality for $p=2$. Hence, it follows from (19)
that
$\|\mathcal{K}_{\overline{z}}\|_{L_{p}(\mathbb{D})}\leqslant\|K\|_{p}\cdot\|\varphi\|_{L_{p}(\mathbb{D})}$,
i.e. $\|\mathcal{K}_{\overline{z}}\|_{p}\leqslant\|K\|_{p}$ for $p\geqslant 2$
with equality for $p=2$. The proposition is proved. ∎

###### Proof of Theorem 1.

First we are going to establish the convergence of the series (9) and its
partial derivatives of the first order in the norm of $L_{p}(\mathbb{D})$.
Lemma 1 yields that $F_{0}=\mathcal{P}h$ belongs to the Sobolev space
$W^{1}_{p}(\mathbb{D})$ for $p<(2(1-\alpha))^{-1}$. Suppose that $F_{n-1}\in
L_{p}(\mathbb{D})$, $p>2$, for some number $n$. Then

$\partial
F_{n}=\mathcal{K}_{z}[(\overline{\omega^{\prime}}/\omega^{\prime})\partial
F_{n-1}],\qquad\overline{\partial}F_{n}=(\mathcal{I}+\mathcal{K}_{\overline{z}})[(\overline{\omega^{\prime}}/\omega^{\prime})\partial
F_{n-1}]$

in the sense of distributions (see [11, page 90]). Using these equalities and
applying Lemma 2, we derive from (16) that the following estimates hold

$\|\partial F_{n}\|_{L_{p}(\mathbb{D})}\leqslant\|K\|_{p}\cdot\|\partial
F_{n-1}\|_{L_{p}(\mathbb{D})},\qquad\|\overline{\partial}F_{n}\|_{L_{p}(\mathbb{D})}\leqslant(1+\|K\|_{p})\|\partial
F_{n-1}\|_{L_{p}(\mathbb{D})}.$

This implies

$\|\partial F_{n}\|_{L_{p}(\mathbb{D})}\leqslant\|K\|_{p}^{n}\cdot\|\partial
F_{0}\|_{L_{p}(\mathbb{D})},\quad\|\overline{\partial}F_{n}\|_{L_{p}(\mathbb{D})}\leqslant(1+\|K\|_{p})\|K\|_{p}^{n-1}\cdot\|\partial
F_{0}\|_{L_{p}(\mathbb{D})}.$ (20)

Then (16) yields

$\|F_{n}\|_{L_{p}(\mathbb{D})}\leqslant\|P\|_{p}\cdot\|\partial
F_{n-1}\|_{L_{p}(\mathbb{D})}\leqslant\|P\|_{p}\cdot\|K\|_{p}^{n-1}\cdot\|\partial
F_{0}\|_{L_{p}(\mathbb{D})}.$ (21)

Estimate (21) gives the convergence of the series (9) in $L_{p}(\mathbb{D})$
to its sum $F\in L_{p}(\mathbb{D})$, when $\|K\|_{p}\tau<1$. Moreover,

$\|F\|_{L_{p}(\mathbb{D})}=\left\|\sum_{n=0}^{\infty}F_{n}\tau^{n}\,\right\|_{L_{p}(\mathbb{D})}\leqslant\|F_{0}\|_{L_{p}(\mathbb{D})}+\sum_{n=1}^{\infty}\|F_{n}\|_{L_{p}(\mathbb{D})}\tau^{n}\leqslant\\\
\leqslant\|F_{0}\|_{L_{p}(\mathbb{D})}+\sum_{n=1}^{\infty}\|P\|_{p}\cdot\|K\|_{p}^{n-1}\cdot\|\partial
F_{0}\|_{L_{p}(\mathbb{D})}\cdot\tau^{n}=\\\
=\|F_{0}\|_{L_{p}(\mathbb{D})}+\frac{\tau\|P\|_{p}}{1-\tau\|K\|_{p}}\|\partial
F_{0}\|_{L_{p}(\mathbb{D})}.$ (22)

Estimate (20) gives the convergence of the first partial derivatives of the
series (9) to the corresponding derivatives of $F$ in the same space
$L_{p}(\mathbb{D})$:

$\|\partial F\|_{L_{p}(\mathbb{D})}\leqslant\frac{\|\partial
F_{0}\|_{L_{p}(\mathbb{D})}}{1-\|K\|_{p}\tau},\qquad\|\overline{\partial}F\|_{L_{p}(\mathbb{D})}\leqslant\|\overline{\partial}F_{0}\|_{L_{p}(\mathbb{D})}+\frac{\tau(1+\|K\|_{p})\,\|\partial
F_{0}\|_{L_{p}(\mathbb{D})}}{1-\tau\|K\|_{p}}.$ (23)

The obtained estimates (22) and (23) mean the convergence of the series (9) in
the norm of the Sobolev space $W_{p}^{1}(\mathbb{D})$:

$\lim\limits_{m\to\infty}\|F-S_{m}\|_{W_{p}^{1}(\mathbb{D})}=0.$ (24)

By virtue of Sobolev embedding theorem [13], we have
$W_{p}^{1}(\mathbb{D})\subset C(\overline{\mathbb{D}})$ when $1-2/p>0$, or,
equivalently, $p>2$, and the embedding is compact. As $F\in
W_{p}^{1}(\mathbb{D})$, so $F\in C(\overline{\mathbb{D}})$ and thus
$f=F\circ\omega^{-1}\in C(\overline{\varOmega})$. In view of the compactness
of the embedding, the series (9), and therefore (5), converge uniformly in
$\overline{\mathbb{D}}$ and $\overline{\varOmega}$ respectively.

Now let us prove that the function $f=F\circ\omega^{-1}$ satisfies the
equation $\mathcal{L}f=0$ in $\varOmega$ by establishing it firstly in the
sense of distributions.

Let $\phi$ be arbitrary test function of the class $C^{2}_{0}(\mathbb{D})$ of
twice continuously differentiable functions in $\mathbb{C}$ with compact
support in $\mathbb{D}$ and let

$\langle g|\phi\rangle:=\int_{\mathbb{C}}g(z)\phi(z)d\mu$

be the action of the distribution $g$ on $\phi$. We have

$\langle
F_{n}|\partial\overline{\partial}\phi\rangle=\int_{\mathbb{D}}\partial\overline{\partial}\phi(z)d\mu(z)\int_{\mathbb{D}}\partial_{\zeta}G(\zeta,z)\frac{\overline{\omega^{\prime}(\zeta)}}{\omega^{\prime}(\zeta)}\partial
F_{n-1}(\zeta)d\mu(\zeta)=\\\
=\int_{\mathbb{D}}\frac{\overline{\omega^{\prime}(\zeta)}}{\omega^{\prime}(\zeta)}\partial
F_{n-1}(\zeta)d\mu(\zeta)\partial_{\zeta}\int_{\mathbb{D}}G(\zeta,z)\partial\overline{\partial}\phi(z)d\mu(z)=\\\
=\int_{\mathbb{D}}\frac{\overline{\omega^{\prime}(\zeta)}}{\omega^{\prime}(\zeta)}\partial
F_{n-1}(\zeta)\partial\phi(\zeta)d\mu(\zeta)=\langle(\overline{\omega^{\prime}}/\omega^{\prime})\partial
F_{n-1}|\partial\phi\rangle,$

which means that the equality
$\partial\overline{\partial}F_{n}=-\partial[(\overline{\omega^{\prime}}/\omega^{\prime})\partial
F_{n-1}]$ holds in $\mathbb{D}$ in the sense of distributions. This together
with (8) implies the following chain of equalities that are understood in the
sense of distributions:

$\displaystyle|\omega^{\prime}|^{2}\mathcal{L}s_{m}$
$\displaystyle=\sum_{n=0}^{m}\left(\partial\overline{\partial}F_{n}+\tau\partial\left(\frac{\overline{\omega^{\prime}}}{\omega^{\prime}}\partial
F_{n}\right)\right)\tau^{n}=$
$\displaystyle=\partial\overline{\partial}F_{0}+\sum_{n=1}^{m}\left(\partial\overline{\partial}F_{n}+\partial\left(\frac{\overline{\omega^{\prime}}}{\omega^{\prime}}\partial
F_{n-1}\right)\right)\tau^{n}+\partial\left(\frac{\overline{\omega^{\prime}}}{\omega^{\prime}}\partial
F_{m}\right)\tau^{m+1}=$
$\displaystyle=\partial\left(\frac{\overline{\omega^{\prime}}}{\omega^{\prime}}\partial
F_{m}\right)\tau^{m+1}.$

Thus for each function $\varphi\in C_{0}^{2}(\varOmega)$, one can write

$\langle\mathcal{L}s_{m}\,|\,\varphi\rangle=\langle\partial\left[(\overline{\omega^{\prime}}/\omega^{\prime})\partial
F_{m}\right]\,|\,\phi\rangle\cdot\tau^{m+1}=-\langle\partial
F_{m}\,|\,(\overline{\omega^{\prime}}/\omega^{\prime})\partial\phi\rangle\cdot\tau^{m+1},$

where $\phi:=\varphi\circ\omega\in C_{0}^{2}(\mathbb{D})$. Then

$\langle\mathcal{L}f\,|\,\varphi\rangle:=\langle
f|\mathcal{L}\varphi\rangle=\langle f-s_{m}|\mathcal{L}\varphi\rangle+\langle
s_{m}|\mathcal{L}\varphi\rangle=\\\ =\langle
F-S_{m}|\mathcal{M}\phi\rangle-\langle\partial
F_{m}|(\overline{\omega^{\prime}}/\omega^{\prime})\partial\phi\rangle\cdot\tau^{m+1}.$

Take $1/p+1/q=1$. Applying the Hölder inequality ant taking into account (20)
and (24) we obtain

$\left|\langle\mathcal{L}f\,|\,\varphi\rangle\right|\leqslant\|F-S_{m}\|_{L_{p}(\mathbb{D})}\cdot\|\mathcal{M}\phi\|_{L_{q}(\mathbb{D})}+\|\partial
F_{m}\|_{L_{p}(\mathbb{D})}\cdot\|\partial\phi\|_{L_{q}(\mathbb{D})}\cdot\tau^{m+1}\leqslant\\\
\leqslant\|F-S_{m}\|_{L_{p}(\mathbb{D})}\cdot\|\mathcal{M}\phi\|_{L_{q}(\mathbb{D})}+\|\partial
F_{0}\|_{L_{p}(\mathbb{D})}\cdot\|\partial\phi\|_{L_{q}(\mathbb{D})}\cdot\|K\|_{p}^{m}\tau^{m+1}\to
0$

when $m\to\infty$ and $\|K\|_{p}\tau<1$. Thus
$\langle\mathcal{L}f\,|\,\varphi\rangle=0$, i.e the function $f$ satisfies the
equation $\mathcal{L}f=0$ in $\varOmega$ in the sense of distributions. Since
this equation is of elliptic type, by virtue of Wheil’s lemma it also holds in
the classic sense.

It remains to show that $f|_{\Gamma}=h$. Note that the estimates (20) and (21)
yield $F_{n}\in W_{p}^{1}(\mathbb{D})$. If $p>2$, by the Sobolev embedding
theorem, this implies that $F_{n}\in C(\overline{\mathbb{D}})$. Since $F_{0}$
is the harmonic extension of the boundary function $H\in
C^{\alpha}(\mathbb{T})$, then $F_{0}|_{\mathbb{T}}=H$. The rest of the
functions $F_{n}$ vanish on $\mathbb{T}$: it follows from equations (12) and
from the continuity of $F_{n}$ in $\overline{\mathbb{D}}$. Thus
$S_{m}|_{\mathbb{T}}=H$ for all $m$. The compactness of the embedding
$W_{p}^{1}(\mathbb{D})\subset C(\overline{D})$ and the convergence (24) lead
to uniform convergence $\|F-S_{m}\|_{C(\overline{\mathbb{D}})}\to 0$, so that
$F|_{\mathbb{T}}=S_{m}|_{\mathbb{T}}=H$. Consequently, $f|_{\Gamma}=h$.

The given arguments are valid when $2<p<(1-\alpha)^{-1}$, which is true since
it is assumed that $1/2<\alpha<1$. As $\|K\|_{p}\to 1$ when $p\to 2$, then for
each fixed value of $\tau<1$ one can choose a value of $p>2$ sufficiently
close to $2$ such that $\|K\|_{p}\tau<1$. This finishes the proof of the
theorem. ∎

## Acknowledgement(s)

I am grateful to professors K.Yu. Fedorovskiy and A.P. Soldatov for helpful
discussions.

## References

  * [1] Petrovskiy I.G. On the analiticity of solutions of PDE. Matem. Sbornik. 1939\. V. 5. P. 3–70 (in Russian)
  * [2] Vishik M.I. Strongly elliptic systems of differential equations. Math. Sbornik. 1951\. V. 29. No 3. P. 615–676 (in Russian).
  * [3] Hua L.K., Lin W., Wu C.Q. On the uniqueness of the solution of the Dirichlet problem of the elliptic system of differential equations. Acta Math. Sinica. 1965\. No 15(2).
  * [4] Hua L.K., Lin W., Wu C.Q. Second-order systems of partial differential equations in the plane. Boston, London, Melbourne: Pitman Advanced Publishing Program. 1985.
  * [5] Bagapsh A.O., Fedorovskiy K.Yu. $C^{1}$ approximation of functions by solutions of second-order elliptic systems on compact sets in $\mathbb{R}^{2}$. Proc. Steklov Inst. Mathm. 2017\. V. 298. P. 35–50.
  * [6] Verchota G.C., Vogel A.L. Nonsymmetric systems on nonsmooth planar domains. Trans. Amer. Math. Soc. 1997\. V. 349. No 11. P. 4501–4535.
  * [7] Privalov I.I. On Cauchy type integrals. Dokl. Academ. nauk 1939\. V. 23. No 9. P. 859–862 (in Russian).
  * [8] Duren P. Theory of $H^{p}$ spaces. Academic Press, New York. 1970\.
  * [9] Pommerenke Ch. Boundary behavior of conformal maps. Springer–Verlag Berlin Heidelberg. 1992\.
  * [10] Calderon A., Zigmund A. On the existence of certain singular integrals. Acta Math. 1952\. V. 88. P. 85–139.
  * [11] Ahlfors L.V. Lectures on quasiconformal mappings. Van Nostrand, Princeton, NJ, 1966.
  * [12] Christ M. Lectures on singular integral operators. CBMS Regional Conference Series in Mathematics, Vol.77, Amer. Math. Soc., Providence, RI, 1990.
  * [13] Sobolev S.L. Some applications of functional analysis in mathematical physics. Amer. Math. Soc., Providence, RI, 1991.


# Investigation of the Nonlinear Optical Frequency Conversion in Ultrathin
Franckeite Heterostructures

Alisson R. Cadore  Mackenzie Engineering School, Mackenzie Presbyterian
University, Sao Paulo, 01302907, SP, Brazil.  Alexandre S. M. V. Ore
Mackenzie Engineering School, Mackenzie Presbyterian University, Sao Paulo,
01302907, SP, Brazil.  David Steinberg  MackGraphe, Mackenzie Presbyterian
Institute, Sao Paulo, 01302907, SP, Brazil.  Juan D. Zapata  Faculty of
Engineering, Universidad del Antioquia UdeA, Medellín, Colombia  Eunézio A. T.
de Souza  Mackenzie Engineering School, Mackenzie Presbyterian University, Sao
Paulo, 01302907, SP, Brazil.  Dario A. Bahamon  Mackenzie Engineering School,
Mackenzie Presbyterian University, Sao Paulo, 01302907, SP, Brazil.
Christiano J. S. de Matos  Mackenzie Engineering School, Mackenzie
Presbyterian University, Sao Paulo, 01302907, SP, Brazil.
<EMAIL_ADDRESS>
###### Abstract

Layered franckeite is a natural superlattice composed of two alternating
layers of different compositions, SnS2\- and PbS-like. This creates
incommensurability between the two species along the planes of the layers,
resulting in spontaneous symmetry-break periodic ripples in the a-axis
orientation. Nevertheless, natural franckeite heterostructure has shown
potential for optoelectronic applications mostly because it is a semiconductor
with 0.7 eV bandgap, air-stable, and can be easily exfoliated down to
ultrathin thicknesses. Here, we demonstrate that few-layer franckeite shows a
highly anisotropic nonlinear optical response due to its lattice structure,
which allow for the identification of the ripple axis. Moreover, we find that
the highly anisotropic third-harmonic emission strongly varies with material
thickness. These features are further corroborated by a theoretical nonlinear
susceptibility model and the nonlinear transfer matrix method. Overall, our
findings help to understand this material and propose a characterization
method that could be used in other layered materials and heterostructures to
assign their characteristic axes.

Brazilian Nanotechnology National Laboratory, Brazilian Center for Research in
Energy and Materials, Campinas, 13083-100, SP, Brazil. These authors
contributed equally to this work. MackGraphe, Mackenzie Presbyterian
Institute, Sao Paulo, 01302907, SP, Brazil. These authors contributed equally
to this work. MackGraphe, Mackenzie Presbyterian Institute, Sao Paulo,
01302907, SP, Brazil. MackGraphe, Mackenzie Presbyterian Institute, Sao Paulo,
01302907, SP, Brazil.

##  1 Introduction

Layered materials (LMs) have shown interesting electrical, optical, and
mechanical properties1, 2, 3, 4 and are of great technological importance5, 6.
Many of these LMs are naturally occurring7, 8, 9, 10, 11, 12, 13, while others
can be synthesized14, 15, 16, 17. Two or more LMs can, nowadays, be assembled
together by deterministic stacking or directly grown by chemical vapour
deposition (CVD) and epitaxial growth techniques, creating the so-called
vertical van der Waals heterostructures (vdWHs)18, 19. Such vdWHs offer a
platform to obtain unprecedented functionalities and properties that do not
exist in the individual LMs20, 21. Albeit these vdWHs are well-studied and
applied22, 23, 24, 25, 19, 26, the fabrication process to create these
nanostructures suffers from some drawbacks such as the need for expensive
systems for the growth, or the presence of contaminants between the layers,
besides the difficulty of depositing the LMs with high-precision, in the
manual stacking case21, 27, 28, 18, 29. To solve these issues, mother nature
has given us many minerals that are naturally occurring vdWHs. These vdWHs are
formed with alternating layers of two different LM due to the phase
segregation process during their formation, and thus are free of the
aforementioned problems. So far, several exfoliated and naturally occurring
vdWHs have been identified, such as cylindrite30, 31, 32, lengenbachite33,
levyclaudite34, nagayagile33, teallite35, 36, to name a few, and the most
notably studied, franckeite37, 38, 39, 40, 41, 42, 43, 44, 45, 46.

Franckeite is a structure composed of alternating stacks of PbS-type
pseudotetragonal (Q) layers and SnS2-type pseudohexagonal (H) layers, as shown
in Fig. 2a37, 47, 48. So far, it has been shown that the alignment of each
franckeite monolayer (Q-H layers) is preserved in the exfoliated flakes (i.e.
Q-H/Q-H/…/Q-H)42, 43. Nevertheless, due to its natural occurrence, the
material may have extra impurities in the crystal lattice which can cause Q
and H layer stacking faults (i.e. Q-H/H/Q-H or Q-H/Q/Q-H…) which are random in
the material, consequently, modifying the electronic properties of the
material on a specific flake location40. However, to the best of our
knowledge, the direct correlation between stacking faults and material
properties was not assigned, and this could explain some data deviations
between research groups working on naturally obtained LMs. Franckeite is also
a p-type semiconductor with a bandgap of approximately 0.7 eV45, 42, 43 and it
has gained notoriety for its potential for optoelectronic applications49, 50,
42, 45, 51. Other works have been carried out to understand its structural and
optical properties. For instance, studies have characterized its vibrational
properties 38, 45, its optical absorption46, and its refractive index, by
means of optical contrast44, as well as its third-order nonlinear optical
properties, characterized by Z-scan50, 52 and third harmonic generation
(THG)53 techniques. In particular, Tripathi and co-workers53 have reported on
the polarization and thickness dependence of THG. However, few-layer
thicknesses were not reached and, thus, they did not notice any polarization-
dependence changes with thickness.

It is also known that franckeite exhibits spontaneous symmetry break due to
the incommensurability of its lattice network43, 42, 37, 47, 48. This process
introduces a strain-induced structural deformation, which is known to develop
symmetry-break periodic ripples in the $a$ axis orientation. The rippling
arises from the competition between inhomogeneous van der Waals interlayer
interactions and inhomogeneous elastic deformations45, 40. As the crystal
relaxes mechanically to minimize the total energy, it develops both strong in-
plane strains and a small out-of-plane rippling, modulated along the armchair
direction of the H layer 45, 40. These ripples, as well as the $a-b$ plane of
franckeite have been deeply investigated by means of electron microscope
techniques41, 45, 43, 53, 42, 37, 47, 48, 41 and, more recently, by
polarization-dependent Raman spectroscopy45, 53. However, these techniques
commonly require long operational times. Therefore, it is still desired to
demonstrate a tool that is non-destructive, and that is able to quickly assign
the ripples´ axis. Thus, in this work, we demonstrate a fast and simple way to
assign the franckeite ripples’ axis. It consists of determining the
polarization-dependent THG in the natural heterostructure. Thus, it is
similar, requiring an alike experimental setup, to that widely used for
determining the crystal axes of 2H-phase transition metal dichalcogenides (in
the latter case, however, second-harmonic generation (SHG) needs to be
used)54, 55. The method also determines the anisotropic response of
franckeite’s third-order nonlinear susceptibility. We show that not only the
rippling direction can be identified, but also that the polarization-dependent
THG response of franckeite shows noticeable dependence on the sample
thickness, which is well explained by our theoretical nonlinear susceptibility
model. Most importantly, we determine the elements of third-order
susceptibility tensor ($\chi^{(3)}$) and find that they undergo noticeable
changes for very thin flakes. We envisage that this experimental procedure and
results will contribute to the understanding of naturally occurring vdWHs and
could potentialize the use of these low-cost LMs in future photonic and
optoelectronic applications.

##  2 Methods

###  2.1 Sample Preparation

A bulk-franckeite sample from San Jose mine, Oruro - Bolivia, was mechanically
cleaved by scotch tape and then exfoliated again on a polydimethylsiloxane
(PDMS) stamp placed on a glass slide for inspection under an optical
microscope56. Optical contrast was optimized to identify few-layer franckeite
flakes prior to transfer44. The selected samples were then stamped on a 500
$\mu$m glass thick substrate at 40°C using xyz micro-manipulators under an
optical microscope, before increasing the temperature to 60°C, so the flakes
detach from the PDMS and adhere preferably to the glass surface.

###  2.2 Sample Characterization

After the flake transfer, Atomic Force Microscopy (AFM - Bruker Dimension Icon
system) was performed in tapping mode to characterize the sample height.
Franckeite flakes with thicknesses ranging from 8 nm to 130 nm thick were
transferred and analyzed. These thicknesses were identified by AFM
measurements in two staircase-like franckeite samples that were investigated
in this work. For the Raman measurements, a confocal Raman microscope (WITec
Alpha 300R system) was used in the backscattering configuration at 532 nm
excitation energy with a 100X objective which has a numerical aperture NA =
0.9. Each Raman spectrum was collected with three accumulations of 30 seconds
with a grating of 1800 grooves/mm. The laser power was kept constant at
$\sim$1 mW during all measurements and the laser was previously aligned on a
prime silicon wafer to use the Si peak at $\sim$521 cm-1 as reference57.
Polarized Raman measurements were also performed on selected flakes. For this,
an optical linear polarizer (analyzer) was added to the reflected (collection)
path of the microscope before the signal reaches the spectrometer58. We fixed
the laser’s incident polarization as parallel to the analyzer axis. To collect
the polarized Raman spectrum as a function of the angle, we stepwise rotated
the sample stage using a goniometer, acquiring one spectrum per step.

The THG experimental setup is illustrated in Fig. 1. It is composed of a mode-
locked Erbium-doped fiber laser (at a wavelength of 1560 nm, yielding 150-fs
pulses at a repetition rate of 89 MHz) that is focused on the sample at normal
incidence through a 20X objective lens (NA = 0.5). The sample is placed on an
$xyz$ translation stage that is controlled by piezo-scanners. The average
laser power used in the experiments and measured before the objective is 10
mW, except in the THG power dependence measurements. The $1/e^{2}$ beam
diameter was measured to be 3.6$\mu$m59. The THG signal at 520 nm is collected
in transmission mode, passing through another objective lens, which collimates
the beam, and directs it to a spectrometer (Andor Kymera coupled to a silicon
CCD iDus 416 camera)54. A variable optical attenuator, at the laser output,
controls the laser power while a quarter-wave plate (QWP) followed by a linear
polarizer in a motorized rotation stage are used to ensure a practically
constant pump laser intensity (with power variation below 2%) independently of
the polarization angle. The franckeite flakes, as well as the THG spot, were
imaged in reflection by a microscope setup, consisting of a white LED, a CCD
camera, and some simple optics (Fig. 1) so that the laser position in the
sample could be identified (Fig. 2b-d).

Figure 1: Experimental setup for THG characterization of the samples in
transmission mode. QWP: quarter-wave plate; DM: dichroic mirror; Obj.:
microscope objective. The red box indicates that the analyzer can be used or
removed from the setup, depending on the desired experiment.

For polarization-resolved measurements, the motorized stage containing the
polarizer was rotated in $5^{\circ}$ steps from $0^{\circ}$ to $360^{\circ}$.
Two types of polarization-dependent THG measurements were performed: i) an
analyzer placed after the second objective was inserted and remained at a
(chosen) fixed angle, while the polarizer was rotated; ii) no analyzer was
used while the polarizer was angle scanned.

To obtain the THG signal as a function of pump power, we selected a
representative franckeite region, with homogeneous thickness. Before
collecting the spectra, we maximized the THG signal by changing the sample
focus on the flake surface. Then, the pump power was tuned from 1 mW to 11 mW
and each spectrum was taken with two accumulations of 1 second acquisition
time. In this case, the polarizer and analyzer were kept parallel to each
other, at a fixed angle. The results presented for each pump power were
averaged from five spectra.

To acquire the THG signal as a function of thickness, we carried out the
following steps, for each region: i) the THG signal was maximized by adjusting
the sample-objective distance; ii) full polarization-dependent THG
measurements were conducted; iii) we selected the angle with maximum THG
intensity and averaged the THG signal over five spectra. All THG points
presented in the manuscript, even those that were taken from different
samples, were collected using these steps. Finally, note that the pump photon
energy is beyond franckeite bandgap energy, leading to linear absorption, in
addition to THG. Despite the fact that pump absorption was found to somewhat
reduce transmittance in the characterized flakes ($\sim$90% and $\sim$65%
transmittance for the thicker and thinner flakes measured, respectively), no
detectable damage was observed throughout the experiments.

###  2.3 THG power as function of thickness: theoretical approach

To theoretically investigate the THG power as function of the thickness of the
sample, we use the transfer matrix method for harmonic generation60. Within
this approach, the sample is modeled as a nonlinear slab of thickness $d$
between two linear half-spaces. The THG fields that radiate to air, thus
counter-propagating with the pump (region 1), and to glass, co-propagating
with the pump (region 3), are then easily obtained as:

$\displaystyle\begin{pmatrix}E^{+}_{3}\\\
E^{-}_{1}\end{pmatrix}=\frac{1}{1+\phi^{2}_{2}r_{12}r_{23}}\begin{pmatrix}\phi_{2}t_{23}[S^{+}+r_{12}S^{-}]\\\
t_{21}[\phi^{2}_{2}r_{23}S^{+}-S^{-}]\end{pmatrix},$ (1)

where the refractive index $N_{i}$ of slab $i$, the terms
$r_{ij}=\frac{N_{i}-N_{j}}{N_{i}+N_{j}}$, $t_{ij}=\frac{2N_{i}}{N_{i}+N_{j}}$
($j=1,2,3$) and $\phi_{2}=e^{iN_{2}(3\omega/c)d_{2}}$ are evaluated at
frequency $3\omega$. The source terms

$\displaystyle\begin{pmatrix}S^{+}\\\
S^{-}\end{pmatrix}=\frac{1}{t_{2s}}\begin{pmatrix}\phi^{-1}_{2}\phi_{s}-1&[\phi^{-1}_{2}\phi^{-1}_{s}-1]r_{s}\\\
[\phi_{2}\phi_{s}-1]r_{s}&\phi_{2}\phi^{-1}_{s}-1\end{pmatrix}\vec{E_{2s}}$
(2)

are related to the pump field in the nonlinear slab (layer 2) by:

$\displaystyle\vec{E_{2s}}=\frac{1}{t^{\omega}_{21}}\begin{pmatrix}1+r^{\omega}_{21}r^{\omega}\\\
r^{\omega}_{21}+r^{\omega}\end{pmatrix}\left[\frac{4\pi}{\varepsilon_{s}}\right]\vec{P}_{NL},$
(3)

where the superscript ${\omega}$ means that the quantity is evaluated at the
pump frequency. $r^{\omega}$ is the overall reflection coefficient and can be
calculated from the well-known Airy’s formula for three media
$r^{\omega}_{123}$ or by the recursion formula for multiple reflections.
$\vec{P}_{NL}=\chi^{(3)]}.\vec{E^{\omega}}.\vec{E^{\omega}}.\vec{E^{\omega}}$
is the nonlinear polarization. Note that when the tensorial product in the
nonlinear polarization is expanded, four types of waves appear with the
$3\omega$ frequency: $e^{\pm i3N_{2}k_{0}z}$ and $e^{\pm iN_{2}k_{0}z}$. Each
one of these waves can be related to components of the non-linear polarization
source term with different $\vec{k}$ vectors and effective dielectric
constants $\varepsilon_{s}=\varepsilon_{2}(3\omega)-\varepsilon_{2}(\omega)$
for the waves $e^{\pm i3N_{2}k_{0}z}$ and
$\varepsilon_{s}=\varepsilon_{2}(3\omega)-\varepsilon_{2}(\omega)/9$ for the
terms $e^{\pm iN_{2}k_{0}z}$60. We treat independently each one of these
contributions through the subscript $s$ in equation 3; to obtain the total
electric field in region 3 all contributions are summed up.

##  3 Results and Discussion

Figure 2: a) Crystal structure of franckeite where the two different stacked
layers, the SnS2-type (H) and the PbS-type (Q), can be seen.b) Transmission
optical microscope image of an exfoliated staircase-like franckeite flake.
Flake orientation ($a-b$) is indicated, in addition to the area where the few-
layer franckeite was measured by AFM (white rectangle). c) AFM topography
image collected from the area marked in (b). d) Lower resolution/contrast
reflection image showing the THG emission (bright green spot) obtained with
laser incidence onto a representative region. e) Measured THG spectrum with a
peak wavelength at 520 nm from a 34 nm-thick franckeite flake. f) Double log-
scale plot of the measured THG intensity as a function of the laser power with
a power law fit ($\alpha=3.10\pm 0.02$) to the data obtained from the region
shown in (b).The error bars shown were calculated as the standard deviation of
5 measurements taken at the same laser power.

Fig. 2b) shows a transmission optical microscope image of one of the samples
studied in this work, with the franckeite orientation ($a$ and $b$ axes)
indicated on it. Regions of various colors are visible, which are caused by
thin-film optical interference, indicating different flake thicknesses 44, as
confirmed by AFM measurements (Fig. 2c). Fig. 2d) shows the THG emission
(green spot indicated by an arrow) captured by the CCD camera of the
experimental setup, from a 34 nm-thick franckeite region during the
measurements. A representative THG spectrum is plotted in Fig. 2e),
demonstrating the peak wavelength at 520 nm which is exactly one-third of the
excitation wavelength. Further confirmation of THG can be seen in Fig. 2f),
which shows the pump power dependence of the THG intensity, with a cubic power
law fit ($\alpha=3.10\pm 0.02$) to the experimental data.

Figure 3: a) Representative angle-resolved polarized Raman spectra for a 34
nm-thick franckeite flake collected in the parallel polarization configuration
under 532 nm laser excitation. The dashed black line indicates the most
sensitive Raman mode to the laser polarization. b-c) THG intensity dependence
on the incident linear polarization angle for a representative franckeite
flake with 34 nm thickness, collected with (b) and without (c) the analyzer.
0∘ corresponds to the a-axis (ripples direction). The red
($I{{}_{a}}^{3\omega}$) and black ($I{{}_{b}}^{3\omega}$) symbols in (b)
correspond to the measurements carried out with the analyzer aligned parallel
and perpendicular to the rippling axis, respectively. Part (c) also shows the
trace obtained by the sum of the curves in (b) multiplied by 0.25 (i.e.,
$(I{{}_{a}}^{3\omega}+I{{}_{b}}^{3\omega})\times 0.25$), which accounts for
the 25% attenuation of the analyzer, as stated in its specifications.

As mentioned before, franckeite is a natural vdWH with lattice
incommensurability between the two consecutively stacked PbS-like and
SnS2-like layers43, 42, 37, 47, 48. This arrangement introduces a strain-
induced structural deformation, creating the known periodic ripples in the $a$
axis orientation. Considering that this feature can be identified through
polarized Raman spectroscopy45, 53, we show in Fig. 3a a set of Raman spectra
at various polarization angles collected from a 34 nm-thick flake. A typical
Raman spectrum of franckeite exhibits several distinct Raman modes within the
75–400 cm-1 window, which varies differently with polarization45, 53. These
previous works demonstrated that the peak at $\sim$265 cm-1 (highlighted by
the vertical dashed line in Fig. 3a) is the most polarization-sensitive mode
and can be used to assign the rippling orientation of franckeite. In
particular, the maximum (minimum) mode intensity is observed when the laser
polarization is parallel (orthogonal) to the rippling axis ($a$ axis). Hence,
the 0∘ (90∘) polarization angle is defined as parallel (perpendicular) to the
$a$ axis (ripples). Note that we have used this approach to identify the
franckeite orientation, as shown in Fig. 2b, prior to further experiments. It
is important to stress that polarized Raman measurements taken on different
franckeite thicknesses (Fig. 2c) of the same flake exhibit a similar trend for
the Raman modes as a function of the polarization angle, suggesting that if a
flake is composed of regions with different layer numbers (as the one in Fig.
2b-d), the orientation of the ripples along the entire flake is maintained.
Therefore, by finding the rippling orientation in one region of the flake, one
can assume the same alignment for all regions.

Since the linear optical response of franckeite shows strong in-plane
structural anisotropy with the incident light’s polarization45, it can be
expected that the nonlinear response would show similar behavior. We have then
investigated the effect of this in-plane anisotropy of franckeite through
polarization-dependent THG (see Methods for more details). We start the
polarization-dependent THG characterization by fixing the analyzer (Fig. 1)
either at the 0∘ or the 90∘ angle, as previously identified by Raman
spectroscopy (Fig. 3a). Fig. 3b displays the angular dependence of the THG
intensity on the incident laser polarization angle relative for a 34 nm-thick
franckeite flake. With the analyzer parallel to the $a$ axis, a symmetric two-
lobe pattern is observed (black symbols), with the maximum THG occurring at 0∘
(parallel to the rippling direction). When the analyzer is orthogonal to the
$a$ axis, a second two-lobed curve (red symbols) is observed. As discussed
later, this curve exhibits asymmetry at low THG intensity angles (see Fig.
4a). The emission pattern shows the maximum THG occurring at 0∘ (parallel to
the rippling direction) and the secondary maximum occurring close to 90∘.

The anisotropy ratio of the THG signal
I${{}_{a}}^{3\omega}$($\theta=0^{\circ}$)/I${{}_{b}}^{3\omega}$($\theta=90^{\circ}$)
is around 1.25. Note that the overall maximum THG is conveniently obtained
along the rippling axis (with the analyzer along the $a$ axis), which means it
can potentially be used for franckeite orientation determination. However, in
the set of data shown so far, this axis orientation needs to be known a
priori, so that the analyzer was positioned. Nevertheless, we show in Fig. 3c
(dark blue) that by performing a single polarization-dependent THG
characterization without the analyzer, one gets directly the same lobe
structure and, consequently, the rippling orientation. Indeed, the obtained
angular profile matches that obtained by the sum of the THG profiles measured
with the analyzer along the $a$ and $b$ axes, also shown in Fig. 3c (cyan).
Additionally, the anisotropy ratio of the THG signal
I${{}_{a}}^{3\omega}$($\theta=0^{\circ}$)/I${{}_{b}}^{3\omega}$($\theta=90^{\circ}$)
$\sim$1.26 is recovered as well. It is worth mentioning that we have
previously shown that the THG signal of franckeite can be $\sim$3 orders of
magnitude higher than graphene 61; hence, one can take $<$1 min for
identifying the franckeite axes without further characterization. Therefore,
we demonstrate that this approach can be used to quickly identify the
franckeite orientation. Additionally, it may potentially be applied to other
anisotropic LMs and heterostructures53, 31, 62, 33, 63, 64.

We now focus on the THG dependence on the franckeite thickness. Figs. 4a-c
show representative THG profiles as a function of polarization, for three
flake thicknesses: 10 nm, 34 nm, and 124 nm, respectively. In the
polarization-dependent THG measurements, we observe differences in the lobe
shape for angles close to the THG minima (Figs. 4d-f). Note that, when the
analyzer is parallel to the $b$ axis, the thinner the flake, the more
pronounced is the shoulder-like feature close to the base of the lobes (lower
THG intensities); for thicker samples, a more symmetric lobed structure is
obtained. These changes with flake thickness are reproducible and beyond
experimental errors, and seem to indicate a change in the $\chi^{(3)}$ tensor,
as franckeite is thinned down. To quantify the $\chi^{(3)}$ variation, we use
a nonlinear susceptibility model53, and perform the fitting of these three
representative experimental THG traces (solid lines in Fig. 4a-f). Even though
this procedure is common practice in the literature, we stress that accurate
determination of the nonlinear tensor components is tricky and cannot be done
univocally solely from THG measurements. Thus, the purpose here is just to
provide fitting to the experimental curves, as well as to indicate which
tensor components might dominantly change with thickness.

Figure 4: Angular dependence of the THG intensity on the incident linear
polarization angle for franckeite flakes with thicknesses of (a) 10 nm, (b) 34
nm, and (c) 124 nm. d-f) are zoomed graphs close to the lobe bases of the
respective samples. Red and black symbols correspond to the measurements
carried out with the analyzer aligned parallel and perpendicular to the
rippling axis, respectively. The corresponding theoretical fittings are
plotted as solid curves.

In the triclinic crystal system all elements of the third-order susceptibility
are independent and non-zero. Despite this, the tensor can be contracted for a
simplified computation of the THG fields. 65 The element $\chi^{(3)}_{ijkl}$
of the fourth-rank tensor is expressed in the contracted notation as
$\chi^{(3)}_{mn}$. Here, the index $m=1,2,3$ denotes the spatial component
along the directions $\hat{a},~{}\hat{b}~{}\text{or}~{}\hat{c}$ of the THG
field. The contracted index $n$ represents permutations of the three spatial
components of the incident wave
($\hat{e}_{x},\hat{e}_{y}~{}\text{and}~{}\hat{e}_{z}$) and varies from 1 to
10. 65 With these considerations, the THG field is

$\left(\begin{matrix}E_{a}^{3\omega}\\\ E_{b}^{3\omega}\\\
E_{c}^{3\omega}\end{matrix}\right)\propto\begin{pmatrix}\chi_{11}&\chi_{12}&\chi_{13}&\chi_{14}&\chi_{15}&\chi_{16}&\chi_{17}&\chi_{18}&\chi_{19}&\chi_{110}\\\
\chi_{21}&\chi_{22}&\chi_{23}&\chi_{24}&\chi_{25}&\chi_{26}&\chi_{27}&\chi_{28}&\chi_{29}&\chi_{210}\\\
\chi_{31}&\chi_{32}&\chi_{33}&\chi_{34}&\chi_{35}&\chi_{36}&\chi_{37}&\chi_{38}&\chi_{39}&\chi_{310}\end{pmatrix}\cdot\vec{L}.$
(4)

Where the rows of the column vector $\vec{L}$ are the combinations of the
three components of the pump field, as shown below. For the analysis, we take
into account that the pump laser has a fundamental frequency of 1560 nm and
that the beam is linearly polarized. Moreover, as the experiment is carried
out at normal incidence, the pump laser electric field at the focal point of
the objective only presents components in the $\hat{e}_{x}-\hat{e}_{y}$ plane.
Consequently, we can disregard all terms where $\hat{e}_{z}$ appears. That is,
to focus on the rows of the column vector $\vec{L}$ with index
$n=1:\hat{e}_{x}\hat{e}_{x}\hat{e}_{x},2:\hat{e}_{y}\hat{e}_{y}\hat{e}_{y},8:(\hat{e}_{x}\hat{e}_{y}\hat{e}_{y}+\hat{e}_{y}\hat{e}_{x}\hat{e}_{y}+\hat{e}_{y}\hat{e}_{y}\hat{e}_{x})~{}\text{and}~{}9:(\hat{e}_{x}\hat{e}_{x}\hat{e}_{y}+\hat{e}_{x}\hat{e}_{y}\hat{e}_{x}+\hat{e}_{y}\hat{e}_{x}\hat{e}_{x})$.
If the incident electric field is defined as
$\vec{E}=|E|(\cos\theta\hat{e}_{x}+\sin\theta\hat{e}_{y})$ after substitution
in equation 4 the intensities of the THG are obtained:

$\displaystyle
I^{3\omega}_{a}\propto|E^{3\omega}_{a}|^{2}\propto\left(\chi_{11}\cos^{3}\theta+\chi_{12}\sin^{3}\theta+3\chi_{18}\cos\theta\sin^{2}\theta+3\chi_{19}\cos^{2}\theta\sin\theta\right)^{2},$
(5) $\displaystyle
I^{3\omega}_{b}\propto|E^{3\omega}_{b}|^{2}\propto\left(\chi_{21}\cos^{3}\theta+\chi_{22}\sin^{3}\theta+3\chi_{28}\cos\theta\sin^{2}\theta+3\chi_{29}\cos^{2}\theta\sin\theta\right)^{2}.$
(6)

The contracted tensor components are then used as fitting parameters to fit
the experimental data in Fig. 4. Note that, particularly in the zoomed regions
(Figs. 4d,e,f), the analytical curves fit the experimental data
satisfactorily. The $\chi_{mn}$ components, normalized by $\chi_{11}$,
obtained after the least square fitting process is carried out, are presented
in Table 1.

It is clearly observed that the elements of the $\chi^{(3)}$ tensor change
with the sample thickness, explaining the observed differences in the lobes
shapes. In particular, note that, as mentioned, the shoulder-like structure
appears in the curve associated with $I_{b}^{3\omega}$ (equation 6), for which
the tensor index $m=2$. Furthermore, the shoulder is obtained when the pump is
practically aligned with the $x$ axis. Thus, we can, to a first approximation,
neglect indices $n$ corresponding to two, or more, $y$ pump components,
relative to indices corresponding to two, or more, $x$ pump components. This
means that the most relevant $n$ indices are $n=1$ and $n=9$. Indeed,
inspection of Table 1 shows that $mn$ components 21 and 29, in fact, present
appreciable values (up to $\sim 30\%$ and $\sim 10\%$ of the value of
$\chi_{11}$, respectively) and a monotonic variation as the thickness
increases ($\chi_{21}$ drops and $\chi_{29}$ increases). Table 1 also includes
a comparison with the results reported by Tripathi et al.53. Note that, in
this case, the authors have only measured flakes down to a thickness of 25 nm
and noticed no significant changes in the tensor components. The values
provided are, therefore, thickness-averaged values. Indeed, they exhibit a
fair agreement with the values we found, especially for the case of our two
thicker flakes, for which case $\chi_{18}$, $\chi_{19}$, $\chi_{21}$, and
$\chi_{28}$ show lower values than for the 10-nm-thick sample. Note that
changes in the amplitude (and phase) of the $\chi^{(3)}$ tensor components
ultimately trace back to changes in the spectral position of resonances and/or
to the oscillator strength between different electronic states in franckeite.
Therefore, the thickness-dependent changes reported here are like to be a
consequence of changes in the electronic structure and, consequently can also
impact the performance of photonic and optoelectronic devices (e.g.
photodetectors and optical modulators 42, 66, 67) based on franckeite.
However, quantification of such an impact is beyond the scope of this paper.

Sample | $\chi_{11}$ | $\chi_{12}$ | $\chi_{18}$ | $\chi_{19}$ | $\chi_{21}$ | $\chi_{22}$ | $\chi_{28}$ | $\chi_{29}$
---|---|---|---|---|---|---|---|---
10 nm | 1 | 0 | 0.0719 | 0.0424 | 0.3454 | 0.8692 | 0.0748 | 0.0370
34 nm | 1 | 0.0052 | 0.0476 | 0.0018 | 0.2819 | 0.7796 | 0.2139 | 0.0638
124 nm | 1 | 0 | 0.0633 | 0 | 0 | 0.9491 | 0 | 0.1084
Tripathi et al.53 | 1 | 0.012 | 0.043 | 0.015 | 0.017 | 0.815 | 0.009 | 0.037
Table 1: Values of $\chi_{mn}$ normalized by $\chi_{11}$ obtained by the
least-square fitting process, and comparison with Ref. 53

We now analyze the influence of the material thickness on the maximum THG
intensity (obtained, as discussed, for a pump polarization along the $a$
axis). Fig. 5 plots the dependence of the THG intensity on the franckeite
thickness for the incident linear polarization along the rippling direction.
In order to maintain consistency between all the points collected, even those
that were in different flakes, the measurements were collected using the
careful calibration criteria described in Methods. Color symbols express
experimental data with error bars, and analytical fittings (obtained as
discussed in Methods) are plotted with solid lines. Compared to previously
reported measurements on thicker samples, where a wide Gaussian-like shape
peak is seen around $\sim 100$ nm,53 in our samples, the THG dependence on
thickness shows a more complex structure. Surprisingly, the peak (maximum THG
value) appears for much thinner samples, being at around 25 nm, with a
shoulder-like feature located at about 70 nm.

Figure 5: THG intensity plotted versus franckeite thickness (symbols).
Thicknesses were extracted from AFM measurements in two staircase-like
franckeite samples (black circles and gray star symbols represent the
different samples measured). The solid lines correspond to the best fit
analytical curves obtained with various models: (purple) fully coherent model
with infinite number of reflections; (blue) coherent model with no reflections
(the intensity was multiplied by 7.97 to match the maximum experimental
intensity); (green) coherent model with $l=3$ reflections; and (red) partially
coherent model. All models are explained in the text. The error bars shown
were calculated as the standard deviation of 5 measurements taken at the same
flake thickness.

To understand the observed trends, we model the sample as a layered structure
and use the transfer matrix method for harmonic generation outlined in the
Methods section. At normal incidence and with the polarization of the pump
parallel to the rippling axis, we can, with no loss of generality, treat the
THG fields as a one-dimensional problem. First, it is clear from equation 1
that the source field determines the THG fields in the air and the silica
glass regions. With this in mind, it seems reasonable to assume that multiple
pump reflections within the flake cannot be neglected. We initially consider
an infinite number of pump reflections within the nonlinear slab; that is,
using for $r^{\omega}$ in equation 3 the Airy reflection formula for three
media
($r^{\omega}=\frac{r^{\omega}_{12}+r^{\omega}_{23}e^{i2N_{2}(\omega)(\omega/c)z}}{1+r^{\omega}_{12}r^{\omega}_{23}e^{i2N_{2}(\omega)(\omega/c)z}}$).
Second, it is equally clear that the process of calculating the THG fields can
be reversed and the optical constants of the material can be determined from
the measurement of the THG fields. Thus, reversing the process we obtain the
refractive index values $N_{2}(\omega)=4.1159+i0.5$ and
$N_{2}(3\omega)=3.5+i1.5$, which satisfactorily reproduces the THG thickness
dependence for thicknesses smaller than 30 nm (purple line in Fig. 5). Indeed,
these refractive index values are very similar to those previously reported53,
44. The peak position of THG versus thickness curve may be attributed to phase
mismatching between the pump and the third-harmonic wave. Indeed, the
coherence length ($L_{c}$) of nonlinear process is given by68
$L_{c}=\lambda/6(N_{2}(\omega)+N_{2}(3\omega))=33.920$ nm for backward
propagating waves and $L_{c}=\lambda/6(N_{2}(\omega)-N_{2}(3\omega))=419.44$
nm for the forward propagating waves. As multiple reflections take place,
forward and backward-generated waves become co-propagating, meaning that the
shorter $L_{c}$ dominates. Indeed, the length for backward THG closely matches
the THG peak position.

In an effort to improve the fit to the experimental data in the thicker sample
range, we evaluate the effect of limiting the number of reflections inside the
nonlinear layer that effectively contributes to THG. This is done by the
recursion formula
$r^{\omega}=r^{\omega}_{12}+t^{\omega}_{21}r^{\omega}_{23}t^{\omega}_{12}\phi^{\omega}_{2}\sum_{j=0}^{l}q^{j}$,
where $q=r^{\omega}_{21}r^{\omega}_{23}\phi^{\omega}_{2}$ and
$\phi^{\omega}_{2}=e^{-i2N_{2}(\omega/c)d}$. For $l=0$, i.e., completely
disregarding multiple reflections of the pump field, it is impossible to fit
the experimental data (blue curve in Fig. 5). This model gives a maximum THG
intensity at $\sim$85 nm. It is important to mention that the calculated
intensities for this case were scaled (x 7.97) in an effort to better fit the
data. For $l=3$ (green curve in Fig. 5), the peak THG value increases and
moves towards smaller thicknesses, rapidly approaching the reaching the value
for $l=\infty$. For $l>5$ the curve for infinite reflections is recovered.
However, the shoulder feature that is observed in the experiments can never be
reproduced.

Note that the entire nonlinear optical treatment provided here assumes the
pump to be monochromatic, consisting of one single frequency $\omega$. In
practice, being a femtosecond pulse source, the pump is composed of a band of
frequencies around $\omega$, which significantly scales the complexity of the
problem. For example, sum-frequency generation involving different pump
frequencies occurs, which can further reduce the coherence of the process. To
emulate these features without resorting to a more complex model, we
empirically add dephasing factors that affect the pump, by multiplying its
Fresnel reflection (transmission) coefficients69, 70. The modified
coefficients are now rewritten as $r_{ij}=\alpha r_{ij}^{0}$ and $t_{ij}=\beta
t_{ij}^{0}$, where the superscript $0$ indicates the Fresnel coefficient for
the monochromatic case. We find that with
$r^{\omega}_{12}=0.5\left(r^{\omega}_{12}\right)^{0}$,
$t^{\omega}_{21}=0.5\left(t^{\omega}_{21}\right)^{0}$ and
$r^{\omega}_{23}=0.9\left(r^{\omega}_{23}\right)^{0}$ the calculated
intensities satisfactorily fit the experimental data for thin and thick
samples, as can be seen for the red curve in Fig. 5. Finally, note that if
dephasing increases (with the factors tending to zero), the curve for $l=0$
(i.e., single pass THG) is recovered. We notice that the single-pass model
was, indeed, used to fit the data reported by Tripathi et al.53, with a good
experimental agreement. We speculate that this is so because those authors
used a pump consisting of significantly shorter pulses (90 fs), for which case
dephasing increases, approaching the single-pass picture.

###  3.1 Conclusion

In summary, we report on an experimental and theoretical investigation of the
nonlinear optical properties of layered franckeite. The third-order nonlinear
optical response (third-harmonic generation) was found to be highly
anisotropic, varying in a nontrivial way with the incident polarization. We
demonstrated that THG measurements are able to quickly identify the franckeite
rippling axis with high accuracy and with no need for additional measurements.
Additionally, we found the THG to be strongly dependent on the material
thickness, not only in terms of absolute achieved harmonic intensity, but also
in terms of pump polarization dependence. We attribute the latter to changes
in the third-order nonlinear susceptibility tensor. A maximum THG intensity
was observed at a franckeite thickness of about 25 nm, which is mainly
attributed to multiple pump reflections within the franckeite flake and to
phase matching. We believe that our findings contribute to a more
comprehensive understanding of anisotropic layered materials and can
potentialize their use in future photonic and optoelectronic applications.

All authors acknowledge the Coordenação de Aperfeiçoamento de Pessoal de Nível
Superior (CAPES), Conselho Nacional de Desenvolvimento de Científico e
Tecnológico (CNPq) and the São Paulo State Foundation (FAPESP; grant numbers:
2020/04374-6, 2018/25339-4, and 2015/11779-4), and the Brazilian Nanocarbon
Institute of Science and Technology (INCT/Nanocarbono). A.R.C., A.S.M.V.O.,
D.A.B., and C.J.S.d.M. acknowledge the support from Fundo Mackenzie de
Pesquisa e Inovação (MackPesquisa grant number 221017).

## References

  * Choi et al. 2017 Choi, W.; Choudhary, N.; Han, G. H.; Park, J.; Akinwande, D.; Lee, Y. H. Recent development of two-dimensional transition metal dichalcogenides and their applications. _Materials Today_ 2017, _20_ , 116–130
  * Backes et al. 2020 Backes, C.; Abdelkader, A. M.; Alonso, C.; Andrieux-Ledier, A.; Arenal, R.; Azpeitia, J.; Balakrishnan, N.; Banszerus, L.; Barjon, J.; Bartali, R.; others Production and processing of graphene and related materials. _2D Materials_ 2020, _7_ , 022001
  * Barcelos et al. 2020 Barcelos, I. D.; Bechtel, H. A.; de Matos, C. J. S.; Bahamon, D. A.; Kaestner, B.; Maia, F. C. B.; Freitas, R. O. Probing Polaritons in 2D Materials with Synchrotron Infrared Nanospectroscopy. _Advanced Optical Materials_ 2020, _8_ , 1901091
  * de Matos et al. 2023 de Matos, C. J. S.; Rosa, H. G.; Zapata, J. D.; Steinberg, D.; Maldonado, M.; de Souza, E. A. T.; de Paula, A. M.; Malard, L. M.; Gomes, A. S. L. Nonlinear optics in 2D materials: focus on the contributions from Latin America. _J. Opt. Soc. Am. B_ 2023, _40_ , C111–C132
  * Akinwande et al. 2014 Akinwande, D.; Petrone, N.; Hone, J. Two-dimensional flexible nanoelectronics. _Nature communications_ 2014, _5_ , 1–12
  * Liao et al. 2019 Liao, W.; Huang, Y.; Wang, H.; Zhang, H. Van der Waals heterostructures for optoelectronics: progress and prospects. _Applied Materials Today_ 2019, _16_ , 435–455
  * Riccardo et al. 2020 Riccardo, F.; Yue, N.; Gant, P.; Muñoz, M.; Castellanos-Gomez, A. Naturally occurring van der Waals materials. _npj 2D Materials and Applications_ 2020, _4_
  * Santos et al. 2019 Santos, J. C.; Barboza, A. P. M.; Matos, M. J.; Barcelos, I. D.; Fernandes, T. F.; Soares, E. A.; Moreira, R. L.; Neves, B. R. Exfoliation and characterization of a two-dimensional serpentine-based material. _Nanotechnology_ 2019, _30_ , 445705
  * de Oliveira et al. 2022 de Oliveira, R.; Guallichico, L. A.; Policarpo, E.; Cadore, A. R.; Freitas, R. O.; da Silva, F. M.; de C. Teixeira, V.; Paniago, R. M.; Chacham, H.; Matos, M. J.; Malachias, A.; Krambrock, K.; Barcelos, I. D. High throughput investigation of an emergent and naturally abundant 2D material: Clinochlore. _Applied Surface Science_ 2022, _599_ , 153959
  * Cadore et al. 2022 Cadore, A. R.; de Oliveira, R.; Lobato, R. L. M.; Teixeira, V.; Nagaoka, D. A.; Alvarenga, V. T.; Soares, J. R.; Watanabe, K.; Taniguchi, T.; Paniago, R.; others Exploring the structural and optoelectronic properties of natural insulating phlogopite in van der Waals heterostructures. _2D Materials_ 2022, _9_ , 035007
  * de Oliveira et al. 2023 de Oliveira, R.; Cadore, A. R.; Freitas, R. O.; Barcelos, I. D. Review on infrared nanospectroscopy of natural 2D phyllosilicates. _J. Opt. Soc. Am. A_ 2023, _40_ , C157–C168
  * Longuinhos et al. 2023 Longuinhos, R.; Cadore, A. R.; Bechtel, H. A.; J S De Matos, C.; Freitas, R. O.; Ribeiro-Soares, J.; Barcelos, I. D. Raman and Far-Infrared Synchrotron Nanospectroscopy of Layered Crystalline Talc: Vibrational Properties, Interlayer Coupling, and Symmetry Crossover. _The Journal of Physical Chemistry C_ 2023, _127_ , 5876–5885
  * Barcelos et al. 2023 Barcelos, I. D.; de Oliveira, R.; Schleder, G. R.; Matos, M. J. S.; Longuinhos, R.; Ribeiro-Soares, J.; Barboza, A. P. M.; Prado, M. C.; Pinto, E. S.; Gobato, Y. G.; Chacham, H.; Neves, B. R. A.; Cadore, A. R. Phyllosilicates as earth-abundant layered materials for electronics and optoelectronics: Prospects and challenges in their ultrathin limit. _Journal of Applied Physics_ 2023, _134_ , 090902
  * Mannix et al. 2017 Mannix, A. J.; Kiraly, B.; Hersam, M. C.; Guisinger, N. P. Synthesis and chemistry of elemental 2D materials. _Nature Reviews Chemistry_ 2017, _1_ , 1–14
  * Shivayogimath et al. 2019 Shivayogimath, A.; Thomsen, J. D.; Mackenzie, D.; Geisler, M.; Stan, R.-M.; Holt, A. J.; Bianchi, M.; Crovetto, A.; Whelan, P. R.; Carvalho, A.; others A universal approach for the synthesis of two-dimensional binary compounds. _Nature communications_ 2019, _10_ , 1–7
  * Hamza Safeer et al. 2022 Hamza Safeer, S.; Ore, A. S. M. V.; Cadore, A. R.; Gordo, V. O.; Vianna, P. G.; Carvalho, I. C. S.; Carozo, V.; de Matos, C. J. S. CVD growth and optical characterization of homo and heterobilayer TMDs. _Journal of Applied Physics_ 2022, _132_ , 024301
  * Nagaoka et al. 2023 Nagaoka, D. A.; Grasseschi, D.; Cadore, A. R.; Fonsaca, J. E. S.; Jawaid, A. M.; Vaia, R. A.; de Matos, C. J. S. Redox exfoliated NbS2: characterization, stability, and oxidation. _Phys. Chem. Chem. Phys._ 2023, _25_ , 9559–9568
  * Geim and Grigorieva 2013 Geim, A. K.; Grigorieva, I. V. Van der Waals heterostructures. _Nature_ 2013, _499_ , 419–425
  * Liu et al. 2016 Liu, Y.; Weiss, N. O.; Duan, X.; Cheng, H.-C.; Huang, Y.; Duan, X. Van der Waals heterostructures and devices. _Nature Reviews Materials_ 2016, _1_ , 1–17
  * Andrei et al. 2021 Andrei, E. Y.; Efetov, D. K.; Jarillo-Herrero, P.; MacDonald, A. H.; Mak, K. F.; Senthil, T.; Tutuc, E.; Yazdani, A.; Young, A. F. The marvels of moiré materials. _Nature Reviews Materials_ 2021, _6_ , 201–206
  * Novoselov et al. 2016 Novoselov, K. S.; Mishchenko, A.; Carvalho, A.; Neto, A. H. C. 2D materials and van der Waals heterostructures. _Science_ 2016, _353_ , aac9439
  * Mak and Shan 2016 Mak, K. F.; Shan, J. Photonics and optoelectronics of 2D semiconductor transition metal dichalcogenides. _Nature Photonics_ 2016, _10_ , 216–226
  * Sun et al. 2016 Sun, Z.; Martinez, A.; Wang, F. Optical modulators with 2D layered materials. _Nature Photonics_ 2016, _10_ , 227–238
  * Feres et al. 2023 Feres, F. H.; Barcelos, I. D.; Cadore, A. R.; Wehmeier, L.; Nörenberg, T.; Mayer, R. A.; Freitas, R. O.; Eng, L. M.; Kehr, S. C.; Maia, F. C. B. Graphene Nano-Optics in the Terahertz Gap. _Nano Letters_ 2023, _23_ , 3913–3920, PMID: 37126430
  * Nutting et al. 2021 Nutting, D.; Prando, G. A.; Severijnen, M.; Barcelos, I. D.; Guo, S.; Christianen, P. C. M.; Zeitler, U.; Galvão Gobato, Y.; Withers, F. Electrical and optical properties of transition metal dichalcogenides on talc dielectrics. _Nanoscale_ 2021, _13_ , 15853–15858
  * Gadelha et al. 2019 Gadelha, A. C.; Cadore, A. R.; Watanabe, K.; Taniguchi, T.; de Paula, A. M.; Malard, L. M.; Lacerda, R. G.; Campos, L. C. Gate-tunable non-volatile photomemory effect in MoS2 transistors. _2D Materials_ 2019, _6_ , 025036
  * Chakraborty et al. 2022 Chakraborty, S. K.; Kundu, B.; Nayak, B.; Dash, S. P.; Sahoo, P. K. Challenges and opportunities in 2D heterostructures for electronic and optoelectronic devices. _Iscience_ 2022, _25_ , 103942
  * Montblanch et al. 2021 Montblanch, A. R.-P.; Kara, D. M.; Paradisanos, I.; Purser, C. M.; Feuer, M. S.; Alexeev, E. M.; Stefan, L.; Qin, Y.; Blei, M.; Wang, G.; others Confinement of long-lived interlayer excitons in WS2/WSe2 heterostructures. _Communications Physics_ 2021, _4_ , 1–8
  * Román et al. 2020 Román, R. J. P.; Auad, Y.; Grasso, L.; Alvarez, F.; Barcelos, I. D.; Zagonel, L. F. Tunneling-current-induced local excitonic luminescence in p-doped WSe 2 monolayers. _Nanoscale_ 2020, _12_ , 13460–13470
  * Niu et al. 2019 Niu, Y.; Villalva, J.; Frisenda, R.; Sanchez-Santolino, G.; Ruiz-González, L.; Pérez, E. M.; García-Hernández, M.; Burzurí, E.; Castellanos-Gomez, A. Mechanical and liquid phase exfoliation of cylindrite: a natural van der Waals superlattice with intrinsic magnetic interactions. _2D Materials_ 2019, _6_ , 035023
  * Dasgupta et al. 2021 Dasgupta, A.; Gao, J.; Yang, X. Natural van der Waals heterostructure cylindrite with highly anisotropic optical responses. _npj 2D Materials and Applications_ 2021, _5_ , 1–11
  * Martín-Pérez and Burzurí 2021 Martín-Pérez, L.; Burzurí, E. Optimized Liquid-Phase Exfoliation of Magnetic van der Waals Heterostructures: Towards the Single Layer and Deterministic Fabrication of Devices. _Molecules_ 2021, _26_ , 7371
  * Dasgupta et al. 2021 Dasgupta, A.; Yang, X.; Gao, J. Naturally occurring van der Waals heterostructure lengenbachite with strong in-plane structural and optical anisotropy. _npj 2D Materials and Applications_ 2021, _5_ , 1–13
  * Chen et al. 2019 Chen, H.; He, J.; Malliakas, C. D.; Stoumpos, C. C.; Rettie, A. J.; Bao, J.-K.; Chung, D. Y.; Kwok, W.-K.; Wolverton, C.; Kanatzidis, M. G. A Natural 2D Heterostructure [Pb3. 1Sb0. 9S4][Au x Te2–x] with Large Transverse Nonsaturating Negative Magnetoresistance and High Electron Mobility. _Journal of the American Chemical Society_ 2019, _141_ , 7544–7553
  * Gusmão et al. 2018 Gusmão, R.; Sofer, Z.; Luxa, J.; Pumera, M. Layered franckeite and teallite intrinsic heterostructures: shear exfoliation and electrocatalysis. _Journal of Materials Chemistry A_ 2018, _6_ , 16590–16599
  * Shu et al. 2020 Shu, Z.; Peng, Q.; Huang, P.; Xu, Z.; Suleiman, A. A.; Zhang, X.; Bai, X.; Zhou, X.; Zhai, T. Growth of Ultrathin Ternary Teallite (PbSnS2) Flakes for Highly Anisotropic Optoelectronics. _Matter_ 2020, _2_ , 977–987
  * Williams and Hyde 1988 Williams, T.; Hyde, B. Electron microscopy of cylindrite and franckeite. _Physics and Chemistry of Minerals_ 1988, _15_ , 521–544
  * Costa et al. 2021 Costa, V. Z.; Liang, L.; Vaziri, S.; Miller, A.; Pop, E.; Newaz, A. Vibrational Properties of a Naturally Occurring Semiconducting van der Waals Heterostructure. _The Journal of Physical Chemistry C_ 2021, _125_ , 21607–21613
  * Garcia-Basabe et al. 2021 Garcia-Basabe, Y.; Steinberg, D.; Daminelli, L. M.; Mendoza, C. D.; de Souza, E. A. T.; Vicentin, F. C.; G. Larrudé, D. Charge-transfer dynamics in van der Waals heterojunctions formed by thiophene-based semiconductor polymers and exfoliated franckeite investigated from resonantly core-excited electrons. _Physical Chemistry Chemical Physics_ 2021, _23_ , 16795–16805
  * Zschiesche et al. 2021 Zschiesche, H.; Aygar, A. M.; Langelier, B.; Szkopek, T.; Botton, G. A. Atomic scale chemical ordering in franckeite—a natural van der Waals superlattice. _Journal of Physics: Condensed Matter_ 2021, _34_ , 055403
  * Paz et al. 2021 Paz, W. S.; Menezes, M. G.; Batista, N. N.; Sanchez-Santolino, G.; Velický, M.; Varela, M.; Capaz, R. B.; Palacios, J. J. Franckeite as an exfoliable naturally occurring topological insulator. _Nano Letters_ 2021, _21_ , 7781–7788
  * Molina-Mendoza et al. 2017 Molina-Mendoza, A. J.; Giovanelli, E.; Paz, W. S.; Niño, M. A.; Island, J. O.; Evangeli, C.; Aballe, L.; Foerster, M.; Van Der Zant, H. S.; Rubio-Bollinger, G.; others Franckeite as a naturally occurring van der Waals heterostructure. _Nature communications_ 2017, _8_ , 1–9
  * Velickỳ et al. 2017 Velickỳ, M.; Toth, P. S.; Rakowski, A. M.; Rooney, A. P.; Kozikov, A.; Woods, C. R.; Mishchenko, A.; Fumagalli, L.; Yin, J.; Zólyomi, V.; others Exfoliation of natural van der Waals heterostructures to a single unit cell thickness. _Nature communications_ 2017, _8_ , 1–11
  * Gant et al. 2017 Gant, P.; Ghasemi, F.; Maeso, D.; Munuera, C.; López-Elvira, E.; Frisenda, R.; De Lara, D. P.; Rubio-Bollinger, G.; Garcia-Hernandez, M.; Castellanos-Gomez, A. Optical contrast and refractive index of natural van der Waals heterostructure nanosheets of franckeite. _Beilstein journal of nanotechnology_ 2017, _8_ , 2357–2362
  * Frisenda et al. 2020 Frisenda, R.; Sanchez-Santolino, G.; Papadopoulos, N.; Urban, J.; Baranowski, M.; Surrente, A.; Maude, D. K.; Garcia-Hernandez, M.; van der Zant, H. S.; Plochocka, P.; others Symmetry breakdown in franckeite: spontaneous strain, rippling, and interlayer Moiré. _Nano Letters_ 2020, _20_ , 1141–1147
  * Xu et al. 2023 Xu, Z.; Xu, Z.; Li, N.; Zhao, C. Ultrafast optical nonlinearity in natural van der Waals heterostructure nanosheets of franckeite. _Results in Physics_ 2023, _49_ , 106492
  * Wang and Kuo 1991 Wang, S.; Kuo, K. Crystal lattices and crystal chemistry of cylindrite and franckeite. _Acta Crystallographica Section A: Foundations of Crystallography_ 1991, _47_ , 381–392
  * Henriksen et al. 2002 Henriksen, R. B.; Makovicky, E.; Stipp, S.; Nissen, C.; Eggleston, C. M. Atomic-scale observations of franckeite surface morphology. _American Mineralogist_ 2002, _87_ , 1273–1278
  * Ray et al. 2017 Ray, K.; Yore, A. E.; Mou, T.; Jha, S.; Smithe, K. K. H.; Wang, B.; Pop, E.; Newaz, A. K. M. Photoresponse of Natural van der Waals Heterostructures. _ACS Nano_ 2017, _11_ , 6024–6030, PMID: 28485958
  * Li et al. 2020 Li, J.; Du, L.; Huang, J.; He, Y.; Yi, J.; Miao, L.; Zhao, C.; Wen, S. Passive photonic diodes based on natural van der Waals heterostructures. _Nanophotonics_ 2020, _10_ , 927–935
  * Costa et al. 2022 Costa, V. Z.; Baker, B.; Sinn, H.-L.; Miller, A.; Watanabe, K.; Taniguchi, T.; Newaz, A. Observation of photoluminescence from a natural van der Waals heterostructure. _Applied Physics Letters_ 2022, _120_ , 253101
  * Li et al. 2020 Li, J.; Yang, K.; Du, L.; Yi, J.; Huang, J.; Zhang, J.; He, Y.; Huang, B.; Miao, L.; Zhao, C.; Wen, S. Nonlinear Optical Response in Natural van der Waals Heterostructures. _Advanced Optical Materials_ 2020, _8_ , 2000382
  * Tripathi et al. 2021 Tripathi, R. P.; Gao, J.; Yang, X. Naturally occurring layered mineral franckeite with anisotropic Raman scattering and third-harmonic generation responses. _Scientific Reports_ 2021, _11_ , 1–10
  * Vianna et al. 2021 Vianna, P. G.; Almeida, A. d. S.; Gerosa, R. M.; Bahamon, D. A.; de Matos, C. J. Second-harmonic generation enhancement in monolayer transition-metal dichalcogenides by using an epsilon-near-zero substrate. _Nanoscale Advances_ 2021, _3_ , 272–278
  * Sousa et al. 2021 Sousa, F. B.; Lafeta, L.; Cadore, A. R.; Sahoo, P. K.; Malard, L. M. Revealing atomically sharp interfaces of two-dimensional lateral heterostructures by second harmonic generation. _2D Materials_ 2021, _8_ , 035051
  * Castellanos-Gomez et al. 2014 Castellanos-Gomez, A.; Buscema, M.; Molenaar, R.; Singh, V.; Janssen, L.; Van Der Zant, H. S.; Steele, G. A. Deterministic transfer of two-dimensional materials by all-dry viscoelastic stamping. _2D Materials_ 2014, _1_ , 011002
  * Ardito et al. 2018 Ardito, F.; Mendes-de Sá, T.; Cadore, A.; Gomes, P.; Mafra, D.; Barcelos, I.; Lacerda, R.; Iikawa, F.; Granado, E. Damping of Landau levels in neutral graphene at low magnetic fields: A phonon Raman scattering study. _Physical Review B_ 2018, _97_ , 035419
  * Marangoni et al. 2021 Marangoni, V. S.; Cadore, A. R.; Ribeiro, H. B.; Hostert, L.; de Matos, C. J.; Silva, C. C.; Seixas, L.; Maroneze, C. M. Long-term environmental stability of nitrogen-healed black phosphorus. _Applied Surface Science_ 2021, _564_ , 150450
  * Woodward et al. 2016 Woodward, R. I.; Murray, R. T.; Phelan, C. F.; de Oliveira, R. E. P.; Runcorn, T. H.; Kelleher, E. J. R.; Li, S.; de Oliveira, E. C.; Fechine, G. J. M.; Eda, G.; de Matos, C. J. S. Characterization of the second- and third-order nonlinear optical susceptibilities of monolayer MoS2 using multiphoton microscopy. _2D Materials_ 2016, _4_ , 011006
  * Bethune 1989 Bethune, D. Optical harmonic generation and mixing in multilayer media: analysis using optical transfer matrix techniques. _Journal of the Optical Society of America B_ 1989, _6_ , 910–916
  * Steinberg et al. 2018 Steinberg, D.; Oré, A. S. M. V.; Zapata, J. D.; de Souza, E. A. T.; de Matos, C. J. S. Third harmonic generation in mechanically exfoliated franckeite. Latin America Optics and Photonics Conference. 2018; p Th3E.3
  * Sar et al. 2020 Sar, H.; Gao, J.; Yang, X. In-plane anisotropic third-harmonic generation from germanium arsenide thin flakes. _Scientific Reports_ 2020, _10_ , 1–10
  * Dasgupta et al. 2020 Dasgupta, A.; Gao, J.; Yang, X. Anisotropic Third-Harmonic Generation in Layered Germanium Selenide. _Laser & Photonics Reviews_ 2020, _14_ , 1900416
  * Youngblood et al. 2017 Youngblood, N.; Peng, R.; Nemilentsau, A.; Low, T.; Li, M. Layer-tunable third-harmonic generation in multilayer black phosphorus. _Acs Photonics_ 2017, _4_ , 8–14
  * Yang and Xie 1995 Yang, X.-L.; Xie, S.-W. Expression of third-order effective nonlinear susceptibility for third-harmonic generation in crystals. _Applied Optics_ 1995, _34_ , 6130–6135
  * Ray et al. 2017 Ray, K.; Yore, A. E.; Mou, T.; Jha, S.; Smithe, K. K. H.; Wang, B.; Pop, E.; Newaz, A. K. M. Photoresponse of Natural van der Waals Heterostructures. _ACS Nano_ 2017, _11_ , 6024–6030, PMID: 28485958
  * Li et al. 2021 Li, J.; Du, L.; Huang, J.; He, Y.; Yi, J.; Miao, L.; Zhao, C.; Wen, S. Passive photonic diodes based on natural van der Waals heterostructures. _Nanophotonics_ 2021, _10_ , 927–935
  * Youngblood et al. 2017 Youngblood, N.; Peng, R.; Nemilentsau, A.; Low, T.; Li, M. Layer-Tunable Third-Harmonic Generation in Multilayer Black Phosphorus. _ACS Photonics_ 2017, _4_ , 8–14
  * Mitsas and Siapkas 1995 Mitsas, C. L.; Siapkas, D. I. Generalized matrix method for analysis of coherent and incoherent reflectance and transmittance of multilayer structures with rough surfaces, interfaces, and finite substrates. _Applied Optics_ 1995, _34_ , 1678–1683
  * Katsidis and Siapkas 2002 Katsidis, C. C.; Siapkas, D. I. General transfer-matrix method for optical multilayer systems with coherent, partially coherent, and incoherent interference. _Applied Optics_ 2002, _41_ , 3978–3987