of 12
PHYSICAL REVIEW E
111
, 015505 (2025)
Controlling the roughening of growing electrochemical interfaces using temperature gradients
Asghar Aryanfar
,
1
,
*
Ali Tayyar,
2
and William A. Goddard, III
3
1
Bo ̆
gaziçi University
, Mechanical Engineering, Bebek, Istanbul, Turkey 34342
2
American University of Beirut
, Mechanical Engineering, Riad El Solh, Beirut, Lebanon 1107 2020
3
California Institute of Technology
, 1200 E California Blvd, Pasadena, California 91125, USA
(Received 3 November 2024; accepted 2 January 2025; published 24 January 2025)
The excessive dendritic development during the electrochemical evolution of the microstructures in recharge-
able batteries can ultimately cause a short circuit, thermal instability, or runaway, and loss of active material.
We initially develop a computational framework to quantify the bias of the electrodeposition on the roughened
interface favoring the convex zones. Subsequently, we impose a countering temperature effect to enhance the
diffusion on the trailing concave zones. Consequently, we establish a stability criterion for controlling surface
roughening where the visualized space of parameters establishes a relationship between the geometry of the
interface, the physical properties of the electrolyte, and the charging conditions. The developed framework could
be useful for controlling the propagation of the microstructures and the prevention of runaway, during prolonged
cycles, particularly when the surface roughness gets pronounced in the later stage of cycle life.
DOI:
10.1103/PhysRevE.111.015505
I. INTRODUCTION
The increasing necessity for energy in past years has re-
quired the innovation of higher-energy density energy storage
systems with enhanced performance [
1
]. Among many tech-
nological innovations, batteries have shown great potential as
a clean source of electrical energy. More recently, lithium-ion
batteries (LIBs) have been widely used in portable electronics,
manufacturing, and the electric cars and scooters industry
[
2
,
3
] as well as the green energy sector [
4
,
5
], which exceeds
billions of units manufactured throughout the world [
6
].
In a more specific way, LIBs have been widely investigated
due to lithium’s extremely high theoretical specific capac-
ity (3860 mAhg
1
), low density (0
.
59 gcm
3
), the lowest
negative electrochemical potential (
3
.
04
V
vs the standard
hydrogen electrode) [
7
], and design flexibility [
8
]. In this
regard, the growth of dendrites, a treelike microstructure in-
duced by the expansion of conductive filaments on the lithium
electrode has been a critical issue since they grow with ex-
treme porosity, leading to excessive reach and causing short
circuits [
9
,
10
]. As well, it may detach from the thinner part of
the branches during the dissolution (i.e., discharge), and create
the so-called dead lithium [
11
], causing thermal instability
[
12
,
13
] and capacity decay [
14
,
15
].
The precise process of dendritic growth remains not en-
tirely grasped, yet it is generally recognized to be impacted
by several factors such as current density [
16
18
], electrolyte
composition [
12
,
19
,
20
], solid electrolyte interphase (SEI)
[
21
,
22
], and surface defects (i.e., kinks) [
23
]. The physical
characterizations comprise porosity [
24
], tortuosity [
25
], and
MacMullin number [
26
], the geometry and microscale cur-
vature [
27
,
28
], and the sand time effect on growth transition
[
29
].
*
Contact author: aryanfar@caltech.edu
The larger scale consequence of dendritic development can
be attributed to a roughened interface where the developed
models for stabilizing such an electrodeposited surface could
be categorized into two realms of continuum and atomistics.
From a continuum perspective, a small enough wavelength of
perturbations in overlimiting currents leads to stable growth
[
27
], imposing prestress can increase the stability span of
the undulations [
30
], using viscous electrolyte can level the
growing interface via limiting the electroconvection [
31
,
32
],
and simultaneous evolution of curved interface with the SEI
can alter the interface geometry due to thickness dependent
variation in the resistance of SEI, below its critical value [
33
].
As well, the artificial design of the strong solid-electrolyte
interface can suppress the formation of roughened surfaces
(i.e., dendrites) [
34
].
From an atomistic perspective, the convex electrode curva-
ture (bent into the electrolyte) has shown a positive effect on
stabilizing the growing morphology [
28
], curvature-induced
local straining can favorably alter lithophilicity on an elec-
trode surface [
35
], and formation of small lithium aqueous
and nonaqueous electrolytes could provide additional elec-
trode or electrolyte passivation and control [
36
]. As well,
the fractal dimension of the microstructures has been shown
to be curvature-independent [
37
]. Other works on the elec-
trochemical interfaces include phase-field modeling [
38
40
],
coarse-grained modeling in time [
41
], sharp interface mod-
eling via an extending space-charge region, and considering
nonlinear electrode reactions [
42
].
Several engineering ideas have tried to control the dendritic
propagation, including shielding techniques [
43
], curvature-
induced structural modification [
44
], imposing external mag-
netic fields [
45
], guiding scaffolds [
46
], using nanofiber arrays
in the polymer electrolyte [
47
], applying prestress to stabilize
planar growth and reduce surface roughening, [
30
], and using
an-isotropic electrolytes to controlling electric-field gradients
offering potential strategies to improve battery stability and
2470-0045/2025/111(1)/015505(12)
015505-1
©2025 American Physical Society
ARYANFAR, TAYYAR, AND GODDARD III
PHYSICAL REVIEW E
111
, 015505 (2025)
FIG. 1. The convex and concave regions in the perturbed interface are compared in terms of (a) electric field and (b) growth stability;
(a) Surface heterogeneity model: Black dashed arrows illustrate the average electric field for the equivalent flat geometry.
h
,
l
: domain height
and length, ̄
y
: average elevation of interface,
x
,
y
: coordinates. (b) The illustration of evolution dynamics, where the reference interface (i.e.,
black) can either grow stable (i.e., blue) or unstable (i.e., red). In stable growth, the convex (top) and concave (bottom) regions get closer to
each other (blue arrows), while during unstable growth they get further away from each other (red arrows).
potentially suppressing dendrite formation [
48
]. Additionally,
they have ascribed that the growth rate could be controlled
when lowering the charge transfer coefficients [
49
], applying
forced advection and cross flow [
50
,
51
]. The higher surface
roughness and thinner SEI layer could exacerbate the den-
dritic propagation [
52
54
] and describe the kinetics of rough-
ening interfaces via incorporating surface tension and elastic
and viscous effects [
55
]. As well, for the engineered undu-
lations, superfills have been used as an additive chemical to
enable preferential metal deposition on the bottom regions and
avoiding the formation of cavities [
56
], and their efficacy of
correlates directly with the curvature of the undulations [
57
].
Another factor for the electrodeposition, which has long
been studied in piezoelectric materials [
58
,
59
] and can alter
the fate of the formed interface, is mechanical factors such
as surface tension [
60
,
61
] and imposing external pressure
and stress [
62
], as an additional term in the electrochemical
potential [
63
,
64
].
In the meantime, the temperature is a key parameter, the
effect of which has proven to be instrumental in the experi-
mental surface film growth [
65
], as well as in the molecular
dynamics simulations of Li
+
diffusional transport in the SEI
(i.e., 250–400
K
)[
66
]. More specifically, temperature has been
explored to shape the surface morphology and dendrite for-
mation of microstructures [
67
69
]. It has been shown that
elevated temperatures (i.e., 60
80
) transform the SEI layer
into a uniform compact layer which can delay the instigation
of dendrites [
70
,
71
]. However, later, lowering the system tem-
perature increases the diffusion resistance and decreases the
surface film thickness, which is favorable for suppressing the
dendrites [
72
,
73
].
While the presence of surface roughness, consisting of the
convex (i.e., peak) and concave (i.e., valley) regions, favors
the electrodeposition in the convex zones, in this paper we
counter such effects via imposing a temperature field in the
opposite direction, which forms a bias for favoring the elec-
trodeposition in the concave (i.e., valleys) sites. Obtaining
the voltage
V
and temperature field
T
, we characterize the
space of physical parameters for the efficacy of the imposed
temperature gradient. The obtained framework could be useful
for controlling the growing dendrites in electrochemical sys-
tems, particularly those consuming larger power, where the
microstructures could grow with runaway behavior.
II. METHODOLOGY
From the atomic perspective, dendritic propagation and
formation of rough interfaces are inevitable, and fast pro-
cesses occur under nonequilibrium conditions. Such a
kinetics-dominated event rapidly accumulates the ions with
random walks (i.e., Brownian dynamics), which is distin-
guishable from typical thermodynamics-dominated bonding
under slow, close-to-equilibrium conditions.
From the continuum perspective, the development of the
interfacial morphology depends on the macroscale electro-
chemical parameters which include the electric field and
electrolyte diffusivity. Such variables control the accessibility
of ions to the electrodeposition sites and the rate of interface
growth.
As shown in Fig.
1(a)
, such stochastic interactions should
yield a fast-growing asperity (i.e., convex:

), followed by a
slow-growing hollow region (i.e., concave:

). The difference
in the growth rate is due to the peaks having better proximity
for the upcoming ions as well as possessing geometry-induced
higher directed electric field (
E

>
E

). Hence, the evolved
microstructures will acquire a certain roughness, which grows
with runaway behavior [
29
,
74
]. Therefore, we aim to counter
such instability by imposing a reaction-rate bias in favor of
015505-2
CONTROLLING THE ROUGHENING OF GROWING ...
PHYSICAL REVIEW E
111
, 015505 (2025)
the concave region, which is controlled by the temperature
field
T
.
A. Electric field
The perturbation in the growing interface can be modeled
in a sinusoidal form, where concave (i.e., valley:

) and
convex (i.e., peak:

) regions are connected as shown in
Fig.
1(b)
. The elevation of the interfacial roughness
y
I
could
be expressed in terms of the horizontal variable
x
as
y
I
=
̄
y
+
A
sin

2
πω
x
l

,
(1)
where ̄
y
is the average elevation of the growing mi-
crostructure,
A
is the amplitude of the perturbation,
ω
is the
frequency of the repetition, and
l
is the domain length. Thus,
the maximum curvature
κ
occurs in the peaks and valleys
(
|
sin(
...
)
|=
1
x
l
/
4
ω
)as
κ
=
y

I


1
+
y

2
I

3





x
=
l
4
ω
=−
A

2
πω
l

2
.
(2)
Therefore, one can calculate the ratio of the maximum-
to-minimum electric field
ˆ
E
:
=
E

/
E

, which is a measure
of how critically the roughening interface. In this regard, the
maximum and minimum electric field depend on both the
radius of curvature
r
(i.e.,
r
=
1
), the applied voltage of
V
+
and
V
between two electrodes as well as the remaining dis-
tance of the interface from the counter electrode (i.e.,
h
y
I
).
Thus, the electric field in any position can be described as
E
=
f
(
r
)
V
+
V
h
y
I
,
(3)
where
f
(
r
) is the curvature dependents augmentation (or
diminishing) factor based on the interface geometry. Partic-
ularly, the ratio of the electric field
ˆ
E
between the convex
E

and concave
E

surfaces (i.e.,
ˆ
E
=
E

/
E

) is obtained as
ˆ
E
=
f
max
(
ˆ
r
)
f
min
(
ˆ
r
)
h
y
min
h
y
max
,
where
y
min
and
y
max
are the respective elevations at concave
and convex surfaces and ˆ
r
:
=
r
/
h
is the normalized radius of
curvature. In fact, for the convex interface, the electric field
is augmented [
f
max
r
)
>
1], and vice versa for the concave
interface it gets diminished [
f
min
r
)
<
1]. Hence, in terms of
the curvature effect, one gets the following order:
E

>
̄
E
>
E

,
(4)
where the indexes


,
̄

, and


represent the convex (i.e.,
peak), flat, and concave (i.e., valley) regions, respectively.
B. Temperature gradient
The temperature has a deterministic role in the rate of the
electrodepositing reactions, which is generally expressed via
the Arrhenius relationship. The diffusivity of the ions
D
within
the electrolyte is governed by the following relationship [
75
]:
D
=
D
0
exp

Q
k
B
T

,
(5)
FIG. 2. The temperature map for the roughened surface
is characterized by convex (
peaks
) and concave (
valleys
)re-
gions. Lines show the dimension-less temperature map [i.e.,
ˆ
T
=
(
T
T
+
)
/
(
T
T
+
)], with the bottom (i.e.,
T
) and top (i.e.,
T
+
) electrode temperatures. The growing interface remains isother-
mal with the connected electrode (i.e.,
T
). The temperature
distribution is obtained later in Sec.
III
.
where
D
0
is a prefactor and,
Q
is the activation energy for the
diffusivity, and
k
B
is the Boltzmann’s constant. Figure
2
illus-
trates a sample of a formed temperature field in the vicinity of
the convex (blue circle) and concave (red) sites when a higher
temperature is imposed in the negative electrode (
T
>
T
+
).
Since the concave region is surrounded by a higher temper-
ature interface, it will get a higher temperature. Conversely,
the sticking-out convex peak will have less of the higher
temperature interface in its surroundings and hence gets lower
temperature values. Consequently, the concave regions pos-
sess warmer electrolytes than their convex counterpart. In this
context, while the electric field forms a bias for more growth
in the convex region, imposing the higher temperature at the
negative electrode could cause a countering bias for more
electrodeposition in the concave sites. This could be explained
in terms of the ionic flux
J
, which is described as [
76
]:
1
J
=−
D
C
zF
RT
DC
V
(6)
where
D
is the diffusivity of the ions in the electrolytic
medium,
C
is the local concentration of the electrolyte,
z
is
the valence number of the charge carriers,
F
is the Faraday’s
constant,
R
is the gas constant,
T
is the local temperature, and
V
is the local voltage. Therefore, while the electric field
E
,
as a nonmaterial-independent property, favors higher flux in
the convex interface, the diffusivity as a material-dependent
property could favor higher ionic flux at the concave sites and
balance out the electrodeposition for forming more uniform
interface, and one gets
D

>
̄
D
>
D

,
(7)
which tends to counter the the electric field effect obtained
earlier in Eq. (
4
).
1
In the absence of convection.
015505-3
ARYANFAR, TAYYAR, AND GODDARD III
PHYSICAL REVIEW E
111
, 015505 (2025)
TABLE I. Physical parameters.
Variable
DC
0
lR TV
0
zF
Value
2
.
3
×
10
10
10
3
25
8.3
298
3.6
1
96.4
Unit
m
2
s
1
mol m
3
μm
J mol
1
K
1
KV
[]
kC
mol
1
Ref.
[
77
]
Assumed
[
78
][
79
]
Ambient
[
80
]forLi
+
[
81
][
79
]
C. Stability condition
The comparative scale of the convex and concave zones,
illustrated in Fig.
2
, falls in the continuum realms, which
extends significantly larger than the smaller-scale electrode-
position sites in the interface (i.e., Helmholtz layer, double
layer, space-charge, etc.). Such a large-scale zone is in the
electroneutral condition [
82
], where the continuum relation-
ship for ionic flux is in effect. In fact, the ionic flux [Eq. (
6
)]
is composed of mainly diffusion (
D
C
) and electromigration
[(
zF
/
RT
)
DC
V
]. Herein, we compare these two terms using
the typical value in Table
I
, as below:
J
diff
=
D
C
2
DC
0
l
10
2
mol
m
2
s
diffusion
J
mig
=
zF
RT
DC
V
De
k
B
T
zFC
0
V
0
l
1
mol
m
2
s
electromigration,
which means that for the typical value
J
m

J
diff
; not to men-
tion that for the voltage
V
, the effective zone of variation is
the space-charge region (i.e.,
x
I
)[
82
], which would make
J
mig
even larger. In fact, the migration efficacy for ionic flux occurs
when the applied voltage is high enough. In such conditions,
the drift velocity
v
d
of the charge carriers will be effective,
which feeds the electrodeposition sites and determines the rate
of interface growth, and is defined as [
81
]
v
d
=
μ
E
,
(8)
where
μ
and
E
are the local ionic mobility and electric
field, respectively. As well, considering the correlation of
the electrical mobility
μ
with the diffusivity
D
(i.e., Einstein
relationship where
μ
D
/
T
), one would require faster elec-
trodeposition in the concave zone and the convex counterpart
for achieving stability (i.e.,
v
d
,
>
v
d
,
). Hence, assuming
similar grown porosity in the convex and concave zones, one
gets
D

E

T

>
D

E

T

,
(9)
where considering the correlation of the diffusivity
D
with the
temperature
T
[Eq. (
5
)], one achieves
T

T

exp

Q
k
B

1
T

1
T


>
E

E

,
(10)
and the required activation energy of the electrolyte
Q
, such
that the concave zone (

) could catch up with the convex zone
(

), is obtained as
Q
>
T

T

T

T

k
B
ln

E

E

T

T


,
(11)
which provides the stability criterion for the stable evolution
of the interfacial perturbation when the temperature gradient
is imposed.
III. NUMERICAL COMPUTATION
A. Voltage distribution
V
To calculate the maximum
E

and minimum
E

electric
fields, one needs to compute the voltage distribution in the
vicinity of these sites. In this regard, Poisson’s equation could
describe the continuum voltage mapping, given as [
83
]
2
V
=
ρ

,
(12)
where
ρ
is the local charge density and

is the permittivity of
the electrolyte medium. Due to the large scale of the interface
perturbation (
μm), a fair assumption is that a significant
portion of the medium falls out of the double layer region
(
nm). Therefore, a simplified version of Eq. (
12
)forthe
voltage distribution in 2D would be [
84
]
2
V
x
2
+
2
V
y
2
0
.
(13)
Therefore, the domain could get partitioned into a 2D grid
of
n
x
×
n
y
points where
n
x
and
n
y
are the number of grid
points in
x
and
y
directions, respectively, which are selected
large enough to be able to capture the voltage variations in the
convex (

) and concave (

) sites.
Hence,
V
i
,
j
represents the voltage at nodes
x
j
and
y
i
. Per-
forming a finite difference scheme on Eq.
13
, one has

V
i
,
j
+
1
2
V
i
,
j
+
V
i
,
j
1

δ
x
2
+
V
i
+
1
,
j
2
V
i
,
j
+
V
i
1
,
j
δ
y
2
0
.
(14)
Rearrangement gives
V
i
,
j
V
i
+
1
,
j
+
V
i
1
,
j
+
δ
y
2
δ
x
2
(
V
i
,
j
+
1
+
V
i
,
j
1
)
2

1
+
δ
y
2
δ
x
2

.
(15)
Regarding the boundary conditions, from Fig.
1(a)
, since
the metallic interface is physically connected to the negative
015505-4
CONTROLLING THE ROUGHENING OF GROWING ...
PHYSICAL REVIEW E
111
, 015505 (2025)
FIG. 3. The distribution in the normalized voltage
ˆ
V
for various normalized curvatures ˆ
κ
. Gray area represents the growing (i.e., dendritic)
interface.
pole, one has
V
(
x
,
0
)
=
V
y
=
0
V
=
V
y

y
I
V
(
x
,
h
)
=
V
+
y
=
h
,
(16)
which are the Dirichlet boundary conditions for the dis-
cretized Eq. (
15
), and
V
and
V
+
are the voltage values for
the negative electrode and positive electrode.
The initial voltage value
V
0
is set as the uniform dis-
tribution in the domain between the negative and positive
electrode, except the interface which is physically touching
the negative electrode (
V
)as
V
0
(
x
,
y
)
=
V
y

y
I
V
+

V
+
V
h

yy
>
y
I
.
(17)
Subsequently, the voltage of the grid points
V
i
,
j
was up-
dated from the neighboring grid points (top
, bottom
,
left
, and right
), through the successive steps based on
Eq. (
15
). During each iteration, the boundary conditions, in
Eq. (
16
) were enforced from the electrodes and the interfaces.
As well, regarding the left and right boundaries, the periodic
boundary condition (PBC)was applied, which indicates that
the voltage values between the left and right boundaries are
transferred to each other.
Throughout the successive runs, the convergence criterion
is defined by a maximum difference of the voltage val-
ues
V
i
,
j
varying of next iteration (
k
+
1) from the previous
one (
k
), which is tracked as the error. The iterations were
stopped when the error value got smaller (or equal) to the
assigned threshold value. Figure
3
visualizes the converged
final voltage distribution for both low- (i.e., small
ω
) and
high-curvature (i.e., large
ω
) perturbations [
85
].
Finally, the average electric field across the highlighted
regions of the convex (i.e.,
E

) and concave (i.e.,
E

) regions
in Fig.
2
are calculated via approximating
E
=−∇
V
as
E
V
top
V
bot
δ
,
(18)
where
V
top
and
V
bot
are the converged voltage values in the top
and bottom ends of the focused zones in both convex (

) and
concave (

) zones and
δ
is their respective diameter.
B. Temperature distribution
T
The relationship for transition and dispersal of the temper-
ature is expressed as [
86
]
T
t
=
α

2
T
x
2
+
2
T
y
2

,
(19)
since the rate of growth in the moving boundary
y
I
/∂
t
is far less than the dynamics of the convergence of the tem-
perature (
T
/∂
t
)[(
y
I
/∂
t
)
(
T
/∂
t
)]; during the growth
process, the boundary could get assumed as quasistationary
with respect to the transition of the temperature field, which
means that the temperature profile can be treated as the qua-
sisteady state
T
/∂
t
0, and the temperature profile is given
via
2
T
x
2
+
2
T
y
2
0
(20)
regarding the boundary conditions; the temperature values
are constant in the electrodes after applying the temperature
gradient (
T
>
T
+
). As well, since the thermal diffusivity of
the growing metallic interface
α
I
is significantly larger than
the nonmetallic metallic boundary
α
El
(
α
I

α
El
), the rate of
temperature variation will also be significantly larger and it
remains relatively isothermal with the electrode (i.e.,
T
). This
means that the interface temperature could be assumed as the
quasisteady state compared to the electrolytic domain, and the
015505-5
ARYANFAR, TAYYAR, AND GODDARD III
PHYSICAL REVIEW E
111
, 015505 (2025)
FIG. 4. The normalized temperature profile
ˆ
T
versus the normalized vertical distance ˆ
y
ˆ
y
I
, above (a) the convex (i.e., peak) and (b) the
concave (i.e., valley) regions for various normalized curvature values ˆ
κ
.
temperature boundary condition will be established as
T
(
x
,
0
)
=
T
T
=
T
y

y
I
T
(
x
,
h
)
=
T
+
.
(21)
It is noticeable that the relationships and boundary con-
ditions for the temperature in Eqs. (
20
) and (
21
) are similar
to the relationships for the voltage distribution in Eqs. (
13
)
and (
16
). Hence, one can avoid repeated computations by
projecting the temperature profile
T
from the voltage profile
V
. In this regard, we define the dimensionless parameters for
the voltage
ˆ
V
and temperature
ˆ
T
as
ˆ
V
=
V
V
V
+
V
0

ˆ
V

1
ˆ
T
=
T
T
+
T
T
+
0

ˆ
T

1
.
Noting that the gradient of the temperature field is imposed
as opposite to the voltage field, the temperature field is ob-
tained as
ˆ
T
=
1
ˆ
V
.
(22)
Figures
4(a)
and
4(b)
illustrate the variation of the temper-
ature versus the normalized vertical distance above the convex
(i.e., peak) and concave (i.e., valley) sites (ˆ
y
ˆ
y
I
), respec-
tively. Here, one can infer the effect of the curved boundary
on the temperature profile as
2
T
ˆ
y
2






<
0
,
2
T
ˆ
y
2






>
0
,
which means that the average temperature in the convex
̄
T

and concave
̄
T

sites (blue and red circles in Fig.
2
)differ
considerably, which can be obtained via integrating in the
region as
̄
T
=
1
ˆ
δ
ˆ
δ
0
Tdy
,
and one gets
̄
T

>
̄
T

,
(23)
which are shown as dashed lines in these two figures.
IV. RESULTS AND DISCUSSIONS
A. Stability analysis
The temperature difference, given in Eq. (
23
), provides
the possibility for higher diffusion of the ions in the concave
(

) sites relative to the convex (

) counterpart regions. One
effective parameter for the imposed temperature gradient
T
is the sensitivity of the resulted diffusivity
D
, which can be
obtained as
D
T
=−
Q
k
B
T
D
.
Here it is evident that one effective parameter for such
sensitivity is the activation energy
Q
. In other words, for an
electrolyte more responsive to the temperature (i.e.,
Q
),
less temperature is needed (
T
) to form stable growth.
Considering the halfway progress of the interface ( ̄
y
=
h
/
2),
Figures
5(a)–5(c)
visualize the range of required tempera-
ture gradient
T
versus the diffusion activation energy
Q
of the electrolyte for three quartiles of destabilization
ˆ
A
=
{
1
/
4
,
1
/
2
,
3
/
4
}
h
/
2, which is additionally characterized in
terms of interface curvature ˆ
κ
=
κ/
h
[
85
].
In this regard, the earlier slower growth of the stable zone
(i.e.,
ˆ
A
={
1
/
4
1
/
2
}
h
/
2) and its later fast diminishing (i.e.,
ˆ
A
={
1
/
2
3
/
4
}
h
/
2) could be explained via the nonlinear
overall effect of the electric field
E
versus the temperature
field
T
[Eq. (
10
)]. Labeling the temperature and electric field
roles as
f
T
:
=
T

/
T

exp(
Q
/
k
B
(1
/
T

1
/
T

)) and
f
E
:
=
E

/
E

, respectively, the stable growth could get obtained
when
f
T
>
f
E
. Therefore, during the earlier transition from
015505-6
CONTROLLING THE ROUGHENING OF GROWING ...
PHYSICAL REVIEW E
111
, 015505 (2025)
FIG. 5. Visualized zones of the required temperature difference
T
versus the activation energy
Q
, for halfway progress of the interface
(i.e., ̄
y
=
h
/
2) for three-quarters of the destabilized growth (
ˆ
A
={
1
/
4
,
1
/
2
,
3
/
4
}
). The stable zones (blue) signify diffusion domination, and
vice versa, the unstable zones (red) illustrate electromigration domination.
the first to second quartile, it is obvious from Eq. (
11
) that
while the competitive ratios of
E

/
E

and
T

/
T

carry
smaller weight in the logarithm term, the primary impact is
derived from the temperature difference as
Q
>
1
T

T

.
In this regard, since the dimensionless voltage
ˆ
V
and tem-
perature
ˆ
T
fields have been found to be complementary of
each other in Eq. (
22
), one could ascribe the following corre-
lation:
ˆ
T
=
1
δ
0
ˆ
Edy
,
which means that the temperature field would be averaging
the values across the focus region of scale
δ
(Fig.
2
). In this
regard, while considering solely the temperature ratio
T

/
T

between the concave and convex sites would carry a smaller
effect than the electric field ratio
ˆ
E
, the exponential term
exp(
Q
/
k
B
(1
/
T

1
/
T

)) empowers overall temperature ef-
fect
f
T
, which makes it ultimately dominate.
However, during the later transition from the second to
third quartile, since the difference in the remaining gap
between the convex and concave sites grows rapidly, the over-
all electric field effect accelerates further and the required
temperature gradient needs to be much higher. Such a desta-
bilizing trend is additionally observed versus the curvature
increase in all these figures.
B. Curvature effect formulation
While the role of the curvature
f
r
) has been qualita-
tively explored in the previous sections, such an augmentation
(diminishing) coefficient for convex (concave) zones could
get addressed quantitatively. In this regard, Eq. (
3
) can get
normalized as
ˆ
E
=
f
(
ˆ
r
)
1
ˆ
y
I
,
where ˆ
y
I
=
y
I
/
h
is the normalized average interface elevation.
Considering a single semicircle with the prescribed normal-
ized radius of curvature ˆ
r
=
r
/
h
, set in the center of the
domain, one could recalculate the voltage field, repeating
the procedure in Sec.
III
. Figures
6(a)
and
6(b)
illustrate
the resulting voltage distributions for a convex and concave
interface, respectively.
From Fig.
6(a)
, the augmentation coefficient
f
max
of the
electric field in the convex sites due to the radius of curvature
could be formulated. As an initial step, the limits can be
placed, where for the minimal radius of curvature this radius
gets indefinitely large, and vice versa for the very large radius
of curvature it becomes ineffective. Hence,
lim
ˆ
r
0
f
max
=∞
lim
ˆ
r
→∞
f
max
=
1
.
Therefore, one could assign an interpolating function sat-
isfying these boundary conditions, as below:
f
max
=
exp
(
α
1
ˆ
r
)
+
β
1
exp
(
α
1
ˆ
r
)
1
,
(24)
where
α
1
and
β
1
are the relaxation coefficients. Similarly,
for Fig.
6(b)
, the diminishing coefficient
f
min
of the electric
field in the concave regions due to the radius of curvature
could be interpolated. In this regard, for the minimal radius
of curvature, this coefficient gets indefinitely small, and vice
versa for a very large radius of curvature it yields to unity (i.e.,
becomes ineffective). Hence,
lim
ˆ
r
0
f
min
=
0
lim
ˆ
r
→∞
f
min
=
1
,
and one can define an interpolating function form to satisfy
the above boundary condition as
f
min
=
exp
(
α
2
ˆ
r
)
1
exp
(
α
2
ˆ
r
)
+
β
2
,
(25)
where
α
2
and
β
2
are the relaxation coefficients. For either
augmentation or diminishing ratios, the coefficients of
α
i
and
015505-7