Extended landslide velocity and analytical drag

Shiva P. Pudasaini
Technical University of Munich

School of Engineering and Design, Chair of Landslide Research

Arcisstrasse 21, D-80333, Munich, Germany

E-mail<EMAIL_ADDRESS>
Abstract: The landslide velocity plays a dominant role in estimating the
impact force and devastated area. Here, based on Pudasaini and Krautblatter
(2022), I develop a novel extended landslide velocity model that includes the
force induced by the hydraulic pressure gradient, which was neglected by all
the existing analytical landslide velocity models. By a rigorous conversion
between this force and inertia, which facilitates constructing exact
analytical solutions for velocity, I develop two peer systems expecting to
produce the same results. However, this contradicts with our conventional
wisdom. This raises a legitimate question of whether we should develop some
new balance equations such that these phenomena can be better explained
naturally. I compare the two velocity models that neglects and includes the
force induced by the hydraulic pressure gradient. Analytical solutions
produced by the two systems are fundamentally different. The new model is
comprehensive, elegant, and yet an extraordinary development as it reveals
serendipitous circumstances resulting in a pressure-inertia-paradox.
Surprisingly, the mass first moves upstream for quite a while, then it winds
back and continues accelerating down slope. The difference between the
extended and simple solution is significant, and widens strongly as the force
associated with the hydraulic pressure gradient increases, demonstrating the
importance of this force in the landslide velocity. The viscous drag is an
essential dissipative force mechanism and plays an important role in
controlling the landslide dynamics. However, no explicit mechanical and
analytical model exists to date for this. The careful sagacity of the graceful
form of new velocity equation results in a plain, yet mechanically extensive,
analytical model for viscous drag, the first of this kind. It contains several
physical and geometrical parameters, and evolves dynamically as it varies
inversely with the flow depth. A dimensionless drag number is constructed
characterizing the drag dynamics. Importance of the drag model is explained.
In contrast to the prevailing practices, I have proven that drags are
essentially different for the expanding and contracting motions. This is an
entirely novel perception. Drag coefficients are close to the empirical or
numerical values often used in practice. But, now, I offer an innovative,
physically-founded analytical model for the drag that can be instantly applied
in mass flow simulations.

##  1 Introduction

Mountain torrents and valleys are often hit by landslides and debris flows
threatening the lives, devastating civil properties, and eternally swapping
the landscape (Evans et al., 2009; Mergili et al., 2020, Shugar et al., 2021).
Adequate understanding of the dynamics, runout and deposition processes of
mass transports still remains a great challenge to the mountain engineers,
practitioners, flow modelers and the hazard assessment professionals (Cuomo et
al., 2016; Dietrich and Krautblatter, 2019). The dynamics of landslide and the
associated hazard is directly related to, and primarily controlled by, its
velocity. So, the landslide velocity plays a fundamental role for the
assessment of landslide hazards, mitigation measures, design of protective
structures and landuse planning (Johannesson et al., 2009; Dowling and Santi,
2014; Faug, 2015). However, our understanding of the mechanical controls of
evolving velocity, runout and impact energy of the landslide is still poor.
Appropriate modelling of landslide impact force demands for a comprehensive
knowledge of its velocity.
There are three major ways of studding the landslide dynamics: field or
laboratory investigations and experiments, physical-mathematical modelling and
numerical simulations, and analytical solutions (Pudasaini and Hutter, 2007).
The available field data on landslides play important role in acquiring basic
information on their dynamics. However, such data are insufficient and cannot
cover the broad spectrum. The properties of the landslide material and its
interaction with the boundary are difficult to determine to a high degree of
accuracy. The dynamic field data are scarce and the aftermath static data are
often only available for single locations. Much of the low resolution
measurements are locally or discretely based on points in time and space
(McCoy et al., 2012; Theule et al., 2015; de Haas et al., 2020). So, the
proper understanding and interpretation of the field data pose challenges. The
well controlled and desirably obtained data from laboratory or field
experiments (Iverson and Ouyang, 2015; de Haas and van Woerkom, 2016; Baselt
et al., 2021) can overcame many of the problems associated with the field
events. But, still such data cannot cover the whole picture of the dynamics
and deposition because of their largely reduced scales, substantially
simplified materials and flow paths, and neglection of many of the important
interactions between the involved phases. Other fundamental problems are also
directly associated with the complex interactions of the natural events with
their boundaries, which cannot be fully understood by controlled experiments.
Due to its overall nature and capabilities of covering the complex materials
and their evolutions, interactions with the environment, applicabilities to
any natural slopes and torrents, the theoretical modelling (Le and Pitman,
2009; Pudasaini, 2012; Pudasaini and Mergili, 2019) have gain a tremendous
pace in the past decades, and have become the major solutions of the problems
associated with the mass flow dynamics. A number of comprehensive numerical
modelling for mass transports are available (McDougall and Hungr, 2005;
Christen et al., 2010; Frank et al., 2015; Iverson and Ouyang, 2015; Cuomo et
al., 2016; Mergili et al., 2020; Liu et al. 2021). Yet, numerical simulations
are approximations of the physical-mathematical model equations and their
validity is often evaluated empirically (Mergili et al., 2020). In contrast,
exact, analytical solutions can provide better insights into the complex flow
behaviors, and provide the fastest, the most cost-effective and the best
rigorous answer to the problem (Gauer, 2018; Faraoni, 2022; Pudasaini and
Krautblatter, 2022). Moreover, exact solutions to non-linear model equations
are necessary to elevate the accuracy of numerical solution methods based on
complex numerical schemes (Chalfen and Niemiec, 1986; Pudasaini and
Krautblatter, 2022). This is very useful to interpret complicated simulations
and avoid mistakes in numerical simulations. However, the numerical solutions
(Mergili et al., 2020; Shugar et al., 2021) can cover the widest spectrum of
the complex flow dynamics described by advanced mass flow models (Pudasaini
and Mergili, 2019). Once such simulations are tested and validated against the
analytical solutions, they may provide even more accurate results than the
simplified analytical solutions.
Since the middle of the 20th century, a number of analytical solutions have
been presented for landslides, avalanches and debris flows (Voellmy, 1955;
Salm, 1966; Perla et al., 1980; McClung, 1983). The existing analytical
velocity models and their solutions are effectively valid to the mass point or
center of mass motion, and did not consider internal deformation. Pudasaini
and Krautblatter (2022) overcame those deficiencies by introducing a simple,
physics-based general analytical landslide velocity model that simultaneously
incorporates the internal deformation and externally applied forces,
consisting of the net driving force and the viscous resistant. They showed
that the non-linear advection and external forcing fundamentally regulate the
state of motion and deformation. They constructed several general exact
analytical solutions. Those solutions cover the wider spectrum of landslide
velocity and directly reduce to the mass point motion.
Here, following Pudasaini and Krautblatter (2022), I develop a novel extended
landslide velocity model that includes the lateral force induced by the
hydraulic pressure gradient, which, however, was neglected by all the existing
analytical landslide velocity models. By a fully consistent conversion between
this force and inertia, I develop two peer systems. Surprisingly, the
solutions produced by these two systems appear to be fundamentally different.
The difference widens strongly as the force associated with the hydraulic
pressure gradient increases. The new extended landslide velocity model
describes some unexpected situations emerging for the first time. The
mechanism of the new solution is explained based on the contribution from the
hydraulic pressure gradient. I also construct a simple, physically-founded
distinctive analytical model for the viscous drag. A new dimensionless drag
number is formulated. I prove that drags are essentially different for the
expanding and contracting motions. The exact, analytical solutions derived in
this paper may find applications in environmental, engineering and industrial
mass transports down slopes and channels where the force induced by the
hydraulic pressure gradient is significant which is the realistic situation.
Similarly, the analytical drag model may serve the purpose of the
practitioners and engineers involved in mass flow simulations and hazard
mitigation.

##  2 The model development

A geometrically two-dimensional motion down a slope is considered. Let $t$ be
time, $(x,z)$ be the coordinates and $\left(g^{x},g^{z}\right)$ the gravity
accelerations along and perpendicular to the slope, respectively. Let, $h$ and
$u$ be the flow depth and the mean flow velocity of the landslide along the
slope. Similarly, $\gamma,\alpha_{s},\mu$ be the density ratio between the
fluid and the solid particles $\left(\gamma=\rho_{f}/\rho_{s}\right)$, volume
fraction of the solid particles (coarse and fine solid particles), and the
basal friction coefficient $\left(\mu=\tan\delta\right)$, where $\delta$ is
the basal friction angle of the solid particles, in the mixture material.
Furthermore, $K$ is the earth pressure coefficient, and $\beta$ is the viscous
drag coefficient. Then, consider the mass and momentum balance equations
describing the motion of a landslide (Pudasaini and Krautblatter, 2022):

$\frac{\partial h}{\partial t}+\frac{\partial}{\partial x}\left(hu\right)=0,$
(1) $\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
x}+\lambda\frac{\partial h}{\partial x}=\alpha-\beta u^{2},$ (2)

where

$\alpha:=g^{x}–(1-\gamma)\alpha_{s}g^{z}\mu,$ (3)

is the net driving force (this includes the forces due to gravity, Coulomb
friction, lubrication, and liquefaction) and $\beta$ is the viscous drag
coefficient. So, $\alpha$ and $\beta$ constitute the net driving and the
resisting forces in the system. Moreover,

$\lambda:=g^{z}\left\\{\left(\left(1-\gamma\right)K+\gamma\right)\alpha_{s}+\left(1-\alpha_{s}\right)\right\\},$
(4)

is the effective hydraulic pressure parameter. This is (associated with) the
lateral load, and combined with the free-surface gradient ($\partial
h/\partial x$), constitutes the force due to the hydraulic pressure gradient
($\lambda\partial h/\partial x$).
Assuming the steady-state ($\partial u/\partial t=0$), equations (1) and (2)
can, respectively, be written as

$\frac{\partial h}{\partial x}=-\frac{h}{u}\frac{\partial u}{\partial x},$ (5)
$u\frac{\partial u}{\partial x}=-\lambda\frac{\partial h}{\partial
x}+\alpha-\beta u^{2}.$ (6)

As we will see below, it is crucial to carefully analyze (6) in relation with
(5). For the landslide accelerating downslope, both the parameters $\alpha$
and $\lambda$ are positive, so is the total system force, $\alpha-\beta
u^{2}$. As the landslide accelerates down the slope, the velocity increases
and its surface gradient decreases, i.e., $\partial u/\partial x>0$, $\partial
h/\partial x<0$. Moreover, the velocity and depth of the landslide are
positive quantities. Then, the term $-\lambda\partial h/\partial x$ in (6),
which stems from the hydraulic pressure gradient, is a positive force that
enhances the motion. This means, all the force components on the right hand
side of (6) are positive, and consequently, $\left(-\lambda\partial h/\partial
x+\alpha-\beta u^{2}\right)>0$. The left hand side of (6) is the inertia of
the system. The advective term, $u\partial u/\partial x$, tells that the
landslide propagates downslope with the advection velocity $u$. And, this
advection is caused by the overall composite forcing $\left(-\lambda\partial
h/\partial x+\alpha-\beta u^{2}\right)$. So, the positive force pulls the
landslide along the slope in the positive $x$-direction, as expected.
Dynamically, the left and right sides of (5) are equivalent as it satisfies
the mass balance. Now, combining the momentum balance equation (6) with the
mass balance equation (5), I obtain:

$\displaystyle u\frac{\partial u}{\partial x}-\lambda\frac{h}{u}\frac{\partial
u}{\partial x}=\alpha-\beta u^{2}.$ (7)

This can be written as

$\displaystyle\left(u-\frac{\lambda h}{u}\right)\frac{\partial u}{\partial
x}=\alpha-\beta u^{2},$ (8)

which describes the landslide motion with the advection velocity
$\displaystyle{\left(u-\frac{\lambda h}{u}\right)}$.
The pressure-inertia-paradox: As convenient, I analyze the landslide dynamics
with respect to (7) or (8). It is clear that as compared with the inertia
$u\partial u/\partial x$ in (6), the inertia $\left(u-\lambda
h/u\right)\partial u/\partial x$ in (8) is decreased, because $\lambda h/u>0$.
In principle, the systems represented by (6) and (7) should be equivalent,
because there is a rigorous and fully consistent conversion between the force
and inertia. In (6), $\left(-\lambda\partial h/\partial x\right)$ acts as a
force contribution induced by the hydraulic pressure gradient, which is
positive for the accelerating landslide, so it enhances the motion. This
means, the landslide described by the system (6) accelerates downslope with
the higher rate when $-\lambda\partial h/\partial x$ is present than when it
is absent. In (7), $-\left(\lambda h/u\right)\partial u/\partial x$ acts as an
inertial contribution, and due to its negative sign, this also enhances the
motion and the mass should accelerate downslope with the higher rate when
$-\left(\lambda h/u\right)\partial u/\partial x$ is present than when it is
absent. Dynamically, increased force is equivalent to consistently decreased
inertia. Consequently, the two peer systems (6) and (7) are expected to
produce the same results. However, this does not appear to be happening here.
To prove this, consider an initially relatively small velocity $u_{0}>0$. This
is a legitimate choice, as often landslide is triggered with zero or small
positive initial velocity. For the accelerating system under consideration,
the total system force in (8) is positive, i.e., $\left(\alpha-\beta
u^{2}\right)>0$. Which says that the mass must slide and accelerate down the
slope. However, the inertial term on the left hand side of (8) tells that, as
long as $u_{0}<\sqrt{\lambda h}$, the mass advects in the negative direction,
means, it accelerates upslope for a certain time and distance. This
contradicts with our conventional wisdom that the mass must accelerate
downslope. This is a paradox, which I call the pressure-inertia-paradox,
$P_{pi}$. This is a novel understanding. Nevertheless, if we begin with the
initial velocity such that $u_{0}>\sqrt{\lambda h}$, then the landslide
described by the system (8) also accelerates downslope as that described by
the system (6). Yet, the solutions produced by these two systems are
fundamentally different. The results are presented in Section 3.1.
As combining momentum balance with the mass balance is a common practice in
fluid mechanics, I have proven that, one must very carefully analyze the
dynamics of the resulting system, and check if it is consistent with the
common perception. This is important. Note that, in principle, the solutions
and results presented here are valid to any dynamical systems that can be
structurally cast in the form (8), irrespective of their specific physical
parameters and meanings.

##  3 Extended landslide velocity model

For ease of analysis, I write (8) as

$\displaystyle\left(u^{2}-\lambda h\right)\frac{1}{u}\frac{\partial
u}{\partial x}=\alpha-\beta u^{2}.$ (9)

This is a simple, yet very interesting equation, that characterizes some key
and abnormal aspects of the flow dynamics. This stems from the elegant form of
the new model. It is important to perceive that (9) implicitly includes the
dynamics associated with the hydraulic pressure gradient via the mass balance.
It also explicitly includes the depth of the landslide, $h$. The model (9) can
be analyzed differently. I explore this in detail.
Equation (9) extends the Pudasaini and Krautblatter (2022) velocity model in
which the hydraulic pressure gradient $\partial h/\partial x=h_{g}$ was
assumed, and the term $\lambda\partial h/\partial x=\lambda h_{g}$ was
combined with $\alpha$ to define the effective driving force
$\alpha^{e}:=\alpha-\lambda h_{g}$. This means, Pudasaini and Krautblatter
(2022) ignored the dynamic contribution of the hydraulic pressure gradient,
but only parameterized it. Moreover, all the existing analytical landslide
velocity models do not consider the dynamics of the hydraulic pressure
gradient. For this reason, I call (9) the extended landslide velocity model.
Assume that the flow depth $h$ can be parameterized with some representative
value, say, typically the initial depth of the landslide, and set
$\Lambda=\lambda h$, which has the dimension of the square of velocity. Then,
(9) takes the form:

$\displaystyle\frac{\partial u}{\partial x}=\left[\frac{u\left(\alpha-\beta
u^{2}\right)}{u^{2}-\Lambda}\right].$ (10)

This equation can be solved analytically to obtain the exact solution:

$\displaystyle\displaystyle{2\alpha\beta\left(x-x_{0}\right)=\ln\left[\left(\frac{u_{0}}{u}\right)^{\left(2\beta\Lambda\right)}\left(\frac{\alpha-\beta
u_{0}^{2}}{\alpha-\beta
u^{2}}\right)^{\left(\alpha-\beta\Lambda\right)}\right]},$ (11)

where, $u_{0}=u(x_{0})$ is the initial (or boundary) condition at a given
channel location $x=x_{0}$. The model (9), and (10) as well as its solution
(11) are new. As $\Lambda\to 0$ (this means: very thin flow depth, or almost
vanishing free-surface gradient, or neglection of the hydraulic pressure
gradient), (11) reduces to the simpler solution to the landslide velocity
model by Pudasaini and Krautblatter (2022):

$u=\sqrt{\frac{\alpha}{\beta}\left[1-\left(1-\frac{\beta}{\alpha}u_{0}^{2}\right)\frac{1}{\exp(2\beta(x-x_{0}))}\right]}.$
(12)

When $\Lambda\neq 0$, the solution (11) becomes physically more interesting
and intricate. Depending on the values of $\Lambda$ and $u_{0}$, it describes
different dynamics. These aspects can now be explained by a single unified
exact analytical solution in (11). Whereas the extended landslide velocity
model and its solution (11) include the dynamically important force associated
with the free-surface gradient, or the hydraulic pressure gradient, and the
flow depth itself, the previous simple model and its solution (12) ignored
this crucial feature. This is the novel aspect of the landslide velocity model
(9).

###  3.1 Results on extended landslide velocity

I consider plausible physical parameters (Pudasaini and Krautblatter, 2022)
as: the slope angle $\zeta=50^{\circ}$, basal friction angle
$\delta=20^{\circ}$, the fluid to solid density ratio $\gamma=1100/2700$,
solid volume fraction $\alpha_{s}=0.65$, and the earth pressure coefficient
$K=1.0$ (extensional motion, Pudasaini and Hutter, 2007). This implies
$\alpha=6.6$, and $\lambda=6.3$. Similarly, consider $\beta=0.0019$. The
results are shown in Fig. 1 for initial conditions $x_{0}=0$, $u_{0}=0.2$, and
different landslide depths $h=2,4,6$. As $\Lambda=\lambda h$, the
corresponding values for $\Lambda$ are $\Lambda=12.61,25.22,37.80$. And,
$\Lambda=0$ indicates the simple solution given by (12), which does not
include the effect of the flow depth $h$ and the hydraulic pressure gradient.
Figure 1 reveals very special aspects of landslide dynamics. The important
observation is that for any (positive) value of the flow depth $h$, the mass
first moves backwards for quite a while, then it winds back and accelerates
down slope. This feature intensifies rapidly as $h$ increases. The model (9)
and its solution (11) describe two very peculiar situations, which I examine
below.

####  3.1.1 The upwind solution

For relatively small values of $u$ such that $\left(u-\Lambda/u\right)<0$, the
mass first advects in the upstream (negative) direction. This means, when the
upstream dynamic wave celerity is negative $u-\sqrt{\Lambda}<0$, or the flow
is sub-critical, $Fr=u/\sqrt{\Lambda}<1$, where $Fr$ is the Froude number,
mass propagates upstream. However, due to the consideration of the mixture as
composed of the solid particles and the fluid, the Froude number here is more
general than the usual Froude number defined for the fluid or the dry grain
flows (Pudasaini, 2012).

Figure 1: Evolution of the landslide velocity as given by the extended
landslide velocity model (9) and its analytical solution (11). The solutions
for $\Lambda>0$ include both the upwind branch and the downwind branch. The
symbols $\bullet$ indicate the terminus (turning points) of the upwind
branches and beginning of the downwind branches at (0,0), (-4.52, 3.55),
(-10.37, 5.02), (-16.70, 6.15) as $\Lambda$ increases, whose values are
indicated in the legend. For $\Lambda=0$, there exists only the downwind
branch.

In order to facilitate the discussion on the landslide dynamics, consider
$\Lambda=37.80$ corresponding to $h=6$. For this, the landslide accelerates
upstream as far as $x_{t}=-16.7$ m with its maximum velocity of about
$u_{t}=6.15$ ms-1, where t stands for terminus as shown in Fig. 1 with the
symbol $\bullet$. I call this the upwind branch of the solution. This is the
situation of the pressure-inertia-paradox, $P_{pi}$. Although mathematically
rigorously derived solution, we still have to find likely scenarios for this
branch of solution to be representative of a physical or natural situation.
However, there can be something extra than our current understanding. We might
even need to develop more general balance equations such that these phenomena
can be better explained naturally.

####  3.1.2 The downwind solution

In another situation, for relatively large values of $u$ such that
$\left(u-\Lambda/u\right)>0$, the mass advects in the downstream (positive)
direction. This means, when the dynamic wave celerity is positive
($u-\sqrt{\Lambda}>0$), or the flow is super-critical
($Fr=u/\sqrt{\Lambda}>1$), the mass propagates downstream. With the above
considered parameters, as the landslide attains its upstream maximum value of
$u_{t}=6.15$ ms-1 at $x_{t}=-16.7$ m, then it winds down (come back), as
indicated by terminus in Fig. 1, to the downstream and continues to accelerate
along the down slope. I call this the downwind branch of the solution, which
can be physically easily explained as usual for any landslide that is
triggered with the initial velocity of $u_{0}\geq\sqrt{\Lambda}$. This can be
a scenario for landslides triggered by a strong seismic shaking or a volcanic
eruption. Or, this solution can be applied by considering the initial velocity
at a location such that the condition $u_{0}\geq\sqrt{\Lambda}$ is satisfied.

####  3.1.3 The mechanism of solution (11)

A closer look at the mathematical structure of the extended solution (11)
reveals its intrinsic mechanism and the associated dynamics. There are two
groups of terms internal to (11) inside the logarithm function $\ln$. These
are: ($i$) $\displaystyle{\left(\frac{u_{0}}{u}\right)^{2\beta\Lambda}}$, and
($ii$) $\displaystyle{\left(\frac{\alpha-\beta u_{0}^{2}}{\alpha-\beta
u^{2}}\right)^{\left(\alpha-\beta\Lambda\right)}}$. Before analyzing these
terms, note that, the $\ln$ function is negative when its argument is less
than unity, positive when its argument is more than unity, and zero otherwise.
Now, as the velocity increases from its initial value, that is when $u>u_{0}$,
($i$) becomes less than unity, and $(ii)$ becomes more than unity. However, as
long as $u<\sqrt{\Lambda}$, the product of these terms remains less than unity
and steadily decreases as $u$ increases. For this, the entire $\ln$ function
takes the negative value, resulting in the continuous decrease in $x$, so is
the winding point. This is exactly the reason why for $u<\sqrt{\Lambda}$ the
velocity solution is upwinding as discussed in Section 3.1.1. At the terminus
$\left(u=\sqrt{\Lambda}\right)$, $x$ is minimum with its negative value. After
that, as soon as $u>\sqrt{\Lambda}$, the solution starts downwinding. Yet, it
takes quite a while until the winding point comes back to the positive domain
(or, $x\geq x_{0}$). How long it takes, depends on the value of ${\Lambda}$,
including that of $\alpha,\beta,x_{0}$ and $u_{0}$, such that the condition
$\displaystyle{\left(\frac{u_{0}}{u}\right)^{2\beta\Lambda}\left(\frac{\alpha-\beta
u_{0}^{2}}{\alpha-\beta u^{2}}\right)^{\left(\alpha-\beta\Lambda\right)}}>1$
is satisfied. The higher the value of ${\Lambda}$ the farther the minimum of
$x$ is pushed in the negative $x$ direction. So, the extended solution (11)
has a lot of interesting and important physical aspects to explore.

###  3.2 Essence of the new landslide velocity model

Below, I illustrate some physical situations that are better explained by the
new model than the existing ones. Figure 1 shows that there are significant
differences between the extended solution (11) and the simple solution (12).
The difference widens strongly as $\Lambda$ increases.

Figure 2: Evolution of the landslide velocity as given by the extended
landslide velocity model (9) and its analytical solution (11) for the initial
conditions $(x_{0},u_{0})=(-16.70,6.15)$, and $(x_{0},u_{0})=(0.0,0.0)$. Only
the downwind branch of the extended solution with $\Lambda=37.80$ is
considered. The left end of the extended solution can be extrapolated such
that it hits at some location of the upstream sector of the channel, then
practically this location may be taken as the landslide trigger point.

Comparison between the solutions (11) and (12) should be done carefully. One
way to achieve this is by looking at the physically intuitive downwind
dynamics starting from the terminus. This has been presented in Fig. 2 for
$\Lambda=0$ and $\Lambda=37.80$. The velocity solution that includes the force
induced by the hydraulic pressure gradient (indicated by $\Lambda\neq 0$) and
the depth of the landslide, which is physically important, is substantially
higher than the one that neglects this force contribution. The difference is
huge in the left and large in the right side of the figure. Even at $x=100$ m,
the difference between these solutions is 3.8 ms-1. The difference increases
steadily as the depth of the landslide increases. This is crystal clear from
the physical perspective, because the higher depth of landslide induces the
higher overburden load, that in turn increases the force induced by the
hydraulic pressure gradient resulting in the enhanced motion. The new model
and its solution explicitly demonstrate the essence of this basic physical
fact which was ignored by all previous simple landslide models and their
solutions.
However, as seen in Fig. 2, the left end of the extended solution does not
continuously decrease and hit the upstream sector of the slope. Yet, given the
physical essence and superiority of the extended solution as compared with the
simple solution, this problem can be technically eliminated. There are two
possible ways. First, the left end point of the extended solution can be
reasonably extrapolated such that it hits at some location of the upstream
sector of the slope. Practically, this location can be taken as the triggering
point of the landslide. Second, one might also seek to consider the velocity
of the landslide at a slope position such that $u\geq\sqrt{\Lambda}$. Then,
from that point, the extended solution can be applied customarily.

Figure 3: Evolution of the landslide velocity as given by the extended
landslide velocity model (9) and its analytical solution (11) for the
transformed initial condition $(x_{0},u_{0})=(0.00,6.15)$. Only the downwind
branch of the extended solution is considered.

Another probably practically more relevant comparison is presented in Fig. 3,
where the physically more pertinent downwind solutions are contrasted. This
was made possible by transforming the general solution in Fig. 2 (to the
positive domain) with the initial condition $(x_{0},u_{0})=(0.00,6.15)$, where
the initial velocity is obtained from the relation $u=\sqrt{\Lambda}$. This
allows a more closer look between the two solutions. However, note that, this
transformation reduces the gap between the two solutions. Even for this, the
difference in the velocity is significant. The difference is about 3.3 ms-1 at
position $x=20$ m, which is about 22% larger than the simple solution that
disregards the force induced by the hydraulic pressure gradient that is
present in landslide.
Except for some paradoxical situations (for $u<\Lambda$), solutions presented
in Fig. 1 $-$ Fig. 3 are deemed natural, and clearly manifest the importance
of the inclusion of the dynamically important force induced by the hydraulic
pressure gradient and the depth of the landslide. Moreover, the difference
between these solutions increases steadily with the increase in the landslide
depth (because $\Lambda$ is proportional to the flow depth). This can also be
explained. If we carefully look at the term associated with the hydraulic
pressure gradient, it contains $\partial h/\partial x$ and the overburden
pressure ($\lambda h$, per unit density) that includes the material depth,
$h$. It is more convenient now to realize this in (8) with the inertial term
$-(\lambda h/u)\partial u/\partial x$ that explicitly accommodates the
landslide depth $h$. It tells that, as $h$ increases the inertia decreases,
resulting in the increase of the landslide velocity. This has exactly been
demonstrated in Fig. 1 $-$ Fig. 3 in their admissible domains.
Because, the landslide velocity plays a dominant role in estimating the impact
force exerted by the landslide, which is proportional to the square of its
velocity, the difference between the two velocities mentioned above, in fact,
demonstrates the importance of the inclusion of the force induced by the
hydraulic pressure gradient and the depth in determining the landslide
velocity. However, all the existing analytical landslide velocity models
ignore this important aspect, including the model by Pudasaini and
Krautblatter (2022). This signifies the essence of the new landslide velocity
model (9) and its solution (11).

##  4 A simple mechanical model for drag

###  4.1 The drag model

The viscous drag plays a very important role in controlling and describing the
dynamics of mass flow (Pudasaini and Krautblatter, 2022). However, no clear
and explicit mechanical and analytical model exists to date for the drag
coefficient, and in practical and engineering applications, drag coefficients
are exclusively numerically selected based on the model validation (Frank et
al., 2015; Mergili et al., 2020; Shugar et al., 2021). Next, I utilize (9) and
construct an analytical expression for the viscous drag coefficient $\beta$.
The structure of (9) manifests that as the particle propagates closely with
the wave speed $u\approx\sqrt{\lambda h}$, the force balance $\alpha-\beta
u^{2}\approx 0$ could be attained. This is a novel understanding resulting
from the very graceful form of (9). These two relations can be combined to
eliminate $u$, and obtain a simple analytical, mechanical model for the
viscous drag coefficient:

$\displaystyle{\displaystyle\beta\approx\frac{\alpha}{\lambda}\frac{1}{h}},$
(13)

where, as defined in (3) and (4), $\alpha$ is the net driving force (per unit
mass) along the direction of motion composed of the gravitational acceleration
and the frictional dissipation, and $\lambda$ is the coefficient of the
hydraulic pressure gradient, that constitutes the lateral force (here, in the
main downslope direction). Moreover, $h$ is the flow depth. Structurally, (13)
is simple, yet, a mechanically extensive, innovative model. It is important to
note that (13) is a dynamically evolving drag coefficient which may vanish as
the net driving force nears zero, e.g., as the slope angle approaches the
friction angle.
For dry and extensional ($K\approx 1$) landslides, $\alpha$ and $\lambda$ are
simplified. Then, for the slope angle $\zeta$ and the friction angle $\delta$,
we have $\alpha=g\sin\zeta-g\tan\delta\cos\zeta$, $\lambda=g\cos\zeta$. For
this, (13) reduces to

$\displaystyle{\displaystyle\beta=\frac{1}{h}}\tan\zeta\left[1-\frac{\tan\delta}{\tan\zeta}\right].$
(14)

There is nothing unknown in (13) or (14), because the flow depth $h$ can
either be assumed (known), or be considered as a dynamical quantity, that is
already in the balance equations, and all other physical, mechanical and
geometrical parameters in $\alpha$ and $\lambda$ are known or given. So, (13)
and (14) are closed, and can directly be applied (in the momentum equation) in
modeling and simulation of mass transport.

###  4.2 Importance of the new drag model

The more general drag coefficient (13) should be used in modelling mass flows
by replacing $\beta$ in the momentum balance equation. There are several
important aspects associated with (13).

  * ($i$)

As all the involved physical parameters are known, (13) offers a new physical-
mathematical model for the drag coefficient for mass flow simulation, covering
dense to dilute flows (due to the presence of the solid fraction,
$\alpha_{s}$, and the lubrication factor, $1-\gamma$, in $\alpha$ and
$\lambda$). The Manning-Strickler formula (Ackers and White, 1973) for open
channel hydraulics for the friction coefficient explicitly depends inversely
on the one-third of the flow depth (i.e., $1/h^{1/3}$). However, the drag
coefficient in (13) depends inversely on the flow depth (i.e., $1/h$).
Analogous to Manning’s or Strickler’s number, I call ${\alpha}/{\lambda}$ the
drag number (parameter), and write as
$\displaystyle{\mathcal{P}_{d}={\alpha}/{\lambda}}$, which is a dimensionless
number, as a ratio of two forces; the net driving force along the slope and
the lateral force associated with the hydraulic pressure gradient. Thus, the
drag coefficient $\beta$ is the ratio between the drag number
$\mathcal{P}_{d}$ and the flow depth $h$: $\beta=\mathcal{P}_{d}/h$.

  * ($ii$)

The drag coefficient increases with decreasing flow height. This means, thin
flows experience higher drag which is in agreement with the fluid-type drag
often used in simulation of mass flows (Mergili et al., 2020) to acquire
observed results. But now, I have presented a physically-founded analytical
drag model that serves the purpose of practitioners and engineers.

  * ($iii$)

Drag increases with increasing slope angle and decreases with increasing basal
friction angle. This is consistent with the flow dynamics, because the higher
slope and lower friction result in higher motion, which in turn, increases the
drag force. In other words, the lower slope angle and higher friction reduces
motion and the drag too.

  * ($iv$)

As the slope angle decreases and approaches the friction angle (equivalently,
as the friction angle approaches the slope angle), both the net driving force
$\alpha$ and the viscous drag coefficient $\beta$ tend to vanish. This is the
transition from motion to deposition or the flow stopping. This is how
collectively and synchronously the entire system force ceases and the mass
comes to a halt.

  * ($v$)

The drag is also influenced by the particle concentration and density ratio
between the fluid and the solid in the mixture as well as the extensional and
compactional motions.

For the first time, these realistic dynamically evolving situations could be
described by my modelling approach which was not possible by any existing drag
models as all those models use empirical drag coefficients.

###  4.3 Dynamics of the drag coefficient: the drag curves

Here, I show, how the dynamics of drag coefficient $\beta$ are influenced by
the flow depth ($h$), the slope angle ($\zeta$), the friction angle
($\delta$), and the earth pressure coefficient ($K$) for a dry landslide.
These are called the drag curves.

####  4.3.1 Influence of the flow depth

Figure 4: The dynamics of the drag coefficient $\beta$ as a function of the
flow depth $h$ showing the exponential increase of the drag coefficient as the
flow depth decreases. Also shown are the physical parameters, the slope angle
$\zeta$, the basal friction angle $\delta$, and the earth pressure coefficient
$K$.

Figure 4 shows the dynamics of the viscous drag coefficient $\beta$ as a
function of the flow depth $h$ for given parameters for the slope angle,
friction angle, and the earth pressure coefficient, representing a landslide.
This figure shows that the drag coefficient increases rapidly as the flow
depth decreases, which is consistent with our intuition and the practical
experience. Interestingly, as the flow depth becomes very high, the drop in
the drag is very small. This means, in practice, we may set a realistic lower
limit to the drag coefficient. The influence of the earth pressure coefficient
$K$ on the drag is discussed separately in Section 4.3.4.

####  4.3.2 Influence of the friction angle

Figure 5: The dynamics of the drag coefficient $\beta$ as a function of the
basal friction angle $\delta$ indicating the almost linear decrease of the
drag coefficient as the friction angle increases. Also shown are the physical
parameters, the slope angle $\zeta$, the landslide depth $h$, and the earth
pressure coefficient $K$.

Figure 6: The dynamics of the drag coefficient $\beta$ as a function of the
slope angle $\zeta$. The drag coefficient increases non-linearly as the slope
angle increases. Also shown are the physical parameters, the basal friction
angle $\delta$, the landslide depth $h$, and the earth pressure coefficient
$K$.

Figure 6 displays the dynamics of the viscous drag coefficient $\beta$ as a
function of the basal friction angle $\delta$ for given parameters for a
landslide. The drag coefficient decreases almost linearly as the basal
friction angle increases. Interestingly, $\beta$ is eminently small as the
fiction angle approaches the slope angle. This is as expected, because, as
explained above, this is the state closer to the flow halting.

####  4.3.3 Influence of the slope angle

The dynamics of the viscous drag coefficient $\beta$ as a function of the
slope angle $\zeta$ for considered parameters is presented in Figure 6. Now,
the drag coefficient increases non-linearly as the slope angle increases. The
drag coefficient is nominally small as the slope angle reduces significantly
and approaches the basal friction angle. This is the near cease situation of
motion.

####  4.3.4 Expanding and contracting motions

As it slides down, usually the landslide accelerates, it expands and the
corresponding earth pressure coefficient $K$ takes the values in the range
$K\leq 1$. However, as it enters the transition zone, propagates in the run-
out area, and begins deposition, the motion switches from accelerating to
decelerating states. Consequently, the landslide mass contracts. In this
situation, the earth pressure coefficient can become much higher than unity.
The earth pressure coefficient for the expanding or contracting motions can be
determined dynamically as (Pudasaini and Hutter, 2007):

$\displaystyle{\displaystyle
K=2\sec^{2}\phi\left(1\mp\sqrt{1-\cos^{2}\phi\sec^{2}\delta}\right)-1},$ (15)

where, $\phi$ is the friction angle between the grains, the $-$ sign
corresponds to the expanding motion ($\partial u/\partial x\geq 0$), and the
$+$ sign corresponds to the contracting motion ($\partial u/\partial x<0$).
For the substantially expanding and contracting motions, I respectively took
$K=0.9$, and $K=1.5$.
The results in Fig. 4 $-$ Fig. 6 reveal that, as landslide changes its state
from extensional to compactional motion, the drag coefficient $\beta$ behaves
completely differently. It decreases substantially with increasing values of
$K$. This is physically meaningful, because as the mass starts depositing,
being in the state of contraction, naturally, the landslide should experience
less drag than in extending region. Otherwise, the process of mass
deceleration and transition in to deposition is unphysically rapid if the
higher drag is applied as in the accelerating region, resulting in
unrealistically reduced runout. This means, drag coefficients are essentially
different for expanding and contracting motions. This is an entirely novel
understanding. In applications, often, uniform numerical drag coefficient is
applied globally for the entire flow path $-$ from triggering to the main
track to transition to run-out and deposition without making any distinction
between the extensional and compactional motions. The present analysis proves
that this is unphysical. So, different drag coefficients must be applied
consistently in the mass expanding and contracting regions. These regions are
dynamically determined by the accelerating and decelerating motions and their
corresponding $K$ values, less than unity or greater than unity (Pudasaini and
Hutter, 2007).
Importantly, the values of the drag coefficients displayed in Fig. 4 $-$ Fig.
6 are close to the empirical or the numerical drag coefficients often used in
practice in mass flow simulations (Mergili et al., 2020; Shugar et al., 2021).
But, previously, those values were taken without any physical foundation. We
have now the physically described, fully mechanical and analytical drag model
that can be instantly applied in mass flow simulations.

##  5 Summary

The landslide velocity plays a decisive role in estimating the impact force
and havoc. Here, I developed a novel extended landslide velocity model that
includes the force induced by the hydraulic pressure gradient, which, however,
was neglected by all the existing analytical landslide velocity models. By a
rigorous and fully consistent conversion between this force and inertia, I
developed two peer systems that are expected to produce the same results.
However, this contradicts with our conventional wisdom that the mass must
accelerate downslope. This is a new understanding. The solutions produced by
these two systems are fundamentally different. As combining the momentum
balance with the mass balance is a common practice in fluid mechanics, I have
proven that, one must very carefully analyze the dynamics of the resulting
system, and check its functionality.
The new extended landslide velocity model is simple, elegant, and yet an
extraordinary development that characterizes some key features of the flow
dynamics, and leads to the pressure-inertia-paradox. It describes some
unexpected, but important situations emanating for the first time, the upwind
and downwind solutions. The mass first moves upstream, then it winds back and
accelerates down slope. This feature intensifies rapidly as the flow depth
increases. The mechanism of the new solution is explained based on the
contribution from the hydraulic pressure gradient. I compare the two landslide
velocity models that neglects and includes the force induced by the hydraulic
pressure gradient. There are significant to huge differences between the
extended and the simple solutions. The difference widens strongly as the
lateral force associated with the hydraulic pressure gradient increases. This
signifies the essence of the extended landslide velocity model and its
solution.
The viscous drag plays an important role in controlling the dynamics of mass
flow. However, no explicit mechanical and analytical model exists to date for
drag. In practice, drag coefficients are exclusively numerically selected.
Here, the very graceful form of the new velocity model provided the unique
opportunity to construct a simple, mechanically extensive, analytical model
for the viscous drag. It offers an innovative, physically-founded model for
the drag, covering dense to dilute flows. It varies inversely with the flow
depth and contains several physical, mechanical and geometrical parameters
associated with the flow. The model features that as both the net driving
force and the viscous drag coefficient tend to vanish, concomitantly the
entire system force ceases resulting in the mass halting. I constructed a
dimensionless drag number as the ratio between the net driving force and the
lateral force. Importance of the drag model and its practical applicability is
explained. The drag coefficient increases rapidly as the flow depth decreases,
it decreases almost linearly as the basal friction angle increases, and it
increases non-linearly as the slope angle increases. This is consistent with
our experience. In contrast to the prevailing practices, where drags are
chosen uniformly for the entire flow path, I have proven that drags are
essentially different for the expanding and contracting motions. So, different
drags must be applied consistently in different regions of the landslide. This
is an entirely novel perception. The drag coefficients are close to the
empirical or the numerical values often used in practice. So, the new drag
model serves the purpose of practitioners and engineers dealing with mass flow
simulations.

## References

  * [1] Ackers, P., White, W.R., 1973. Sediment Transport: New Approach and Analysis. Journal of the Hydraulics Division, ASCE, No. HY11.
  * [2] Baselt, I., de Oliveira, G.Q., Fischer, J.-T., Pudasaini, S.P., 2021. Evolution of stony debris flows in laboratory experiments. Geomorphology 372, 107431. https://doi.org/10.1016/j.geomorph.2020.107431.
  * [3] Chalfen, M., Niemiec, A., 1986. Analytical and numerical solution of Saint-Venant equations. Journal of Hydrology 86(1-2), 1-13.
  * [4] Christen, M., Kowalski, J., Bartelt, P., 2010. RAMMS: Numerical simulation of dense snow avalanches in three-dimensional terrain. Cold Regions Science and Technology 63, 1-14.
  * [5] Cuomo, S., Pastor, M., Capobianco, V., Cascini, L., 2016. Modelling the space time evolution of bed entrainment for flow-like landslides. Engineering Geology 212, 10-20.
  * [6] de Haas, T., Nijland, W., de Jong, S.M., McArdell, B.W., 2020. How memory effects, check dams, and channel geometry control erosion and deposition by debris flows. Scientific Reports. 10, 14024. https://doi.org/10.1038/s41598-020-71016-8.
  * [7] de Haas, T., van Woerkom, T., 2016. Bed scour by debris flows: experimental investigation of effects of debris flow composition. Earth Surf. Process. Landforms 41, 1951-1966.
  * [8] Dietrich, A., Krautblatter, M., 2019. Deciphering controls for debris-flow erosion derived from a liDAR-recorded extreme event and a calibrated numerical model (Rossbichelbach, Germany). Earth Surf. Process. Landform 44, 1346-1361.
  * [9] Dowling, C.A., Santi, P.M., 2014. Debris flows and their toll on human life: a global analysis of debris-flow fatalities from 1950 to 2011. Nat. Hazards 71(1), 203-227.
  * [10] Evans, S.G., Bishop, N.F., Smoll, L.F., Murillo, P.V., Delaney, K.B., Oliver-Smith, A., 2009. A re-examination of the mechanism and human impact of catastrophic mass flows originating on Nevado Huascaran, Cordillera Blanca, Peru in 1962 and 1970. Eng. Geol. 108, 96-118
  * [11] Faraoni, V., 2022. Helmoltz problem for the Riccati equation from an analogous Friedmann equation. Eur. Phys. J. C 82. https://doi.org/10.1140/eplc/s10052-021-09966-0.
  * [12] Faug, T., 2015. Depth-averaged analytic solutions for free-surface granular flows impacting rigid walls down inclines. Phys. Rev. E 92\. http://dx.doi.org/10.1103/PhysRevE.92.062310.
  * [13] Frank, F., McArdell, B.W., Huggel, C., Vieli, A., 2015. The importance of entrainment and bulking on debris flow runout modeling: examples from the Swiss Alps. Nat. Hazards Earth Syst. Sci. 15, 2569-2583.
  * [14] Gauer, P., 2018. Considerations on scaling behavior in avalanche flow along cycloidal and parabolic tracks. Cold Reg. Sci. Technol. 151, 34-36.
  * [15] Iverson, R. M., Ouyang, C., 2015. Entrainment of bed material by earth-surface mass flows: review and reformulation of depth-integrated theory. Rev. Geophys. 53(1), 27-58.
  * [16] Johannesson, T., Gauer, P., Issler, D., Lied, K., 2009. In: Barbolini, M., Domaas, U., Harbitz, C.B., Johannesson, T., Gauer, P., Issler, D., Lied, K., Faug, T., Naaim, M. (Eds.), The Design of Avalanche Protection Dams. Recent Practical and Theoretical Developments. European Commision. Directorate General for Research.
  * [17] Le, L., Pitman, E.B., 2009. A model for granular flows over an erodible surface. SIAM J. Appl. Math. 70, 1407-1427.
  * [18] Liu, W., Yang, Z., He, S., 2021. Modeling the landslide-generated debris flow from formation to propagation and run-out by considering the effect of vegetation. Landslides 18, 43–58.
  * [19] McClung, D.M., 1983. Derivation of Voellmy’s Maximum Speed and Run-Out Estimates from a Centre-of-Mass Model. Journal of Glaciology 29(102), 350-352.
  * [20] McCoy, S.W., Kean, J.W., Coe, J.A., Tucker, G.E., Staley, D.M., Wasklewicz, T.A., 2012. Sediment entrainment by debris flows: In situ measurements from the headwaters of a steep catchment. J. Geophys. Res. 117, F03016. doi:10.1029/2011JF002278.
  * [21] McDougall, S., Hungr, O., 2005. Dynamic modelling of entrainment in rapid landslides. Can. Geotech. J. 42, 1437-1448.
  * [22] Mergili, M., Pudasaini, S.P., Emmer, A., Fischer, J.-T., Cochachin, A., Frey, H., 2020. Reconstruction of the 1941 GLOF process chain at lake Palcacocha (Cordillera Blanca, Peru). Hydrol. Earth Syst. Sci. 24, 93-114.
  * [23] Perla, R., Cheng, T.T., McClung, D.M., 1980. A two-parameter model for snow-avalanche motion. J. Glaciology 26(94), 197-207.
  * [24] Pudasaini, S.P., Krautblatter, M., 2022. The Landslide Velocity. Earth Surf. Dynam. 10, 165-189. https://doi.org/10.5194/esurf-2021-81.
  * [25] Pudasaini, S.P., Mergili, M., 2019. A multi-phase mass flow model. Journal of Geophysical Research: Earth Surface, 124, 2920-2942.
  * [26] Pudasaini, S.P., 2012. A general two-phase debris flow model. J. Geophysics. Res. 117, F03010. doi:10.1029/2011JF002186.
  * [27] Pudasaini, S.P., Hutter, K., 2007. Avalanche Dynamics: Dynamics of Rapid Flows of Dense Granular Avalanches. Springer, Berlin, New York.
  * [28] Salm, B., 1966. Contribution to avalanche dynamics. International Symposium on Scientific Aspects of Snow and Ice Avalanches, 1965, Davos; pp. 199-214: IAHS Publ. No. 69.
  * [29] Shugar, D.H., Jacquemart, M., Shean, D., et al., 2021. A massive rock and ice avalanche caused the 2021 disaster at Chamoli, Indian Himalaya. Science 373, 300-306.
  * [30] Theule, J.I., Liebault, F., Laigle, D., Loye, A., Jaboyedoff, M., 2015. Channel scour and fill by debris flows and bedload transport. Geomorphology 243, 92-105.
  * [31] Voellmy. A., 1955. Über die Zerstörungskraft von Lawinen. Schweizerische Bauzeitung. Jahrg. 73. Ht. 12., 159-162; Ht. 15, 212-217; Ht. 17, 246-249: Ht. 19, 280-285. On the destructive force of avalanches, Translation No. 2. Alta. Avalanche Study Center, USDA, Forest Service, 1964\.


# Lasing modes in ZnO nanowires coupled to planar metals

Daniel Repp  Institute of Applied Physics, Abbe Center of Photonics, Friedrich
Schiller University Jena, Albert-Einstein-Straße 15, 07745 Jena, Germany
Angela Barreda  Institute of Applied Physics, Abbe Center of Photonics,
Friedrich Schiller University Jena, Albert-Einstein-Straße 15, 07745 Jena,
Germany  Institute of Solid State Physics, Friedrich-Schiller-Universität
Jena, Max-Wien-Platz 1, 07743 Jena, Germany  Francesco Vitale  Institute of
Solid State Physics, Friedrich-Schiller-Universität Jena, Max-Wien-Platz 1,
07743 Jena, Germany  Isabelle Staude  Institute of Applied Physics, Abbe
Center of Photonics, Friedrich Schiller University Jena, Albert-Einstein-
Straße 15, 07745 Jena, Germany  Institute of Solid State Physics, Friedrich-
Schiller-Universität Jena, Max-Wien-Platz 1, 07743 Jena, Germany  Ulf Peschel
Institute of Condensed Matter Theory and Solid State Optics, Friedrich-
Schiller-Universität, Max-Wien-Platz 1, Jena 07743, Germany  Carsten Ronning
Institute of Solid State Physics, Friedrich-Schiller-Universität Jena, Max-
Wien-Platz 1, 07743 Jena, Germany  Thomas Pertsch  Institute of Applied
Physics, Abbe Center of Photonics, Friedrich Schiller University Jena, Albert-
Einstein-Straße 15, 07745 Jena, Germany  Fraunhofer Institute for Applied
Optics and Precision Engineering, Albert-Einstein-Straße 7, 07745 Jena,
Germany

###### Abstract

Semiconductor nanowire lasers can be subject to modifications of their lasing
threshold resulting from a variation of their environment. A promising choice
is to use metallic substrates to gain access to low-volume Surface-Plasmon-
Polariton (SPP) modes. We introduce a simple, yet quantitatively precise model
that can serve to describe mode competition in nanowire lasers on metallic
substrates. We show that an aluminum substrate can decrease the lasing
threshold for ZnO nanowire lasers while for a silver substrate, the threshold
increases compared with a dielectric substrate. Generalizing from these
findings, we make predictions describing the interaction between planar metals
and semiconductor nanowires, which allow to guide future improvements of
highly-integrated laser sources.

##  1 Introduction

Semiconductor nanowires have attracted a broad interest in recent years and
have been studied in fields as diverse as nanosensing [1], nanoscale
optoelectronics [2] and thermoelectrics [3]. They are of particular interest
in the field of nanophotonics, as they are wavelength-scale nanolasers [4, 5].
ZnO nanowires are interesting for UV-applications due to their direct bandgap
in this spectral range, corresponding to an emission wavelength
$\lambda\approx$380\text{\,}\mathrm{nm}$$, and high exciton binding energy,
supporting excitonic emission and electron-hole-plasma formation even at room
temperature [6].
Nanowire-based lasers form a cylindrical, quasi-one-dimensional cavity and
have transversal geometrical extents close to or below optical wavelengths.
Since the fabrication of Bragg-reflectors during nanowire growth is
challenging at the nanoscale [7], the cavity feedback is often restricted to
the back-reflection of the end facets of the nanowires. Furthermore, the
lasing properties of nanolasers can be changed by interaction with the
substrate. For example, the radiation emitted by the nanowire can interact
with the free electron plasma of metallic substrates, forming SPPs [8, 9],
which can facilitate a reduction of the mode volume [10]. Metals known to
provide this mode confinement for ZnO nanowires have plasma frequencies in the
ultraviolet, like aluminum and silver [11].
The reduced mode volume ($V$) leads to a faster emission dynamics via the
Purcell effect [12], which is quantitatively described by the Purcell factor
$F$ [13]:

$F=\frac{\gamma}{\gamma_{0}}=\frac{3\lambda^{3}}{4\pi^{2}}\frac{Q}{V},$ (1)

where $\gamma$ is the emission rate of a recombining electron-hole pair in the
structure of interest (here the nanowire on top of the metal), $\gamma_{0}$ is
the emission rate of the same electron-hole pair in bulk ZnO, $\lambda$ is the
emission wavelength, $Q$ is the quality factor of the resonance and $V$ is the
mode volume. Usually, a small dielectric spacer is sandwiched between the
nanowire and the underlying metal substrate to protect the nanowire emission
from quenching by the metal [14]. An illustration of the entire structure can
be seen in Fig. 1.

Figure 1: Scheme of the system under investigation in this article (not to
scale). The coordinate system used throughout this article as well as typical
geometrical scales are indicated.

The Purcell effect, which is applicable to emission of radiation in general,
is discussed in many different lasing applications [15, 16]. In fact, the
system described above (nanowire on a metallic substrate) is often called a
spaser, which is an acronym for “Surface Plasmon Amplification via Stimulated
Emission of Radiaton”, in analogy to the laser [17]. The physical principle of
a spaser is equivalent to that of a laser, only the amplified bosonic mode is
a surface plasmon polariton instead of a photon [18].
In this article, a thorough investigation of the influence of nanowire size
and substrate material on the Purcell factor is presented. This information is
combined with detailed mode calculations and used in a lasing model to
calculate mode-specific lasing thresholds and thereby quantify mode
competition in the system. As a result, we gain fundamental, non-trivial
insights into the mode competition processes for ZnO nanowire lasers coupled
to planar metals and highlight trade-offs that have to be considered in the
device design. Furthermore, we provide a versatile, quantitative tool for the
description, investigation and optimization of semiconductor lasers in the
spaser geometry. Also, we define two coupling regimes between SPPs and
nanowire emission, having distinct advantages and disadvantages.

##  2 Numerical methods

###  2.1 Mode calculation

An important part of the quantitative description of the lasing process is the
proper characterization of the lasing modes regarding critical parameters like
propagation losses and gain overlap. The MODE package of the Lumerical
software suite was used in order to obtain the transversal eigenmodes of the
simulated structures [19]. This package implements an EigenMode Expansion
(EME) on a Yee-grid. The rectangular Yee grid is known to lose accuracy for
metal interfaces and curved surfaces, both of which are present in our system.
Treating the system therefore often requires very fine grids. In order to
obtain highly accurate results, a technique called “Conformal Meshing” is
used, which calculates effective permittivities close to interfaces and
greatly reduces the computational costs.
To correctly connect the Purcell factor and the different mode profiles, a
normalization procedure of the modes has to be chosen. The Purcell factor
describes the modification of emission rates as a result of changes in the
environment. In quantum mechanics, Fermi’s Golden rule describes the rate of
decay and, correspondingly, photon emission. The matrix element responsible
for the description of decay process in Fermi’s Golden rule, however, is
proportional to the field amplitude. Therefore, a normalization with respect
to the field strength has to be chosen, which is defined as

$\int_{A}\sum_{j}|E_{i,j}(\textbf{r})|^{2}=1,$ (2)

where $A$ is a cross section through the mode perpendicular to the nanowire
axis and $E_{i,j}(\textbf{r})$ is the $j$ component of the electric field
strength of mode $i$ at position r.

Figure 2: Example plots of $|E|^{2}$ of the considered modes for a ZnO
nanowire of a diameter of $200\text{\,}\mathrm{nm}$ on top of a silver
substrate with an interposed silica nanometric spacer. Plotted are the a)
photonic mode ($n_{eff}=2.1+0.002i$), b) plasmonic mode ($n_{eff}=3.78+1.76i$)
and c) hybrid mode ($n_{eff}=2.07+0.02i$).

The profiles of some of the transversal modes are plotted in Fig. 2. They are
similar to the ones described in ref. [20]. The plots include the nanowire,
the metallic substrate and a dielectric spacer in between. We call the
nanowire bound mode with a strong maximum in the nanowire center, and almost
no visible interaction with the metal, a “photonic mode”, as plotted in Fig.
2a.
The main interest of this investigation is however on highly confined
“plasmonic” modes. These modes exhibit a very small mode volume which leads to
a large Purcell factor. The fields of plasmonic modes are mainly confined
within the spacer, as shown in Fig. 2b.
A third “hybrid” mode can consists of a nanowire waveguide mode, which is
dominantly polarized orthogonally to the metallic substrate. Compared to the
plasmonic mode, the overlap with the nanowire material is increased. An
example is plotted in Fig. 2c.

###  2.2 Purcell factor determination

The Purcell factor was determined by performing Finite-Difference-Time-Domain
(FDTD) simulations [19], utilizing a dipole source placed at points inside the
nanowire, which maximize the coupling to the modes described in Sec. 2.1. To
characterize the photonic mode, the dipole was placed at the nanowire center.
For the plasmonic- and hybrid modes, the dipole was placed closely to the
metal. Its emission wavelength was set to $\lambda=$380\text{\,}\mathrm{nm}$$,
which corresponds to the bandgap of ZnO.
The E and H fields emitted by the dipole were recorded on the six sides of a
small cube (side length of 5 nm) placed around the emitter. The Fourier-
transforms of the recorded fields were subsequently used to obtain the
spectrum of Poynting vectors, which was integrated over all sides of the cube,
yielding the spectral power distribution [21]:

$P(\omega)=\frac{1}{2}\int_{\text{Box}}Re(\textbf{E}(\omega)\times\textbf{H}^{*}(\omega))dA.$
(3)

This distribution was subsequently normalized to the injected source power,
yielding the Purcell factors $F_{j}$ by arguments presented in the
supplementary document Sec. 1A, where $j$ indicates the dipole orientation.
To attribute the Purcel factor calculated at position r to a specific mode
$i$, the normalized mode profiles $\textbf{E}_{i}(\textbf{r})$ are used. At
position $r$ and with orientation $j$, the Purcell factor attributed to mode
$i$ was calculated as

$F_{i,j}(\textbf{r})=F_{j}(\textbf{r})\cdot\frac{|E_{i,j}(\textbf{r})|^{2}}{\Sigma_{l}|E_{l,j}(\textbf{r})|^{2}},$
(4)

where the index $l$ runs over all modes.
The Purcell factor, as defined above, is a vectorial quantity. To describe the
Purcell factor attributable to mode $i$ as a scalar value in a rate-equation
model, the stochastic distribution of dipole orientations has to be taken into
account. Assuming the electron-hole plasma responsible for gain in ZnO
nanowire lasers to be isotropic, the dipole orientation was considered equally
distributed, leading to a simple averaging procedure [22]:

$F_{i}(\textbf{r})=\frac{\sum_{j}F_{i,j}(\textbf{r})}{3}$ (5)

where $j$ indicates the dipole orientation.
The nanowire was modeled as an isotropic, non-dispersive dielectric material
with a refractive index of $2.5$ to conform to the conditions at the emission
wavelength of $\lambda=$380\text{\,}\mathrm{nm}$$. The chosen refractive index
value corresponds approximately to the value of the refractive index of
crystalline ZnO from ref. [23] and is supported by a calculation from ref.
[24] at the Mott density $N_{\text{Mott}}=2\cdot 10^{24}${\mathrm{m}}^{-3}$$
(value from [25]), as can be seen also in ref. [26]. The actual value is seen
to vary within a 10-20% range, which is an effect that was neglected to keep
the model simple. The possible errors caused by this particular choice are
elaborated on in Sec. 2.4 and in the supplementary document Sec. 2D. The
imaginary part of the refractive index was set to zero, as losses and gain are
incorporated into the gain profile in the lasing model (see Sec. 2.3).
The nanowire was placed on top of a silica spacer-layer with a thickness of 5
nm, which separates the nanowire from a 100 nm thick metal layer, with the
metal being either aluminum or silver. The simulation domain was discretized
with a uniform mesh of 1 nm to guarantee the convergence of the results.
Perfectly Matched Layer (PML) boundary conditions were applied to avoid
reflections on the domain boundaries. Both nanowire end facets as well as the
metal were extended into the PML to disable end facet reflections, which
effectively corresponds to an infinitely long nanowire. Reflections at the
nanowire end facets were addressed with a separate FDTD calculation (see
supplementary document).

###  2.3 Lasing model

A multimode rate-equation lasing model, adapted from ref. [27], was used to
describes the lasing process. Here, a pumping rate $P$ describes optical
excitation of charge carriers. After the relaxation of charge carriers, they
recombine and emit photons, which couple to specific modes $i$ as described by
the Purcell factors $F_{i}$ calculated above. The rate equations governing the
evolution of the charge carrier density $n$ as well as the photon densities
$s_{i}$ in the three competing modes as shown in Fig. 2 are:

$\displaystyle\derivative{n}{t}$ $\displaystyle=\eta\frac{P}{\hbar\omega
V_{\text{a}}}-\frac{n}{\tau_{\text{nr}}}-(F_{1}\beta_{1}+F_{2}\beta_{2}+F_{3}\beta_{3})\frac{n}{\tau_{\text{sp}}}$
$\displaystyle\quad-\nu_{\text{g,a}}g(n)(\Gamma_{1}F_{1}s_{1}+\Gamma_{2}F_{2}s_{2}+\Gamma_{3}F_{3}s_{3})$
(6) $\displaystyle\derivative{s_{1}}{t}$
$\displaystyle=F_{1}\beta_{1}\frac{n}{\tau_{\text{sp}}}+\nu_{\text{g,a}}g(n)\Gamma_{1}F_{1}s_{1}-\frac{s_{1}}{\tau_{\text{p,1}}}$
(7) $\displaystyle\derivative{s_{2}}{t}$
$\displaystyle=F_{2}\beta_{2}\frac{n}{\tau_{\text{sp}}}+\nu_{\text{g,a}}g(n)\Gamma_{2}F_{2}s_{2}-\frac{s_{2}}{\tau_{\text{p,2}}}$
(8) $\displaystyle\derivative{s_{3}}{t}$
$\displaystyle=F_{3}\beta_{3}\frac{n}{\tau_{\text{sp}}}+\nu_{\text{g,a}}g(n)\Gamma_{3}F_{3}s_{3}-\frac{s_{3}}{\tau_{\text{p,3}}}.$
(9)

In this model, $\eta$ is the absorption efficiency of the optical pump, which
is set to 1 to model efficient one-photon above-bandgap excitation, $P$ is the
pumping rate, $\hbar$ is the reduced Planck constant, $\omega$ is the incoming
photon angular frequency and $V_{\text{a}}$ is the volume of the active
medium. $1/\tau_{\text{nr}}$ is the nonradiative decay rate. This rate is kept
constant, which neglects higher order nonradiative recombination processes
like Auger recombination, to simplify the model. $1/\tau_{\text{sp}}$ is the
rate of spontaneous emission. The exciton lifetime $\tau_{total}$ in ZnO
nanowire is in the range of 200-300 ps [28, 29], combining radiative
($\tau_{sp}$) and non-radiative ($\tau_{nr}$) lifetime. Measurements of the
Internal Quantum Efficiency (IQE) yield values close to 50% [30]. Therefore,
we assume that $\tau_{sp}=\tau_{nr}$, resulting in

$\displaystyle\frac{1}{\tau_{total}}$
$\displaystyle=\frac{1}{\tau_{sp}}+\frac{1}{\tau_{nr}}$ (10)
$\displaystyle=\frac{2}{\tau_{sp}},$ (11)

which yields $\tau_{sp}=\tau_{nr}=2\tau_{total}=$700\text{\,}\mathrm{ps}$$.
The confinement factors $\Gamma_{i}$ were defined as

$\displaystyle\Gamma_{i}$
$\displaystyle=\frac{\int_{V_{\text{a}}}\epsilon(\textbf{r})|E_{i}(\textbf{r})|^{2}}{\int_{V}\epsilon(\textbf{r})|E_{i}(\textbf{r})|^{2}}$
(12)
$\displaystyle\approx\frac{\int_{A_{\text{a}}}\epsilon(\textbf{r})|E_{i}(\textbf{r})|^{2}}{\int_{A}\epsilon(\textbf{r})|E_{i}(\textbf{r})|^{2}}.$
(13)

This expression, although not accurate for open systems, leads to reliable
results. Special care is necessary for open systems, which can be treated by
quasi-normal modes [31]. In Eq. 12, $V$ is the total volume of the simulation
domain, $\epsilon(\textbf{r})$ is the permittivity distribution and
$\textbf{E}_{i}(\textbf{r})$ the field profile of mode $i$. The integration
along the nanowire axis cancels out because of translation symmetry, leaving
the actual integration to be performed on the nanowire cross section $A_{a}$
and the cross section of the simulation domain $A$. There is a slight error
caused by the termination of the nanowire which decreases with increasing
nanowire length.
Note that, while $\Gamma_{i}$ is calculated for the entire simulation domain,
the value is dominated by the maximum of the specific mode. Even more so
because of the high value of the nanowire permittivity index. This justifies
the connection of a confinement factor $\Gamma_{i}$ and a Purcell factor
$F_{i}$ calculated for only one dipole position.
Also, note that the fraction of spontaneous emission into mode $i$ was chosen
as $\beta_{i}=\frac{\Gamma_{i}}{\Gamma_{1}+\Gamma_{2}+\Gamma_{3}}$. This
choice is a sum over all longitudinal cavity modes in the nanowire having the
same transversal mode profile. It neglects all other modes except for the ones
considered, as well as emission in free space, and weights the amount of
spontaneous emission into mode $i$ by its overlap with the gain material. Note
that this is a heuristic approach motivated by ref. [32].
Furthermore, the group velocity is defined as

$\nu_{\text{g,a}}=\frac{c}{\derivative{(\omega n^{\prime})}{\omega}},$ (14)

with $n^{\prime}$ being the active medium refractive index, and $g(n)$ is the
active medium gain as in [27], given by

$g(n)=g_{0}\ln(\frac{n}{n_{\text{tr}}}),$ (15)

with $n_{\text{tr}}=n_{\text{mott}}=2\cdot 10^{24}${\mathrm{m}}^{3}$$ the
transparency density, which we set equal to the Mott density for formation of
an electron-hole plasma [24]. The transparency density changes with
wavelength, but the Mott density gives a good order-of-magnitude estimate.
Furthermore, we set $g_{\text{0}}=$0.7\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{-1}$$ [33].
The modal losses are modelled via a loss rate, given by

$\frac{1}{\tau_{\text{P},i}}=\nu_{\text{gz},i}\left(\alpha_{i}+\frac{1}{2L}\ln(\frac{1}{R_{i}^{2}}),\right)$
(16)

where $\nu_{\text{gz},i}=c/n^{\prime}_{i}$, with $n^{\prime}_{i}$ being the
mode group index, is the mode group velocity, $\alpha_{i}$ is the modal
propagation loss, $L$ the length of the nanowire (we set
$L=$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$) and $R_{i}$ the end facet
reflectivity for mode $i$.

###  2.4 Thresholds

To identify the lasing threshold, the initial value problem defined by the
lasing model was solved using an implicit multi-step Runge-Kutta time-stepping
algorithm, implemented in the Python SciPy package. The initial charge carrier
density was set to $1${\mathrm{m}}^{-3}$$ to treat an ideal, defect-free
semiconductor without intrinsic doping while avoiding a diverging logarithm.
The initial photon densities were set to $0\text{\,}{\mathrm{m}}^{-3}$. The
simulation time was set to $10\text{\,}\mathrm{ns}$, but the results do not
depend on the simulation time in a window of $100\text{\,}\mathrm{ps}$ to
$10\text{\,}\mathrm{ns}$.
Sweeps of the pumping rate $P$ were performed until characteristic relaxation
oscillations of the photon number are observed, which are identified by
comparing the maximum and steady-state values of the photon numbers in the
time trace. If they differed by more than 10%, the corresponding power is
defined as the threshold power, and subsequently divided by the area of a
typical, circular laser spot size of $10\text{\,}\mathrm{\SIUnitSymbolMicro
m}$ to obtain the intensity.
The nanowire material model has an impact on the accuracy of the threshold
calculation, since the chosen refractive index is on the upper end of the
likely values and changes during the build-up of the electron-hole plasma.
First of all, it overestimates the end facet reflectivity, leading to an
underestimation of the thresholds, particularly of the photonic modes. On the
other hand, it also underestimates the Purcell factor of plasmonic/hybrid
modes, leading to an overestimation of the threshold. These effects cause a
correction of the results, but can be incorporated into the model by a
careful, piecewise modeling of the loss channels. Detailed plots showing these
error sources can be found in the supplementary document.

###  2.5 FDTD lasing calculations

Additionally to the model presented in Sec. 2.3, full-wave FDTD calculations
were performed for selected nanowire diameters in order to corroborate the
developed approximative lasing model.
The gain model used was the 4-level-2-electron material with a pumping level
at 350 nm (FWHM 0.8 nm) and a broad lasing level at 380 nm (FWHM 8 nm). The
nanowire was pumped with a Total-Field-Scattered-Field (TFSF) source [19],
which is essentially a uniform source, encompassing the whole nanowire, with
an amplitude of $5\cdot 10^{7}\cdot$\mathrm{V}\text{\,}{\mathrm{m}}^{-1}$$ for
an aluminum and $10^{8}\cdot$\mathrm{V}\text{\,}{\mathrm{m}}^{-1}$$ for a
silver substrate and a wavelength resonant with the pumping transition of the
4-level system. The pump was polarized in the direction of the nanowire axis
and the duration was set to 500 fs for the aluminum and 1000 fs for the silver
substrate.
The nanowire end facets were not extended into the PMLs. The nanowire length
was chosen as 3 µm. A non-uniform mesh was used to decrease computational
requirements because of the large simulation area and the simulation duration.
When the pump is shut off, broadband temporal Fourier-components are excited,
some of which are resonant with the lasing wavelength, seeding the emission of
a pulse. The spatial and spectral field distribution of the emitted pulse is
recorded via a power monitor transversal to the nanowire axis. At the
wavelength showing maximum power transmitted through the monitor, the spatial
field distribution is exported.

##  3 Results and discussion

###  3.1 Lasing of nanowires on a dielectric substrate

The lasing model described above was first applied to a nanowire placed on top
of a silica substrate. In Sec. 3.3., these results were compared with those
obtained considering the nanowire to be placed on metallic substrates,
separated by a thin dielectric spacer.

##### Mode characterization:

The nanowire diameter is varied between $120\text{\,}\mathrm{nm}$ and
$190\text{\,}\mathrm{nm}$ in steps of $10\text{\,}\mathrm{nm}$ to characterize
size-dependent effects. This diameter range was chosen because the photonic
mode is expected to be mainly responsible for lasing in this regime for a
silica substrate, since the modes shown in 2b,c cannot exist because of the
lack of plasmonic confinement. It also includes the range at which the
photonic mode loses support via leakage into free space and thereby is
interesting in terms of mode competition with plasmonic modes, which lack this
cut-off. Since for diameters below $120\text{\,}\mathrm{nm}$, no nanowire-
bound modes could be found, the analysis starts at $120\text{\,}\mathrm{nm}$.
As can be seen in Fig. 3, the overlap with the gain material decreases with
decreasing diameter, as the mode loses confinement at small diameters. The
modal loss rate, on the other hand, increases, as the end facet reflectivity
decreases as a consequence of mode leakage into the environment and free
space.

Figure 3: a) Gain overlap and b) modal loss rate for the photonic mode of a
ZnO nanowire placed on a silica substrate as a function of the nanowire
diameter.

##### Purcell factor:

The results of the Purcell factor calculation as described in Sec. 2.2 can be
found in Fig. 4a. The dipole was placed in the nanowire center, where the
maximum of the photonic mode is located. The Purcell factor of the photonic
mode drops with decreasing diameter as the mode increasingly leaks into the
environment, thereby increasing the mode area.

Figure 4: a) Purcell factor for a dipole decaying at the nanowire center as a
function of the nanowire diameter. b) Threshold intensity as a function of
nanowire diameter.

##### Lasing Threshold calculation:

The threshold intensity, which is plotted in Fig. 4b, depends strongly on the
nanowire diameter. The threshold first decreases with decreasing diameter as a
result of the fact that with smaller nanowires, less material loss needs to be
compensated. The threshold reaches a minimum around 150 nm and increases again
for further decreasing diameters since the photonic mode increasingly leaks
into the surrounding medium.

###  3.2 Dipole position for metallic substrates

Turning towards metallic substrates, aluminum and silver are identified as the
most interesting substrates as they feature plasma frequencies in the UV
spectral range. While the former has a plasma frequency at
$80\text{\,}\mathrm{nm}$, the latter has a plasma frequency of
$330\text{\,}\mathrm{nm}$. One would therefore expect a much stronger
plasmonic enhancement for silver, since its plasma frequency is closer to the
emission wavelength ($\lambda=$380\text{\,}\mathrm{nm}$$) of ZnO.
An increase of the Purcell factor when moving the dipole closer to the metal
is expected, as the interaction with SPPs, which are exponentially confined to
the interface between metal and dielectric, becomes stronger. In other words,
the SPPs add a contribution to the LDOS only accessible for emitters close to
the metal. The Purcell factor was analyzed as a function of position and
orientation for the two metals for a nanowire with a diameter of 200 nm.
Choosing a large diameter has the advantages of small mode leakage into free
space and a large spatial separation between the maxima of the considered
modes. The dipole position is parameterized in terms of the “Dipole offset”,
defined as the distance between dipole and the top of the dielectric spacer.
The requirement of having at least one mesh cell between the dipole box walls
and the emitter as well as having no material change inside the box leads to
the minimum achievable distance between spacer and dipole of
$5\text{\,}\mathrm{nm}$, given a uniform mesh with a 1 nm discretization.
Figure 5 shows an increase of the Purcell factors for a polarization
perpendicular to the metal surface (here $y$) when approaching the aluminum
substrate and for all polarizations on the silver substrate. The
$y$-polarization is dominant with respect to the plasmonic enhancement because
plasmons are mainly polarized in the $y$-direction.
The fact that the Purcell factor increases for all polarization directions on
silver as compared to only one on aluminum might be related to a closer match
to the plasma frequency for silver, leading to a larger degree of confinement
in the metal plane for the silver SPPs.
In the following, when calculating the Purcell factors of the plasmonic as
well as the hybrid mode, the dipole is placed as close as possible to the
metal. To be more precise, the dipole is placed 5 nm above the spacer-nanowire
interface. Within numerical accuracy, this corresponds to the maximum of the
plasmonic and hybrid mode within the nanowire. For calculating the Purcell
factor of the photonic mode, the dipole is still put in the nanowire center.

\begin{overpic}[width=433.62pt]{figures/P_Al_xyz_d200.png} \put(45.0,95.0){a)}
\end{overpic}

\begin{overpic}[width=433.62pt]{figures/P_Ag_xyz_d200.png} \put(45.0,95.0){b)}
\end{overpic}

Figure 5: Purcell factor as a function of dipole position relative to
nanowire-spacer interface (offset) for an a) aluminum and a b) silver
substrate.

###  3.3 Lasing of nanowires on plasmonic substrates

####  3.3.1 Aluminum substrate

##### Mode characterization:

The gain overlap $\Gamma$ and the modal loss rate $1/\tau_{p}$ were analyzed,
as can be seen in Fig. 6 and compared with the case of a silica substrate. The
gain overlap (Fig. 6a) is reduced for the plasmonic mode, but it is still
significant. This is due to the high permittivity of ZnO, which is used for
“weighting” the field distribution (see Eq. 13). The photonic mode on aluminum
shows a slightly higher gain overlap than the one on silica, since leakage
into the substrate is reduced.
For large diameters, the hybrid mode has overlap values between the photonic
and the plasmonic mode, as is expected, given its hybrid character. For
smaller diameters, the hybrid mode cuts-off faster than the photonic mode, and
therefore has a reduced gain overlap.
The modal loss rate (Fig. 6b) of the plasmonic mode and the hybrid mode are
increased strongly with respect to the photonic modes, which is mainly a
result of the non-negligible propagation losses, and decreases with decreasing
diameter. Interestingly, the modal loss rate of the hybrid mode reaches a
value as high as ca. 115 THz at about 180 nm. This is caused by high
propagation losses and a high modal group velocity around this point (see
supplementary document Sec. 2C).

\begin{overpic}[width=433.62pt]{figures/GainOverlap_Al.png}
\put(25.0,85.0){a)} \end{overpic}

\begin{overpic}[width=433.62pt]{figures/Modal_loss_rate_Al.png}
\put(22.0,85.0){b)} \end{overpic}

Figure 6: a) Gain overlap and b) modal loss rate for the three modes of a
nanowire placed on an aluminum substrate and the mode on a silica substrate as
a function of nanowire diameter.

##### Purcell factor:

Another influence of the aluminum substrate is a modification of the Purcell
factor. To characterize the Purcell factors of the three modes, in separate
simulations, one dipole is placed in the center of the nanowire and one at the
bottom. The Purcell factors $F_{x,y,z}$ were calculated and averaged as in
Sec. 3.1. In Fig. 7a, it is shown that the Purcell factor of the photonic mode
is reduced for the aluminum substrate, which supports more modes that compete
for gain.
Interestingly, the hybrid mode as well as the plasmonic mode have an increased
Purcell factor. But while the hybrid mode dominates at large diameters, the
plasmonic mode dominates at small diameters, whereas the hybrid mode cuts off
at small diameters.

Figure 7: a) Purcell factor for a dipole decaying at the nanowire bottom
(plasmonic & hybrid modes) and the nanowire center (photonic modes for both
substrates) as a function of the nanowire diameter. b) Threshold intensity of
the modes under consideration as a function of nanowire diameter.

##### Lasing Threshold calculation:

The threshold intensities of all three modes were calculated and compared with
the threshold intensity on a silica substrate in Fig. 7b. The plasmonic mode
and the hybrid mode on the aluminum substrate have the highest thresholds.
Their threshold depends on the nanowire diameter similarly as the photonic
mode on a silica substrate, as expected from their weak degree of confinement
to the spacer.
The photonic modes have similar thresholds in the small diameter range
independently of the substrate. Remarkably, for larger diameters, the aluminum
substrate leads to a threshold reduction, since leakage into the substrate is
reduced, but effectively no coupling to highly dissipative SPPs is taking
place.
We conclude that on aluminum, the photonic mode dominates the lasing process
under steady-state conditions. The losses of the plasmonic and the hybrid mode
cannot be compensated by the increased Purcell factor.

####  3.3.2 Silver substrate

In contrast to the case of the aluminum substrate, the high losses of silver
in the UV spectral range of interest prevent the possibility of pure plasmonic
lasing with field distributions as shown in Fig. 2b. For example, for a
nanowire diameter of $120\text{\,}\mathrm{nm}$, the photonic mode effective
refractive index has an imaginary part of 0.009 and the hybrid mode 0.3,
whereas the plasmonic mode effective refractive index has an imaginary part of
1.67. Therefore, we only consider the photonic and the hybrid modes. Note,
however, that hybrid modes, as sketched in Fig. 2c, also lead to the
excitation of plasmons, but suffer smaller losses.

##### Mode characterization:

The modes could only be clearly distinguished starting at
$130\text{\,}\mathrm{nm}$. The photonic modes have a similar overlap with the
gain material, as can be seen in Fig. 8a. The hybrid mode, on the other hand,
has an increased overlap with the gain medium caused by the field enhancement
close to the metal. The modal loss rates (Fig. 8b) for both modes are
increased significantly compared to the photonic mode on a silica substrate.
Furthermore, the losses of both modes increase with decreasing diameter
because the mode fraction overlapping with the gain material leaks into free
space.

Figure 8: a) Gain overlap and b) modal loss rate for the two modes of a
nanowire placed on a silver substrate and the mode on a silica substrate as a
function of the nanowire diameter.

##### Purcell factor calculation:

The Purcell factors are much smaller than expected from Fig. 5b, since most of
the high Purcell factor is attributed to the highly damped plasmonic mode. The
Purcell factor of the hybrid mode is small for high diameters, but increases
at lower diameters, as can be seen in Fig. 9a.

Figure 9: a) Purcell factor for a dipole decaying at the nanowire bottom
(hybrid mode) and the nanowire center (photonic modes for both substrates) as
a function of the nanowire diameter for a silver substrate. b) Threshold
intensities of the modes under consideration as a function of nanowire
diameter.

##### Lasing Threshold calculation:

The threshold intensities plotted in Fig. 9b show a large increase overall and
with respect to the nanowire diameter. For almost all but very high diameters,
the hybrid mode will dominate the lasing process, since the increases in gain
overlap and Purcell factor can compensate the higher propagation losses
observed with respect to the photonic mode on aluminum.

####  3.3.3 FDTD calculations of field distribution at the lasing wavelength

As explained in Sec. 2.5, FDTD calculations were performed to corroborate the
findings presented above. The results of these FDTD calculations for two
nanowires with comparable diameters and excitation conditions but two
different substrates are presented in Fig. 10.

Figure 10: Distribution of $|E|^{2}$, normalized to the maximum value, in the
lasing regime for a ZnO nanowire of diameter a) 200 nm on top of an aluminum
substrate and b) 160 nm on top of a silver substrate. A 5 nm silica spacer was
included in both cases. The nanowire outline as well as the spacer are
indicated in white.

Evidently, the field distribution for the case of an aluminum substrate is
concentrated inside the nanowire, whereas for a silver substrate, excitation
of SPPs close to the metal surface can clearly be observed, while still a
strong overlap with the nanowire is present.
The field distributions look different to the ones shown in Fig. 2. Figure 10a
shows three maxima, whereas 2a only shows one. The reason is that our choice
of seeding the lasing modes via shutting of the pump excites a set of modes
that is polarized perpendicular to the waveguide mode in 2a. Nonetheless, the
calculations show that very little interaction with surface plasmons takes
place.
The mode shown in Fig. 10b, on the other hand, more closely resembles 2c
because of the closeness of the lasing wavelength to the plasma frequency. The
main difference lies in the fact that the plasmonic part in 10b shows 2 lobes,
whereas 2a only shows one. This is due to a slight difference in the fitting
of the material properties between both configurations, which has a large
influence on the field distribution.
Together, both calculations show that lasing for a silver substrate is
dominated by hybrid modes, whereas for aluminum, little interaction with
surface plasmons is taking place.

##  4 Conclusions

We have set out to quantify the effect of metallic substrate materials on the
lasing thresholds of ZnO nanowire lasers. In the process, we have derived an
intuitive method to quantify the effect of plasmonically-induced Purcell
factor enhancement on silver and aluminum substrates. We presented detailed
mode calculations and a rate-equation model to study mode competition.
Considering silver, a plasmonic material with a plasma frequency being
resonant with the bandgap of ZnO, we have shown that plasmonic lasing is a
double-edged sword: On the one hand, the low mode volume of SPP-modes offers
interesting possibilities with respect to the strong increase of the Purcell
factor. On the other hand, the high propagation losses make an excitation of
these modes in a steady-state or quasi-steady-state regime very difficult.
Note that these results offer a new perspective for the interpretation of
previously published results [5]. A careful analysis of the mode profiles
shows that the modes found to be responsible for lasing by Oulton et. al. on
the plasmonic substrate are in fact the “hybrid” modes in our terminology. The
order-of-magnitude increase below a 100 nm nanowire diameter fits nicely with
the cut-off expected for nanowire-bound modes (at around 120 nm diameter, the
photonic modes cannot be found in our solver) and the losses of the purely
plasmonic modes (see introductory text for Sec. 3.3.2).
As expected, silver shows a strong effect on lasing properties and even offers
the possibility to access hybrid modes in the steady-state by compensating
increased losses by increased Purcell factors and gain overlap. Aluminum, on
the other hand, does not offer these advantages. But curiously, an aluminum
substrate decreases the threshold for photonic lasing as compared to a
photonic substrate, since the small penetration depth suppresses mode leakage
into the substrate.
These findings should also be valid for other combinations of nanowire
materials and metals. We expect for metals with plasma frequenciess far below
the nanowire emission wavelength $\lambda$ and with diameters comparable to
$\lambda/n$, where $n$ is the nanowire refractive index at the emission
wavelength, a threshold reduction without plasmon excitation. This is caused
by a reduction of leakage into the substrate or, in other words, a small
penetration depth into the metal. On the other hand, for resonant
combinations, we expect a threshold increase caused by plasmon excitation via
Purcell enhancement.
Another interesting finding is related to the polarization properties of
emitters close to the substrate. Highly-resonant combinations of nanowire
lasing wavelength and plasma frequency seem to offer an increase of the
emission rate independently of the emitter polarization, whereas in an off-
resonant configuration, mainly emitters polarized normal to the metal
experience enhancement.

## Acknowledgement

We acknowledge the financial support via the project ID 398816777 in the CRC
1375 “NOA–Nonlinear optics down to atomic scales” funded by the Deutsche
Forschungsgemeinschaft (DFG).
A.B. thanks the financial support of the German Research Foundation (DFG)
under the framework of the International Research Training Group (IRTG) 2675
“META-ACTIVE”.

## References

  * [1] Robert Röder et al.  “Mode Switching and Filtering in Nanowire Lasers”  In _Nano Letters_ 16.4 American Chemical Society, 2016, pp. 2878–2884  DOI: 10.1021/acs.nanolett.6b00811
  * [2] Pai-Chun Chang et al.  “High-performance ZnO nanowire field effect transistors”  In _Applied Physics Letters_ 89.13, 2006, pp. 133113  DOI: 10.1063/1.2357013
  * [3] Jaeyun Moon et al.  “Gate-Modulated Thermoelectric Power Factor of Hole Gas in Ge-Si Core-Shell Nanowires”  In _Nano Letters_ 13.3, 2013, pp. 1196–1202  DOI: 10.1021/nl304619u
  * [4] Shaimaa I. Azzam et al.  “Ten years of spasers and plasmonic nanolasers”  In _Light: Science & Applications_ 9.1, 2020, pp. 90  DOI: 10.1038/s41377-020-0319-7
  * [5] Rupert F. Oulton et al.  “Plasmon lasers at deep subwavelength scale”  In _Nature_ 461.7264, 2009, pp. 629–632  DOI: 10.1038/nature08364
  * [6] Daniël Vanmaekelbergh and Lambert K. Vugt  “ZnO nanowire lasers”  In _Nanoscale_ 3.7 The Royal Society of Chemistry, 2011, pp. 2783–2800  DOI: 10.1039/C1NR00013F
  * [7] Jeremy B. Wright and Igal Brener  “Distributed feedback gallium nitride nanowire lasers”  In _Applied Physics Letters_ 104.4, 2014, pp. 041107  DOI: 10.1063/1.4862193
  * [8] Koichi Okamoto et al.  “Surface-plasmon-enhanced light emitters based on InGaN quantum wells”  In _Nature Materials_ 3.9, 2004, pp. 601–605  DOI: 10.1038/nmat1198
  * [9] Na Gao et al.  “Surface-plasmon-enhanced deep-UV light emitting diodes based on AlGaN multi-quantum wells”  In _Scientific Reports_ 2.1, 2012, pp. 816  DOI: 10.1038/srep00816
  * [10] David J. Bergman and Mark I. Stockman  “Surface plasmon amplification by stimulated emission of radiation: quantum generation of coherent surface plasmons in nanosystems.”  In _Physical review letters_ 90.2, 2003, pp. 027402
  * [11] Davy Gérard and Stephen K Gray  “Aluminium plasmonics”  In _Journal of Physics D: Applied Physics_ 48.18 IOP Publishing, 2014, pp. 184001  DOI: 10.1088/0022-3727/48/18/184001
  * [12] Ying Gu et al.  “Surface-Plasmon-Induced Modification on the Spontaneous Emission Spectrum via Subwavelength-Confined Anisotropic Purcell Factor”  In _Nano Letters_ 12.5, 2012, pp. 2488–2493
  * [13] Sergey V. Gaponenko  “Introduction to Nanophotonics”  Cambridge University Press, 2010
  * [14] Qing-Hua Ren et al.  “Surface-plasmon mediated photoluminescence enhancement of Pt-coated ZnO nanowires by inserting an atomic-layer-deposited Al2O3spacer layer”  In _Nanotechnology_ 27.16 IOP Publishing, 2016, pp. 165705  DOI: 10.1088/0957-4484/27/16/165705
  * [15] Christoph Walther, Giacomo Scalari, Mattias Beck and J. Faist  “Purcell effect in the inductor-capacitor laser”  In _Opt. Lett._ 36.14 OSA, 2011, pp. 2623–2625  DOI: 10.1364/OL.36.002623
  * [16] Hatice Altug, Dirk Englund and Jelena Vučković  “Ultrafast photonic crystal nanocavity laser”  In _Nature Physics_ 2.7, 2006, pp. 484–488  DOI: 10.1038/nphys343
  * [17] Mariam M. Tohari, Andreas Lyras and Mohamad S.  “A Novel Metal Nanoparticles-Graphene Nanodisks-Quantum Dots Hybrid-System-Based Spaser”  In _Nanomaterials_ 10.3, 2020  DOI: 10.3390/nano10030416
  * [18] M.. Noginov et al.  “Stimulated Emission of Surface Plasmon Polaritons”  In _Phys. Rev. Lett._ 101.22 American Physical Society, 2008, pp. 226806  DOI: 10.1103/PhysRevLett.101.226806
  * [19] “Lumerical FDTD Solutions” Lumerical Version 2019aR2, www.lumerical.com
  * [20] Lin Zhu  “Modal Properties of Hybrid Plasmonic Waveguides for Nanolaser Applications”  In _IEEE Photonics Technology Letters_ 22.8, 2010, pp. 535–537  DOI: 10.1109/LPT.2010.2041923
  * [21] “Integrating the Poynting Vector” Version 2019aR2, https://support.lumerical.com/hc/en-us/articles/1500007397442-Integrating-the-Poynting-Vector
  * [22] Tobias D. Schmidt et al.  “Non-isotropic emitter orientation and its implications for efficiency analysis of organic light-emitting diodes”  In _Organic Photonics V_ 8435 SPIE, 2012, pp. 843513  DOI: 10.1117/12.921537
  * [23] Hisashi Yoshikawa and Sadao Adachi  “Optical Constants of ZnO”  In _Japanese Journal of Applied Physics_ 36.Part 1, No. 10 IOP Publishing, 1997, pp. 6237–6243  DOI: 10.1143/jjap.36.6237
  * [24] Marijn A.. Versteegh  “Ultrafast screening and carrier dynamics in ZnO: Theory and experiment”  In _Physical review. B, Condensed matter and materials physics_ 84.3 Published by the American Physical Society through the American Institute of Physics, 2011
  * [25] Marijn A.M. Versteegh, Daniël Vanmaekelbergh and Jaap I. Dijkhuis  “Lasing in ZnO nanowires is electron-hole plasma lasing”  In _CLEO: 2011 - Laser Science to Photonic Applications_ , 2011, pp. 1–2  DOI: 10.1364/QELS.2011.QWH3
  * [26] Marcel Wille et al.  “Carrier density driven lasing dynamics in ZnO nanowires”  In _Nanotechnology_ 27.22 IOP Publishing, 2016, pp. 225702  DOI: 10.1088/0957-4484/27/22/225702
  * [27] Esteban Bermúdez-Ureña et al.  “Plasmonic Waveguide-Integrated Nanowire Laser”  In _Nano Letters_ 17.2 American Chemical Society, 2017, pp. 747–754
  * [28] J.. Reparaz et al.  “Spatial mapping of exciton lifetimes in single ZnO nanowires”  In _APL Materials_ 1.1, 2013, pp. 012103  DOI: 10.1063/1.4808441
  * [29] Sangsu Hong et al.  “Time-resolved photoluminescence of the size-controlled ZnO nanorods”  In _Applied Physics Letters_ 83.20, 2003, pp. 4157–4159  DOI: 10.1063/1.1627472
  * [30] Daniel J. Gargas, Hanwei Gao, Hungta Wang and Peidong Yang  “High Quantum Efficiency of Band-Edge Emission from ZnO Nanowires”  In _Nano Letters_ 11.9, 2011, pp. 3792–3796  DOI: 10.1021/nl201850k
  * [31] A.. Koenderink  “On the use of Purcell factors for plasmon antennas”  In _OPTICS LETTERS_ 35.24, 2010, pp. 4208–4210  DOI: 10.1364/OL.35.004208
  * [32] Jelena Vuckovic et al.  “FDTD calculation of the spontaneous emission coupling factor in optical microcavities”  In _Micro- and Nano-photonic Materials and Devices_ 3937 SPIE, 2000, pp. 2–11  International Society for OpticsPhotonics  DOI: 10.1117/12.382802
  * [33] K. Bohnert, G. Schmieder and C. Klingshirn  “Gain and Reflection Spectroscopy and the Present Understanding of the Electron–Hole Plasma in II–VI Compounds”  In _physica status solidi (b)_ 98.1, 1980, pp. 175–188  DOI: https://doi.org/10.1002/pssb.2220980117