# Dynamics of fluid-driven fractures in the viscous-dominated regime

Sri Savya Tanikella1, Marie C. Sigallon1,2, Emilie Dressaire1
1\. University of California Santa Barbara
2\. Ecole Polytechnique, CNRS, Institut Polytechnique de Paris

###### Abstract

During hydraulic fracturing, the injection of a pressurized fluid in a brittle
elastic medium leads to the formation and growth of fluid-filled fractures. A
disc-like or penny-shaped fracture grows radially from a point source during
the injection of a viscous fluid at a constant flow rate. We report an
experimental study on the dynamics of fractures propagating in the viscous
regime. We measure the fracture aperture and radius over time for varying
mechanical properties of the medium and fluid and different injection
parameters. Our experiments show that the fracture continues to expand in an
impermeable brittle matrix, even after the injection stops. In the viscous
regime, the fracture radius scales as $t^{4/9}$ during the injection. Post
shut-in, the crack continues to propagate at a slower rate, which agrees well
with the predictions of the scaling arguments, as the radius scales as
$t^{1/9}$. The fracture finally reaches an equilibrium set by the toughness of
the material. The results provide insights into the propagation of hydraulic
fractures in rocks.

##  1 Introduction

Hydraulic fracturing is a well-stimulation technique used to recover natural
gas and oil from reservoirs with low permeability, such as shale formations.
The US Environmental Protection Agency reports that the natural gas production
from hydraulically fractured wells in the United States saw a 10-fold increase
between 2000 and 2015 [1]. The formation of fractures in rocks has other
applications, including carbon sequestration and geothermal energy extraction
[2, 3].

It is estimated that almost 1 million wells have been hydraulically fractured
since the 1940s. As hydraulic fracturing has become more prevalent, so has the
need to characterize the associated risk to the local environment and
populations. Over the past two decades, groundwater contamination and induced
seismicity have been linked to hydraulic fracturing operations. Between 2000
and 2013, there was at least one hydraulically fractured well system within 1
mile of the water sources of 3900 public water systems in the continental
United States [4]. The water from these systems was distributed to more than
8.6 million people year-round in 2013. Another major concern is the induced
earthquakes associated with hydraulic fracturing [5]. Earthquakes of magnitude
4.0 can be caused by the disposal of wastewater in fractured reservoirs [6, 7,
8, 9]. The associated risks increase as the distance between stimulated wells
and groundwater wells or fault lines decreases. Therefore, understanding the
dynamics of a fracture during and after the fluid injection is critical for
risk assessment.

When a pressurized Newtonian fluid is injected from a point source into a
uniform impermeable brittle matrix, a disk-like hydraulic fracture forms and
propagates. This penny-shaped fracture results from the coupling of three
mechanisms: (1) the elastic deformation of the fracture surfaces, (2) the
propagation of the fracture at the rim of the fluid-filled region, and (3) the
flow of the fluid in the fracture. Theoretical and numerical modeling of the
penny-shaped crack has been developed since the seminal work of Sneddon $\&$
Mott [10]. Yet the fracture dynamics remain complicated to model owing to the
multi-scale nature of the problem, [11, 12], with viscous dissipation
associated with fluid transport through the volume and the stress
concentration at the tip of the fracture. In limiting regimes, in which the
viscous dissipation or the fracture opening controls the dynamics of the
fracture, tip asymptotes can be defined [13, 14, 15, 16, 17, 18].

In the viscous-dominated or zero-toughness regime, the elastic stresses drive
the radial fracture propagation, which is limited by the viscous stresses. In
the toughness-dominated regime, the stress concentration or stress intensity
factor controls the fracture expansion, and the viscous stresses are
negligible. To study the two-propagation regimes, laboratory-scale experiments
use hydrogels, whose brittle elastic properties are analogous to those of
rocks. For example, gelatin is a clear material that allows fracture
visualization for a wide range of mechanical properties [19, 20]. Upon
injection of an aqueous solution or oil from a needle, a penny-shaped fracture
forms at the injection point and expands radially through the gelatin,
perpendicularly to the needle. Radius and aperture measurements in
experimental model systems agree with the theoretical predictions during the
fluid injection. [20, 21, 22, 23].

The dynamics of a fracture after the injection, i.e., post-shut-in, differs
depending on the propagation regime [24]. Here we assume that the fracture
propagates in a mobile equilibrium with no fluid lag. The influence of gravity
on fracture propagation is negligible [25]. In the toughness-dominated regime,
the fluid-filled fracture propagates when the stress intensity factor is equal
to the toughness of the matrix. When the fluid injection stops, the elastic
pressure is no longer sufficient to sustain the propagation of the fracture.
In the toughness regime, the final geometry of the fracture is reached when
the injection stops. In the viscous regime, however, the material’s toughness
does not limit the propagation of the fracture. The viscous dissipation in the
fluid balances the elastic stresses during and after the injection. This study
focuses on the dynamics of fractures that propagate in the viscous-dominated
regime during and post-injection. Previous experimental studies have
investigated the fracture dynamics during the injection of a viscous fluid in
a hydrogel matrix [20, 23] and the closure of the fracture post-injection in
porous materials like cement and plaster with a typical Young’s modulus of the
order of 1 GPa [26]. Those studies characterized the decline in pressure
inside the fracture post-shut-in, where the pressure decrease is due to
various phenomena such as leak-off. Numerical studies predicted fracture
growth and pressure reduction after shut-in in non-permeable matrices [27,
28]. Here, we conduct injection experiments in gelatin blocks to study the
post shut-in dynamics of a fracture in the viscous regime. The fracture
initially forms and propagates during the injection of a viscous Newtonian
liquid. The injection is stopped, and we record the time dependence of the
fracture radius and aperture until both properties become constant, indicating
that the fracture has reached its equilibrium configuration. We observe that
the fracture created in the viscous regime continues to propagate even after
the injection stops. Consequently, we identify three different regimes of
propagation in our experiments: (1) propagation during injection, (2)
propagation post shut-in or at constant volume, and (3) saturation. To the
best of our knowledge, the experiments presented here are the first
observations for the three regimes in a hydrogel matrix, and the data agree
well with the corresponding scaling laws [24, 27, 28].

This paper is structured as follows. In §2, we discuss the experimental set-up
and methods and our observations. In §3, we summarize the scaling arguments
and derivation of the dimensionless parameters. The experimental results and
theoretical predictions are compared in §4. Our conclusions are summarized in
§5.

##  2 Experiments and observations

To study the fracture dynamics in the viscous-dominated regime, we inject a
viscous liquid in a high Young’s modulus and low toughness gelatin block.
Using dyed silicone oil, we can characterize the geometry of the fracture,
i.e., its radius and aperture, as a function of the radial distance from the
injection point during and after the injection.

###  2.1 Experimental methods

Figure 1: Experimental set-up: (a) schematic and (b) penny-shaped fracture
formed by injecting silicone oil of viscosity $\mu=10$ Pa.s at a flow rate
$Q_{0}=10$ ml.min-1) into a block of gelatin with Young’s modulus $E=88$ kPa.
The oil is dyed with red oil-soluble food color.

The gelatin is prepared by heating ultra-pure water to $60^{o}$C and slowly
adding gelatin powder (Gelatin type A; Sigma-Aldrich, USA) while mixing. The
gelatin is then allowed to be set over 24 hours at room temperature in a cubic
clear container (12.5 cm $\times$ 12.5 cm $\times$ 12.5 cm) around a blunt
needle as represented in figure 1. The Young’s modulus of the gelatin $E$ is
measured with cylindrical samples of height and diameter equal to $2.5$ cm.
The cylinders are tested under compression using a custom-built displacement-
controlled load frame. The Young’s moduli range between $88$ to $144$ KPa
${\pm 10\%}$ for mass fractions of gelatin powder in water between $20$ \-
$25\%$. The fracture energy and Poisson’s ratio of the gelatin are assumed
constant with $\gamma_{S}\approx 1$ J.m-2 and $\nu\approx 0.5$ respectively
[29]. The inner diameter of the blunt needle is equal to $2.15$ mm. A plastic
washer of diameter of about 6 mm is placed at the tip of the needle to ensure
the propagation of the fracture in a plane perpendicular to the needle and
parallel to the lens of the camera. We inject silicone oils of different
viscosity $\mu$ to fracture the gelatin. Viscosity measurements are conducted
using an MCR 92 Anton Parr rheometer with a parallel plate measuring system.
The values obtained at $20^{o}$C are listed in table 1 and have an error of
$\pm 1\,\%$. We use a syringe pump (KDS Legato 200 series infusion syringe
pump) to inject the fracturing fluid at a controlled flow rate $Q_{0}$ ranging
from $5$ to $28$ ml.min-1. The injection stops when a volume $V_{0}$ of fluid
has been injected.
To ensure that the experiments are in the viscous regime, we estimate the
ratio of the toughness-related pressure and the viscous pressure. This ratio
is called the dimensionless toughness ($\mathcal{K}_{s}$), and its maximum
value is reached at the time of shut-in [24]

$\mathcal{K}_{s}=K^{\prime}\frac{t_{s}^{1/9}}{E^{\prime 13/18}\mu^{\prime
5/18}Q_{o}^{1/6}},$ (1)

with the effective viscosity $\mu^{\prime}=12\mu$, the effective toughness
$K^{\prime}=4\left(\frac{2}{\pi}\right)^{\frac{1}{2}}K_{IC}$, the toughness
$K_{IC}=\sqrt{2\gamma_{S}E^{\prime}}$ and the effective Young’s modulus
$E^{\prime}=E/\left(1-\nu^{2}\right)$. The dimensionless toughness is of order
1 for all experiments as presented in table 1. We, therefore, expect the
fracture propagation to be limited by the viscous dissipation associated with
the fluid flow in the fracture.
The silicone oil is dyed using oil-based food color to help visualize the
propagation of the fracture. The list of experiments and the corresponding
parameters are summarized in table 1. The propagation of the fracture is
recorded using a Nikon D5300 camera with a Phlox® LED panel ensuring uniform
backlighting. The images are processed using a custom-made MATLAB code to
determine the radius of the fracture $R$.

Exp. |  | Markers |  | E (kPa) |  |  $\mu$ (Pa.s) |  |  $Q_{0}$ (ml.min-1) |  |  $V_{0}$ (ml) |  | $\mathcal{K}_{s}$
---|---|---|---|---|---|---|---|---|---|---|---|---
1 |  |  $\color[rgb]{0.9290,0.6940,0.125}{{{\bullet}}}\color[rgb]{0,0,0}{{{{\circ}}}}$ / $\color[rgb]{0,0,0}{{{{\circ}}}}$ |  | 88 |  | 10.3 |  | 15 |  | 6 |  | 1.6
2 |  |  $\blacktriangle$ $\bigtriangleup$ /  $\bigtriangleup$ |  | 88 |  | 10.3 |  | 5 |  | 6 |  | 2.17
3 |  |  ${\color[rgb]{0,0,0}{{{\blacktriangle}}}}$ $\bigtriangleup$ / $\bigtriangleup$ |  | 88 |  | 10.3 |  | 10 |  | 6 |  | 1.79
4 |  |  $\blacktriangle$ $\bigtriangleup$ /  $\bigtriangleup$ |  | 88 |  | 10.3 |  | 25 |  | 4.3 |  | 1.34
5 |  |  $\blacktriangle$ $\bigtriangleup$ /  $\bigtriangleup$ |  | 88 |  | 10.3 |  | 28 |  | 4.3 |  | 1.30
6 |  |  $\color[rgb]{0.3010,0.7450,0.9330}{{\blacklozenge}}\color[rgb]{0,0,0}$ $\lozenge$ / $\lozenge$ |  | 88 |  | 10.3 |  | 15 |  | 8 |  | 1.65
7 |  |  ${\color[rgb]{0.6350,0.0780,0.1840}{{{{{\times}}}}}}\color[rgb]{0,0,0}\square$ / ${\color[rgb]{0.6350,0.0780,0.1840}{{{{{\times}}}}}}$ |  | 88 |  | 10.3 |  | 15 |  | 4 |  | 1.53
8 |  |  $\color[rgb]{0,0,0}\blacksquare\color[rgb]{0,0,0}\square$ / $\square$ |  | 88 |  | 20 |  | 25 |  | 3 |  | 1.06
9 |  |  ${\color[rgb]{0.6050,0.4023,0.2344}{{{{{+}}}}}}\color[rgb]{0,0,0}\square$ / ${\color[rgb]{0.6050,0.4023,0.2344}{{{{{+}}}}}}$ |  | 116 |  | 10.3 |  | 10 |  | 4 |  | 1.61
10 |  |  $\color[rgb]{0.3010,0.7450,0.9330}\blacksquare\color[rgb]{0.3010,0.7450,0.9330}\square$ / $\color[rgb]{0.3010,0.7450,0.9330}\square$ |  | 144 |  | 10.3 |  | 20 |  | 3 |  | 1.23
11 |  |  $\color[rgb]{0.89,0.141,0.169}\blacksquare\color[rgb]{0.89,0.141,0.169}\square$ / $\color[rgb]{0.89,0.141,0.169}\square$ |  | 144 |  | 30 |  | 15 |  | 3.5 |  | 0.996
Table 1:  List of Experiments. The markers on the left and right of the $/$
symbol correspond to the data recorded during and after the injection
respectively.

###  2.2 Thickness measurements

We use the light absorption technique pioneered by Bunger [30] to measure the
fracture aperture using a soluble dye in the injected fluid. In our system,
white light illuminates the sample and a filter is placed on the camera to
measure the light intensity at a single wavelength. The filtered wavelength
corresponds to the maximum absorbance of the dye. At this wavelength, the
absorbance $A_{\lambda}$ follows Beer’s law:

$A_{\lambda}=-\log_{10}\left(\frac{I_{\lambda}}{I_{\lambda,0}}\right)=\epsilon_{\lambda}c\,h$
(2)

where $I_{\lambda,0}$ is the background intensity and $I_{\lambda}$ is the
intensity of light after it passes through a fluid layer of thickness $h$,
with a dye concentration $c$. The fitting parameter $\epsilon_{\lambda}$
depends on the dye-fluid combination and the concentration and is obtained
through calibration. Here, the fracturing fluid is dyed with Nile Red at a
concentration of 0.2 g.L-1, and the wavelength of the optical filter used was
632 nm. The value of the fitting parameter $\epsilon_{\lambda}c=3.54\times
10^{-3}$ mm-1 is obtained through the calibration process using liquid layers
whose thickness ranges from $0.14$ mm to $2$ mm.

###  2.3 Observations

During the injection process, the fracture forms at the tip of the needle and
propagates along the washer and beyond, expanding radially. Figure 2 presents
a time series of the fracture propagation (see also electronic supplementary
material, Movie S1). As the fracture grows radially, its thickness increases
as its color becomes darker. When the injection stops, the fracture expands
radially at a slower pace. As the fracture grows, its color fades, indicating
that the aperture decreases with time. Since the amount of fluid in the
fracture needs to be conserved, as the radius increases, the width of the
fracture decreases. Finally, the fracture stops growing. To ensure that the
finite volume of material and the bounding container walls are not affecting
fracture growth, we experimented with a larger volume of gelatin (see Appendix
A for more details). The maximum fracture radius is independent of the size of
the gelatin block, indicating that the viscous-dominated fracture stops
expanding when it reaches equilibrium.

##  3 Scaling arguments

The fracture results from the injection of a high-viscosity Newtonian fluid in
a brittle elastic matrix that is impermeable. We assume there is no lag
between the fracture tip and the fluid front. The fracture is initially driven
by an incompressible fluid of viscosity of $\mu$ pumped at a constant flow
rate $Q_{0}$. The elastic medium is characterized by Young’s modulus $E$,
Poisson’s ratio $\nu$, and toughness $K_{IC}$. The injection through a point
source leads to the formation and radial propagation of a penny-shaped
fracture with no lag between the fluid and the fracture tip, as represented in
figure 3. The experiments conducted in this study are in the viscous-dominated
regime during the injection. The fracture grows as the elastic stresses in the
fracture boundaries drive the fluid outward. The viscous dissipation
associated with fluid transport in the fracture limits growth. The material
toughness is negligible and does not contribute to the fracture dynamics
during the injection.

Figure 2: Fracture growth during and after the injection. The fluid is dyed
with red food color to enhance the contrast. The grayscale images presented
here are obtained by filtering the red channel from the color images. The
first images are recorded during the injection: (a) $t=0$ s, (b) $t=18$ s, (c)
$t=36$ s. The fracture continues to grow after the injection stops (d) $t=54$
s, (c) $t=126$ s, (f) $t=198$ s, (g) $t=270$ s. The recording ends when the
fracture has reached its equilibrium configuration (h) $t=396$ s. The
experimental parameters are: Young’s modulus E = 88 KPa, flow rate $Q_{0}$ =
10 ml/min, volume injected $V_{0}$ = 6 ml and fluid viscosity $\mu$ = 10 Pa.s.

The evolution of the radius and width of the fracture in the viscous dominated
regime can be modeled using scaling arguments, derived for both the injection
and post ”shut-in” stages of the propagation [15, 24] We make the following
assumptions: (i) the stresses in the matrix are well described by Linear
Elastic Fracture Mechanics (LEFM), (ii) the lubrication theory can be used to
model the flow; (iii) the fracture propagates continuously in a mobile
equilibrium; and (iv) the matrix is impermeable (no leak-off). To describe the
fracture aperture $w(r,t)$, radius $R(t)$, and pressure $p(r,t)$, we solve the
coupled equations that describe (a) the viscous flow of the fracturing fluid
in the time-dependent fracture, see equation (3), (b) the elastic deformation
of the solid material, see equation (4), (c) the fracture propagation criteria
based on LEFM, see equation (5), and the global mass balance or fluid mass
conservation, see equations (6) or (7). The net pressure couples these
equations. The non-dimensional forms of the equations are summarised below,
and a detailed derivation is provided in Appendix B. The equations are non-
dimensionalized with $R=R_{o}\hat{R}$, $w=w_{o}\hat{w}$, $P=P_{o}\hat{P}$, and
$t=t_{o}\hat{t}$, where $R_{o}$, $w_{o}$, $P_{o}$ and $t_{o}$ represent the
characteristic radius, aperture, pressure of the fracture and timescale of the
propagation, respectively. As proposed by Savitski and Detournay [15], we
define the effective viscosity $\mu^{\prime}=12\mu$, the effective toughness
$K^{\prime}=4\left(\frac{2}{\pi}\right)^{\frac{1}{2}}K_{IC}$ and the effective
Young’s modulus $E^{\prime}=E/\left(1-\nu^{2}\right)$. The dimensionless
lubrication equation writes [23]:

$\frac{\partial\hat{w}}{\partial\hat{t}}=\frac{t_{0}w_{0}^{2}p_{0}}{\mu^{\prime}R_{0}^{2}}\frac{1}{\hat{r}}\frac{\partial}{\partial\hat{r}}\left(\hat{r}\hat{w}^{3}\frac{\partial\hat{p}}{\partial\hat{r}}\right).$
(3)

The elastic deformation leads to the following:

$\hat{w}=\frac{8}{\pi}\frac{p_{0}R_{0}\hat{R}}{w_{0}E^{\prime}}\int_{\hat{r}/\hat{R}}^{1}\frac{\xi}{\sqrt{\xi^{2}-(\hat{r}/\hat{R})^{2}}}\int_{0}^{1}\frac{x\hat{p}}{\sqrt{1-x^{2}}}\mathrm{d}x\mathrm{d}\xi.$
(4)

The fracture propagation requires:

$\frac{K^{\prime}}{p_{0}R_{0}^{1/2}}=\frac{2^{7/2}}{\pi\sqrt{\hat{R}}}\int_{0}^{\hat{R}}\frac{\hat{p}}{\sqrt{\hat{R}^{2}-\hat{r}^{2}}}\hat{r}\mathrm{~{}d}\hat{r}.$
(5)

Since the gelatin is considered impermeable, the volume of the fracture is
equal to the volume of fluid injected. During the injection, the volume of the
fracture is equal to

$\hat{Q}\hat{t}=2\pi\frac{R_{0}^{2}w_{0}}{Q_{0}t_{0}}\int_{0}^{\hat{R}}\hat{r}\hat{w}\mathrm{~{}d}\hat{r}.$
(6)

The injection stops when the fracture volume is equal to $V_{0}=Q_{0}t_{0}$.
After the injection, the volume of the fracture writes:

$\hat{V}=2\pi\frac{R_{0}^{2}w_{0}}{V_{0}}\int_{0}^{\hat{R}}\hat{r}\hat{w}\mathrm{~{}d}\hat{r}.$
(7)

To determine the fracture dynamics during and after the injection, we now
solve the sub-set of equations relevant to each stage of fracture propagation.

Figure 3: Diagram showing the radial fracture geometry.

###  3.1 Viscous regime: constant flow rate propagation

To obtain the scaling relations that describe the fracture dynamics during the
injection at a constant flow rate, we assume that the major form of
dissipation of energy in the fracture is the viscous dissipation of the fluid
flow. We set the dimensionless parameters in equations (3), (4) and (6) equal
to 1. We recover the scaling relations for the radius $R_{o}$ and aperture
$w_{o}$ of the fracture originally derived by Savitski and Detournay [15]:

$\displaystyle R_{0}$ $\displaystyle\approx$
$\displaystyle\left(\frac{Q_{0}^{3}E^{\prime}t_{0}^{4}}{\mu^{\prime}}\right)^{1/9}$
(8) $\displaystyle w_{0}$ $\displaystyle\approx$
$\displaystyle\left(\frac{\mu^{\prime 2}Q_{0}^{3}t_{0}}{E^{\prime
2}}\right)^{1/9}.$ (9)

During the injection, both the radius and the aperture are increasing
functions of time, which is consistent with our observations.

###  3.2 Viscous regime: propagation at constant volume

Once the injection is complete, the volume inside of the fracture is constant
[24]. The viscous dissipation of the fluid is assumed to be the limiting
factor in the propagation of the fracture. Setting the dimensionless
parameters of equations (3), (4) to 1 and (7) and solving for the
characteristic radius $R_{0}$ and aperture $w_{0}$, we get:

$\displaystyle R_{0}$ $\displaystyle\approx$
$\displaystyle\left(\frac{V_{0}^{3}E^{\prime}t_{0}}{\mu^{\prime}}\right)^{1/9}$
(10) $\displaystyle w_{0}$ $\displaystyle\approx$
$\displaystyle\left(\frac{\mu^{\prime 2}V_{0}^{3}}{E^{\prime
2}t_{0}^{2}}\right)^{1/9}.$ (11)

The radius is still increasing with time, at a slower rate than during the
injection. The aperture of the fracture is now a decreasing function of time.
This is consistent with the fact that the total volume of fluid in the
fracture ($\propto R_{0}^{2}w_{0}$) needs to be conserved, independent of
time.

###  3.3 Toughness regime: saturation

Based on the scaling relations derived for the viscous regime, the fracture
radius is an increasing function of time at constant fracture volume. Yet, the
elastic pressure in the fracture is a decreasing function of time, as the wall
deformation decreases. Eventually, the material toughness is no longer
negligible. Indeed, if the pressure in the fracture is too low, the material
no longer fractures and the propagation of the fracture stops. We assume that
saturation, unlike the previous two stages of propagation, is controlled by
the fracture opening or material toughness. The fracture opening criterion as
a function of the pressure and stress intensity factor is given by equation
(5). This equation coupled with the elastic stress in the gelatin matrix,
equation (4) and the volume conservation of the fracture, equation (7) leads
to the following scaling [24]:

$\displaystyle R_{0}$ $\displaystyle\approx$
$\displaystyle\left(\frac{V_{0}E^{\prime}}{K^{\prime}}\right)^{2/5}$ (12)
$\displaystyle w_{0}$ $\displaystyle\approx$
$\displaystyle\left(\frac{K^{\prime 4}V_{0}}{E^{\prime 4}}\right)^{1/5}.$ (13)

The radius and the aperture of the fracture are now functions of the volume of
fluid and mechanical properties of the matrix. They are independent of time.

##  4 Results and observations

In the following section, we report quantitative experimental results, and
their comparison with the scaling arguments derived above for both the radius
and width of the fracture.

Figure 4: Radius of the fracture formed during experiments 1-9, see table 1
for corresponding parameters. Data collected during and after the injection
are displayed with solid and open symbols, respectively. (a) The radius of the
fracture over time. (b) Rescaled radius using equation (8) for the regime (i)
as a function of time. (c) Rescaled radius using equation (10) for the regime
(ii) as a function of time. (d) Rescaled radius using equation (12) for the
regime (iii) as a function of time. For each regime, the best fit line is
represented by a solid black line.

###  4.1 Radius measurements

We first measure the radius of the fracture as a function of time for the
experimental parameters summarized in table 1. The raw experimental data are
shown in figure 4(a). The solid symbols correspond to the values recorded
during the injection. In contrast, the open markers indicate that the data
were recorded after the injection. For some of the experiments, such as
experiment $8$, we see a $50\,\%$ increase in the radius of the fracture after
the injection stops.

We successively rescale the data using the scaling laws derived above and
obtain figures 4(b-d), where the fit line is shown in black. On the log-log
plot, the slope of the line is set to the value of the power-law derived using
scaling arguments. The y-intercept is obtained by minimizing the mean squared
error and corresponds to the prefactor, which can be predicted theoretically
[17]. The propagation dynamics exhibit three stages described by (i) the
viscous-dominated propagation during the injection, (ii) the viscous-dominated
propagation at constant volume after the injection stops, and (iii) a
toughness-controlled saturation. For the fracture dynamics during the
injection [regime (i) and figure 4(b)], we rescale the data using the viscous
scaling. After the early times that correspond to the fracture formation and
its propagation over the washer, the data points collapse onto a best-fit line
which has a prefactor, $k_{1}=0.28$ and an exponent, $\alpha_{1}=4/9$ where

$R(t)=k_{1}\left(\frac{E^{\prime}Q_{0}^{3}}{\mu^{\prime}}\right)^{1/9}t^{\alpha_{1}}.$
(14)

The theoretical prefactor derived by Savitski and Detournay [17] is equal to
0.7. The low value of the prefactor obtained in our experiments is consistent
with measurements previously reported for the viscous-dominated regime [20,
23].

In the next regime [regime (ii) and figure 4(c)], the fracture fluid continues
to propagate after the injection stops, and the open symbols collapse on the
best fit line

$R(t)=k_{2}\left(\frac{E^{\prime}V_{0}^{3}}{\mu^{\prime}}\right)^{1/9}t^{\alpha_{2}},$
(15)

which has a prefactor $k_{2}=0.39$ and an exponent, $\alpha_{2}=1/9$. The
results indicate that the propagation continues to be dominated by viscous
dissipation. For time values above 100 s, the increase in radius no longer
follows the best fit line: the radius reaches its maximum or equilibrium
value.

In the saturation regime, regime (iii), and figure 4(d)], we rescale the data
using the volume of the fracture and the mechanical properties of the gelatin
matrix. All data collapse on an average value represented with a horizontal
line whose y-intercept $k_{3}\approx 0.58$ where,

$R(t)=k_{3}\left(\frac{E^{\prime}V_{0}}{K^{\prime}}\right)^{2/5}.$ (16)

The experimental prefactor of 0.58 is comparable to the expected theoretical
prefactor of 0.85 [24]. This result indicates that the propagation controlled
by the viscous dissipation stops when the stress intensity factor at the tip
of the fracture is too low to sustain the formation of the fracture.

Figure 5: (a) Fracture profiles measured during and after the injection for
$t=[9,21,33,81,105,117,762]$ s with time increasing from clear red to dark
blue. (b) Rescaled fracture profiles using the scaling laws derived for the
injection in §2.1. (c) Rescaled fracture profiles using the scaling laws
derived for the constant volume propagation in §2.2. Experimental parameters:
$E=88$ KPa, $Q_{0}=10$ ml.min-1, $V_{0}=6$ ml, and $\mu=10$ Pa.s

The time dependence of the radius demonstrates the succession of three regimes
of fracture propagation: growth at constant flow rate, growth at constant
volume and saturation. The growth of the fracture at constant volume is
characterized by its duration and the relative change in fracture radius. The
duration and the relative change in radius decrease as the dimensionless
toughness of the injection increases [24]. Indeed the larger the value of
$\kappa_{s}$, the closer to the toughness regime the injection is. The
fracture is arrested immediately after the injection stops for a cut-off value
of $\mathcal{K}_{s}=2.5$. For the experiments conducted in this study, the
$\mathcal{K}_{s}$ values have been listed in table 1 and range between
0.996-2.17, which explains why the growth at constant volume lasts a few
minutes for a relative change in radius is about 50$\%$.

###  4.2 Thickness measurements

To further characterize the three propagation regimes, we measure the fracture
aperture with the dye absorption method and plot the profiles of the fractures
in figure 5. During the injection, both the radius and aperture of the
fracture increase. Upon rescaling, the data collapse on a self-similar
fracture profile, after an initial transient regime, as shown in figure 5(b).
After the injection, the radius of the fracture increases as the aperture
decreases. The rate of propagation is slower than it was during the injection.
Upon rescaling, all data collapse on a second self-similar profile which
corresponds to the viscous propagation of a fracture of constant volume (see
figure 5(c)).

##  5 Conclusion

As a pressurized fluid is injected in an elastic brittle material, a penny-
shaped fracture forms and propagates. The complex fracture dynamics depend on
the matrix and fluid properties and the injection parameters. During the fluid
injection, modeling and experimental studies have demonstrated two asymptotic
regimes. The fracture expansion is either controlled by the viscous
dissipation in the fluid, in the viscous-dominated regime or the toughness of
the material, in the toughness-dominated regime. Upon shut-in, the continued
propagation of the fracture at constant volume has been observed in porous
materials and predicted for impermeable matrices. In this study, we
experimentally study the propagation regimes of a viscous-dominated fracture
in a hydrogel matrix. We demonstrate the existence of three propagation
regimes: injection growth, post shut-in propagation, and saturation. For each
regime, we show a good agreement with the scaling laws derived for the growth
of fractures in the viscous-dominated regime during injection at a constant
flow rate and post shut-in. The saturation values of the fracture radius and
aperture are reached when the stress intensity factor at the tip of the
fracture becomes lower than the material’s toughness. The experimental results
and corresponding model allow for predicting the relative growth of the
fracture after the injection stops, i.e., the final size of the fracture and
how long it takes for the fracture to reach this equilibrium geometry once the
injection stops.

The datasets supporting this article have been uploaded as part of the
supplementary material [31].

## Appendix A. Effect of box dimensions

All experiments are conducted in a 12.5 cm $\times$ 12.5 cm $\times$ 12.5 cm
block of gelatin set in an acrylic box of the same dimensions. To test the
influence of the box size or edge effects on the fracture dynamics and the
saturation radius, we performed experiment 3 (see table 1 for experimental
parameters in a larger volume of gelatin of dimensions 15 cm $\times$ 15 cm
$\times$ 15 cm.

Figure 6: Influence of matrix dimensions on fracture propagation. Radius vs
time for experiment 3 (see table 1 for experimental parameters) conducted in a
box of size 12.5 cm $\times$ 12.5 cm $\times$ 12.5 cm
(${\color[rgb]{0,0,0}{{{\blacktriangle}}}}$ $\bigtriangleup$) and a box of
size 15 cm $\times$ 15 cm $\times$ 15 cm
(${\color[rgb]{.5,.5,.5}{{{\blacktriangle}}}}$ $\bigtriangleup$).

The radius of the fracture is recorded over time and plotted in fig.6, for the
two block sizes. The results demonstrate that the evolution of the radius does
not depend on the size of the block of gelatin. The saturation radius is not
set by edge effects.

## Appendix B. Governing equations

We review the mathematical derivations that define the radius and aperture of
penny-shaped fracture driven by a fluid [23].

The mechanical deformation in the elastic matrix associated with the fracture
thickness w relates to the pressure in the fracture and the fracture radius
through the following integral relation, initially derived by Sneddon $\&$
Lowengrub [32].

$w(r,t)=\frac{8R}{\pi
E^{\prime}}\int_{r/R}^{1}\frac{\xi}{\sqrt{\xi^{2}-(r/R)^{2}}}\int_{0}^{1}\frac{xp(x\xi
R,t)}{\sqrt{1-x^{2}}}\mathrm{~{}d}x\mathrm{~{}d}\xi.$ (17)

To describe the fluid flow in the fracture, we use the non-linear lubrication
equation called the Reynolds equations [33] which relates the aperture of the
fracture to the pressure and the fracture radius.

$\frac{\partial w(r,t)}{\partial
t}=\frac{1}{12\mu}\frac{1}{r}\frac{\partial}{\partial
r}\left(rw^{3}(r,t)\frac{\partial p}{\partial r}\right)$ (18)

The stress intensity factor $K_{I}$ defines the stress concentration at the
tip of the fracture. The fracture propagates if the stress intensity factor
$K_{I}$ is equal to $K_{I}C$, i.e., the material toughness. For a penny-shaped
fracture, the stress intensity is equal to [34]:

$K_{I}=\frac{2}{\sqrt{\pi
R}}\int_{0}^{R(t)}\frac{p(r,t)}{\sqrt{R^{2}-r^{2}}}r\mathrm{~{}d}r.$ (19)

The boundary conditions are set by the fracture geometry. The integral
representation of the fluid mass conservation in the fracture is in equation
(20). The fracture thickness at the tip is 0, equation (21). There is no flow
through the tip of the fracture in the elastic medium, equation (22).

$\displaystyle Qt=2\pi\int_{0}^{R(t)}rw(r,t)\mathrm{d}r$ (20) $\displaystyle
w=0,\quad r=R(t)$ (21) $\displaystyle w^{3}(r,t)\frac{\partial
p(r,t)}{\partial r}=0,\quad r=R(t)$ (22)

In this study, this set of coupled equations and boundary conditions are used
to derive scaling laws for the fracture radius and aperture.

## References

  * [1] U.S. EPA (U.S. Environmental Protection Agency). 2016 Hydraulic Fracturing for Oil and Gas: Impacts from the Hydraulic Fracturing Water Cycle on Drinking Water Resources in the United States. Technical report. Office of Research and Development, Washington, DC.
  * [2] Huppert HE, Neufeld JA. 2013. The fluid mechanics of carbon dioxide sequestration. Annu. Rev. Fluid Mech. 46 (1), 255–272.
  * [3] Caulk RA, Ghazanfari E, Perdrial JN, PERDRIAL N. 2016. Experimental investigation of fracture aperture and permeability change within enhanced geothermal systems. Geothermics 62, 12–21.
  * [4] Jasechko S, Perrone D. 2017. Hydraulic fracturing near domestic groundwater wells. Proc. National Acad. Sci. 114, 13138–13143.
  * [5] Yu H, Harrington RM, Kao H, Liu Y, Wang B. 2021. Fluid-injection-induced earthquakes characterized by hybrid-frequency waveforms manifest the transition from aseismic to seismic slip. Nat. Commun. 12, 6862.
  * [6] Folger P, Tienmann M. 2016 Human-Induced Earthquakes from Deep-Well Injection: A Brief Overview. Congr. es. Serv. 42.
  * [7] Schultz R, Atkinson G, Eaton DW, Gu YJ, Kao H. 2018 Hydraulic fracturing volume is associated with induced earthquake productivity in the Duvernay play. Science 359, 304-308.
  * [8] U.S. EPA (U.S. Environmental Protection Agency). 2016 Minimizing and Managing Potential Impacts of Injection Induced Seismicity from Class II Disposal Wells: Practical Approaches. Technical report. Underground Injection Control National Technical Workgroup. Washington, DC.
  * [9] Ellsworth WL. 2013 Injection-Induced Earthquakes. Science 341, 1225942.
  * [10] Sneddon IN, Mott NF. 1946 The distribution of stress in the neighbourhood of a crack in an elastic solid. Proc. R. Soc. Lond. A 187, 229-260.
  * [11] Detournay E. 2016 Mechanics of Hydraulic Fractures. Annu. Rev. Fluid Mech. 48, 311-339.
  * [12] Osiptsov AA. 2017 Fluid Mechanics of Hydraulic Fracturing: a Review. J. Petrol. Sci. Eng. 156, 513–535.
  * [13] Spence DA, Sharp P. 1985 Self-similar solutions for elastohydrodynamic cavity flow. Proc. R. Soc. Lond. A 400, 289 – 313.
  * [14] Savitski AA, Detournay E. 2001 Similarity solution of a penny-shaped fluid-driven fracture in a zero-toughness linear elastic solid. C. R. Mec. 329, 255 – 262.
  * [15] Savitski AA, Detournay E. 2002 Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: asymptotic solutions. Int. J. Solids Struct. 39, 6311 – 6337.
  * [16] Garagash DI, Detournay E. 2000 The tip region of a fluid-driven fracture in an elastic medium. J. Appl. Mech. 67, 183 – 192.
  * [17] Garagash DI, Detournay E. 2005 Plane-Strain Propagation of a Fluid-Driven Fracture: Small Toughness Solution. J. Appl. Mech. 72, 916 - 928.
  * [18] Garagash DI, Detournay E, Adachi JI. 2011 Multiscale tip asymptotics in hydraulic fracture with leak-off. J. Fluid Mech. 669, 260 – 297.
  * [19] Giuseppe ED, Funiciello F, Corbi F, Ranalli G, Mojoli G. 2009 Gelatins as rock analogs: A systematic study of their rheological and physical properties. Tectonophysics 473, 391–403.
  * [20] Lai CY, Zheng Z, Dressaire E, Wexler JS, Stone HA. 2015 Experimental study on penny-shaped fluid-driven cracks in an elastic matrix. Proc. R. Soc. Lond. A 471, 20150255 – 10.
  * [21] Bunger AP, Detournay E. 2008 Experimental validation of the tip asymptotics for a fluid-driven crack. J. Mech. Phys. Solids 56, 3101–3115.
  * [22] Lai CY, Zheng Z, Dressaire E, Stone HA. 2016 Fluid-driven cracks in an elastic matrix in the toughness-dominated limit. Phil. Trans. R. Soc. Lond. A 374, 20150425.
  * [23] O’Keeffe NJ, Huppert HE, Linden PF. 2018 Experimental exploration of fluid-driven cracks in brittle hydrogels. J. Fluid Mech. 844, 435-458.
  * [24] Möri A, and Brice L. 2021. Arrest of a radial hydraulic fracture upon shut-in of the injection. Int. J. of Solids Struct.,219-220, p.151-165.
  * [25] Lister JR, Kerr RC. 1991 Fluid-mechanical models of crack propagation and their application to magma transport in dykes. J. Geophys. Res. 96, 10049–10077.
  * [26] Van Dam D, de P, and Romĳn R. 2000 Analysis of Hydraulic Fracture Closure in Laboratory Experiments. SPE Production & Facilities, 15(03), p.151-158.
  * [27] Garagash DI. 2006 Plane-strain propagation of a fluid-driven fracture during injection and shut-in: Asymptotics of large toughness. Eng. Fract. Mech., 73-4, p.456-481.
  * [28] Lakirouhani A, Detournay E, Bunger AP. 2016. A reassessment of in situ stress determination by hydraulic fracturing. Geophysical Journal International, 205-3, p.1859–1873
  * [29] Menand T, Tait SR. 2002 The propagation of a buoyant liquid-filled fissure from a source under constant pressure: An experimental approach. J. Geophys. Res. 107, ECV 16–1–ECV 16–14..
  * [30] Bunger AP. 2006 A photometry method for measuring the opening of fluid-filled fractures. Meas. Sci. Technol.  17, 3237.
  * [31] Tanikella SS, Sigallon MC, Dressaire E. 2022 Data from: Dynamics of fluid-driven fractures in the viscous-dominated regime
  * [32] Sneddon IN, Lowengrub M. 1969 Crack Problems in the Classical Theory of Elasticity. Wiley.
  * [33] Batchelor GK. 1967 An introduction to fluid dynamics. Cambridge University Press.
  * [34] Rice JR. 1968 Mathematical analysis in the mechanics of fracture. Academic Press, N. Y.


# On Dancer’s conjecture for stable solutions with sign-changing nonlinearity

Yong Liu  Department of Mathematics, University of Science and Technology of
China, Hefei, China, 230026<EMAIL_ADDRESS>,  Kelei Wang  School of
Mathematics and Statistics, Wuhan University, Wuhan, China, 430072
<EMAIL_ADDRESS>,  Juncheng Wei  Department of Mathematics, University of
British Columbia, Vancouver, B.C., Canada, V6T 1Z2<EMAIL_ADDRESS>and  Ke
Wu  School of Mathematics and Statistics, Wuhan University, Wuhan, China,
430072<EMAIL_ADDRESS>
###### Abstract.

We establish a Liouville type result for stable solutions for a wide class of
second order semilinear elliptic equations in $\mathbb{R}^{n}$ with sign-
changing nonlinearity $f$. Under the hypothesis that the equation does not
have any nonconstant one dimensional stable solution, and a further
nondegeneracy condition of $f$ at its zero points, we show that in any
dimension, stable solutions of the equation must be constant. This partially
answers a question raised by Dancer.

##  1\. Introduction

We are interested in the classification of stable solutions for general
elliptic equation of the form

$-\Delta u=f(u),\quad\text{in}\quad\mathbb{R}^{n},$ (1.1)

where the nonlinearity $f$ is assumed to be a $C^{1}$ function. Recall that a
$C^{2}$ solution $u$ of this equation is said to be stable, if and only if

$\int_{\mathbb{R}^{n}}(|\nabla\psi|^{2}-f^{\prime}(u)\psi^{2})dx\geq
0,\quad\forall\psi\in C_{0}^{\infty}(\mathbb{R}^{n}).$ (1.2)

Stability is an important notion, because many solutions obtained from
variational method have finite Morse index, and therefore the classification
of stable solutions is useful in the analysis of these solutions. In [8],
Dancer proposed the following stable solution conjecture:

######  Conjecture 1.1.

Assume $n\leq 8$ and $u$ is a bounded stable solution of the equation (1.1).
Then $u$ is one dimensional.

A solution is called one dimensional if it only depends on one direction. Note
that the constant solutions are also one dimensional. In the Allen-Cahn case,
that is, $f(u)=u-u^{3}$, Conjecture 1.1 is closely related to the well-known
De Giorgi conjecture concerning the classification of bounded monotone
solutions, which are automatically stable:

######  Conjecture 1.2.

Let $u$ be a bounded solution of the equation

$\Delta u+u-u^{3}=0\quad\text{in}\quad\mathbb{R}^{n}$

such that $\frac{\partial u}{\partial x_{n}}>0$. Then $u$ is one dimensional,
at least if $n\leq 8$.

De Giorgi conjecture was proved in dimension $n=2$ by Ghoussoub and Gui in
[15]. For $n=3$, this is proved by Ambrosio and Cabré in [1]. Savin proved in
[21] that for $4\leq n\leq 8$, Conjecture 1.2 is true under the additional
limit condition that

$u(x_{1},...,x_{n})\rightarrow\pm 1\quad\text{as}\quad
x_{n}\rightarrow\pm\infty.$ (1.3)

For $n\geq 9$, counterexamples have been constructed in [9]. We emphasize that
without any additional assumption like (1.3), the De Giorgi conjecture is
still open for dimensions $4\leq n\leq 8$.

The properties of stable solutions actually have a delicate dependence on the
sign of $f$. In the case that $f$ is nonnegative, in [7], Cabrè, Figalli, Ros-
Oton and Serra obtained several a priori estimates on stable solutions. Based
on these results, they are able to solve the Brezis’s conjecture (see [5]) on
the $L^{\infty}$ regularity of “extremal solutions”. Later on, in [10],
Dupaigne and Farina used the estimates in [7] to prove the following Liouville
theorem.

######  Theorem 1.3.

Assume that $u\in C^{2}(\mathbb{R}^{n})$ is a stable solution of (1.1) bounded
from below, where $f$ is assumed to be nonnegative. If $n\leq 10$, then $u$
must be a constant.

In view of Dupaigne and Farina’s result, Dancer’s conjecture is completely
solved when the nonlinearity $f$ is nonnegative. However, if the nonlinearity
$f$ changes sign, to the best of our knowledge, there are few results on
Dancer’s conjecture. When the solution is radially symmetric, Cabré and
Capella [6] proved that if $n\leq 10$, any bounded stable solution of (1.1) is
a constant solution, regardless of the type of the nonlinearity $f$. In [11]
Dupaigne and Farina established Liouville type results for some convex
nonlinearities. We refer to this paper and its references for more results in
this direction. In recent papers [3, 4], Berestycki and Graham proved
Liouville type and half-space theorems for monostable, ignition and bistable
nonlineariries. Let us also mention that for equations of the special form

$-\Delta u+|u|^{q-1}u=|u|^{p-1}u,$

where $1\leq q<p$, there are some works on Liouville type results on stable
solutions, see [12], [19], [22].

In this paper, we will consider Dancer’s conjecture for sign-changing
nonlinearities. To state our result in a more precise way, let us introduce
the following three conditions.

(H1). Isolated zeros:

The zeros of $f$ are isolated, or more precisely, there exists an integer $m$
such that $f$ vanishes exactly at $m$ points $a_{1}<a_{2}<\cdots<a_{m}$.

(H2). Nondegeneracy at stable zeros:

There exist constants $c_{0},\varepsilon_{0}>0$ and

$1\leq p<p^{*}(n),$ (1.4)

such that for any $a_{j}$ satisfying $f^{\prime}(a_{j})\leq 0$, we have

$c_{0}|t|^{p}\geq-f(a_{j}+t)sgn(t)\geq c_{0}^{-1}|t|^{p},\quad\text{for}\quad
t\in(-\varepsilon_{0},\varepsilon_{0}).$

Here $sgn$ is the sign function and the exponent $p^{*}(n)$ is defined by

$p^{*}(n)=\left\\{\begin{array}[]{lll}+\infty,&\text{if}\quad 2\leq n\leq
3,\\\ \frac{n(n-2)+4}{n(n-2)-4},&\text{if}\quad 4\leq n.\end{array}\right.$

(H3). Non-existence of stable one dimensional profile:

The ODE

$g^{\prime\prime}+f(g)=0\quad\text{in}\quad\mathbb{R}$ (1.5)

does not have any bounded nonconstant stable solution.

Note that monostable and bistable nonlinearies considered in [3, 4] all
satisfy (H1)-(H3).

Before proceeding, let us explain these conditions in more details. First of
all, unstable zeros of $f$, which by definition are the $a_{j}$ with
$f^{\prime}(a_{j})>0$, are actually not important in our analysis. The reason
is that they are unstable constant solutions of (1.1).

Secondly, (H2) includes the special case when all zeros of $f$ are
nondegenerate, i.e. $f^{\prime}(a_{j})\neq 0$ for every $a_{j}$. In general,
(H2) can not be removed. This nondegeneracy condition implies that if $a$ is a
zero of $f$ such that $f^{\prime}(a)=0$, then $f$ will change sign around $a$
and has to be negative at the right of $a$ in a small neighbourhood. Now
consider the special case $f(u)=u^{p}$, then the equation

$u^{\prime\prime}+u^{p}=0,\quad\text{in}\quad\mathbb{R}$

does not have one dimensional stable solution. However, if
$p>p_{\text{JL}}(n)$, then the equation

$\Delta u+u^{p}=0,\quad\text{in}\quad\mathbb{R}^{n}$

do have a nonconstant radially symmetric stable solution, where

$p_{\text{JL}}(n)=\left\\{\begin{array}[]{lll}+\infty,&\text{if}\quad 2\leq
n\leq 10,\\\ \frac{(n-2)^{2}-4n+8\sqrt{n-1}}{(n-2)(n-10)},&\text{if}\quad
n\geq 11.\end{array}\right.$

Note that this case does not satisfy the nondegeneracy condition.

On the other hand, we would like to point out that the upper bound on $p$ in
(1.4) may not be optimal. However, it is still not clear to us what should be
the optimal upper bound.

Thirdly, the Hypothesis (H3) will play an important role in our analysis.
Roughly speaking, it holds for “generic” nonlinearities. To have a better
understanding of this condition, we recall that for any solution of (1.5), if
it is stable, then it is monotone. This fact has already been observed by
Dancer.111This follows by applying Sturm-Liouville comparison theorem to
$g^{\prime}$, which is always a solution to the linearized equation of (1.5),
that is, an eigenfuction for the quadratic form (1.2) with zero eigenvalue.
Hence for any bounded stable solution $g$ of (1.5), the quantity

$g_{\pm}:=\lim_{t\to\pm\infty}g(t)$

are well-defined. It follows that

$\frac{1}{2}g^{\prime}(t)^{2}\equiv C-F(g(t)),$

where $F^{\prime}=f$ and $C$ is a constant determined by $g_{\pm}$. We then
see that $W(s):=C-F(s)$ is a double-well type potential in the interval
$[g_{-},g_{+}]$:

$0=W(g_{\pm})=\min_{s\in[g_{-},g_{+}]}W(s).$

The above reasoning tells us that if there is a stable solution to (1.5), then
the equation is essentially of Allen-Cahn type in the corresponding interval.
In this sense, our hypothesis exclude the classical Allen-Cahn equation. It is
worth pointing out that at this stage, the classification of stable solutions
for the Allen-Cahn equation is still open except in dimension $2$.

Now we can state our main result, which is the following

######  Theorem 1.4.

Under (H1-H3), any bounded stable solution of (1.1) is a constant.

As a direct consequence of Theorem 1.4, we have the following nonexistence
result.

######  Corollary 1.5.

Suppose $f$ has a unique zero $a$, and $f^{\prime}(a)>0$. Then (1.1) does not
have bounded stable solution.

Another consequence of Theorem 1.4 is the following result on a half space
problem.

######  Corollary 1.6.

Suppose the nonlinearity $f$ satisfies (H1-H3) on $(0,+\infty)$. If $f(0)\geq
0$ and $u$ is a bounded positive stable solution of (1.1) in
$\mathbb{R}^{n+1}_{+}:=\\{x_{n+1}>0\\}$, with the boundary condition $u=0$ on
$\partial\mathbb{R}^{n+1}_{+}$. Then $u$ is one dimensional.

######  Remark 1.7.

It is conjectured by Berestycki, Caffarelli and Nirenberg in [2] that every
bounded positive solution of the above half space problem is one dimensional.
For results on the conjecture we refer to [14, 13] and the references therein.

The proof of Theorem 1.4 uses the Sternberg-Zumbrun (see [23] and [24])
inequality, which is an equivalent formulation of the stability condition
(1.2). This inequality requires that for any $\phi\in
C_{0}^{\infty}(\mathbb{R}^{n})$,

$\int_{\mathbb{R}^{n}}\left(|\nabla^{2}u|^{2}-|\nabla|\nabla
u||^{2}\right)\phi^{2}dx\leq\int_{\mathbb{R}^{n}}|\nabla
u|^{2}|\nabla\phi|^{2}dx.$ (1.6)

Here

$|B(u)|^{2}:=|\nabla^{2}u|^{2}-|\nabla|\nabla u||^{2}$

can be regarded as a curvature term. In fact, if the level set $\\{u=u(x)\\}$
is smooth near $x$, then

$|B(u)|^{2}=|A|^{2}|\nabla u|^{2}+|\nabla_{T}|\nabla u||^{2},$

where $A$ is the second fundamental form of this level set, and $\nabla_{T}$
denotes the tangential derivative.

In our proof, the Sternberg-Zumbrun inequality will be used in the form of a
Caccioppoli inequality. It implies that at most places, $u$ is almost one
dimensional (i.e. $|B(u)|^{2}$ is small). Then we take a cube decomposition of
$\mathbb{R}^{n}$, and divide cubes into bad ones and good ones, depending
whether the integral of $|B(u)|^{2}$ on the cube is small. The integral in the
left hand side of (1.6) can be used to bound the number of bad cubes, while in
good cubes, by the nondegeneracy hypothesis (H2), $|\nabla u|$ has a decay
away from the bad set. These two estimates are combined to perform an
iteration, which eventually leads to a quadratic bound on the energy growth.
We then conclude the proof by a standard log cut-off function technique.

Notation. For every point
$\bar{x}=(\bar{x}_{1},\bar{x}_{2},\cdots,\bar{x}_{n})\in\mathbb{R}^{n}$, let

$Q_{r}(\bar{x})=\left\\{x\in\mathbb{R}^{n}:|x_{i}-\bar{x}_{i}|\leq\frac{r}{2},~{}i=1,2,\cdots,n\right\\}$

be the closed cube centered at $\bar{x}$ with side length $r$ .

We assume the zeros of $f$ satisfy, for any $j$, $a_{j+1}-a_{j}\geq
3\varepsilon_{0}$, where $\varepsilon_{0}$ is the constant in (H2).

Acknowledgement Y. Liu is supported by the National Key R&D Program of China
2022YFA1005400 and NSFC 11971026, NSFC 12141105. K. Wang is supported by
National Key R&D Program of China (No. 2022YFA1005602) and NSFC (No. 12131017
and No. 12221001). J. Wei is supported by NSERC of Canada.

##  2\. Proof of Theorem 1.4

The key in our proof is to use an iteration scheme to establish the quadratic
energy growth estimate

$\int_{Q_{R}(0)}|\nabla u|^{2}\leq CR^{2}.$

Once this is proved, a standard argument will lead to the desired theorem.
Hence in the rest of this section, we assume without loss of generality that
$n\geq 3$, since in the dimension two case, the energy automatically has at
most quadratic growth in terms of $R$. At this stage, it is worth pointing out
that for the one dimensional heteroclinic solution (which is stable) of the
Allen-Cahn equation, its energy grows like $O(R^{n-1})$.

We begin with the following

######  Lemma 2.1.

There exists a positive constant $C$ such that for any $x\in\mathbb{R}^{n}$
and every $R>1$,

$\int_{Q_{R}(x)}\left(|\nabla^{2}u|^{2}-|\nabla|\nabla u||^{2}\right)\leq
CR^{-2}\int_{Q_{2R}(x)}|\nabla u|^{2}.$ (2.1)

###### Proof.

This follows from substituting a standard cut-off function into Sternberg-
Zumbrun inequality (1.6). More precisely, let $\eta$ be a cutoff function such
that $\eta(s)=1$ for $|s|<1$, and $\eta(s)=0$ for $|s|>2$. Then we simply
choose $\phi(\cdot):=\eta(|\cdot-x|/R)$ in (1.6). ∎

Since $u$ is bounded, standard regularity theories (see [16]) imply that there
exists a constant $C$ such that

$|\nabla u|\leq C\quad\text{in}\quad\mathbb{R}^{n}.$ (2.2)

Plugging (2.2) into (2.1), we get

$\int_{Q_{R}(x)}\left(|\nabla^{2}u|^{2}-|\nabla|\nabla u||^{2}\right)\leq
CR^{n-2}.$ (2.3)

This bound will be the starting point of our iteration procedure, where the
exponent $n-2$ in the right hand side will eventually be decreased to a
negative one.

######  Lemma 2.2.

For any $\varepsilon>0$, there exists a $\delta>0$ such that, for any
$x\in\mathbb{R}^{n}$, if

$\int_{Q_{1}(x)}(|\nabla^{2}u|^{2}-|\nabla|\nabla u||^{2})<\delta,$ (2.4)

then there exists a stable zero $a_{j}$ such that

$\|u-a_{j}\|_{L^{\infty}(Q_{1}(x))}<\varepsilon.$ (2.5)

###### Proof.

Assume to the contrary that there is a sequence of points $x^{(k)}$ such that

$\int_{Q_{1}(x^{(k)})}(|\nabla^{2}u|^{2}-|\nabla|\nabla u||^{2})\to 0,$ (2.6)

but for any stable zeros $a_{j}$,

$\|u-a_{j}\|_{L^{\infty}(Q_{1}(x^{(k)}))}\geq\varepsilon.$ (2.7)

Let $u_{k}(x)=u(x+x^{(k)})$, which is a bounded stable solution of (1.1)
satisfying

$\int_{Q_{1}(0)}(|\nabla^{2}u_{k}|^{2}-|\nabla|\nabla u_{k}||^{2})\to 0.$
(2.8)

Since $\\{u_{k}\\}$ are uniformly bounded, combining standard elliptic
estimates with Arzelà-Ascoli theorem, we know that there exists a function
$u_{\infty}\in C^{2}(\mathbb{R}^{n})$ such that up to a suitable subsequence,
$\\{u_{k}\\}$ converge to $u_{\infty}$ locally uniformly in $\mathbb{R}^{n}$.
Passing to the limit in (2.8) gives

$\int_{Q_{1}(0)}(|\nabla^{2}u_{\infty}|^{2}-|\nabla|\nabla
u_{\infty}||^{2})=0.$ (2.9)

Hence $u_{\infty}$ is one dimensional in $Q_{1}(0)$. By unique continuation,
it is one dimensional in the whole $\mathbb{R}^{n}$. By (H3), there exists a
stable zero $a_{j}$ such that $u_{\infty}\equiv a_{j}$. As a consequence,

$\lim_{k\to\infty}\|u-a_{j}\|_{L^{\infty}(Q_{1}(x^{(k)}))}=\lim_{k\to\infty}\|u_{k}-a_{j}\|_{L^{\infty}(Q_{1}(0))}=0.$

This is a contradiction with (2.7). ∎

Let

$\mathbb{Z}^{n}=\\{(k_{1},k_{2},\cdots,k_{n})\in\mathbb{R}^{n}:k_{i}\in\mathbb{Z},~{}i=1,2,\cdots
n\\}$

be the integer lattice in $\mathbb{R}^{n}$. Points in $\mathbb{Z}^{n}$ are
denoted as $\textbf{k}=(k_{1},k_{2},\cdots,k_{n})$. The standard distance
function on $\mathbb{Z}^{n}$ is denoted by dist.

Construction of good and bad cubes. Take $\delta_{0}$ to be the small constant
in Lemma 2.2 with $\varepsilon=\varepsilon_{0}$. Let

$\displaystyle\mathcal{B}$
$\displaystyle=\left\\{\textbf{k}\in\mathbb{Z}^{n}:~{}~{}\int_{Q_{1}(\textbf{k})}(|\nabla^{2}u|^{2}-|\nabla|\nabla
u||^{2})dx\geq\delta_{0}\right\\},$ $\displaystyle\mathcal{G}$
$\displaystyle=\left\\{\textbf{k}\in\mathbb{Z}^{n}:~{}~{}\int_{Q_{1}(\textbf{k})}(|\nabla^{2}u|^{2}-|\nabla|\nabla
u||^{2})dx<\delta_{0}\right\\}.$

The corresponding bad and good cubes form the sets

$\widehat{\mathcal{B}}=\mathop{\cup}\limits_{\textbf{k}\in\mathcal{B}}Q_{1}(\textbf{k}),\quad\widehat{\mathcal{G}}=\mathop{\cup}\limits_{\textbf{k}\in\mathcal{G}}Q_{1}(\textbf{k}).$

We first assume that $\mathcal{B}\neq\emptyset$.

For each $R\in\mathbb{N}$, define

$\displaystyle N(R)$ $\displaystyle:=\sharp\left(\mathcal{B}\cap
Q_{R}(0)\right),$ $\displaystyle D(R)$
$\displaystyle:=\int_{Q_{R}(0)}\left(|\nabla^{2}u|^{2}-|\nabla|\nabla
u||^{2}\right),$ $\displaystyle E(R)$ $\displaystyle:=\int_{Q_{R}(0)}|\nabla
u|^{2}.$

By Lemma 2.1, we have the following estimates for these three quantities as
the starting point:

$\left\\{\begin{aligned} &N(R)\leq CR^{n-2},\\\ &D(R)\leq CR^{n-2},\\\
&E(R)\leq CR^{n}.\end{aligned}\right.$ (2.10)

These three quantities are actually intimately related to each other. Indeed,
Lemma 2.1 gives

$D(R)\leq CR^{-2}E(2R).$ (2.11)

On the other hand, by a counting of bad cubes, we have

$N(R)\leq CD(R).$ (2.12)

To close the loop, we need one more recursive relation to control $E$ by $N$ .
This is the content of the following

######  Proposition 2.3.

Let $p$ be the constant defined in (H2). For any $R$ large, $10<L\ll R$, we
have

$E(R)\leq CL^{n}N(2R)+CL^{-2\frac{p+1}{p-1}}R^{n},$ (2.13)

where the constant $C$ does not depend on $L$ and $R$.

To prove this proposition, we need a decay estimate of $|\nabla u|$ in the
good set $\widehat{\mathcal{G}}$. This is the content of the following

######  Lemma 2.4.

There exists a constant $C$ such that for any $\textbf{k}\in\mathcal{G}$, if
$\text{dist}\left(\textbf{k},\mathcal{B}\right)>10$, then

$\sup_{Q_{1}(\textbf{k})}|\nabla u|\leq
C\text{dist}\left(\textbf{k},\mathcal{B}\right)^{-\frac{p+1}{p-1}}.$ (2.14)

###### Proof.

Let us set $D:=\text{dist}\left(\textbf{k},\mathcal{B}\right)/2$. By Lemma
2.2, for any $x\in Q_{D}(\textbf{k})$, there exists a stable zero $a_{x}$ of
$f$ such that

$\sup_{Q_{1}(x)}|u-a_{x}|\leq\varepsilon_{0}.$

Because $u$ is continuous and zeros of $f$ are isolated, this zero is the same
one for all these $x$. In other words, there exists a stable zero $a$ of $f$
such that

$\sup_{Q_{D}(\textbf{k})}|u-a|\leq\varepsilon_{0}.$

Now let us define

$v(x)=u(x)-a,\quad\text{in}\quad Q_{D}(\textbf{k}).$

For any $x\in Q_{D}(\textbf{k})$, if $u(x)>a$, then by (H2),

$\Delta v(x)=\Delta u(x)=-f(v(x)+a)\geq c_{1}|v(x)|^{p}.$

In the same way, if $u(x)<a$, then

$-\Delta v(x)=-\Delta u(x)=f(v(x)+a)\geq c_{1}|v(x)|^{p}.$

Summing up, we conclude that $v$ satisfies

$\Delta|v|\geq c_{1}|v|^{p},\quad\text{in}\quad Q_{D}(\textbf{k}).$ (2.15)

Then by the Keller-Osserman theory (see [17, 18]) and a standard blow up
argument using the doubling lemma (see [20] ), we can find a constant $C$,
which depends only on $p$, such that

$\sup_{Q_{D/2}(\textbf{k})}|v|\leq CD^{-\frac{2}{p-1}}.$ (2.16)

The estimate (2.14) then follows from applying standard interior gradient
estimate to $u$ in $Q_{D/2}(\textbf{k})$. ∎

######  Remark 2.5.

If $p=1$, (2.14) should be replaced by

$\sup_{Q_{1}(\textbf{k})}|\nabla u|\leq
Ce^{-c\text{dist}\left(\textbf{k},\mathcal{B}\right)}.$ (2.17)

This follows from a similar a priori estimate for the $p=1$ case of (2.15).

###### Proof of Proposition 2.3.

We consider the good and bad set separately.

First of all, let us use $\mathcal{B}_{L}$ to denote the $L$ neighborhood of
$\mathcal{B}$ in the lattice $\mathbb{Z}^{n}$, and $\widehat{\mathcal{B}}_{L}$
will represent the corresponding domain in $\mathbb{R}^{n}$. We have

$\left|\widehat{\mathcal{B}}_{L}\cap
Q_{2R}(0)\right|\leq\sharp\left(\mathcal{B}_{L}\cap Q_{2R}(0)\right)\leq
CL^{n}N(2R).$ (2.18)

Then combining (2.18) with the gradient bound (2.2), we get

$\int_{\widehat{\mathcal{B}}_{L}\cap Q_{R}(0)}|\nabla u|^{2}\leq CL^{n}N(2R).$
(2.19)

Next, by Lemma 2.4, for any $\textbf{k}\notin\mathcal{B}_{L}$,

$\sup_{Q_{1}(\textbf{k})}|\nabla u|\leq CL^{-\frac{p+1}{p-1}}.$

Therefore

$\int_{Q_{R}(0)\setminus\widehat{\mathcal{B}}_{L}}|\nabla u|^{2}\leq
CL^{-2\frac{p+1}{p-1}}|Q_{R}(0)|\leq CL^{-2\frac{p+1}{p-1}}R^{n}.$ (2.20)

Summing up (2.19) and (2.20), we obtain (2.13). ∎

With all these preparations, we are ready to prove our main result.

###### Proof of Theorem 1.4.

If $\mathcal{B}=\emptyset$, a direct application of Lemma 2.4 shows that $u$
is a constant. Hence we assume this is not the case, and we want to establish
a quadratic growth bound on $E(R)$ by an iteration with (2.10) as a starting
point.

We will assume without loss of generality that $p>1$. Combining (2.11), (2.12)
and (2.13), we get

$E(R)\leq CL^{n}R^{-2}E(4R)+CL^{-2\frac{p+1}{p-1}}R^{n}.$ (2.21)

Let us define $\alpha_{0}=0$ and

$\alpha_{k+1}=\left(\alpha_{k}+2\right)\frac{2(p+1)}{n(p-1)+2(p+1)}.$ (2.22)

Since $p\geq 1$, a simple induction argument tells us that the sequence
$\\{\alpha_{k}\\}$ is increasing and

$\alpha_{k}\rightarrow\frac{4(p+1)}{n(p-1)}\quad\text{as}\quad
k\rightarrow+\infty.$ (2.23)

We claim that for any $k=0,\dots$, there exist $C_{k}$ such that

$E(R)\leq C_{k}R^{n-\alpha_{k}},\quad\forall R>1.$ (2.24)

The case of $k=0$ is just(2.10). To see that (2.24) is indeed true for all
$k$, we simply choose

$L=R^{\frac{(2+\alpha_{k})(p-1)}{n(p-1)+2(p+1)}}$

in (2.21) and use an induction argument again.

Now if $p<p^{*}(n)$, then we can check directly that

$\frac{4(p+1)}{n(p-1)}>n-2.$

Therefore, in view of (2.23), there exists an index $k_{0}$ such that
$\alpha_{k_{0}}\geq n-2$. It then follows from (2.24) that there exists a
$\beta\leq 2$ such that

$E(R)\leq CR^{\beta},\quad\forall R>1.$

With this estimate at hand, taking a standard log cut-off function in the
Sternberg-Zumbrun inequality (as in [1] or [15]), we deduce that $u$ is one
dimensional. But since we have assumed (H3), $u$ must be a constant. ∎

With Theorem 1.4 at hand, we can prove Corollary 1.6 in the same spirit as in
[10].

###### Proof of Corollary 1.6.

By the same proof of [10, Theorem 4], we know that the function
$v:=\lim_{x_{n}\to+\infty}u$ is a positive bounded stable solution of the
equation

$-\Delta v=f(v),\quad\text{in}\quad\mathbb{R}^{n-1}.$

Applying Theorem 1.4, we find that $v$ is a constant. Then still following the
proof of [10, Theorem 4], we conclude that $u$ is one dimensional. ∎

## References

  * [1] L. Ambrosio and X. Cabré.  Entire solutions of semilinear elliptic equations in $\mathbb{R}^{3}$ and a conjecture of De Giorgi.  J. Amer. Math. Soc., 13:725–739, 2000.
  * [2] H. Berestycki, L. Caffarelli, and L. Nirenberg.  Symmetry for elliptic equations in a half space, In: Boundary value problems for partial differential equations and applications.  pages 27–42, 1993.
  * [3] Henri Berestycki and Cole Graham.  Reaction-diffusion equations in the half-space.  Ann. Inst. H. Poincaré C Anal. Non Linéaire, 39(5):1053–1095, 2022.
  * [4] Henri Berestycki and Cole Graham.  The steady states of positive reaction-diffusion equations with dirichlet boundary conditions.  arXiv:2309.16642, 2023.
  * [5] H. Brezis.  Is there failure of the inverse function theorem? Morse theory, minimax theory and their applications to nonlinear differential equations.  Int. Press, Somerville, 2003.
  * [6] X. Cabré and C. Antonio.  Regularity of radial minimizers and extremal solutions of semilinear elliptic equations.  J. Funct. Anal., 238(2):709–733, 2005.
  * [7] X. Cabré, A. Figalli, Ros-Oton, and J. Serra.  Stable solutions to semilinear elliptic equations are smooth up to dimension 9.  Acta Math., 224:187–252, 2019.
  * [8] E. Dancer.  Finite Morse index and linearized stable solutions on bounded and unbounded domains,In Proceedings of the International Congress of Mathematicians. Volume III.  Hindustan Book Agency, New Delhi, 2010.
  * [9] M. del Pino, M. Kowalczyk, and J. Wei.  On De Giorgi’s conjecture in dimension $n\geq 9$.  Ann. of Math., 174(3):1485–1569, 2011.
  * [10] L. Dupaigne and A. Farina.  Classification and Liouville-type theorems for semilinear elliptic equations in unbounded domains.  Anal. PDE, 15(2):551–566, 2022.
  * [11] A. Farina and L. Dupaigne.  Liouville theorems for stable solutions of semilinear elliptic equations with convex nonlinearities.  Nonlinear Anal., 70(8):2882–2888, 2009.
  * [12] A. Farina, Y. Sire, and E. Valdinoci.  Stable solutions of elliptic equations on Riemannian manifolds.  J. Geom. Anal., 23:1158–1172, 2008.
  * [13] Alberto Farina.  Some results about semilinear elliptic problems on half-spaces.  Math. Eng., 2(4):709–721, 2020.
  * [14] Alberto Farina and Berardino Sciunzi.  Monotonicity and symmetry of nonnegative solutions to $-\Delta u=f(u)$ in half-planes and strips.  Adv. Nonlinear Stud., 17(2):297–310, 2017.
  * [15] N.A. Ghoussoub and C. Gui.  On a conjecture of De Giorgi and some related problems.  Math. Ann., 311(3):481–491, 1998.
  * [16] D. Gilbarg and N. S. Trudinger.  Elliptic Partial Differential Equations of Second Order, 2nd ed. Springer-Verlag, New York, 1983.
  * [17] J. Keller.  On solutions of ${\Delta}u=f(u)$.  Comm Pure Appl Math., 10(4):503–510, 1957.
  * [18] Robert Osserman.  On the inequality ${\Delta}u{\geq}f(u)$.  Pacific J. Math., 7(4):1641–1647, 1957.
  * [19] Phuong P. Le and B. Rahal.  On stable and finite morse index solutions to quasilinear Schrödinger equations.  NoDEA Nonlinear Differential Equations Appl., 29(4):19pp, 2022.
  * [20] P. Quittner P. Poláčik and P. Souplet.  Singularity and decay estimates in superlinear problems via Liouville-type theorems, I: Elliptic equations and systems.  Duke Math. J., 139(3):555–579, 2015.
  * [21] O. Savin.  Regularity of flat level sets in phase transitions.  Ann. of Math., 169(1):41–78, 2009.
  * [22] A. Selmi, A. Harrabi, and C. Zaidi.  Nonexistence results on the space or the half space of $-{\Delta}u+{\lambda}u={|u|^{p-1}}u$ via the Morse index.  Commun. Pure Appl. Anal., 19(5):2839–2852, 2020.
  * [23] P. Sternberg and K. Zumbrun.  On the connectivity of boundaries of sets minimizing perimeter subject to a volume constraint.  Comm. anal. geom., 7(1):199–220, 1999.
  * [24] K. Zumbrun and P. Sternberg.  A Poincaré inequality with applications to volume-constrained area-minimizing surfaces.  J. Reine Angew. Math., 1998(503), 1998.