# NeUDF: Learning Unsigned Distance Fields from Multi-view Images for
Reconstructing Non-watertight Models

Fei Hou
Institute of Software, CAS
<EMAIL_ADDRESS>Equal contributions. Junkai Deng11footnotemark: 1
Institute of Software, CAS
<EMAIL_ADDRESS>Xuhui Chen
Institute of Software, CAS
<EMAIL_ADDRESS>Wencheng Wang
Institute of Software, CAS
<EMAIL_ADDRESS>Ying He
Nanyang Technological University
<EMAIL_ADDRESS>
###### Abstract

Volume rendering-based 3D reconstruction from multi-view images has gained
popularity in recent years, largely due to the success of neural radiance
fields (NeRF). A number of methods have been developed that build upon NeRF
and use neural volume rendering to learn signed distance fields (SDFs) for
reconstructing 3D models. However, SDF-based methods cannot represent non-
watertight models and, therefore, cannot capture open boundaries. This paper
proposes a new algorithm for learning an accurate unsigned distance field
(UDF) from multi-view images, which is specifically designed for
reconstructing non-watertight, textureless models. The proposed method, called
NeUDF, addresses the limitations of existing UDF-based methods by introducing
a simple and approximately unbiased and occlusion-aware density function. In
addition, a smooth and differentiable UDF representation is presented to make
the learning process easier and more efficient. Experiments on both texture-
rich and textureless models demonstrate the robustness and effectiveness of
the proposed approach, making it a promising solution for reconstructing
challenging 3D models from multi-view images.

##  1 Introduction

|  |
---|---|---
Reference Image | Ours | NeuS [33]
Figure 1: We show a reconstruction result using our method and NeuS [33] on a
mostly blue short trouser with little texture. Our method successfully
reconstructs the structure of the trouser, while NeuS fails to learn the
correct structure of the object, and instead represents it as a closed model.

In recent years, volume rendering-based 3D model and scene reconstruction from
multi-view images has gained popularity, largely due to the success of neural
radiance fields (NeRF) [27]. A number of methods have been developed that
build upon NeRF and use neural volume rendering to learn signed distance
fields (SDFs) for reconstructing 3D models [37, 33, 35, 7]. While these SDF-
based methods are effective in reconstructing watertight models, they cannot
represent non-watertight models, as SDFs distinguish between the interior and
exterior of a model, and therefore cannot capture open boundaries.

Recent research has attempted to address the limitation of signed distance
fields in reconstructing non-watertight models by using unsigned distance
fields (UDFs). For instance, NeuralUDF [22] extends NeuS [33] to learn UDFs
for reconstructing open models. However, UDF-based methods have their own
limitations. For example, as shown in Figure 1, NeuralUDF struggles to
reconstruct textureless models or models with few distinguishable features, as
it relies on texture information to improve the accuracy of reconstruction.
Also, UDF is not differentiable at the zero level-set, making it difficult to
learn. Reconstructing open, textureless models from multi-view images using
volume rendering remains a challenging problem.

We propose a new algorithm for learning an accurate UDF from multi-view
images, which we call NeUDF. Our method is specifically designed for
reconstructing non-watertight, textureless models, which are particularly
challenging to reconstruct using existing UDF-based methods. The key to
learning an accurate UDF is to design an unbiased and occlusion-aware density
function of the UDF. Since the inside and outside of the model are not
distinguished, the S-shape density function used in SDF-based methods [33] is
not appropriate for UDFs. Long et al. [22] proposed a piecewise density
function adapted from NeuS [33], but this function is rather complicated and
difficult to learn, limiting its applicability to texture-rich models.

To overcome these limitations, we propose a new UDF density function for
learning UDFs, which can tolerate a small amount of bias to ease learning. Our
density function is dense enough to be almost opaque, making it occlusion-
aware and effective for reconstructing non-watertight models. Additionally, we
present a smoothed UDF representation that makes the UDF differentiable at the
zero level set, which is crucial for the learning process. We evaluate our
proposed approach on both texture-rich and textureless models from the
DeepFashion3D [41] and DTU [17] datasets. Our experiments demonstrate that our
method is robust and can effectively learn UDFs for both types of inputs.

Our work makes the following contributions:

  1. 1.

We propose a method for reconstructing non-watertight models based on NeRF
that is capable of reconstructing both texture-rich and challenging
textureless models.

  2. 2.

We introduce a theoretically sound density function for UDFs and in practice
adopt a variant that allows a small bias and is approximately occlusion-aware.

  3. 3.

We present a smooth and differentiable UDF representation that makes the
learning process easier and more efficient.

  4. 4.

We present a simple yet effective method for extracting the target surface
from the learned UDFs, which can reduce the bias.

##  2 Related Work

###  2.1 3D Reconstruction from Multi-View Images

Surface reconstruction from multi-view images has been a subject of study for
several decades, and can generally be classified into two categories: voxel-
based and point-based methods. Voxel-based methods [2, 3, 18, 19, 32] divide
the 3D space into voxels and determine which ones belong to the object. These
methods can be computationally expensive and may not be suitable for
reconstructing complex surfaces. Point-based methods [12, 30, 36] use
structure-from-motion [15] to calibrate the images and generate a dense point
cloud using multi-view stereo [11]. Finally, surface reconstruction methods
(e.g., [1, 20, 16]) are used to generate a mesh. Since multi-view stereo
requires dense correspondences to generate a dense point cloud, which are
often difficult to compute, its results often contain various types of
artifacts, such as noise, holes and incomplete structures.

Neural network-based 3D surface reconstruction has received attention in
recent years with the emergence of neural rendering [27]. Several methods have
been proposed for volume rendering and surface reconstruction using neural
networks. VolSDF [37] uses the cumulative distribution function of Laplacian
distribution to evaluate the density function from SDF for volume rendering
and surface reconstruction. NeuS [33] adopts an unbiased density function for
SDFs for more accurate reconstruction. SparseNeuS [23] extends NeuS to use
fewer images for reconstruction. HF-NeuS [35] improves NeuS by proposing a
simplified and unbiased density function and using hierarchical MLPs for
detail reconstruction. Geo-NeuS [9] incorporates structure-from-motion to add
more constraints. NeuralWarp [7] improves the accuracy by optimizing
consistency between warped views of different images. All of these methods
learn SDFs, which can only reconstruct watertighjt models.

More recently, Long _et al_. proposed NeuralUDF [22] for learning UDF for
reconstructing open models. It adapts the density function of NeuS to UDFs by
introducing an indicator function. However, this method can only learn
texture-rich models due to the complex density function used in training. In
contrast, our proposed UDF learning method is capable of reconstructing both
texture-rich and textureless models without the need for masks.

###  2.2 3D Reconstruction from Point Clouds

There has been recent interest in surface representation using signed distance
fields (SDFs) and occupation fields. Several methods have been proposed for
learning SDFs [28, 4, 31, 24, 34], while occupation fields have been used in
methods such as [26, 5]. However, both SDFs and occupation fields can only
represent watertight models.

To represent non-watertight models, some methods are proposed to learn UDF
from 3D point clouds [6, 39, 40]. Our proposed method also uses UDF for non-
watertight models representation, but we learn it directly from multi-view
images, which is a challenging problem.

##  3 Method

Our goal is to extract a set of surface points from a set of input images
along with their corresponding known camera poses and intrinsic parameters. To
achieve this, we first learn an implicit unsigned distance field, whose zero
level set represents the 3D surface. Unlike signed distance fields, UDFs have
positive distance values, making them suitable for representing non-watertight
models. However, learning a UDF from multi-view images is a challenging task
that requires a proper density function balancing occlusion-awareness, ease of
learning, and unbiasedness. Another technical issue is that UDFs are not
differentiable at the zero level-set, making it difficult to train. Our
proposed method is designed to address the above challenges of learning UDFs
from multi-view images.

In this section, we first review the key concepts of volume rendering and
neural radiance fields. Next, we propose an unbiased density function and use
its variant to balance the density properties for UDF learning. We then
introduce a differentiable form of UDF that can be used for stable learning.
Finally, we detail our loss function design.

###  3.1 Review of Volume Rendering

Volume rendering is a crucial component of neural radiance fields [27]. During
volume rendering, a ray $\mathbf{r}$ is cast from the camera position
$\mathbf{o}$ to each pixel of the virtual canvas. The direction of the ray is
denoted by $\mathbf{d}$, and any point $\mathbf{r}(t)$ along the ray can be
expressed as $\mathbf{o}+t\mathbf{d}$, where $t$ is the arc-length parameter.
In the following, we refer to the point $\mathbf{r}(t)$ using the parameter
$t$ if there is no confusion.

The color of the corresponding pixel is determined by the color of the ray,
which is the weighted integral of all colors $c(t)$ along the ray $\bf r$,
$c(\mathbf{r})=\int_{0}^{\infty}w(t)c(t)\,\mathrm{d}t.$

The integral is often computed numerically using quadrature, which is a sum of
color values at a set of sampled points $t_{1},t_{2},\cdots,t_{n}$ on the ray,
i.e., $c(\mathbf{r})=\sum_{i=1}^{n}w(t_{i})c(t_{i}).$

The weight $w(t)$ of each color value is the product between the volume
density $\sigma$ and transparency $T$ [25] $w(t)=T(t)\sigma(t)$, which
measures the accumulated transmittance and occlusion of the ray up to that
point. The transparency $T$, which is defined as
$T(t)=\exp\left(-\int_{0}^{t}\sigma(u)\,\mathrm{d}u\right),$ is the
probability that the ray travels from $0$ to $t$ without hitting any other
particle. $T(t)$ is a monotonic decreasing function with a starting value of
$T(0)=1$.

In NeRF [27], the volume density $\sigma$ is directly computed by a multi-
layer perceptron (MLP), which takes the 3D coordinates of a point as input and
outputs the corresponding density value. This allows for a flexible
representation of the volume density field and enables the synthesis of novel
views from arbitrary viewpoints. However, since the volume density field does
not provide explicit information about the surface geometry, extracting a
high-quality surface from it is difficult. Moreover, the density field may not
correspond to a valid surface at all due to the ill-posed nature of the
reconstruction problem [38].

To overcome this limitation, NeuS [33] utilizes signed distance functions to
represent 3D geometry. Specifically, it uses a scalar function $f$ that maps
each point in 3D space to the signed distance from that point to the surface,
where the zero level set of the function represents the target surface. It
then maps the signed distance function $f$ to an S-density function
$\phi_{s}(f)$ using a logistic density distribution
$\phi_{s}(x)=se^{-sx}/(1+e^{-sx})^{2}$, which assigns a normalized density
value to each point in space based on its signed distance from the surface.
The normalized S-density function is then used to define the weight function
$w(t)$ for volume rendering.

Wang et al. [33] proved that the weight function $w(t)$ in NeuS is both
unbiased and occlusion-aware. Being unbiased means that the weight function
attains a locally maximal value at a surface intersection point, while an
occlusion-aware weight function implies that when two points have the same SDF
value, the point closer to the viewpoint has a larger contribution to the
output color than the other point. This accounts for the occlusion effect,
where points behind the surface are not visible and should not contribute to
the output color.

###  3.2 Density Function

A core part of NeUDF is an unbiased and occlusion-aware weight function.
Inspired by HF-NeuS [35], we propose a bell-shaped weight function that maps
unsigned distance to density. Our weight function is _unbiased_ in that the
weights are maximum on the surface.

(a) (b)

Figure 2: Weight function. (a) A biased weight function can have its maximum
at a point away from the surface, reducing the accuracy of the reconstructed
surfaces. (b) In contrast, an unbiased weight function attains its local
maximum at the point where the ray intersects the surface. The use of an
unbiased weight function is essential for accurate surface reconstruction in
volume rendering.

######  Theorem 1.

A ray $\mathbf{r}=\mathbf{o}+t\mathbf{d}$ hits a planar object $M$. Let
$\theta$ be the angle between the UDF gradient and the ray direction
$\mathbf{d}$. Define the bell-shaped density $\sigma$ as

$\sigma(t)=\frac{se^{-sf(t)}}{1+e^{-sf(t)}}|\cos(\theta)|,$ (1)

where $s>0$ is a learnable parameter. Then the weight function
$w(t)=T(t)\sigma(t)$ is unbiased.

###### Proof.

Let $\mathbf{r}(t_{0})$ be the intersection point.

We first consider the transparency $T(t)$ with $t<t_{0}$, i.e., for any point
along the ray up to $\mathbf{r}(t_{0})$. Since the point $t$ is in front of
$M$, we have $90^{\circ}<\theta\leq 180^{\circ}$ and $\cos(\theta)<0$. The
transparency $T(t)$, $t<t_{0}$, is given by

$\displaystyle T(t)$
$\displaystyle=\exp\left(-\int^{t}_{0}\frac{se^{-sf(u)}}{1+e^{-sf(u)}}(-\cos(\theta))\mathrm{d}u\right)$
$\displaystyle=\exp\left(\int^{t}_{0}\frac{se^{-sf(u)}}{1+e^{-sf(u)}}\mathrm{d}f(u)\right)$
$\displaystyle=\exp\left(\ln(1+e^{-sf(0)})-\ln(1+e^{-sf(t)})\right)$
$\displaystyle=\frac{1+e^{-sf(0)}}{1+e^{-sf(t)}}.$

Next we compute the weight function $w(t)$, which is proportional to

$\displaystyle w(t)=T(t)\sigma(t)$
$\displaystyle\propto\frac{e^{-sf(t)}}{(1+e^{-sf(t)})^{2}}.$

It is straightforward to verify that the sign of
$\frac{\mathrm{d}w}{\mathrm{d}f}$ is the same as $e^{-sf(t)}-1$. Using the
chain rule, we can show that the derivative $\mathrm{d}w/\mathrm{d}t>0$ is
positive for $t<t_{0}$, implying that the weight $w$ increases as the ray
approaches $M$ from the front.

After the ray passes through $M$, $\cos\theta$ becomes positive and both
$T(t)$ and $\sigma(t)$ are monotonically decreasing for $t>t_{0}$. Therefore,
the product $T(t)\sigma(t)$ decreases as the ray leaves $M$.

Since the weight is continuous, $w$ attains its maximum on the plane $M$,
implying that it is unbiased. ∎

The density function $\sigma$ in Equation (1) is theoretically sound,
satisfying the requirements of unbiased sampling. However, in practice, it may
be too “transparent” to use (due to its bell-shaped geometry), leading to poor
sampling efficiency. Moreover, the dependence of the bell-shaped function on
the UDF gradients, as indicated by the $\cos(\theta)$ term, can introduce
instability and sensitivity to noise and oscillations, thereby posing
challenges in learning the function.

To address these issues, we replace the $\cos(\theta)$ term with a constant
$c>1$, resulting in a modified density function

$\hat{\sigma}(t)=\frac{cse^{-sf(t)}}{1+e^{-sf(t)}},\;s>0,\;c>1.$

This modification increases both numerical stability and opacity of the
original density function $\sigma$.

However, unlike the original density $\sigma$, the modified density
$\hat{\sigma}$ introduces bias in the weight function, since the maximum value
of weight occurs at a point $t^{*}$ in front of $M$, which has a distance
value

$f(t^{*})=\frac{1}{s}\ln\frac{-c}{\cos(\theta)}.$ (2)

While the modified density function $\hat{\sigma}$ is not theoretically
unbiased, setting $c$ as a small constant can greatly reduce the bias in
practice. Additionally, we expect $\hat{\sigma}$ to be approximately
occlusion-aware. To further understand the numerical properties of
$\hat{\sigma}$, we consider the extreme case where the incident light ray is
perpendicular to the planar surface $M$. In this case, the unsigned distance
function is $f(t)=1-t$ for points in front of $M$. As $\hat{\sigma}$ is
symmetric for the two sides of $M$, the surface transparency is the square of
the transparency of the front side. Our computation shows

$\displaystyle\left(e^{-\int_{0}^{1}\hat{\sigma}(u)\mathrm{d}u}\right)^{2}$
$\displaystyle=\left[\exp\left(-\int^{1}_{0}\frac{cse^{-s(1-u)}}{1+e^{-s(1-u)}}\mathrm{d}u\right)\right]^{2}=\left(\frac{1+e^{-s}}{2}\right)^{2c},$

and to reduce transparency, we should choose a relatively large $c$.

In our implementation, we set the constant $c=5$ based on the typical value of
the learned parameter $s$, which ranges between $1000$ and $2000$ given that
the models are learned in a unit sphere to become unitless. This enables us to
estimate the upper bound of the bias. For points in front of the surface, the
incident angle $\theta$ between the ray and the surface normal is obtuse, so
we restrict $\theta$ to the range of $[91^{\circ},180^{\circ}]$. By setting
$c=5$, the offset width between 0.00161 and 0.00566 is obtained relative to
the true zero level set, indicating that the maximum relative bias is below
0.5%. This error level is acceptable for most application scenarios. Moreover,
the surface transparency in the extreme case mentioned above is less than
$0.001$. When a ray has a larger incident angle, its transparency becomes even
smaller, resulting in an almost opaque density $\hat{\sigma}$. As a result,
the weight function $w$ is approximately occlusion-aware. Thus, setting the
constant $c=5$ offers a good balance between occlusion-awareness and
unbiasedness.

It is worth mentioning that our modified density function $\hat{\sigma}$ is
much simpler than the density function used in NeuralUDF [22]. Although both
densities introduce biases, we can control the bias of our density function
within a very small range. In contrast, it is unclear what the range of bias
is for the density function used in NeuralUDF. Thanks to the simplicity, our
density function is easier to learn, thereby our method can work for a wider
range of models. This includes textureless models, which pose challenges for
NeuralUDF.

###  3.3 Training

Differentiable UDFs. NeuS uses an MLP network to learn the signed distance
function $f$, which is a differentiable function. In contrast, UDF is not
differentiable at the zero level set, making the network difficult to learn
the values and gradients of the UDF close to the zero level set.

Another crucial requirement is to ensure non-negative values for the computed
distances, which seems a trivial task as one may simply apply absolute value
or normalization such as ReLU to the MLP output. However, applying the
absolute value to the distance is not viable due to its non-differentiability
at zero. Similarly, normalizing the output value using ReLU is not feasible as
it is also non-differentiable at zero, and its gradient vanishes for negative
inputs. This can be particularly problematic for learning UDFs, since when the
MLP returns a negative distance value, the ReLU gradient vanishes, hindering
the update of the distance to a positive value in the subsequent iterations.

We add a softplus function after the output layer of the MLP. The softplus
function [8] is a smooth and differentiable approximation of the ReLU function
that is defined as

$\text{softplus}(x)=\frac{1}{\beta}\ln(1+e^{\beta x}).$

Softplus has the same shape as ReLU, but it is continuous and differentiable
at every point, and its gradients do not vanish anywhere. Using the softplus
function allows us to ensure that the output of the MLP is non-negative and
differentiable, making it suitable for learning the UDF. We set $\beta=100$ in
our experiments.

Loss functions. Following NeuralUDF [22], we adopt an iso-surface regularizer
to penalize the UDF values of the non-surface points from being zero,
therefore encouraging smooth and clean UDFs. The regularization loss is
defined as [22]

$\mathcal{L}_{reg}=\frac{1}{MN}\sum_{i,k}\exp{\left(-\tau\cdot
f(t_{i,k})\right)},$

where $\tau=5.0$ is a constant scalar that scales the learned UDF values, $M$
is the total number of sampled rays per training iteration, and $N$ is the
number of sampled points on a single ray.

The value of $s$, which is learnable in our method, significantly affects the
quality of the reconstruction. When $s$ is small, it introduces a larger bias
and leads to a more blurred output. We observe that $s$ typically converges to
a relatively large value between 1000 and 2000, leading to visually pleasing
results. However, in rare cases when $s$ stops increasing during training, we
apply a penalty to force it to increase. The penalty is defined as follows

$\mathcal{L}_{s}=\frac{1}{M}\sum_{i,k}\frac{1}{s_{i,k}},$

where $M$ is the number of rays during a training epoch. This term
$\mathcal{L}_{s}$ aggregates the reciprocals of all $s$ values used for the
point $t_{i,k}$ on ray $r_{i}$. Intuitively speaking, it encourages a larger
$s$ during the early stage of training. In our implementation, we make this
term optional since $s$ generally increases with a decreasing rate during
training, and the penalty term is only necessary in rare cases when $s$ stops
at a relatively low value.

As in other SDF- and UDF-based methods [33, 35, 22], we adopt color loss and
Eikonal loss in our approach. Specifically, the color loss
$\mathcal{L}_{color}$ is the $L_{1}$ loss between the predicted color and the
ground truth color of a single pixel as used in [33]. The Eikonal loss
$\mathcal{L}_{eik}$ is used to regularize the learned distance field to have a
unit gradient [13]. Putting it all together, we define the combined loss
function as a weighted sum

$\mathcal{L}=\mathcal{L}_{color}+\lambda_{1}\mathcal{L}_{eik}+\lambda_{2}\mathcal{L}_{reg}+\lambda_{3}\mathcal{L}_{s},$

where $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$ are hyperparameters that
control the weight of each loss term. In our experiments, we empirically set
$\lambda_{1}=0.1$, $\lambda_{2}=0.01$, and $\lambda_{3}=0.001$, although
$\lambda_{2}$ is occasionally set to $0.02$, and $\lambda_{3}$ is optional.

Adaptive hierarchical sampling. We propose an adaptive hierarchical sampling
(AHS) strategy for efficiently finding sample points along each ray for color
computation. Existing approaches, such as NeuS [33] and NeuralUDF [22], adopt
a fixed set of pre-defined sampling rates for all models, regardless of their
geometry differences. They first uniformly sample 64 points on the ray and
then iteratively conduct importance sampling on top of the coarse probability
estimation. The sampling rate is doubled after each iteration, resulting in a
maximum 512 samples on each ray.

We argue that an effective sampling strategy should be self-adaptive to
geometry, which is connected to the learnable parameter $s$. Essentially, a
larger value of $s$ leads to a narrower bell-shaped density function,
resulting in higher weights given to points closer to the surface and more
concentrated samples on the surface. Our density function depends on the
learnable parameter $s$, which usually ranges between 1000 and 2000. We use
the learned $s$ to guide more precise sample selection. Specifically, in the
$i$-th iteration of sampling, we set $s$ to $\max\left(32\times
2^{i},\frac{s}{2^{k-i}}\right)$, where $k$ is the maximum number of sampling
iterations, which we set as $k=4$ in our implementation. During the early
stage of training, when $s$ is still relatively small, we use the sampling
rates as in NeuS and NeuralUDF. As the training progresses and the value of
$s$ consistently increases, our adaptive sampling strategy is activated. Our
ablation study confirms that adaptive hierarchical sampling increases the
accuracy of our approach by up to 25%. See Section 4.3.

###  3.4 Surface Extraction

After learning the UDF $f$, extracting the object surface is essential.
Instead of extracting the zero level set of the function $f$, we locate the
sample points with the maximum weights to construct the surface. The reason
for doing so is that extracting the zero level set from UDFs is notoriously
unstable [14], often leading to artifacts such as holes, flipped triangles,
and self-intersections in the output meshes.

The main technical challenge in locating the samples with maximum weights is
that they depend on the view direction, due to the $\cos\theta$ term in
Equation (2). However, since our method penalizes small values of $s$ and
yields a relatively large $s\in[1000,2000]$, the $\ln\frac{-c}{\cos\theta}$
term has little effect to the distance value $f(t^{*})$.

In our implementation, we obtain a subset of rays for each input image by
uniformly sampling rays every $k=5$ pixels in both horizontal and vertical
directions. For each ray $\bf r$, we sample points along it to find the
location $t^{*}$ with the maximum weight. We classify the ray $\bf r$ as a
foreground ray if the sum of weights for the samples along $\bf r$ inside the
region of interest (a unit sphere centered at origin) is greater than 0.5.
Otherwise, we consider $\bf r$ as a background ray that misses the target
object.

After collecting all foreground rays from all input images, we use the sample
with the maximum weight for each ray as a surface point. Computational results
show that our strategy is effective in producing satisfactory results (see
Section 4.3).

##  4 Experiments

_Implementation Details._ The MLP for the UDF network consists of 8 hidden
layers, each with 256 elements. We also use skip connections after every 4
hidden layers. The output of the UDF network is a single value representing
the predicted UDF and a 256-dimensional feature vector used in the color
network.

For the color network, we use another MLP with 4 hidden layers, each having
256 elements. We use the coarse-to-fine strategy proposed by Park _et al_.
[29] for position encoding, setting the maximum number of frequency bands to
16 for the UDF network and 6 for the color network. For background rendering,
we use NeRF++ [38] for background prediction. During training, we use the Adam
optimizer [21] with a global learning rate of 5e-4. We sample 512 rays per
batch and train our model for 300,000 iterations.

_Data Sets._ To evaluate our method, we use two datasets: DeepFashion3D [41]
and DTU [17]. The DeepFashion3D dataset consists of clothing models, which are
open models with boundaries. As only 3D points are available, we render 72
images of resolution $1024\times 1024$ with a white background from different
viewpoints for each model. The DTU dataset consists of models captured in a
studio, and all the models are watertight. We use this dataset to validate
that our method also works well for watertight models. These datasets have
been widely used in previous works such as [37, 33, 35].

_Baselines._ Several methods have been proposed for learning signed distance
functions (SDFs) from multi-view images, which generate watertight models. To
validate the effectiveness of our proposed method, we compare with state-of-
the-art methods, namely VolSDF [37], NeuS [33], and HF-NeuS[35].

To our knowledge, NeuralUDF [22] 111Since the source code of NeuralUDF is
unavailable at the moment of submission, we cannot make quantitative
comparison with it in the paper. is the only method designed for open model
reconstruction, but it is limited to texture-rich models. Our method, on the
other hand, can reconstruct both textureless and texture-rich models thanks to
the new and simpler density function we propose in this paper.

|  |
---|---|---
Reference Image | Ours | HF-NeuS [35]
Figure 3: Qualitative comparison between our method and mask-supervised HF-NeuS [35]. Since the mask does not show the hole of the sleeve, HF-NeuS generates a watertight sleeve for the model. Our method, which is mask-free, correctly learns the structure of the sleeve. | Reference Image | GT | Ours | VolSDF [37] | NeuS [33] | HF-NeuS [35]
---|---|---|---|---|---|---
#1 |  |  |  |  |  |
#5 |  |  |  |  |  |
#6 |  |  |  |  |  |
Figure 4: Qualitative comparisons with VolSDF [37], NeuS [33], and HF-NeuS
[34] on the DeepFashion3D [41] dataset. Row 1 is texture-rich, while rows 2
and 3 don’t contain highly contrasting colors, and thus they are textureless.
The surfaces produced by NeuS and VolSDF are closed watertight models,
resulting in large reconstruction errors near the boundaries. HF-NeuS with the
aid of masks can produce reasonable results. However, the results are still
double-covered surfaces, which are closed from a topology perspective, and
post-processing is required to extract the single-layered surface. In
contrast, our NeUDF can effectively reconstruct non-watertight models, leading
to more faithful reconstruction results without relying on masks. See the
supplementary material for additional results.

###  4.1 Comparisons on Open Models

We evaluate our method and compare it with baselines using the clothes from
DeepFashion3D, where the models have multiple open boundaries. VolSDF, NeuS,
and HF-NeuS always close the boundaries since they learn SDFs. In contrast,
our method learns UDFs, which can generate open models. Table 1 shows the
point-to-point Chamfer distances of the results. Some of the Chamfer distances
of the compared methods are large because the open holes are closed, resulting
in significant errors.

It is worth noting that HF-NeuS requires mask supervision to produce
reasonable results on the DeepFashion3D dataset. With mask supervision, it is
capable of producing double-covered surfaces that capture the geometric
features well. However, these surfaces are still closed from a topological
perspective, and their Chamfer distances are higher than those produced by our
method. Additionally, the ability of HF-NeuS to generate valid results relies
heavily on the mask; if the mask lacks boundary information, HF-NeuS may close
the boundaries. For example, as shown in Figure 3, the sleeve structure is not
labeled in the mask, resulting in a closed sleeve. In contrast, our method is
able to properly reconstruct the open sleeve structure without relying on a
mask.

As demonstrated in Figure 4, we test various types of garments, some of which
have rich textures, while others are nearly a single color. Learning UDFs for
textureless models is more challenging since various regions of a model are
ambiguous without clear color differences. However, our NeUDF generates
satisfactory results even without masks. On the other hand, NeuralUDF [22] is
unable to properly reconstruct textureless models, possibly due to their
complex density function which is difficult to converge.

Method | #1 | #2 | #3 | #4 | #5 | #6 | #7 | #8 | #9 | #10 | #11 | #12 | Mean
---|---|---|---|---|---|---|---|---|---|---|---|---|---
HF-NeuS [35] | 3.28 | 2.93 | 4.11 | 3.22 | 5.41 | 4.18 | 5.51 | 5.51 | 4.61 | 5.46 | 4.30 | 3.37 | 4.32
VolSDF [37] | 6.22 | 5.75 | 9.45 | 8.33 | 8.99 | 16.37 | 17.01 | 5.95 | 12.88 | 9.08 | 17.39 | 11.54 | 10.75
NeuS [33] | 6.75 | 4.60 | 4.35 | 7.95 | 13.52 | 10.74 | 14.54 | 6.23 | 16.69 | 17.07 | 13.21 | 5.13 | 10.07
NeUDF w/ MeshUDF [14]) | 2.68 | 4.33 | 2.90 | 3.81 | 4.34 | 3.32 | 4.53 | 3.57 | 3.65 | 3.78 | 3.26 | 4.73 | 3.74
NeUDF w/ max-weight | 2.09 | 1.77 | 1.75 | 2.00 | 2.34 | 2.16 | 3.54 | 2.60 | 2.81 | 3.63 | 2.64 | 2.61 | 2.50
Table 1: Quantitative evaluation of Chamfer distance (lower is better) on the
DeepFashion3D [41] dataset. HF-NeuS is trained with additional object mask
supervision, while the other methods are not. Our NeUDF results extracts by
maximum weights achieve the lowest Chamfer distance for all cases,
demonstrating its superior reconstruction capability for open models even
without mask supervision. For comparison, our results extracted by MeshUDF
[14] are also quantified, which achieve the second best for most cases.

###  4.2 Comparisons on Watertight Models

We compare the performance of NeUDF and HF-NeuS in reconstructing watertight
models using the DTU dataset, which is known for its rich geometric details.
The comparison focuses on the visual quality of the models’ output,
specifically surface smoothness and the presence of artifacts. The DTU dataset
poses challenges in 3D reconstruction due to the fact that many of the images
only show a part of the object-of-interest. Despite this challenge, NeUDF is
able to reconstruct watertight models with acceptable visual quality,
especially in regions that are only visible in a few images. In comparison,
HF-NeuS produces smoother surfaces than ours but introduces noticeable
artifacts in regions with limited visibility. To illustrate the comparison,
Figure 5 shows two examples of the output from both NeUDF and HF-NeuS on the
DTU dataset.

Reference Image | Ours | HF-NeuS
---|---|---
|  |
|  |
Figure 5: Qualitative comparison on the DTU dataset. While HF-NeuS is able to
produce smooth surfaces, it also introduces noticeable artifacts in regions
with limited visibility. Our method uses UDFs which are generally more
difficult to learn than SDFs, thereby our resulting surfaces are not as smooth
as theirs. However, our results contain very small artifacts.

###  4.3 Ablation Studies

_Regularization loss._ The effectiveness of the iso-surface regularizer
$\mathcal{L}_{reg}$ is verified through an ablation study. Figure 6
demonstrates that when trained with the regularization loss, the network
successfully removes small and unwanted pieces, such as the parts covering the
cuff of the sleeve in the reconstructed model. The iso-surface regularizer
encourages clean and smooth UDFs, which reduces artifacts and produces a high-
quality reconstructed model.

|
---|---
w/ reg. | w/o reg.
|
w/ AHS | w/o AHS
Figure 6: Ablation studies. Left two: The regularization loss
$\mathcal{L}_{reg}$ helps to remove unwanted surfaces on the cuff of the
sleeve. The Chamfer distance drops by 5.85%. Right two: Adaptive hierarchical
sampling leads to more accurate point sampling and helps to learn more
accurate models. The Chamber distance also drops by 25%.

_Adaptive hierarchical sampling._ We train the model with and without AHS and
compare their performance in terms of accuracy and visual quality. We found
that using adaptive hierarchical upsampling helps to sample points more
closely to the surface, especially for the parts close to boundaries,
resulting in more accurate color computation in volume rendering. As shown in
Figure 6, training with AHS produces better reconstruction on the sleeves.

Figure 7: Ablation study on the usage of ReLU (orange) [10] vs softplus (blue)
[8] in training. The former is non-differentiable at $0$ and its gradient
vanishes for negative input, whereas the latter is differentiable everywhere.
Using ReLU after the output layer of the MLP, the network makes progress at
the early stage of training, but collapses after 40K iterations, leading to a
training loss reduction through the rendering of only backgrounds. In
contrast, softplus leads to correct learning of both geometry and color, and
consistently decreases the training loss over iterations.

_Non-negativity._ Ensuring that the computed distances in the proposed method
are non-negative is important, and can be achieved by applying either ReLU or
softplus to the MLP output. However, ReLU is not differentiable at 0 and has
vanishing gradients for negative inputs, which can make the network difficult
to train. An ablation study confirms that training with ReLU only results in
early progress, but fails to learn a valid UDF later on. See Figure 7 for
details.

_Extracting surface points._ Instead of extracting the zero level-set from the
learned UDFs, we construct the object surface by finding the points with the
maximum weights. We compare our surface extraction method with MeshUDF [14],
the state-of-the-art method for extracting the zero level-set from UDFs. Since
our method is tailored to the proposed density function $\hat{\sigma}$, it
produces more accurate results than MeshUDF, which achieves the second best
for most cases. See Table 1.

_Limitations._ Our method introduces a small bias during training, we can
reduce the bias, but cannot eliminate it. Also, if the object has a very
similar color to the background, our method has difficulty differentiating
between the foreground object and the background, leading to incorrect result.

##  5 Conclusions

Overall, NeUDF offers a promising approach to the problem of reconstructing
both open and watertight models from multi-view images. Its advantages over
existing methods lie in the use of a simpler and more accurate density
function, a smooth and differentiable UDF representation, and a simple yet
effective surface reconstruction strategy tailored to the density function,
which greatly improves the learning process. Results from our experiments on
the DeepFashion3D and DTU datasets demonstrate the effectiveness of our
method, particularly in reconstructing textureless models and regions with
limited visibility. Moreover, our method does not rely on object masks, making
it more practical in real-world applications.

In the future, we plan to investigate strictly unbiased and occlusion-aware
density functions to fill the theoretical gap and further improve the accuracy
of our method. We also aim to explore the use of sparse views for UDF
learning.

## References

  * [1] Fausto Bernardini, Joshua Mittleman, Holly Rushmeier, Cláudio Silva, and Gabriel Taubin.  The ball-pivoting algorithm for surface reconstruction.  IEEE Transactions on Visualization and Computer Graphics, 5(4):349–359, oct 1999.
  * [2] Jeremy S. De Bonet.  Poxels: Probabilistic voxelized volume reconstruction.  In Proc. Int. Conf. on Computer Vision (ICCV), 1999.
  * [3] A. Broadhurst, T.W. Drummond, and R. Cipolla.  A probabilistic framework for space carving.  In Proceedings Eighth IEEE International Conference on Computer Vision. ICCV 2001, volume 1, pages 388–393 vol.1, 2001.
  * [4] Rohan Chabra, Jan E. Lenssen, Eddy Ilg, Tanner Schmidt, Julian Straub, Steven Lovegrove, and Richard Newcombe.  Deep local shapes: Learning local sdf priors for detailed 3d reconstruction.  In Computer Vision – ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XXIX, page 608–625, Berlin, Heidelberg, 2020. Springer-Verlag.
  * [5] J. Chibane, T. Alldieck, and G. Pons-Moll.  Implicit functions in feature space for 3d shape reconstruction and completion.  In 2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 6968–6979, Los Alamitos, CA, USA, jun 2020. IEEE Computer Society.
  * [6] Julian Chibane, Aymen Mir, and Gerard Pons-Moll.  Neural unsigned distance fields for implicit function learning.  In Proceedings of the 34th International Conference on Neural Information Processing Systems, NIPS’20, Red Hook, NY, USA, 2020. Curran Associates Inc.
  * [7] François Darmon, Bénédicte Bascle, Jean-Clément Devaux, Pascal Monasse, and Mathieu Aubry.  Improving neural implicit surfaces geometry with patch warping.  In 2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 6250–6259, 2022.
  * [8] Charles Dugas, Yoshua Bengio, François Bélisle, Claude Nadeau, and René Garcia.  Incorporating second-order functional knowledge for better option pricing.  In Proceedings of the 13th International Conference on Neural Information Processing Systems, NIPS’00, page 451–457, Cambridge, MA, USA, 2000\. MIT Press.
  * [9] Qiancheng Fu, Qingshan Xu, Yew-Soon Ong, and Wenbing Tao.  Geo-neus: Geometry-consistent neural implicit surfaces learning for multi-view reconstruction.  Advances in Neural Information Processing Systems (NeurIPS), 2022\.
  * [10] K Fukushima.  Cognitron: a self-organizing multilayered neural network.  Biological cybernetics, 20(3-4):121—136, November 1975.
  * [11] Yasutaka Furukawa and Carlos Hernández.  Multi-view stereo: A tutorial.  Found. Trends. Comput. Graph. Vis., 9(1–2):1–148, jun 2015.
  * [12] Silvano Galliani, Katrin Lasinger, and Konrad Schindler.  Massively parallel multiview stereopsis by surface normal diffusion.  In 2015 IEEE International Conference on Computer Vision (ICCV), pages 873–881, 2015.
  * [13] Amos Gropp, Lior Yariv, Niv Haim, Matan Atzmon, and Yaron Lipman.  Implicit geometric regularization for learning shapes.  In ICML, volume 119, pages 3789–3799, 2020.
  * [14] Benoit Guillard, Federico Stella, and Pascal Fua.  Meshudf: Fast and differentiable meshing of unsigned distance field networks.  In European Conference on Computer Vision, 2022.
  * [15] Richard Hartley and Andrew Zisserman.  Multiple View Geometry in Computer Vision.  Cambridge University Press, 2 edition, 2004.
  * [16] Fei Hou, Chiyu Wang, Wencheng Wang, Hong Qin, Chen Qian, and Ying He.  Iterative poisson surface reconstruction (ipsr) for unoriented points.  ACM Trans. Graph., 41(4), jul 2022.
  * [17] Rasmus Jensen, Anders Dahl, George Vogiatzis, Engil Tola, and Henrik Aanæs.  Large scale multi-view stereopsis evaluation.  In 2014 IEEE Conference on Computer Vision and Pattern Recognition, pages 406–413, 2014.
  * [18] Mengqi Ji, Jinzhi Zhang, Qionghai Dai, and Lu Fang.  Surfacenet+: An end-to-end 3d neural network for very sparse multi-view stereopsis.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 43(11):4078–4093, 2021.
  * [19] Abhishek Kar, Christian Häne, and Jitendra Malik.  Learning a multi-view stereo machine.  In Proceedings of the 31th International Conference on Neural Information Processing Systems, NIPS’17, page 364–375, Red Hook, NY, USA, 2017\. Curran Associates Inc.
  * [20] Michael Kazhdan and Hugues Hoppe.  Screened poisson surface reconstruction.  ACM Trans. Graph., 32(3), July 2013.
  * [21] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In 3rd International Conference on Learning Representations, ICLR, San Diego, CA, USA, May 2015.
  * [22] Xiaoxiao Long, Cheng Lin, Lingjie Liu, Yuan Liu, Peng Wang, Christian Theobalt, Taku Komura, and Wenping Wang.  Neuraludf: Learning unsigned distance fields for multi-view reconstruction of surfaces with arbitrary topologies.  ARXIV, 2022.
  * [23] Xiaoxiao Long, Cheng Lin, Peng Wang, Taku Komura, and Wenping Wang.  Sparseneus: Fast generalizable neural surface reconstruction from sparse views.  In Computer Vision – ECCV 2022: 17th European Conference, Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part XXXII, page 210–227, Berlin, Heidelberg, 2022. Springer-Verlag.
  * [24] Baorui Ma, Zhizhong Han, Yu-Shen Liu, and Matthias Zwicker.  Neural-pull: Learning signed distance function from point clouds by learning to pull space onto surface.  In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, ICML 2021, 18-24 July 2021, Virtual Event, volume 139 of Proceedings of Machine Learning Research, pages 7246–7257. PMLR, 2021.
  * [25] N. Max.  Optical models for direct volume rendering.  IEEE Transactions on Visualization and Computer Graphics, 1(2):99–108, 1995.
  * [26] Lars Mescheder, Michael Oechsle, Michael Niemeyer, Sebastian Nowozin, and Andreas Geiger.  Occupancy networks: Learning 3d reconstruction in function space.  In 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 4455–4465, 2019.
  * [27] Ben Mildenhall, Pratul P. Srinivasan, Matthew Tancik, Jonathan T. Barron, Ravi Ramamoorthi, and Ren Ng.  Nerf: Representing scenes as neural radiance fields for view synthesis.  In European Conference on Computer Vision, 2020.
  * [28] Jeong Joon Park, Peter R. Florence, Julian Straub, Richard A. Newcombe, and S. Lovegrove.  Deepsdf: Learning continuous signed distance functions for shape representation.  2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 165–174, 2019.
  * [29] Keunhong Park, Utkarsh Sinha, Jonathan T. Barron, Sofien Bouaziz, Dan B Goldman, Steven M. Seitz, and Ricardo Martin-Brualla.  Nerfies: Deformable neural radiance fields.  ICCV, 2021.
  * [30] Johannes L. Schönberger, Enliang Zheng, Jan-Michael Frahm, and Marc Pollefeys.  Pixelwise view selection for unstructured multi-view stereo.  In Bastian Leibe, Jiri Matas, Nicu Sebe, and Max Welling, editors, Computer Vision – ECCV 2016, pages 501–518, Cham, 2016. Springer International Publishing.
  * [31] Vincent Sitzmann, Julien N. P. Martel, Alexander W. Bergman, David B. Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions.  In Proceedings of the 34th International Conference on Neural Information Processing Systems, NIPS’20, Red Hook, NY, USA, 2020. Curran Associates Inc.
  * [32] J. Sun, Y. Xie, L. Chen, X. Zhou, and H. Bao.  Neuralrecon: Real-time coherent 3d reconstruction from monocular video.  In 2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 15593–15602, Los Alamitos, CA, USA, jun 2021. IEEE Computer Society.
  * [33] Peng Wang, Lingjie Liu, Yuan Liu, Christian Theobalt, Taku Komura, and Wenping Wang.  Neus: Learning neural implicit surfaces by volume rendering for multi-view reconstruction.  In Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021, December 6-14, 2021, virtual, pages 27171–27183, 2021.
  * [34] Yifan Wang, Lukas Rahmann, and Olga Sorkine-Hornung.  Geometry-consistent neural shape representation with implicit displacement fields.  In The Tenth International Conference on Learning Representations, ICLR 2022, Virtual Event, April 25-29, 2022. OpenReview.net, 2022.
  * [35] Yiqun Wang, Ivan Skorokhodov, and Peter Wonka.  Hf-neus: Improved surface reconstruction using high-frequency details.  In Proceedings of the 36th International Conference on Neural Information Processing Systems, 2022.
  * [36] Y. Yao, Z. Luo, S. Li, T. Shen, T. Fang, and L. Quan.  Recurrent mvsnet for high-resolution multi-view stereo depth inference.  In 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 5520–5529, Los Alamitos, CA, USA, jun 2019. IEEE Computer Society.
  * [37] Lior Yariv, Jiatao Gu, Yoni Kasten, and Yaron Lipman.  Volume rendering of neural implicit surfaces.  In Thirty-Fifth Conference on Neural Information Processing Systems, 2021.
  * [38] Kai Zhang, Gernot Riegler, Noah Snavely, and Vladlen Koltun.  Nerf++: Analyzing and improving neural radiance fields.  arXiv:2010.07492, 2020.
  * [39] F. Zhao, W. Wang, S. Liao, and L. Shao.  Learning anchored unsigned distance functions with gradient direction alignment for single-view garment reconstruction.  In 2021 IEEE/CVF International Conference on Computer Vision (ICCV), pages 12654–12663, Los Alamitos, CA, USA, oct 2021. IEEE Computer Society.
  * [40] Junsheng Zhou, Baorui Ma, Liu Yu-Shen, Fang Yi, and Han Zhizhong.  Learning consistency-aware unsigned distance functions progressively from raw point clouds.  In Advances in Neural Information Processing Systems (NeurIPS), 2022\.
  * [41] Heming Zhu, Yu Cao, Hang Jin, Weikai Chen, Dong Du, Zhangye Wang, Shuguang Cui, and Xiaoguang Han.  Deep fashion3d: A dataset and benchmark for 3d garment reconstruction from single images.  In Andrea Vedaldi, Horst Bischof, Thomas Brox, and Jan-Michael Frahm, editors, Computer Vision – ECCV 2020, pages 512–530, Cham, 2020. Springer International Publishing.

##  6 Supplementary Material

###  6.1 More Results

We present the remaining results on DeepFashion3D [41] dataset in Figure 8.
They have various textures. We consider models having mostly pure color, or
without highly contrasting colors to be textureless. Hence, #7, #8 and #9 are
textureless, models #2 and #11 are with grid textures, models #3 and #4 are
with stripe textures, and models #10 and #12 are with fancy textures.

We also generate three more results of challenge textureless models, as shown
in Figure 9.

| Reference Image | GT | Ours | VolSDF [37] | NeuS [33] | HF-NeuS [35]
---|---|---|---|---|---|---
#2 |  |  |  |  |  |
#3 |  |  |  |  |  |
#4 |  |  |  |  |  |
#7 |  |  |  |  |  |
#8 |  |  |  |  |  |
#9 |  |  |  |  |  |
#10 |  |  |  |  |  |
#11 |  |  |  |  |  |
#12 |  |  |  |  |  |
Figure 8: The remaining qualitative comparisons with VolSDF [37], NeuS [33], and HF-NeuS [34] on the DeepFashion3D [41] dataset. | Reference Image | GT | Ours
---|---|---|---
#13 |  |  |
|  |  | CD$=2.40\times 10^{-3}$
#14 |  |  |
|  |  | CD$=2.03\times 10^{-3}$
#15 |  |  |
|  |  | CD$=2.51\times 10^{-3}$
Figure 9: Additional results on textureless models. CD is the Chamfer
distance.


Contraction in $L^{1}$ and large time behavior for a system arising in
chemical reactions and molecular motors

Michel Chipot 111Angewandte Mathematik, Universitaet Zurich, CH-8057 Zurich
E-mail<EMAIL_ADDRESS>Danielle Hilhorst 222Laboratoire de
Mathématiques, CNRS and Université de Paris-Sud, 91405 Orsay Cédex, France;
E-mail<EMAIL_ADDRESS>
David Kinderlehrer 333Center for Nonlinear Analysis and Department of
Mathematical Sciences, Carnegie Mellon University, Pittsburgh, PA 15213;
E-mail<EMAIL_ADDRESS>Michał Olech 444Instytut Matematyczny
Uniwersytetu Wrocławskiego, pl. Grunwaldzki 2/4, 50-384 Wrocław, Polska;
Laboratoire de Mathématiques, CNRS Université de Paris-Sud, 91405 Orsay Cédex,
France; E-mail<EMAIL_ADDRESS>
The preparation of this article has been partially supported by a KBN/MNiI
grant 1 P03A 008 30 and a Marie Curie Transfer of Knowledge Fellowship of the
European Community’s Sixth Framework Programme under the contracts MTKD-
CT-2004-013389, DMS 0305794, DMS 0405343, and DMS 0806703. We would also like
to mention the support of the Swiss National Science Foundation under the
contract # 20-113287/1.

###### Abstract

We prove a contraction in $L^{1}$ property for the solutions of a nonlinear
reaction–diffusion system whose special cases include intercellular transport
as well as reversible chemical reactions. Assuming the existence of stationary
solutions we show that the solutions stabilize as $t$ tends to infinity.
Moreover, in the special case of linear reaction terms, we prove the existence
and the uniqueness (up to a multiplicative constant) of the stationary
solution.

Key words: weakly coupled system, molecular motor, transport, parabolic
systems, contraction property.

AMS subject classification: 34D23, 35K45, 35K50, 35K55, 35K57, 92C37, 92C45.

##  1 Introduction

We start with two specific reaction-diffusion systems. The first one describes
a reversible reaction and the other one a molecular motor. We first consider
the reversible chemical reaction (see also Bothe [4], Bothe and Hilhorst [5],
Desvillettes and Fellner [10] and Érdi and Tóth [11]). It involves a reaction-
diffusion system of the form

$\begin{split}u_{t}=&\,d_{1}\Delta u-\alpha
k\big{(}r_{A}(u)-r_{B}(v)\big{)}\quad\text{in}\quad\Omega\times(0,T),\ \
\Omega\subset\mathbb{R}^{d},\\\ v_{t}=&\,d_{2}\Delta v+\beta
k\big{(}r_{A}(u)-r_{B}(v)\big{)}\quad\text{in}\quad\Omega\times\ (0,T),\ \
\Omega\subset\mathbb{R}^{d},\end{split}$ (1.1)

together with the homogeneous Neumann boundary conditions, where $d_{1}$,
$d_{2}$, $\alpha$, $\beta$, $k$ and $T$ are positive constants and where
$\Omega$ is a bounded subset of $\mathbb{R}^{d}$ with smooth boundary. Such
systems describe, with a suitable choice of the functions $r_{A}$ and $r_{B}$,
chemical reactions for two mobile species. For example, functions
$r_{A}(u)=u^{k},\ r_{B}(v)=v^{m}$ correspond to a reversible reaction
$kA\rightleftharpoons mB$. Reactions of the type $q_{1}A_{1}+\ldots
q_{k}A_{k}\rightleftharpoons q_{1}B_{1}+\ldots q_{m}B_{m}$ can also be
described by similar systems with more complicated reactions terms.

Another model problem is a system in $d=1$ space dimension and $n$ unknown
variables $u_{1},\dots,u_{n}$, $n>1$, for intercellular transport, namely

$\begin{split}\frac{\partial u_{i}}{\partial t}&=\frac{\partial}{\partial
x}\left(\sigma\frac{\partial u_{i}}{\partial
x}+u_{i}\psi_{i}^{\prime}\right)\\\
&\qquad\qquad+\sum_{j=1}^{n}a_{ij}u_{j}\quad\text{in}\quad
Q_{T}=[0,1]\times(0,T)\\\ &\sigma\frac{\partial u_{i}}{\partial
x}+u_{i}\psi_{i}^{\prime}=0\quad\text{on}\quad\partial
Q_{T}=\\{0,1\\}\times(0,T),\end{split}$

where

$\begin{split}&a_{ii}\leq 0,\ a_{ij}\geq 0\mbox{~{}~{}for~{}~{}all~{}~{}}\ \
i\in\\{1,\ldots,n\\},i\neq j,\\\
&\sum_{i=1}^{n}a_{ij}=0\mbox{~{}~{}for~{}~{}all~{}~{}}i,j\in\\{1,\ldots,n\\}.\end{split}$
(1.2)

It models transport via motor proteins in the eukaryotic cell where chemical
energy is transduced into directed motion. A derivation of the system from a
mass transport viewpoint is given in [7]. For an analysis of the steady state
solutions and for further references we refer to [6], [12], [13], and [20].

In this paper we study the corresponding system in higher space dimension,
namely

$\begin{split}\frac{\partial u_{i}}{\partial
t}=\text{div}\big{(}\sigma_{i}\nabla u_{i}&+u_{i}\nabla\psi_{i}\big{)}\\\
&+\alpha_{i}\bigg{(}\sum_{j=1}^{n}\lambda_{ij}r_{j}\big{(}u_{j}(x,t),x\big{)}\bigg{)}\quad\mbox{~{}~{}in~{}~{}}Q_{T},\end{split}$
(1.3a) where $i\in\\{1,\ldots,n\\}$, and $u_{i}(x,t):Q_{T}\to\mathbb{R}^{+}$,
with $Q_{T}=\Omega\times(0,T)$, $\Omega$ an open bounded subset of
$\mathbb{R}^{d}$ with smooth boundary, and $T$ some positive constant. We
supplement this system with the Robin (no-flux) boundary conditions
$\sigma_{i}\frac{\partial
u_{i}}{\partial\nu}+u_{i}\frac{\partial\psi_{i}}{\partial\nu}=0,\quad
i\in\\{1,\ldots,n\\},\quad\text{on}\quad\partial\Omega\times(0,T),$ (1.3b)
where $\nu$ is the outward normal vector to $\partial\Omega$, and the initial
conditions $u_{1}(x,0)=u_{0,1}(x),\ \ldots\ ,u_{n}(x,0)=u_{0,n}(x),\quad
x\in\Omega.$ (1.3c)

We assume that the following hypotheses hold

  1. 1.

The constants $\sigma_{i}$ and $\alpha_{i}\in\mathbb{R}$, where
$i\in\\{1,\ldots,n\\}$, are strictly positive;

  2. 2.

For $i,j\in\\{1,\ldots,n\\},\ \lambda_{ii}\leq 0,\ \lambda_{ij}\geq 0$ if
$i\neq j$, $\sum_{k=1}^{n}\lambda_{kj}=0$;

  3. 3.

for all $i\in\\{1,\ldots,n\\}$, the smooth functions $r_{i}$ are nondecreasing
with respect to the first variable; $r_{i}(0,x)=0$ and we assume that the
functions $\psi_{i}$ are smooth as well;

  4. 4.

$u_{i}(.,0)=u_{0i}\in C(\overline{\Omega}),\ u_{0i}\geqslant 0$.

In the linear case of the molecular motors, it amounts to choosing

$r_{i}(s,x)=s,\ \lambda_{ij}=a_{ij}\ \text{and}\ \alpha_{i}=1\ \text{for all}\
i,j\in\\{1,\ldots,n\\}.$ (1.4)

We denote by Problem (P) the system (1.3a) together with the boundary and
initial conditions (1.3b), (1.3c), and admit without proof that Problem (P)
possesses a unique smooth and bounded solution on each time interval $(0,T]$.
An essential idea for proving the existence of a solution would be to apply
the Comparison principle Theorem 2.2 below to deduce that any solution of
Problem (P) has to be nonnegative and bounded from above by a stationary
solution.

Finally, we note that because of the boundary conditions (1.3b) the quantity

$\sum_{i=1}^{n}\frac{1}{\alpha_{i}}\int_{\Omega}u_{i}(x,t)\,\text{d}x$ (1.5)

is conserved in time.

The organization of this paper is as follows. In Section 2 we prove a
comparison principle for Problem (P). The main idea, which permits to show
that Problem (P) is cooperative, is a change of functions which transforms the
Robin boundary conditions into the homogeneous Neumann boundary conditions. In
Section 3 we establish a contraction in $L^{1}$ property for the corresponding
semigroup solution. Let us point out the similarity with an old result due to
Crandall and Tartar [8] where they proved in a scalar case that in the
presence of a conservation of the integral property such as (1.5), a
comparison principle such as Theorem 2.2 is equivalent to a contraction in
$L^{1}$ property such as the inequality (3.4) below. As far as we know such an
abstract result is not known in the case of systems.
Section 4 deals with the large time behavior of the solutions. Supposing the
existence of a stationary solution, we construct a continuum of stationary
solutions and prove that the solutions stabilize as $t$ tends to infinity. Let
us mention a result by Perthame [19] who proved the stabilization in the case
of the two component one-dimensional molecular motor problem. Finally in
Section 5, show the existence and uniqueness (up to a multiplicative constant)
of the stationary solution of the molecular motor problem.

Acknowledgment: The authors acknowledge the preliminary master thesis work of
Aude Brisset about the corresponding two component system. They are grateful
to the professors Piotr Biler, Stuart Hastings, Annick Lesne and Hiroshi
Matano for very fruitful discussions.

##  2 Comparison principle

First, we remark that the system of equations (1.3a) is cooperative. However,
since nothing is known about the sign of the coefficients
$\displaystyle{\frac{\partial\psi_{i}}{\partial\nu}}$ in the Robin boundary
conditions (1.3b), we cannot decide whether the Problem (P) is cooperative.
This leads us to perform a change of variables which transforms the Robin
boundary conditions into the homogeneous Neumann boundary conditions.

###  2.1 The change of unknown functions

Performing the change of variables

$w_{i}(x,t)=u_{i}(x,t){\,\text{e}}^{\nicefrac{{\psi_{i}(x)}}{{\sigma_{i}}}},\quad
i\in\\{1,\ldots,n\\},$ (2.1)

we deduce from (1.3) that ${\vec{w}}:=(w_{1},\ldots,w_{n})$ satisfies the
parabolic problem

$\begin{split}&\frac{\partial w_{i}}{\partial
t}=\sigma_{i}{\,\text{e}}^{\nicefrac{{\psi_{i}(x)}}{{\sigma_{i}}}}\text{div}\big{(}{\,\text{e}}^{\nicefrac{{-\psi_{i}(x)}}{{\sigma_{i}}}}\nabla
w_{i}\big{)}\\\
&+\alpha_{i}{\,\text{e}}^{\nicefrac{{\psi_{i}(x)}}{{\sigma_{i}}}}\bigg{(}\sum_{j=1}^{n}\lambda_{ij}r_{j}\big{(}w_{j}(x,t){\,\text{e}}^{\nicefrac{{-\psi_{j}(x)}}{{\sigma_{j}}}},x\big{)}\bigg{)}\quad\mbox{~{}~{}in~{}~{}}Q_{T},\end{split}$
(2.2)

together with the homogeneous Neumann boundary conditions

$\frac{\partial w_{i}}{\partial\nu}=0,\quad
i\in\\{1,\ldots,n\\},\quad\text{on}\quad\partial\Omega,$ (2.3)

and the initial conditions

$w_{i}(x,0)=u_{0,i}(x){\,\text{e}}^{\nicefrac{{\psi_{i}(x)}}{{\sigma_{i}}}},\quad
i\in\\{1,\ldots,n\\},\quad x\in\Omega.$ (2.4)

In the following, we denote by Problem $P_{N}$ — the problem (2.2), (2.3),
(2.4). To begin with we define the operators

$\begin{split}&{\cal L}_{i}(w_{i})=\frac{\partial w_{i}}{\partial
t}-\sigma_{i}{\,\text{e}}^{\nicefrac{{\psi_{i}(x)}}{{\sigma_{i}}}}\text{div}\big{(}{\,\text{e}}^{\nicefrac{{-\psi_{i}(x)}}{{\sigma_{i}}}}\nabla
w_{i}\big{)}\\\
&-\alpha_{i}{\,\text{e}}^{\nicefrac{{\psi_{i}(x)}}{{\sigma_{i}}}}\bigg{(}\sum_{j=1}^{n}\lambda_{ij}r_{j}\big{(}w_{j}(x,t){\,\text{e}}^{\nicefrac{{-\psi_{j}(x)}}{{\sigma_{j}}}},x\big{)}\bigg{)}\quad\mbox{~{}~{}in~{}~{}}Q_{T}.\end{split}$
(2.5)

We say that $\left(\underline{w}_{1},\ldots,\underline{w}_{n}\right)$ is a
subsolution of Problem $P_{N}$ if

$\begin{split}{\cal L}_{i}(\underline{w}_{i})\leqslant&\;0\quad\text{in}\quad
Q_{T},\\\
\frac{\partial\underline{w}_{i}}{\partial\nu}\leqslant&\;0\quad\text{on}\quad\partial\Omega\times(0,T),\\\
\underline{w}_{i}(x,0)\leqslant&\;w_{i}(x,0),\quad x\in\Omega\end{split}$
(2.6)

for all $i\in\\{1,\ldots,n\\}$. We define similarly a supersolution
$\left(\overline{u}_{1},\ldots,\overline{u}_{n}\right)$ of Problem $P_{N}$ by
the inequalities

$\begin{split}{\cal L}_{i}(\overline{w}_{i})\geqslant&\;0\quad\text{in}\quad
Q_{T},\\\
\frac{\partial\overline{w}_{i}}{\partial\nu}\geqslant&\;0\quad\text{on}\quad\partial\Omega\times(0,T),\\\
\overline{w}_{i}(x,0)\geqslant&\;w_{i}(x,0),\quad x\in\Omega.\end{split}$
(2.7)

The following comparison theorem holds ([2], [21]).

######  Theorem 2.1.

Let $(\underline{w}_{1},\ldots,\underline{w}_{n})$ and
$(\overline{w}_{1},\ldots,\overline{w}_{n})$, be a sub - and a super -
solution, respectively, for the operators ${\cal L}_{j}$ defined by (2.5) with
$j\in\\{1,\ldots,n\\}$, which means that (2.6) and (2.7) hold for
$i\in\\{1,\ldots,n\\}$. Then $\underline{w}_{i}\leqslant\overline{w}_{i}$ in
$Q_{T}$. Moreover, for all $i\in\\{1,\ldots,n\\}$ such that
$\underline{w}_{i}\leqslant\overline{w}_{i}$ and
$\underline{w}_{i}\not\equiv\overline{w}_{i}$ on $\\{t=0\\}\times\Omega$ then
$\underline{w}_{i}<\overline{w}_{i}$ in $Q_{T}$. $\blacksquare$

This comparison theorem immediately translates into a comparison theorem for
solutions of the original Problem (P). For all $i\in\\{1,\ldots,n\\}$, we
define the operators

$\begin{split}L_{i}(u_{i})=(u_{i})_{t}-\text{div}\big{(}&\sigma_{i}\nabla
u_{i}+u_{i}\nabla\psi_{i}\big{)}\\\
&-\alpha_{i}\bigg{(}\sum_{j=1}^{n}\lambda_{ij}\,r_{j}\,(u_{j},x)\bigg{)}\quad\text{in}\quad
Q_{T}.\end{split}$ (2.8)

The following result holds.

######  Theorem 2.2.

Let $(\underline{u}_{1},\ldots,\underline{u}_{n})$ and
$(\overline{u}_{1},\ldots,\overline{u}_{n})$, be a sub - and a super -
solution, respectively, for the operators $L_{j}$, defined by (2.8) with
$j\in\\{1,\ldots,n\\}$. Then $\underline{u}_{i}\leqslant\overline{u}_{i}$ in
$Q_{T}$. Moreover, for all $i\in\\{1,\ldots,n\\}$ such that
$\underline{u}_{i}\leqslant\overline{u}_{i}$ and
$\underline{u}_{i}\not\equiv\overline{u}_{i}$ on $\\{t=0\\}\times\Omega$ then
$\underline{u}_{i}<\overline{u}_{i}$ in $Q_{T}$. $\blacksquare$

Next we state two immediate corollaries of Theorem 2.2.

######  Corollary 2.3.

(uniqueness) If $(u_{1}^{1},\ldots,u_{n}^{1})$ and
$(u_{1}^{2},\ldots,u_{n}^{2})$ are solutions of Problem (P) with the same
initial condition
$(u_{0,1},\ldots,u_{0,n})\in\big{(}C(\overline{\Omega}))^{n}$, then for all
$i\in\\{1,\ldots,n\\},\ u_{i}^{1}=u_{i}^{2}$ . $\blacksquare$

######  Corollary 2.4.

(positivity) If $(u_{1},\ldots,u_{n})$ is the solution of Problem (P) with the
nonnegative initial condition
$(u_{0,1},\ldots,u_{0,n})\in\big{(}C(\overline{\Omega})\big{)}^{n}$, then for
all $i\in\\{1,\ldots,n\\},$ $u_{i}\geqslant 0$. Moreover, for all
$i\in\\{1,\ldots,n\\}$, such that $u_{0,i}\geqslant 0$ and $u_{0,i}\not\equiv
0$, $u_{i}>0$ in $\Omega$. $\blacksquare$

##  3 Contraction property

The purpose of this section is to show a contraction in
$\big{(}L^{1}(\Omega)\big{)}^{n}$ property for the solutions of Problem (P)
with the initial conditions belonging to
$\big{(}L^{\infty}(\Omega)\big{)}^{n}$. The main steps of the proof rely upon
arguments due to [3] and [18].

We first introduce some notation. We suppose that the functions
$(u_{1}^{1},\ldots,u_{n}^{1})$ and $(u_{1}^{2},\ldots,u_{n}^{2})$ are the
solutions of Problem (P) with the initial conditions
$(u_{0,1}^{1},\ldots,u_{0,n}^{1})$ and $(u_{0,1}^{2},\ldots,u_{0,n}^{2})$,
respectively. Define

$(U_{1},\ldots,U_{n}):=(u_{1}^{1}-u_{1}^{2},\ldots,u_{n}^{1}-u_{n}^{2}).$
(3.1)

Then

$\begin{split}&(U_{i})_{t}=\text{div}\big{(}\sigma_{i}\nabla
U_{i}+U_{i}\nabla\psi_{i}\big{)}\\\ &+\
\alpha_{i}\sum_{j=1}^{n}\lambda_{ij}\big{(}r_{j}(u_{j}^{1}(x,t),x)-r_{j}(u_{j}^{2}(x,t),x)\big{)}\quad\text{in}\quad
Q_{T},\\\ &\sigma_{i}\frac{\partial
U_{i}}{\partial\nu}+U_{i}\frac{\partial\psi_{i}}{\partial\nu}=0\quad\text{on}\quad\partial\Omega\times(0,T),\\\
&U_{i}(x,0)=U_{0,i}(x)\quad\text{for}\quad x\in\Omega,\end{split}$ (3.2)

together with

$U_{0,i}=u_{0,i}^{1}-u_{0,i}^{2},$ (3.3)

for each $i\in\\{1,\ldots,n\\}$.

Next we prove the following contraction in $L^{1}$ property.

######  Theorem 3.1.

For all $t>0$,

$\frac{1}{\alpha_{1}}\|U_{1}(\cdot,t)\|_{L^{1}(\Omega)}+\ldots+\frac{1}{\alpha_{n}}\|U_{n}(\cdot,t)\|_{L^{1}(\Omega)}\\\
\leqslant\frac{1}{\alpha_{1}}\|U_{0,1}(\cdot)\|_{L^{1}(\Omega)}+\ldots+\frac{1}{\alpha_{n}}\|U_{0,n}(\cdot)\|_{L^{1}(\Omega)},$
(3.4)

where $U_{i}$ and $U_{0,i},\ i\in\\{1,\ldots,n\\}$, are defined by (3.1) and
(3.3), respectively.

Proof  Dividing each partial differential equation of (3.2) by $\alpha_{i}$
and summing them up, we obtain

$\begin{split}\frac{\text{d}}{\text{d}t}\bigg{(}\sum_{i=1}^{n}\frac{1}{\alpha_{i}}&U_{i}\bigg{)}=\sum_{i=1}^{n}\frac{1}{\alpha_{i}}\,\text{div}\left(\sigma_{i}\nabla
U_{i}+U_{i}\nabla\psi_{i}\right)\\\
&+\sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_{ij}\Big{(}r_{j}(u_{j}^{1}(x,t),x)-r_{j}(u_{j}^{2}(x,t),x)\Big{)}\\\
=&\sum_{i=1}^{n}\frac{1}{\alpha_{i}}\,\text{div}\left(\sigma_{i}\nabla
U_{i}+U_{i}\nabla\psi_{i}\right)\\\
&+\sum_{j=1}^{n}\bigg{\\{}\Big{(}r_{j}(u_{j}^{1}(x,t),x)-r_{j}(u_{j}^{2}(x,t),x)\Big{)}\sum_{i=1}^{n}\lambda_{ij}\bigg{\\}}\\\
=&\sum_{i=1}^{n}\frac{1}{\alpha_{i}}\,\text{div}\left(\sigma_{i}\nabla
U_{i}+U_{i}\nabla\psi_{i}\right),\end{split}$

where we have used Hypothesis 2.

This, together with the boundary conditions (1.3b), implies the conservation
in time of the quantity

$\frac{\text{d}}{\text{d}t}\sum_{i=1}^{n}\frac{1}{\alpha_{i}}\int_{\Omega}U_{i}(x,t)\,\text{d}x=0.$
(3.5)

Let us look closer at the nonlinear term in (3.2). We can write, for fixed
index $i$

$\sum_{j=1}^{n}\lambda_{ij}\big{(}r_{j}(u_{j}^{1}(x,t),x)-r_{j}(u_{j}^{2}(x,t),x)\big{)}\\\
=\sum_{j=1}^{n}\lambda_{ij}U_{j}\int_{0}^{1}\frac{\partial}{\partial
u}r_{j}(\theta
u_{j}^{1}+(1-\theta)u_{j}^{2},x)\text{d}\theta=\sum_{j=1}^{n}A_{ij}U_{j}.$

Freezing the functions $u_{i}^{k}$ for $i\in\\{1,\ldots,n\\},\ k\in\\{1,2\\}$,
we deduce that the functions $U_{1},\ldots,U_{n}$ satisfy a system of the form

$(U_{i})_{t}=\text{div}\Big{(}\sigma_{i}\nabla
U_{i}+U_{i}\nabla\psi_{i}\Big{)}+\sum_{j=1}^{n}A_{ij}U_{j}\quad\text{in}\quad
Q_{T},$ (3.6)

with the boundary and initial conditions

$\begin{split}&\sigma_{i}\frac{\partial
U_{i}}{\partial\nu}+U_{i}\frac{\partial\psi_{i}}{\partial\nu}=0\quad\text{on}\quad\partial\Omega\times(0,T),\\\
&U_{i}(x,0)=U_{0,i}(x),\quad x\in\Omega.\end{split}$ (3.7)

for $i\in\\{1,\ldots,n\\}$, where $A_{ij}$ are functions of space and time.

In order to make the notation more concise, we write

$\begin{split}\vec{U}_{0}=&\,\big{(}U_{0,1},\ldots,U_{0,n}\big{)},\\\
\vec{U}=&\,\big{(}U_{1},\ldots,U_{n}\big{)},\\\
\vec{U}_{0}^{\pm}=&\,\big{(}U_{0,1}^{\pm},\ldots,U_{0,n}^{\pm}\big{)},\\\
\vec{U}^{\pm}=&\,\big{(}U_{1}^{\pm},\ldots,U_{n}^{\pm}\big{)},\\\ \end{split}$

where $s^{+}=\max\\{s,0\\},\ s^{-}=\max\\{-s,0\\}$. By (3.6), (3.7) and
Corollary 2.3 we can write $\vec{U}$ in the form

$(\vec{U})(x,t)=\mathcal{S}(t)\,\vec{U}_{0}(x)=\big{(}{\mathcal{S}}_{1}(t)\vec{U}_{0},\ldots,{\mathcal{S}}_{n}(t)\vec{U}_{0}\big{)}(x)$

with some operator ${\mathcal{S}}(t)$. We set

$\big{(}W_{1},\ldots,W_{n}\big{)}=-\big{(}U_{1}{\,\text{e}}^{\nicefrac{{\psi_{1}(x)}}{{\sigma_{1}}}},\ldots,U_{n}{\,\text{e}}^{\nicefrac{{\psi_{n}(x)}}{{\sigma_{n}}}}\big{)},$

and
$\widetilde{A}_{ij}=A_{ij}\,{\,\text{e}}^{\nicefrac{{\psi_{i}(x)}}{{\sigma_{i}}}}\,{\,\text{e}}^{-\nicefrac{{\psi_{j}(x)}}{{\sigma_{j}}}}$.
Then, the system of equations (3.6) can be expressed in the form

$\big{(}W_{i}\big{)}_{t}={\sigma_{i}}{\,\text{e}}^{\nicefrac{{\psi_{i}(x)}}{{\sigma_{i}}}}\text{div}\Big{(}{\,\text{e}}^{-\nicefrac{{\psi_{i}(x)}}{{\sigma_{i}}}}\nabla
W_{i}\Big{)}+\sum_{j=1}^{n}\widetilde{A}_{ij}W_{j}\quad\text{in}\quad Q_{T},$
(3.8)

with the boundary and initial conditions

$\displaystyle\frac{\partial W_{i}}{\partial\nu}$
$\displaystyle=0\quad\text{on}\quad\partial\Omega\times(0,T),$ (3.9)
$\displaystyle W_{i}(x,0)$
$\displaystyle=-U_{0,i}{\,\text{e}}^{\nicefrac{{\psi_{i}(x)}}{{\sigma_{i}}}},\quad
x\in\Omega,$ (3.10)

for $i\in\\{1,\ldots,n\\}$.

Next we show that the solutions $W_{i}$ of the problem (3.8) – (3.10) with
nonpositive initial conditions are nonpositive in $\overline{\Omega}$ for all
$t\in(0,T)$. To that purpose we consider the auxiliary problem

$\displaystyle\big{(}W_{i}\big{)}_{t}-\vartheta_{i}(x)\text{div}\Big{(}\zeta_{i}(x)\nabla
W_{i}\Big{)}-\sum_{j=1}^{n}\gamma_{ij}W_{j}\leqslant 0\quad\text{in}\quad
Q_{T},$ (3.11) $\displaystyle\frac{\partial W_{i}}{\partial\nu}\leqslant
0\quad\text{on}\quad\partial\Omega\times(0,T),$ (3.12) $\displaystyle
W_{i}(x,0)=W_{0,i}(x)\leqslant 0\quad x\in\Omega,$ (3.13)

for $i\in\\{1,\ldots,n\\}$. We assume that $\vartheta_{i}(x)$ and
$\zeta_{i}(x)$ are nonnegative in $\overline{\Omega}$ and that the
coefficients $\gamma_{ij}$ satisfy the same assumptions as the coefficients
$\lambda_{ij}$ in Problem (P). The following result holds.

######  Lemma 3.2.

Let $(W_{1},\ldots,W_{n})$ be a smooth and bounded solution of the problem
(3.11) – (3.13) with nonpositive initial conditions $W_{0,i}$ on a time
interval $[0,T]$. Then $W_{i}(x,t)\leqslant 0$ in
$\overline{\Omega}\times(0,T]$. Moreover, for each $i\in\\{1,\ldots,n\\}$ such
that $W_{0,i}\leqslant 0$ and $W_{0,i}\not\equiv 0,\ W_{i}<0$ in
$\overline{\Omega}\times(0,T]$.

Proof  The result of Lemma 3.2 follows from the fact that the system (3.11),
(3.12), (3.13), with the inequalities $\\{\leqslant\\}$ replaced by the
equalities $\\{=\\}$, is a cooperative system. However, for the sake of
completeness, we present a proof below. We first remark that, in view of [21,
Remark (i), p. 191], one can always satisfy the condition

$\sum_{j=1}^{n}\gamma_{ij}\leqslant 0\text{\ for all\ }i\in\\{1,\ldots,n\\},$
(3.14)

for the matrix of coefficients $\big{(}\gamma_{ij}\big{)}_{i,j=1}^{n}$ by
performing the change of variables $\overline{W}_{i}=W_{i}{\,\text{e}}^{-ct}$
for all $i\in\\{1,\ldots,n\\}$ and $c>0$ large enough.
Thanks to the regularity of each $W_{i}$, we can apply Theorem 15, p. 191 from
[21] to conclude that $W_{i}-M\leqslant 0$ in $\overline{\Omega}\times[0,T]$
for some $M\>0$ and all $i\in\\{1,\ldots,n\\}$. In fact, we can deduce that
$W_{i}-M<0$ in $\overline{\Omega}\times(0,T)$.
Indeed, if for some $k\in\\{1,\ldots,n\\}$, $W_{k}=M$ in an interior point
$(\tilde{x},\tilde{t})\in\Omega\times(0,T)$, then Theorem 15, p. 191 in [21]
implies that $W_{k}\equiv M$ for all $0\leqslant t<\tilde{t}$, which is
impossible since $W_{k}(x,0)\leqslant 0$. If the maximum $M$ of $W_{k}$ is
attained at a boundary point $P\in\partial\Omega\times(0,T)$ then either there
exists an open ball $K\subset\Omega\times(0,T)$ such that $P\in\partial K$ and
$W_{k}-M<0$ in $K$, and the last part of Theorem 15, p. 191 in [21]
contradicts the boundary inequality (3.12), or for all open balls
$K\subset\Omega\times(0,T)$ such that $P\in\partial K$ there exists a point
$(\tilde{x},\tilde{t})\in K$ such that $W_{i}(\tilde{x},\tilde{t})=M$, and we
proceed as in the case before.
Hence, there exists $\widetilde{M}>0$, such that
$W_{i}\leqslant\widetilde{M}<M$ in $\overline{\Omega}\times[0,T]$ for all
$i\in\\{1,\ldots,n\\}$. Then we can repeat the reasoning for all $M>0$ until
$M=0$. Indeed, if this would not be the case, we find the least real number
$\overline{M}>0$, with $W_{i}\leqslant\overline{M}\leqslant\widetilde{M}$ in
$\overline{\Omega}\times[0,T]$, which leads again to the existence of a real
number $0\leqslant\widehat{M}<\overline{M}$ with the same property. This
contradicts the fact that $\overline{M}$ was defined as the least such real
number. $\blacksquare$

Since the functions $u_{i}^{1},\ u_{i}^{2}$ are bounded on
$\overline{\Omega}\times[0,T]$, it follows that the functions $W_{i}$ are
bounded on $\overline{\Omega}\times[0,T]$ for all $i\in\\{1,\ldots,n\\}$.
Then we are in a position to apply Lemma 3.2 with
$\vartheta_{i}(x)={\,\text{e}}^{\nicefrac{{\psi_{i}}}{{\sigma_{i}}}}$,
$\zeta_{i}(x)=\sigma_{i}{\,\text{e}}^{-\nicefrac{{\psi_{i}}}{{\sigma_{i}}}}$
and $\gamma_{ij}=\widetilde{A}_{ij}$ for $i,j\in\\{1,\ldots,n\\}$. We deduce
that the solutions $W_{i}$ of the problem (3.8) – (3.10) with nonpositive
initial conditions are nonpositive in $\overline{\Omega}$ for all $t\in(0,T)$.

Next we remark that the above reasoning can be applied either with
$\vec{U}_{0}$ replaced by $U^{+}_{0}$ or with $\vec{U}_{0}$ replaced by
$U^{-}_{0}$. This permits to show that
${\mathcal{S}}_{i}(t)\vec{U}_{0}^{+},{\mathcal{S}}_{i}(t)\vec{U}_{0}^{-}\geqslant
0$ and that

${\mathcal{S}}_{i}(t)\vec{U}_{0}^{\pm}>0\quad\text{if}\quad\vec{U}_{0}^{\pm}\not\equiv
0.$ (3.15)

We easily compute

$\begin{split}\sum_{i=1}^{n}\frac{1}{\alpha_{i}}&\big{\|}U_{i}(\cdot,t)\big{\|}_{L^{1}(\Omega)}-\sum_{i=1}^{n}\frac{1}{\alpha_{i}}\big{\|}U_{0,i}(\cdot)\big{\|}_{L^{1}(\Omega)}\\\
=&\
\sum_{i=1}^{n}\frac{1}{\alpha_{i}}\big{\|}{\mathcal{S}}_{i}(t)\vec{U}_{0}^{+}-{\mathcal{S}}_{i}(t)\vec{U}_{0}^{-}\big{\|}_{L^{1}(\Omega)}-\sum_{i=1}^{n}\frac{1}{\alpha_{i}}\big{\|}U_{0,i}(\cdot)\big{\|}_{L^{1}(\Omega)}\\\
=&\
\sum_{i=1}^{n}\int_{\Omega}\frac{1}{\alpha_{i}}\Big{\\{}\max\big{\\{}{\mathcal{S}}_{i}(t)\vec{U}_{0}^{+},{\mathcal{S}}_{i}(t)\vec{U}_{0}^{-}\big{\\}}\end{split}$
(3.16)
$\begin{split}&-\frac{1}{\alpha_{i}}\min\big{\\{}{\mathcal{S}}_{i}(t)\vec{U}_{0}^{+},{\mathcal{S}}_{i}(t)\vec{U}_{0}^{-}\big{\\}}\Big{\\}}\,\text{d}x-\sum_{i=1}^{n}\frac{1}{\alpha_{i}}\int_{\Omega}\big{\\{}U_{i,0}^{+}+U_{i,0}^{-}\big{\\}}\,\text{d}x\\\
=&\
\sum_{i=1}^{n}\int_{\Omega}\frac{1}{\alpha_{i}}\big{(}{\mathcal{S}}_{i}(t)\vec{U}_{0}^{+}+{\mathcal{S}}_{i}(t)\vec{U}_{0}^{-}\big{)}\,\text{d}x-\sum_{i=1}^{n}\frac{1}{\alpha_{i}}\int_{\Omega}\big{\\{}U_{i,0}^{+}+U_{i,0}^{-}\big{\\}}\,\text{d}x\\\
&-2\sum_{i=1}^{n}\int_{\Omega}\frac{1}{\alpha_{i}}\min\big{\\{}{\mathcal{S}}_{i}(t)\vec{U}_{0}^{+},{\mathcal{S}}_{i}(t)\vec{U}_{0}^{-}\big{\\}}\,\text{d}x\\\
=&\
-2\sum_{i=1}^{n}\int_{\Omega}\frac{1}{\alpha_{i}}\min\big{\\{}{\mathcal{S}}_{i}(t)\vec{U}_{0}^{+},{\mathcal{S}}_{i}(t)\vec{U}_{0}^{-}\big{\\}}\,\text{d}x\leqslant
0,\end{split}$ (3.17)

which completes the proof of (3.4). $\blacksquare$

######  Corollary 3.3.

Let $(u_{0,1}^{1},\ldots,u_{0,n}^{1}),\
(u_{0,1}^{2},\ldots,u_{0,n}^{2})\in\big{(}C(\overline{\Omega})\big{)}^{n}$ be
as in Theorem 3.1. Moreover, let us assume that for at least one index
$k\in\\{1,\ldots,n\\}$ the difference $u^{1}_{0,k}-u^{2}_{0,k}$ changes the
sign. Then, the inequality (3.4) is strict for all $t>0$, so that solution
satisfies a strict contraction property.

##  4 Large time behavior of solutions

In this section we assume the existence and uniqueness of a positive solution
$\vec{v}=(v_{1},\ldots,v_{n})\in\big{(}C(\overline{\Omega})\cap
C^{2}(\Omega)\big{)}^{n}$ of the elliptic problem

$\displaystyle\text{div}\big{(}\sigma_{i}\nabla
v_{i}+v_{i}\nabla\psi_{i}\big{)}+\alpha_{i}\bigg{(}\sum_{j=1}^{n}\lambda_{ij}r_{j}\big{(}v_{j}(x),x\big{)}\bigg{)}$
$\displaystyle=0\quad\text{in}\quad\Omega,$ (4.1)
$\displaystyle\sigma_{i}\frac{\partial
v_{i}}{\partial\nu}+v_{i}\frac{\partial\psi_{i}}{\partial\nu}$
$\displaystyle=0\quad\text{on}\quad\partial\Omega,$ (4.2)
$\displaystyle\sum_{i=1}^{n}\frac{1}{\alpha_{i}}\int_{\Omega}v_{i}(x)\,\text{d}x$
$\displaystyle=1,$ (4.3)

for $i\in\\{1,\ldots,n\\}$.

######  Definition 4.1.

We say that a vector function
$\vec{v}=(v_{1},\ldots,v_{n})\in\big{(}C(\overline{\Omega})\big{)}^{n}$ is
nonnegative (resp. positive) if $v_{i}(x)\geqslant 0$ (resp. $v_{i}(x)>0$) for
all $x\in\overline{\Omega}$ and all $i\in\\{1,\ldots,n\\}$.

Next we introduce the semigroup notation for the unique solution of Problem
(P), namely

$\vec{u}(t)=\mathcal{T}(t)\,\vec{u}_{0}=\Big{(}\mathcal{T}_{1}(t)\,\vec{u}_{0},\ldots,\mathcal{T}_{n}(t)\,\vec{u}_{0}\Big{)},$

with the initial data $\vec{u}_{0}\in\big{(}C(\overline{\Omega})\big{)}^{n}$.
The method of the proof is based upon an idea of Osher and Ralston [18]. It
mainly exploits the contraction properties for the nonlinear semigroup
$\mathcal{T}(t)\,$ given by Theorem 3.1 and Corollary 3.3. A similar reasoning
was developed in other contexts by Bertsch and Hilhorst [3], Hilhorst and
Hulshof [14] and Hilhorst and Peletier [15].

We suppose there exists a set
$\mathscr{H}\subset\big{(}C(\overline{\Omega})\cap C^{2}(\Omega)\big{)}^{n}$
of positive stationary solutions with the following property which we denote
by $\mathscr{S}$:
For each $\vec{f}=(f_{1},\ldots,f_{n})\in\big{(}C(\overline{\Omega})\cap
C^{2}(\Omega)\big{)}^{n}$ either $\vec{f}\in\mathscr{H}$ or there exists
$(\xi_{1},\ldots,\xi_{n})\in\mathscr{H}$, such that $f_{i}-\xi_{i}$ changes
the sign for at least one index $i\in\\{1,\ldots,n\\}$.

One can prove that a set $\mathscr{H}$ satisfying Property $\mathscr{S}$
exists in at least two cases:

  * i)

In the case of the system (1.1) where the Robin boundary conditions reduce to
the homogeneous Neumann boundary conditions, the set $\mathscr{H}$ is given by

$\rule{20.0pt}{0.0pt}\mathscr{H}=\Big{\\{}(a,b):\ a>0,\
b=r_{B}^{-1}(r_{A}(a))\\\ \text{and}\
\frac{a}{\alpha}+\frac{b}{\beta}=\int_{\Omega}\Big{(}\frac{u}{\alpha}+\frac{v}{\beta}\Big{)}\,\text{d}x\Big{\\}}.\rule{20.0pt}{0.0pt}$

For more details we refer to [5].

  * ii)

In the case of the molecular motor with a linear $n$-component system the set
$\mathscr{H}$ is given by

$\mathscr{H}=\big{\\{}c\vec{v}:\ c\in\mathbb{R}^{+}\big{\\}},$

where $\vec{v}$ is a unique solution of the elliptic problem (4.1) – (4.3).

######  Proposition 4.2.

The continuum $\mathscr{H}$ is such that for each

$\vec{f}=(f_{1},\ldots,f_{n})\in\big{(}C(\overline{\Omega})\cap
C^{2}(\Omega)\big{)}^{n}$

either $\vec{f}\in\mathscr{H}$, or there exists
$(\xi_{1},\ldots,\xi_{n})\in\mathscr{H}$ such that $f_{i}-\xi_{i}$ changes the
sign for at least one index $i\in\\{1,\ldots,n\\}$.

Proof

  * i)

In the case of system (1.1) the proof is rather obvious since the continuum
$\mathscr{H}$ is composed of constant pairs.

  * ii)

In the case of the molecular motor, let us assume that
$\vec{f}\not\in\mathscr{H}$. Then there does not exist any positive constant
$c$ such that $c\,\vec{v}=\vec{f}$. In particular, there exists an index
$i\in\\{1,\ldots,n\\}$ such that $v_{i}$ is not proportional to $f_{i}$, or in
other words $cv_{i}\neq f_{i}$ for all $c>0$. Without loss of generality we
can assume that the first coordinate has this property. Let $x_{0}\in\Omega$
be arbitrary. Since $v_{1}$ is strictly positive in $\overline{\Omega}$, we
can define

$c_{0}=\frac{f_{1}(x_{0})}{v_{1}(x_{0})},$

so that

$\big{(}f_{1}-c_{0}\,v_{1}\big{)}(x_{0})=0.$

Let $\mathcal{Z}=\big{\\{}x\in\overline{\Omega}:\
\big{(}f_{1}-c_{0}\,v_{1}\big{)}(x)=0\big{\\}}$. From the continuity of
$f_{1}$ and $v_{1}$, $\mathcal{Z}$ is closed as a subset of $\Omega$. If there
exist $x_{1},\,x_{2}\in\mathcal{Z}^{c}$, such that
$\big{(}f_{1}-c_{0}\,v_{1}\big{)}(x_{1})$ and
$\big{(}f_{1}-c_{0}\,v_{1}\big{)}(x_{2})$ are of different signs, then the
proof is complete. Now suppose that $\big{(}f_{1}-c_{0}\,v_{1}\big{)}(x)$ is
positive for all $x\in\mathcal{Z}^{c}$. In particular

$\big{(}f_{1}-c_{0}\,v_{1}\big{)}(\tilde{x})=d>0$

for some fixed $\tilde{x}\in\mathcal{Z}^{c}$. Then choosing
$\displaystyle\varepsilon=\frac{d}{2v_{1}(\tilde{x})}$ we see that

$\big{(}f_{1}-(c_{0}+\varepsilon)v_{1}\big{)}(\tilde{x})=\frac{d}{2}>0.$

However

$\big{(}f_{1}-(c_{0}+\varepsilon)v_{1}\big{)}(x_{0})<0.$

We proceed similarly when $\big{(}f_{1}-c_{0}\,v_{1}\big{)}(x)$ is negative
for all $x\in\mathcal{Z}^{c}$. $\blacksquare$

In the sequel we suppose that the initial data
$\vec{u}_{0}=(u_{0,1},\ldots,u_{0,n})$ from
$\big{(}C(\overline{\Omega})\big{)}^{n}$ also satisfy the following property:

$\text{{There exists\ }}\vec{h}\in\mathscr{H}\text{{\ such that\
}}0\leqslant\vec{u}_{0}\leqslant\vec{h}\quad\text{in}\quad\overline{\Omega},$
(4.4)

and remark that this property is satisfied in both the cases (i) and (ii).

######  Proposition 4.3.

Let
$\vec{u}_{0}=(u_{0,1},\ldots,u_{0,n})\in\big{(}C(\overline{\Omega})\big{)}^{n}$
satisfy the property (4.4). Then the solution $(u_{1},\ldots,u_{n})$ of
Problem (P) is such that $0\leqslant\vec{u}(t)\leqslant\vec{h}$ for all $t>0$.

Proof  We remark that $\vec{0}$ is a subsolution of Problem (P) and that
$\vec{h}$ is a supersolution, and apply Theorem 2.2. $\blacksquare$

Next we prove the main result of this section. To that purpose we first define
the norm $\boldsymbol{\big{\|}}\cdot\boldsymbol{\big{\|}_{1}}$ by

$\boldsymbol{\big{\|}}\vec{f}\
\boldsymbol{\big{\|}_{1}}=\sum_{i=1}^{n}\frac{1}{\alpha_{i}}\big{\|}f_{i}\big{\|}_{L^{1}(\Omega)}\hskip
5.0pt.$

Note that this norm is equivalent to the usual product norm in the space
$\big{(}L^{1}(\Omega)\big{)}^{n}$.

######  Theorem 4.4.

For all nonnegative
$\vec{u}_{0}=(u_{0,1},\ldots,u_{0,n})\in\big{(}C(\overline{\Omega})\big{)}^{n}$
there exists $\vec{f}=(f_{1},\ldots,f_{n})\in\mathscr{H}$, such that

$\lim_{t\rightarrow\infty}\boldsymbol{\big{\|}}\
\mathcal{T}(t)\,\vec{u}-\vec{f}\ \boldsymbol{\big{\|}_{1}}=0.$

Proof
The proof consists of several steps. To begin with we define the
$\omega$-limit set

$\omega(\vec{u}_{0})=\Big{\\{}\vec{g}\in\big{(}L^{1}(\Omega)\big{)}^{n}:\
\text{there exists a sequence \ }t_{k}\rightarrow\infty\\\ \text{\ as\
}k\rightarrow\infty,\text{\ such that\
}\lim_{k\rightarrow\infty}\boldsymbol{\big{\|}}\mathcal{T}(t_{k})\,\vec{u}_{0}-\vec{g}\
\boldsymbol{\big{\|}_{1}}=0\Big{\\}},$ (4.5)

The organization of the proof is as follows. First we show that
$\omega(\vec{u}_{0})$ is not empty. In the second step we define the Lyapunov
functional

$\mathcal{V}(\vec{\xi})=\boldsymbol{\big{\|}}\vec{\xi}-\vec{w}\boldsymbol{\big{\|}_{1}},$

where $\vec{w}$ is a stationary solution and check that it is constant on
$\omega(\vec{u}_{0})$. We then deduce that
$\omega(\vec{u}_{0})\subset\mathscr{H}$, and finally prove that
$\omega(\vec{u}_{0})$ consists of exactly one function.

Step 1. $\omega(\vec{u}_{0})$ is not empty.
Let $\varepsilon>0$ be arbitrary. Suppose that
$\Omega^{\prime}\subset\subset\Omega$ satisfy

$\big{|}\Omega\setminus\Omega^{\prime}\big{|}\leqslant\frac{\varepsilon}{2K}.$

and set

$K=\sum_{i=1}^{n}\frac{2}{\alpha_{i}}\|h_{i}\|_{C(\overline{\Omega})},$ (4.6)

where $\vec{h}$ has been introduced in (4.4). We have already proved in
Proposition 4.3 that $\mathcal{T}(t)\,\vec{u}_{0}$ is bounded in
$\big{(}L^{\infty}(\Omega)\big{)}^{n}$. Therefore there exist a vector
function $\vec{g}\in\big{(L^{\infty}}(\Omega))^{n}$ and a sequence
$\\{\vec{u}(t_{k})\\}$ such that

$\vec{u}(t_{k})\rightharpoonup\vec{g}\mbox{~{}~{}weakly in
~{}~{}}(L^{2}(\Omega))^{n},$ (4.7)

as $t_{k}\rightarrow\infty$. Next we deduce from [16, Chap. III, Theorem 10.1]
that there exists a positive constant $C$ such that

$\big{|}u_{i}(x_{1},t)-u_{i}(x_{2},t)\big{|}\leqslant C|x_{1}-x_{2}|^{\alpha}$

for all $x_{1},\ x_{2}\in\Omega^{\prime}$ and all $t>0$. Therefore, it follows
from the Ascoli-Arzelà Theorem (see, e.g., [1, Theorem 1.33]) that
$\vec{u}(t_{k})\rightarrow\vec{g}$ as $t_{k}\rightarrow\infty$, uniformly in
$\overline{\Omega}^{\prime}$. We choose $t_{0}$ large enough such that for all
$t_{k}\geqslant t_{0}$

$\boldsymbol{\big{\|}}\vec{u}(\cdot,t_{k})-\vec{g}(\cdot)\boldsymbol{\big{\|}}_{\boldsymbol{1},\Omega^{\prime}}\leqslant\frac{\varepsilon}{2},$
(4.8)

where
$\boldsymbol{\big{\|}}\cdot\boldsymbol{\big{\|}}_{\boldsymbol{1},\Omega^{\prime}}$
corresponds to the $L^{1}$ norm in $\Omega^{\prime}$. We deduce that, in view
of (4.6) and (4.7) that

$\boldsymbol{\big{\|}}\vec{u}(\cdot,t_{k})-\vec{g}(\cdot)\boldsymbol{\big{\|}}_{\boldsymbol{1},\Omega\setminus\Omega^{\prime}}\leqslant
K\big{|}\Omega\setminus\Omega^{\prime}\big{|}\leqslant\frac{\varepsilon}{2},$

which together with (4.8) yields

$\boldsymbol{\big{\|}}\vec{u}(\cdot,t_{k})-\vec{g}(\cdot)\boldsymbol{\big{\|}_{1}}\leqslant\varepsilon.$

Step 2. $\omega(\vec{u}_{0})\subset\mathscr{H}$.
Indeed, let $\vec{g}\in\omega(\vec{u}_{0})$ and suppose
$\vec{g}\notin\mathscr{H}$. According to Proposition 4.2 we can find a steady
state solution $\vec{w}\in\mathscr{H}$, such that at least one component of
$\vec{w}-\vec{g}$ changes the sign. Without loss of generality we can assume
that it happens for the first component, namely that $f_{1}-w_{1}$ changes the
sign. We remark that, by the contraction property in Theorem 3.1, the
functional

$\mathcal{V}(\vec{\xi})=\boldsymbol{\big{\|}}\vec{\xi}-\vec{w}\boldsymbol{\big{\|}_{1}}$

is a Lyapunov functional for Problem (P), where
$\vec{\xi}\in\big{(}L^{1}(\Omega)\big{)}^{n}$. Next we describe some of its
properties.

Property (a) The functional $\mathcal{V}$ is constant on
$\omega(\vec{u}_{0})$.
Since $\mathcal{T}(t)\,\vec{w}=\vec{w}$ and $\mathcal{T}(t)\,$ has the
contraction property (3.4), the functional $\mathcal{V}$ is nonincreasing in
time along the trajectory $\mathcal{T}(t)\,\vec{u}_{0}$, which yields

$\begin{split}\mathcal{V}\big{(}\mathcal{T}(t)\,\vec{u}_{0}\big{)}&=\boldsymbol{\big{\|}}\mathcal{T}(t)\,\vec{u}_{0}-\vec{w}\boldsymbol{\big{\|}_{1}}\\\
&=\boldsymbol{\big{\|}}\mathcal{T}(t)\,\vec{u}_{0}-\mathcal{T}(t)\,\vec{w}\boldsymbol{\big{\|}_{1}}\leqslant\boldsymbol{\big{\|}}\vec{u}_{0}-\vec{w}\boldsymbol{\big{\|}_{1}}<\infty\hskip
5.0pt.\end{split}$

Thus there exists a finite limit $\mathcal{V}^{\ast}$ of
$\mathcal{V}\big{(}\mathcal{T}(t)\,\vec{u}_{0}\big{)}$ as
$t\rightarrow\infty$. Let $\vec{h}_{1},\ \vec{h}_{2}\in\omega(\vec{u}_{0})$.
We can find a sequence $t_{k}\rightarrow\infty$ as $k\rightarrow\infty$, such
that

$\boldsymbol{\big{\|}}\mathcal{T}(t_{2k})\,\vec{u}_{0}-\vec{h}_{1}\boldsymbol{\big{\|}_{1}}\rightarrow
0\quad\text{and}\quad\boldsymbol{\big{\|}}\mathcal{T}(t_{2k+1})\,\vec{u}_{0}-\vec{h}_{2}\boldsymbol{\big{\|}_{1}}\rightarrow
0,$

as $k$ tends to $\infty$. It follows that
$\mathcal{V}\big{(}\vec{h}_{1}\big{)}=\mathcal{V}\big{(}\vec{h}_{2}\big{)}=\mathcal{V}^{\ast}$.

Property (b) The $\omega$-limit set $\omega(\vec{u}_{0})$ is invariant with
respect to the semigroup $\mathcal{T}(t)\,$, namely if
$\vec{h}\in\omega(\vec{u}_{0})$, then for all $t>0$ also
$\mathcal{T}(t)\,\vec{h}\in\omega(\vec{u}_{0})$.
Let the sequence $t_{k}\rightarrow\infty$ as $k\rightarrow\infty$ be such that
$\boldsymbol{\big{\|}}\mathcal{T}(t_{k})\,\vec{u}_{0}-\vec{h}\boldsymbol{\big{\|}_{1}}\rightarrow
0$. From the contraction property (3.4)

$\displaystyle\boldsymbol{\big{\|}}\mathcal{T}(t_{k}+t)\,\vec{u}_{0}-\mathcal{T}(t)\,\vec{h}\boldsymbol{\big{\|}_{1}}=\boldsymbol{\big{\|}}\mathcal{T}(t)\,\mathcal{T}(t_{k})\,\vec{u}_{0}-\mathcal{T}(t)\,\vec{h}\boldsymbol{\big{\|}_{1}}$
$\displaystyle\leqslant\boldsymbol{\big{\|}}\mathcal{T}(t_{k})\,\vec{u}_{0}-\vec{h}\boldsymbol{\big{\|}_{1}}.$

Since the last term above tends to $0$ as $k$ tends to $\infty$ this shows
that $\mathcal{T}(t)\,\vec{h}\in\omega(\vec{u}_{0})$.

Now, remember that $\vec{g}\in\omega(\vec{u}_{0})$ is such that
$\vec{g}\notin\mathscr{H}$ and $\vec{w}\in\mathscr{H}$ is such that the first
component of $\vec{w}-\vec{g}$ changes the sign in $\Omega$. Then, Corollary
3.3 yields

$\rule{20.0pt}{0.0pt}\mathcal{V}(\mathcal{T}(t)\,\vec{g})=\boldsymbol{\big{\|}}\mathcal{T}(t)\,\vec{g}-\vec{w}\boldsymbol{\big{\|}_{1}}\\\
=\boldsymbol{\big{\|}}\mathcal{T}(t)\,\vec{g}-\mathcal{T}(t)\,\vec{w}\boldsymbol{\big{\|}_{1}}<\boldsymbol{\|}\vec{g}-\vec{w}\boldsymbol{\|_{1}}=\mathcal{V}(\vec{g}),\rule{20.0pt}{0.0pt}$

for all $t>0$, which contradicts Property (a). Therefore
$\vec{g}\in\mathscr{H}$.

Step 3. The set $\omega(\vec{u}_{0})$ contains only one element.
Suppose that $\vec{g}_{1},\vec{g}_{2}\in\omega(\vec{u}_{0})$. Then we can find
two sequences $t_{k},\,s_{k}$ tending to $\infty$ as $k\rightarrow\infty$,
such that $s_{k}\leqslant t_{k}$ and
$\boldsymbol{\big{\|}}\mathcal{T}(t_{k})\,\vec{u}_{0}-\vec{g}_{1}\boldsymbol{\big{\|}_{1}},\,\boldsymbol{\big{\|}}\mathcal{T}(s_{k})\,\vec{u}_{0}-\vec{g}_{2}\boldsymbol{\big{\|}_{1}}\rightarrow
0$ as $t_{k}\rightarrow\infty$. Since $\omega(\vec{u}_{0})\subset\mathscr{H}$,
it follows that

$\displaystyle\boldsymbol{\|}\vec{g}_{1}-\vec{g}_{2}\boldsymbol{\|_{1}}\leqslant\boldsymbol{\big{\|}}\mathcal{T}(t_{k})\,\vec{u}_{0}-\vec{g}_{1}\boldsymbol{\big{\|}_{1}}+\boldsymbol{\big{\|}}\mathcal{T}(t_{k})\,\vec{u}_{0}-\vec{g}_{2}\boldsymbol{\big{\|}_{1}}$
$\displaystyle=\,\boldsymbol{\big{\|}}\mathcal{T}(t_{k})\,\vec{u}_{0}-\vec{g}_{1}\boldsymbol{\big{\|}_{1}}+\boldsymbol{\big{\|}}\mathcal{T}(t_{k}-s_{k})\,\mathcal{T}(s_{k})\,\vec{u}_{0}-\mathcal{T}(t_{k}-s_{k})\,\vec{g}_{2}\boldsymbol{\big{\|}_{1}}$
$\displaystyle\leqslant\boldsymbol{\big{\|}}\mathcal{T}(t_{k})\,\vec{u}_{0}-\vec{g}_{1}\boldsymbol{\big{\|}_{1}}+\boldsymbol{\big{\|}}\mathcal{T}(s_{k})\,\vec{u}_{0}-\vec{g}_{2}\boldsymbol{\big{\|}_{1}}\hskip
3.0pt,$

which tends to $0$ as $k\rightarrow\infty$. $\blacksquare$

##  5 Stationary solutions for the linear molecular motor problem

In this section we show the existence and the uniqueness (up to a
multiplicative constant) of the classical stationary solution of the problem
for the molecular motor. We suppose that $\Omega$ is an open bounded subset of
$\mathbb{R}^{d}$ with smooth boundary $\partial\Omega$.

We consider the linear system

$\text{div}\big{(}\sigma_{i}\nabla
v_{i}(x)+v_{i}(x)\nabla\psi_{i}(x)\big{)}+\sum_{j=1}^{n}a_{ij}v_{j}(x)=0\quad\text{in}\quad\Omega,$
(5.1)

where $i\in\\{1,\ldots,n\\}$, $n>1$. The system (5.1) is supplemented with the
Robin boundary conditions

$\sigma_{i}\frac{\partial
v_{i}}{\partial\nu}+v_{i}\frac{\partial\psi_{i}}{\partial\nu}=0\quad\text{on}\quad\partial\Omega,$
(5.2)

where $i\in\\{1,\ldots,n\\}$. Thus, the problem can be written as

${\cal A}\vec{v}=0,$

with a linear operator $\cal A$ in a suitable Banach space $\cal X$ of
functions on $\Omega$, to be made precise later. Moreover, we impose the
integral constraint

$\sum_{i=1}^{n}\int_{\Omega}v_{i}(x)\,\text{d}x=1.$ (5.3)

The adjoint problem ${\cal A}^{\ast}\vec{\varphi}=0$ to (5.1), in a dual space
${\cal X}^{\ast}$, is now

$\sigma_{i}\Delta\varphi_{i}-\nabla\psi_{i}\cdot\nabla\varphi_{i}+\sum_{j=1}^{n}a_{ji}\varphi_{j}=0,\quad\text{in}\quad\Omega,$
(5.4)

with the Neumann boundary conditions for each $i=1,\,\dots,\,n$

$\frac{\partial\varphi_{i}}{\partial\nu}=0\quad\text{on}\quad\partial\Omega.$
(5.5)

Since $\sum_{j=1}^{n}a_{ji}=0$, the problem (5.4) has the obvious solution

$\vec{\varphi}=(\varphi_{1},\dots,\varphi_{n})=(1,\dots,1).$ (5.6)

We are going to apply the Krein-Rutman theorem on the first eigenvalues and
eigenvectors of positive operators, and this will permit us to conclude that
the problem (5.1)–(5.2) has a one-dimensional space of solutions. Therefore,
under the additional constraint (5.3), the original problem (5.1)–(5.2) has a
unique solution.
Perthame and Souganidis sketched this argument for $n>1$ and $d=1$ in [20].

######  Theorem 5.1.

Under the assumption $\sum_{j=1}^{n}a_{ji}=0$, there exists a unique smooth
solution $\vec{v}$ of the system (5.1)–(5.3).

Before proving Theorem 5.1 we recall some basic definitions as well as the
Krein-Rutman theorem from [9, Ch. VIII, p. 188–191].

######  Definition 5.2 (Reproducing cone).

We say that a closed set $K$ in $\cal X$ is a cone, if it possesses the
following properties:

  * i)

$0\in K$,

  * ii)

$u,\ v\in K\Longrightarrow\alpha u+\beta v\in K$, for all $\alpha,\
\beta\geqslant 0$,

  * iii)

$v\in K$ and $-v\in K\Longrightarrow v=0$.

A cone $K\subset\cal X$ is said to be reproducing if ${\cal
X}=K-K\equiv\big{\\{}k_{1}-k_{2}:\ k_{1},\ k_{2}\in K\big{\\}}$.

######  Definition 5.3 (Dual cone).

If $K$ is a cone in $\cal X$, then the set $K^{\ast}\subset{\cal X}^{\ast}$ is
said to be a dual cone if

$\langle f^{\ast},v\rangle\geqslant 0,$

for every $v\in K$.

######  Definition 5.4 (Strict positivity).

Let $\cal B$ be a linear operator on $\cal X$. Then $\cal B$ is said to be
strongly positive if ${\cal B}v\in K^{o}$ for all $v\in K$ such that $v\neq
0$.

######  Theorem 5.5.

Let $K$ be a reproducing cone in a Banach space $\cal X$, with nonempty
interior $K^{o}\neq\emptyset$, and let $\cal B$ be a strongly positive compact
operator on $K$ in a sense of Definition 5.4. Then the spectral radius of
$\cal B$, $r(\cal B)$, is a simple eigenvalue of $\cal B$ and $\cal B^{\ast}$,
and their associated eigenvectors belong to $K^{o}$ and $(K^{\ast})^{o}$. More
precisely, there exists a unique associated eigenvector in $K^{o}$ (resp.
$(K^{\ast})^{o}$) of norm $1$. Furthermore, all other eigenvalues are strictly
less in absolute value than $r(\cal B)$.

Proof  We will apply Theorem 5.5 to the space ${\cal
X}=\big{(}C(\overline{\Omega})\big{)}^{n}\subset\big{(}L^{1}(\Omega)\big{)}^{n}$
endowed with the usual supremum norm, and the operators

$\displaystyle{\cal B}=(\lambda I-{\cal A})^{-1}:\ {\cal X}\to{\cal X},$
$\displaystyle{\cal B}^{\ast}=(\lambda I-{\cal A}^{\ast})^{-1}:\ {\cal
X}^{\ast}\to{\cal X}^{\ast},$

where $\lambda>0$ is a strictly positive real number to be fixed later.

Let

$K=\big{\\{}\vec{u}\in{\cal X}:\ u_{i}(x)\geqslant 0\ {\rm for\ each\
}x\in\overline{\Omega},\ i=1,\dots,n\big{\\}}.$

We remark that $K$ is a reproducing cone, with nonempty interior

$K^{o}=\big{\\{}\vec{u}\in{\cal X}:\ \inf_{x\in\overline{\Omega}}u_{i}(x)>0,\
i=1,\dots,n\big{\\}}.$

From the standard theory [17, Theorem 2.1 and Theorem 3.1, Ch. 7] for elliptic
partial differential linear systems, the boundary value problem

$\sigma_{i}\Delta\varphi_{i}-\nabla\psi_{i}\cdot\nabla\varphi_{i}+\sum_{i=1}^{n}a_{ji}\varphi_{j}-\lambda\varphi_{i}=f_{i}\
\ {\rm in\ \ }\Omega,$ (5.7)

with the homogeneous Neumann conditions (5.5) on $\partial\Omega$, for
$\lambda=\widetilde{\lambda}>0$ sufficiently large, has a solution
$\vec{\varphi}=(\varphi_{1},\dots,\varphi_{n})\in{\cal X}$ for each
$\vec{f}=(f_{1},\dots,f_{n})\in{\cal X}$. Moreover, if $f_{i}(x)\geqslant 0$
for each $i=1,\dots,n$, and $x\in\overline{\Omega}$, then
$\varphi_{i}(x)\geqslant 0$ (in fact, $\varphi_{i}(x)>0$ in $\Omega$), which
is a consequence of the maximum principle (cf. also Example 3 on p. 196–197 in
[9]). Thus, the operator ${\cal B}^{\ast}=\big{(}\widetilde{\lambda}I-{\cal
A}^{\ast}\big{)}^{-1}$ is a strongly positive and compact operator, and by
Theorem 5.5, the largest eigenvalue $\mu$ of $\cal B$ and ${\cal B}^{\ast}$ is
simple.

Since

$\displaystyle-\sigma_{i}\Delta\varphi_{i}+\nabla\psi_{i}\cdot\nabla\varphi_{i}-\sum_{j=1}^{n}a_{ji}\varphi_{j}+\widetilde{\lambda}\varphi_{i}=\widetilde{\lambda}\varphi_{i}\quad\text{in}\quad\Omega$
$\displaystyle\frac{\partial\varphi_{i}}{\partial\nu}=0\quad\text{on}\quad\partial\Omega,$

for all $i\in\\{1,\ldots,n\\}$, with
$\vec{\varphi}=(\varphi_{1},\ldots,\varphi_{n})=(1,\ldots,1)$, and since
$(1,\ldots,1)\in(K^{\ast})^{o}$, it follows that
$\displaystyle\frac{1}{\widetilde{\lambda}}=r\Big{(}\big{(}\widetilde{\lambda}I-{\cal
A}^{\ast}\big{)}^{-1}\Big{)}$ is a simple eigenvalue of the operator
$\big{(}\widetilde{\lambda}I-{\cal A}^{\ast}\big{)}^{-1}$. Applying again
Theorem 5.5, we deduce that $\displaystyle\frac{1}{\widetilde{\lambda}}$ is
the largest eigenvalue of the operator $\big{(}\widetilde{\lambda}I-{\cal
A}\big{)}^{-1}$ and that it is simple, and that there exists $\vec{v}\in
K^{o}\subset\cal X$ such that

$\Big{(}\widetilde{\lambda}I-{\cal
A}\Big{)}^{-1}\vec{v}=\frac{1}{\widetilde{\lambda}}\,\vec{v},$

which is equivalent to

${\cal A}\vec{v}=0.$

This proves the existence of the solution of the problem (5.1)–(5.3).
$\blacksquare$

## References

  * [1] R. A. Adams, Sobolev Spaces, Pure and Applied Mathematics, vol. 65, Academic Press, New York-London, 1975\.
  * [2] N. D. Alikakos, P. Hess, H. Matano, Discrete order preserving semigroups and stability for periodic parabolic differential equations, J. Differential Equations, 1989, 2, 82, 322–341.
  * [3] M. Bertsch and D. Hilhorst, A density dependent diffusion equation in population dynamics: stabilization to equilibrium, SIAM J. Math. Anal., 1986, 4, 17, 863–883.
  * [4] D. Bothe, Instantaneous limits of reversible chemical reactions in presence of macroscopic convection, J. Differential Equations, 2003, 1, 193, 27–48.
  * [5] D. Bothe and D. Hilhorst, A reaction-diffusion system with fast reversible reaction, J. Math. Anal. Appl., 2003, 1, 268, 125–135.
  * [6] M. Chipot, S. Hastings, and D. Kinderlehrer, Transport in a molecular motor system, M2AN Math. Model. Numer. Anal., 2004, 38, 6, 1011- 1034.
  * [7] M. Chipot, D. Kinderlehrer, and M. Kowalczyk, 2004 A variational principle for molecular motors, Meccanica, 38, 505–518
  * [8] M. G. Crandall and L. Tartar, Some relations between nonexpansive and order preserving mappings, Proc. Amer. Math. Soc., 1980, 3, 78, 385–390.
  * [9] R. Dautray, J.-L. Lions, Mathematical Analysis and Numerical Methods for Science and Technology, Volume 3: Spectral Theory and Applications, Springer-Verlag, Berlin, 1990.
  * [10] L. Desvillettes and K. Fellner, Exponential decay toward equilibrium via entropy methods for reaction-diffusion equations, J. Math. Anal. Appl., 2006, 1, 319, 157–176.
  * [11] P. Érdi and J. Tóth, Mathematical Models of Chemical Reactions, Nonlinear Science: Theory and Applications, Princeton University Press, Princeton, NJ, 1989\.
  * [12] Hastings, S., Kinderlehrer, D. and McLeod, J.B. Diffusion mediated transport in multiple state systems, An elliptic-parabolic problem in combustion theory: convergence to travelling waves, SIAM J. Math. Anal., 2007, 39.4, 1208-1230.
  * [13] Hastings, S., Kinderlehrer, D. and McLeod, J.B. Diffusion mediated transport with a look at motor proteins, Recent advances in nonlinear analysis, (Chipot, M., Lin, C-S,., and Tsai, D-H., eds) (2008) World Scientific, 95 - 112.
  * [14] D. Hilhorst and J. Hulshof, An elliptic-parabolic problem in combustion theory: convergence to travelling waves, Nonlinear Analysis TMA, 1991, 6, 17, 519–546.
  * [15] D. Hilhorst and M. A. Peletier, Convergence to travelling waves in a reaction-diffusion system arising in contaminant transport, J. Differential Equations, 2000, 1, 163, 89–112.
  * [16] O. A. Ladyženskaja and V. A. Solonnikov and N. N. Ural’ceva, Linear and Quasilinear Equations of Parabolic Type, Translations of Mathematical Monographs, Vol. 23, American Mathematical Society, Providence, R.I., 1967\.
  * [17] O. A. Ladyženskaja and N. N. Ural’ceva, Linear and Quasilinear Elliptic Equations, Academic Press, New York, 1968\.
  * [18] S. Osher and J. Ralston, $L^{1}$ stability of travelling waves with applications to convective porous media flow, Comm. Pure Appl. Math., 1982, 6, 35, 737–749.
  * [19] B. Perthame, The general relative entropy principle—applications in Perron-Frobenius and Floquet theories and a parabolic system for biomotors, Rend. Accad. Naz. Sci. XL Mem. Mat. Appl. (5), 2005, 29, 307–325.
  * [20] B. Perthame, P. E. Souganidis, Asymmetric potentials and motor effect: a large deviation approach, 1–16, Arch. Rat. Mech. Anal., to appear, 2008.
  * [21] M. H. Protter and H. F. Weinberger, Maximum Principles in Differential Equations, Springer-Verlag, New York, 1984\.