format: text

text:
# 1 Introduction. A model for contour enhancement

## 1 Introduction. A model for contour enhancement

This paper is devoted to study the behaviour of the solutions of some degenerate parabolic equations with moving boundaries which appear in describing the technique of contour enhancement in image processing. Indeed, computer vision has become in recent decades a mathematical discipline which relies on the differential-geometric approach. More specifically, an appropriate technique of image processing consists of formulating a partial differential equation of evolution type for the image intensity, $`u(x,y)`$. This function, also called the grey level, takes values in the interval $`0u1`$ and is defined on a two-dimensional image domain, $`\mathrm{\Omega }`$. The usual evolution model leads to a nonlinear equation of parabolic type, possibly degenerate or singular. The nonlinearity is created by the law relating the image intensity flux to the image intensity.

It has been observed by Perona and Malik \[PM\], 1990, that for a suitable choice of this fundamental nonlinearity there appears an effect of enhancement of image edges that has a strong interest in the application to processing, denoising and recognition of images. The Perona-Malik Anisotropic Diffusion model has had a deep influence in the field, being at the source of many later developments. The model proposed by Malladi and Sethian \[MS\] leads (after proper scaling) to the following equation for the image intensity:

(1.1)
$$u_t=(1+|Du|^2)^{1/2}\kappa $$

where $`Du`$ denotes the spatial gradient of $`u`$, and $`\kappa `$ denotes the curvature of the surface $`z=u(x,y)`$. The equation represents movement by curvature (curvature flow) and can be written as

(1.2)
$$u_t=\frac{(1+u_y^2)u_{xx}2u_xu_yu_{xy}+(1+u_x^2)u_{yy}}{1+u_x^2+u_y^2}.$$

We consider here the more general flow given by equation $`u_t=(1+|Du|^2)^{(1+p)/2}\kappa ,`$ where $`p`$ is a constant parameter. In other words, we study the equations

(1.3)
$$u_t=\frac{(1+u_y^2)u_{xx}2u_xu_yu_{xy}+(1+u_x^2)u_{yy}}{(1+u_x^2+u_y^2)^{1+\alpha }}$$

with parameter $`\alpha =p/2`$. Along with the former case $`p=\alpha =0`$, the case $`p=2`$, $`\alpha =1`$ has also attracted the attention of researchers (Beltrami flow, cf. Sochen et al. \[SKMS\]).

The asymptotic treatment of these models done in the papers \[MS, SKMS\] shows the enhancement of the intensity contrasts by formation of regions of large intensity gradients, i.e., the normal component of the image intensity gradient becomes quite large. This phenomenon was analyzed in \[B01\], where a further simplification of the model was proposed in order to focus on the boundary layer where large gradients concentrate. Arguing locally around a sharp gradient point and choosing the $`x`$-axis as the direction normal to the boundary layer or front, we may disregard the effect of $`y`$ derivatives with respect to the $`x`$ derivatives in (1.3). In this way we get the reduced equation, which is just the one-dimensional version of (1.3)

(1.4)
$$u_t=\frac{u_{xx}}{(1+u_x^2)^{1+\alpha }},$$

where we have neglected $`u_y,u_{yy}`$. Different dimensional constants appear in the model, but they have been scaled to unity here without loss of generality. The mathematical problem consists in solving this equation with suitable boundary data, namely, $`u=0`$ on the left-hand side of the contour and $`u=1`$ on the right-hand side (be that a finite or an infinite distance), and initial conditions

$$u(x,0)=u_0(x),$$

satisfying $`0<u_0<1`$ and $`u_0^{}>0`$ in an interval $`I=(a,b)`$ and constant values otherwise, zero to the left, 1 to the right. As was pointed out in \[B01\], the phenomenon of gradient enhancement takes place in this model (in a proper setting) for all $`\alpha 0`$: the spatial gradient of the solutions, $`u_x`$, increases with time, and its support shrinks. Indeed, we show below that such a behaviour can be observed in the larger exponent range $`\alpha >1/2`$. The conditions on $`u_0`$ can be relaxed, but then problems arise. Thus, the less stringent size restriction $`0u_01`$ or the lack of monotonicity create interesting alternatives, that will also be briefly discussed.

Let us remark that this is not the only model which uses nonlinear diffusion equations for image processing and enhancement. We refer for alternatives to the works of Alvarez, Lions and Morel \[ALM\] and Caselles et al. \[CMS\], where further references are found. Let us point out that in \[ALM\] a fundamentally different assumption is made, namely that the flux is perpendicular to the gradient, and the basic equation is substantially different from the class (1.3) treated here. On the other hand, we mention that the mathematical difficulties of the original model by Perona and Malik have been further investigated by several authors, like \[Kch\] and \[WB\]. The application of Nonlinear Diffusion to Image Processing is a every active concern with many issues being discussed in the literature, cf. e.g. \[KB, Wck\].

## 2 Asymptotic self-similarity and enhancement

Evolution equations like (1.4) and many other variants have been studied and are known in the literature under the general name of nonlinear parabolic equations of diffusion type, or nonlinear diffusion equations for short. They are typically used in describing processes of mass diffusion or thermal propagation. A quite general one-dimensional form popular among PDE experts is

$$u_t=(a(x,t,u,u_x))_x$$

with suitable conditions to make it parabolic, at least in a formal sense, like $`a/u_x>0`$. However, many practical applications (as in the present case) involve functions $`a`$ for which $`a/u_x0`$, but the values zero or infinity can also be taken, and then the equations are known as degenerate parabolic or singular parabolic resp., cf. \[DB, Ka\]. In any case, a general feature of this wide class of equations is the diffusive character, which roughly means the spreading of the level sets of the solutions as time advances. This property goes squarely against the desired enhancement, therefore an extra mechanism must be present if enhancement is to occur. We recall that in the Perona-Malik model this mechanism was negative diffusion.

In the range of exponents $`\alpha 0`$, equation (1.4) falls into the class of degenerate parabolic equations with degeneracy at $`u_x=\mathrm{}`$, precisely the limit value which is of concern in contour enhancement. An investigation of the phenomenon of gradient enhancement in equation (1.4) is performed in the paper \[B01\] by relating it to the convergence toward self-similar asymptotics of an approximate model. It goes as follows: the observed evolution of the solutions towards a configuration with large gradients makes it plausible to further simplify the expression $`1+u_x^2`$ in (1.4) into $`u_x^2`$, so that the relevant reduced equation becomes

(2.1)
$$u_t=u_x^{2(1+\alpha )}u_{xx}.$$

It is further observed that self-similar solutions for this equation with end-levels $`u=0`$ and $`u=1`$ exist for all $`\alpha >0`$ and exhibit the similarity form

(2.2)
$$u(x,t)=F(\xi ),\xi =(xx_0)(t+t_0)^{1/2\alpha },$$

where $`x_0`$ and $`t_0`$ are parameters to be fixed, and the profile $`F(\xi )`$ is an increasing function joining the levels $`F=0`$ at a finite distance $`\xi =c<0`$ to the level $`F=1`$ at $`\xi =c`$. At these levels, which are taken at a finite distance $`\xi =\pm c`$, the gradients are infinite (actually, the analysis in \[B01\] deals with decreasing profiles joining $`u=1`$ to $`u=0`$ but these two problems are obviously equivalent after a mirror symmetry. We have chosen the increasing option to avoid chasing around many minus signs). The scaling implies that

$$u_x(x,t)=(t+t_0)^{1/2\alpha }F^{}(\xi ),\xi =(xx_0)(t+t_0)^{1/2\alpha },$$

which shows that the solution is concentrated in an increasingly narrower strip

(2.3)
$$S=\{(x,t):|xx_0|c(t+t_0)^{1/2\alpha }\}$$

with gradients that diverge like $`t^{1/2\alpha }`$ as $`t\mathrm{}`$. Let us remark that the asymptotic divergence of the gradients is known to be exponential for $`\alpha =0`$, but the form is not self-similar of the same type, as we will see below. An important feature to be remembered of these solutions, along with the divergence as $`t\mathrm{}`$, is the infinite gradient condition at the end-points of the domain of definition. On the other hand, these self-similar solutions represent the intermediate asymptotics of the problem with more general data. This is demonstrated numerically in \[B01\] when infinite flux data are imposed at moving end-points located at finite distance.

With these preliminaries we are ready to attack the construction and analysis of solutions with steep fronts for a general class of equations which includes (2.1) and (1.4) and is natural for our application, and with the general data mentioned in the Introduction. We use the theory of nonlinear diffusion equations which has been strongly developed in the last decades. It turns out that the standard initial value problem solved in this theory does not suit our model, since it leads to dispersion instead of concentration of the image. A new free boundary problem with singular boundary data is then introduced and solved by means of a number of non-trivial problem transformations available for nonlinear diffusion problems, a subject that is developed in detail in \[V03c\]. The concentration property of these solutions is possible thanks to the degenerate character of the diffusivity for high gradient values.

Our theory below covers existence and uniqueness of solutions of suitable problems, existence and behaviour of the bounding interfaces, and large-time behaviour. In particular, we obtain the rates of convergence of the interfaces, an important question for image enhancement.

## 3 Nonlinear diffusion equations

As we have just said, we are going to construct a mathematical theory for the one-dimensional evolution problem with initial and boundary conditions of the type mentioned above, but replacing the special class of equations (1.4) by the wider class of nonlinear diffusion equations of the form

(3.1)
$$u_t=(\mathrm{\Phi }(u_x))_x.$$

The nonlinear function $`\mathrm{\Phi }`$ expresses the dependence of the image flux on the gradient of the image intensity. It can be called the flux law of our process, or also the constitutive function. In this general setting, $`\mathrm{\Phi }`$ can be any increasing real function defined in a suitable interval of gradients. We will assume for convenience that

(H1)$`\mathrm{\Phi }(s)`$ is $`C^{\mathrm{}}`$ smooth and strictly increasing in the intervals $`0<s<\mathrm{}`$ and $`\mathrm{}<s<0`$, allowing for a quite arbitrary behaviour as $`s0`$ or $`|s|\mathrm{}`$.

In our opinion, this level of generality instead of obscuring the problem makes more apparent the relation between nonlinear diffusion and contour enhancement. Slightly more general assumptions can be made on $`\mathrm{\Phi }`$, but they are unessential for our present purposes, cf. Section 11.

We want to characterize the class of those $`\mathrm{\Phi }`$ for which we can construct solutions with gradient-enhancement. We also want to describe the rates at which the formation of steep profiles takes place, thus justifying and extending the results of \[B01\] to the actual equations like (1.4), and to general data. The mathematical problem we pose consists in solving the nonlinear parabolic equation (3.1) with initial conditions

(3.2)
$$u(x,0)=u_0(x)$$

satisfying $`0u_01`$ and other suitable conditions (see below), and boundary data $`u=0`$ on the left-hand side of the contour and $`u=1`$ on the right-hand side. Regarding the initial data, we remark that we are mainly interested in monotone solutions, i.e., $`u_x0`$, for equation (1.4), and this condition will naturally follow from a similar monotonicity condition on the initial data. As for the boundary conditions, it turns out that, depending on the form of $`\mathrm{\Phi }`$, the boundary can be chosen to be located either at infinity or at a finite distance. This latter case will be the one of interest for us, and then the problem must be properly posed as a free-boundary problem.

Let us recall the simplest examples of $`\mathrm{\Phi }`$ and some of the difficulties we will encounter. Indeed, a quite important and simple example is the power function, that we write as $`\mathrm{\Phi }(s)=(1/m)s^m`$, defined for $`s0`$, so that it agrees with (2.1) with $`m=12\alpha `$ since $`\mathrm{\Phi }^{}(s)=s^{m1}`$. There is in principle no reason to restrict the generality of the exponent $`m`$ in the mathematical treatment to follow, and this will lead to quite different behaviour types inside this family. The case $`m=0`$ is included in the form $`\mathrm{\Phi }(s)=\mathrm{log}(s)`$, i.e., $`\mathrm{\Phi }^{}(s)=1/s`$. Note that in the cases with $`m<0`$ the function $`\mathrm{\Phi }`$ is negative, but the important quantity for the parabolic character of the equation, i.e., $`\mathrm{\Phi }^{}(u_x)`$, is always positive. Finally, we note that equation (1.4) corresponds to $`\mathrm{\Phi }^{}(s)=(1+s^2)^{(1+\alpha )}`$; it degenerates as $`s\mathrm{}`$ in the prescribed range $`\alpha 0`$, even if $`\alpha >1`$, but the equation is is perfectly parabolic in regions of bounded $`u_x`$.

### 3.1 Non-monotone solutions

The above statements are made on the assumption that the solutions are monotone, $`u_x0`$, which is not unjustified in our problem setting but restricts the mathematical generality. When considering non-monotone solutions it is customary in the nonlinear diffusion literature to extend the power nonlinearity to arguments $`s=u_x<0`$ in the simplest symmetric way:

(3.3)
$$\mathrm{\Phi }(s)=\frac{1}{m}|s|^{m1}s,$$

which gives $`\mathrm{\Phi }^{}(s)=|s|^{m1}`$, always nonnegative. This definition poses no problem when $`m>0`$ and solutions with changing sign exist corresponding to initial data with the same property. However, we are interested in exponents $`m<0`$ (so-called very fast diffusion in the literature) where $`\mathrm{\Phi }^{}`$ is singular at $`s=0`$, and the whole function $`\mathrm{\Phi }`$ is no more monotone. The difficulty has been studied in \[RV02\] leading to the consequence that solutions with changing sign do not exist, even in the weak sense when $`m0`$. This is a consequence of the singularity of the equation at the level $`u_x=0`$ and does not affect equation (1.4). Due to this obstruction, we will concentrate here on problems with monotone solutions.

## 4 Second formulation as a nonlinear diffusion equation

If we formally differentiate equation (3.1) with respect to $`x`$ and put $`v=u_x`$, we obtain the equation satisfied by the image intensity gradient:

(4.1)
$$v_t=\mathrm{\Phi }(v)_{xx},$$

which is usually called the nonlinear filtration equation, NLFE, and is the most standard class of nonlinear diffusion equations studied in the literature. We can call it in this context the “differentiated equation”. Inversely, we can recover the $`u`$ formulation from a solution $`v(x,t)`$ for (4.1) by means of the rule

(4.2)
$$u(x,t)=c+_\mathrm{\Gamma }(vdx+\mathrm{\Phi }(v)_xdt),$$

integrated along any curve $`\mathrm{\Gamma }`$ in the domain of definition of $`v`$ which joins a fixed point $`(x_0,t_0)`$ to the generic point $`(x,y)`$, cf. \[RV95, V03c\]. The constant $`c`$ is the value of $`u`$ at $`(x_0,t_0)`$ to be chosen at will in principle. The calculation is justified for smooth solutions $`v`$ and smooth $`\mathrm{\Phi }`$, but holds in a much wider context.

Monotone solutions $`u`$ for equation (1.4) translate into nonnegative solutions $`v`$ for equation (4.1), and conversely. In this context, $`u`$ is usually viewed as the mass function for $`v`$, since when we take $`\mathrm{\Gamma }`$ to be a segment of the line $`t=t_0`$, formula (4.2) becomes

$$u(x,t_0)u(x_0,t_0)=_{x_0}^xv(y,t_0)𝑑y.$$

The phenomenon of gradient enhancement can be then translated into usual diffusion parlance as mass concentration. We will keep in the sequel the denomination intensity or image intensity for the solution $`u`$ of equation (3.1), and we will view the solution $`v=u_x`$ of (4.1) as the intensity gradient. In that context $`\mathrm{\Phi }(u_x)=\mathrm{\Phi }(v)`$ is the intensity flux. Finally, $`\mathrm{\Phi }^{}(u_x)=\mathrm{\Phi }^{}(v)`$ is the diffusivity.

## 5 Solutions of Type I

We now try to apply the more standard theory of nonlinear diffusion equation to solve the evolution problem motivated by the application in the Introduction. We show next that this implementation can be performed in a rather standard way for a wide class of $`\mathrm{\Phi }`$’s, but such a process gives an evolution with no concentration effect, hence useless for our purposes in image processing. The reader already familiar with diffusion theory may check the contents of Theorem 5.1 and skip the rest of this section.

### 5.1 Simplest setting

Generalizing the well-known properties of the heat equation (the choice $`\mathrm{\Phi }(s)=s`$), we consider first the case where $`\mathrm{\Phi }`$ is a $`C^1`$ function and $`\mathrm{\Phi }^{}`$ does not vanish. We pose the Cauchy problem for equation (3.1) in the whole line $`x`$ for $`t>0`$ with bounded initial data $`0u_01`$, and get a unique smooth solution $`u(x,t)`$ defined in $`Q=\times (0,\mathrm{})`$ and such that $`0<u<1`$, $`u`$ is smooth for all $`x`$ and $`t\tau >0`$. We can also work with equation (4.1): then, if $`v_0=u_{0,x}`$ is locally integrable (a Radon measure will also do) and such that

$$_{}v_0(x)𝑑x=1,$$

the solution $`v(x,t)`$ is smooth in $`Q`$ and satisfies the same type of integrability condition, $`_{}v(x,t)𝑑x=1.`$ For solutions $`v0`$ this is called conservation of mass. Using formula (4.2) we obtain a smooth solution $`u(x,t)`$ in $`Q`$ such that the initial condition

$$u(x,0)=u_0(x)=_{\mathrm{}}^xv_0(y)𝑑y$$

holds, as well as the end conditions

$$\underset{x\mathrm{}}{lim}u(x,t)=0,\underset{x\mathrm{}}{lim}u(x,t)=1,$$

hold locally uniformly in time. On the other hand, $`v_00`$ implies that $`u_x=v0`$, which by the Strong Maximum Principle implies $`u_x>0`$. It follows that $`0<u(x,t)<1`$ in $`Q`$. Note that, due to the smoothness condition on $`\mathrm{\Phi }`$, these solutions are smooth in $`(x,t)`$ for all $`x`$ and $`t>0`$, even if the initial data are not.

### 5.2 Dispersion of solutions with time

At first glance, the type of solution we have constructed seems to solve our image problem. However, it lacks a basic ingredient, i.e., the eventual concentration of intensity gradients. On the contrary, if we consider the heat equation $`u_t=u_{xx}`$, the solutions, which can be expressed in terms of error functions, spread in time and its gradients go to zero. for instance, if $`v_0=u_{0,x}`$ is integrable, then $`v(x,t)`$ goes to zero as $`t\mathrm{}`$ at a rate of the order of $`t^{1/2}`$, and takes on a Gaussian space shape,

$$v(x,t)t^{1/2}\text{exp}(x^2/4t).$$

And a similar result (with possibly different rates) applies to more general data $`u_0`$ having definite limits at $`\pm \mathrm{}`$, and to all the functions $`\mathrm{\Phi }`$ of the above class.

### 5.3 Extension of Type I solutions to other $`\mathrm{\Phi }`$

Even if this is not the class of solutions we are looking for, we pursue a bit further the analysis, since it is the standard class found in studies of nonlinear diffusion or thermal propagation. The study will also serve for comparison with the “correct” solutions of Type II. However, the reader may choose to skip this the rest of the section and proceed with Section 6.

The class of constitutive functions $`\mathrm{\Phi }`$ for which there exists a class solutions $`0u1`$ with a dispersive character can be extended to include a quite general choice of $`\mathrm{\Phi }`$ on the condition of allowing for a suitably generalized concept of solution. Thus, it is well-known following Bénilan, Crandall and others \[Be72, BC\] that we may take as $`\mathrm{\Phi }`$ any continuous nondecreasing function and the Cauchy problem is then well-posed in the class of so-called mild solutions with $`L^1()`$ data. Actually, $`\mathrm{\Phi }`$ can be allowed to be discontinuous but this will be of no particular interest here; on the contrary, we will keep the assumption of smoothness and strict monotonicity for $`0<s<\mathrm{}`$ for simplicity of presentation (and since it is satisfied in the application we are dealing with).

Let us comment on the main properties of the solutions for this class of equations. For the case of power nonlinearities $`\mathrm{\Phi }(s)=s^m/m,`$ $`s0`$, mentioned above, rather complete details are known, \[Ar, Ka, V92\]. We can consider that $`m`$ is a real parameter, positive in principle but not necessarily as we will see. The definition is extended to solutions with negative values as in formula (3.3). Thus, for $`m>0`$ equation $`v_t=(v^m)_{xx}`$ (or better, $`v_t=(|v|^{m1}v)_{xx}`$) generates a positive semigroup of contractions in the space $`L^1()`$. In other words, for every $`v_0L^1()`$ there exists a unique function $`vC([0,\mathrm{}:L^1())`$, $`v0`$ with $`v^mL_{loc}^1(Q)`$ such that the equation is satisfied in the sense of distributions in $`Q=\times (0,\mathrm{})`$, the initial data are taken in the $`L^1()`$-sense, and the map $`v_0v(,t)`$ is an $`L^1`$-contraction. Moreover, the total mass is conserved

$$_0^{\mathrm{}}v(x,t)𝑑x=_0^{\mathrm{}}v_0(x)𝑑x.$$

In our application we still have to impose the extra condition of total mass 1, and then we recover the image intensity by means of the formula

$$u(x,t)=_{\mathrm{}}^xv(x^{},t)𝑑x^{},$$

and the intensity level goes from $`u=0`$ at minus infinity to $`u=1`$ at infinity. The asymptotic behaviour of these solutions has been carefully calculated in the literature, cf. e.g. \[An2, FK, V03a\]: solutions $`v(,t)`$ with finite mass go zero uniformly as

(5.1)
$$v(x,t)t^\gamma ,\gamma =\frac{1}{m+1},$$

which in the notation of the Introduction means $`\gamma =1/2\alpha `$. Actually, the asymptotic rate comes from comparison with the source-type self-similar solutions \[Bbk\] which take the form

$$v(x,t)=t^\gamma G(xt^\gamma )$$

for a certain symmetric profile function $`GC_b()`$ such that $`G(\xi )0`$ as $`\xi \mathrm{}`$. This formula suggests that the behaviour will be maintained as long as $`\gamma >0`$, hence, as long as $`m>1`$, and this turns out to be true (for the study when $`1<m0`$ cf. \[ERV\]). Indeed, this exponent range is optimal since there are no solutions with finite integral for $`m1`$.

Let us look a bit closer at the kind of initial data that we want to consider in the outmost generality. We recall that we want $`u_0`$ to be increasing and bounded between 0 and 1. Hence, $`v_0=u_0^{}`$ has to be defined and locally integrable, or at most be a measure, in an interval $`(a,b)`$ with possible divergence at the end-points, but anyway with finite integral. Existence of solutions (with $`v_0`$ extended all of $``$ with value 0 otherwise) offers no problem in all the range $`m>1`$ and the solutions are bounded for all $`t>0`$. But large values will be the origin of the new class of solutions to discuss in the following sections.

By integration we obtain a solution of the problem

$$(\mathrm{P}_I)\{\begin{array}{cc}u_t=\mathrm{\Phi }(u_x)_x\hfill & \text{ in }Q=\times (0,\mathrm{})\hfill \\ u(x,0)=u_0(x)\hfill & \text{ for }x\hfill \\ u(x,t)0,\hfill & \end{array}$$

where $`u_0`$ is any nondecreasing continuous real function joining the levels $`u_0=0`$ at $`x=\mathrm{}`$ to $`u_0=1`$ at $`x=\mathrm{}`$. This is what we term as Solutions of Type I. Note that integration of the self-similar profile gives for $`u`$ the form

$$u(x,t)=F(xt^\gamma )$$

where $`F`$ is a primitive of $`G`$ with $`F(\mathrm{})=0`$, $`F(\mathrm{})=M`$ that we want to normalize to $`M=1`$.

###### Theorem 5.1

Let us consider the initial-value for equation (3.1) posed in $`Q=\times (0,\mathrm{})`$ with power function $`\mathrm{\Phi }(s)=s^m/m`$. If $`m>1`$ then for every nondecreasing $`u_0`$ with $`u_0(\mathrm{})=0`$, $`u_0(\mathrm{})=1`$, there exists a unique continuous weak solution in the sense $`u(x,t)0`$ such that $`u(,t)`$ jumps from $`0`$ to $`1`$ as $`x`$ ranges over the line $`x`$. If $`m0`$ the last condition is essential to ensure uniqueness. This class of solutions has bounded gradients for strictly positive times ($`t\tau >0`$), and spreads in space as time advances and the asymptotic decay formula (5.1) holds for $`v=u_x`$.

On the contrary, if $`m1`$ solutions for this initial-value problem with bounded data do not exist.

The results of the power case reflected in this theorem can be generalized to equations with more general $`\mathrm{\Phi }`$. Only the behaviour of $`\mathrm{\Phi }(s)`$ at $`s=0`$ and $`s=\mathrm{}`$ will determine the different behaviour types. In order to tackle the first, we assume that the initial data (and hence the solutions) $`v`$ are bounded. Then it is known that the condition for existence with finite mass is

(5.2)
$$_0^s\mathrm{\Phi }^{}(s)s𝑑s\mathrm{}.$$

The question of large arguments is similar to the power case. These issues will be discussed at length in \[Vbk\]. For a study of these questions in the class of self-similar solutions cf. \[FV\].

### 5.4 Existence of sharp interfaces

There is an interesting subclass of equations (3.1), or (4.1), where sharp interfaces appear. Let us look first at power nonlinearities. Indeed, for exponent $`m>1`$ solutions of the Cauchy problem for equation (3.1) with initial data having compact support will keep this property for all times, while infinite propagation occurs whenever $`m1`$ cf. \[BV, Ka\]. In the first case, given an integrable function $`u_00`$ with integral 1 the solution of the problem can be seen as the a classical solution equation (3.1) in a domain

$$\mathrm{\Omega }=\{(x,t):l(t)<x<r(t)\},$$

with initial conditions $`u(x,0)=u_0(x)`$ and boundary conditions

$$\{\begin{array}{c}u(x,t)=0,\mathrm{\Phi }(u_x)=0\text{ for }x=l(t),\hfill \\ u(x,t)=1,\mathrm{\Phi }(u_x)=0\text{ for }x=r(t).\hfill \end{array}$$

The lines $`x=l(t)`$ and $`x=r(t)`$ are called interfaces or moving boundaries and are completely determined by the over-specified conditions. They are known to be smooth (analytic) functions of $`t`$ \[An1, AV\], and diverge as $`t\mathrm{}`$ like $`O(t^\gamma )`$, thus giving a quantitative estimate of the dispersion effect. We remark in passing that the presence of interfaces means also that the equation is not uniformly parabolic at those points, and consequently the solutions have limited regularity.

On the other hand, when $`m1`$ the same Cauchy problem leads to positive solutions with $`l(t)=\mathrm{}`$ and $`r(t)=\mathrm{}`$. The property of null flux is equivalent to imposing $`u_x=0`$ at $`\pm \mathrm{}`$, a quite natural condition in view of the values $`u=0,1`$ at $`\pm \mathrm{}`$. This condition is automatic for $`m>0`$. However, for $`m0`$ we can have new solutions with decreasing total mass, i.e., such that

$$\frac{d}{dt}v(x,t)𝑑x<0,$$

and they can be even forced to vanish identically in finite time by controlling the outgoing flux at $`x=\pm \mathrm{}`$. We refer to \[RV93, RV95\] for a detailed analysis. In the integrated version they would lead to solutions with a restricted grey range. For the more general class of $`\mathrm{\Phi }`$ mentioned above the condition to have finite interfaces is

(5.3)
$$_0^s\frac{\mathrm{\Phi }^{}(s)}{s}𝑑s\mathrm{}.$$

The fact that this condition is necessary and sufficient can easily seen on the family of traveling waves. The property was first pointed out in \[OKC\].

## 6  Solutions of Type II. Conjugate formulations

We now address the main question of this paper, the construction of solutions with large gradients, appropriate for the contour enhancement problem.

### 6.1  Basic Free Boundary Problem. It is formulated as follows:

Given an increasing function $`u_0(x)`$ defined in an interval $`(a,b)`$ with end values $`u(a+)=0`$, $`u(b)=1`$, to find a continuous function $`u(x,t)`$ and continuous curves $`x=l(t)`$ and $`x=r(t)`$ such that

(i) $`l(0)=a`$, $`r(0)=b`$, and $`l(t)<r(t)`$ for some time interval $`t(0,T)`$.

(ii) $`u`$ solves the following problem in $`\mathrm{\Omega }=\{(x,t):0<t<T,l(t)<x<r(t)\}`$:

$$(\mathrm{P}_{\mathrm{II}})\{\begin{array}{cc}u_t=\mathrm{\Phi }(u_x)_x\hfill & \text{ in }\mathrm{\Omega }\hfill \\ u(x,0)=u_0(x)\hfill & \text{ for }axb\hfill \\ u(l(t),t)=0,u_x(l(t),t)=+\mathrm{}\hfill & \text{ for }0<t<T,\hfill \\ u(r(t),t)=1,u_x(r(t),t)=+\mathrm{}\hfill & \text{ for }0<t<T.\hfill \end{array}$$

Such a triple $`(u,l,r)`$ is called a solution of Type II. The regularity required from $`u`$ as a solution of the problem will depend on the generality of the data. At least $`u`$ will be continuous in the closure of $`\mathrm{\Omega }`$. Furthermore, to avoid unnecessary generality we will ask $`u`$ to be smooth in the interior of $`\mathrm{\Omega }`$. Finally, the requirement of monotonicity is not intrinsic from the mathematical point of view, but it suits the application and allows for the use of the powerful conjugate formulations. Here is the existence result that we are going to prove:

###### Theorem 6.1

Let $`\mathrm{\Phi }`$ be a flux function that is defined, smooth and with $`\mathrm{\Phi }^{}(s)>0`$ for all $`s>0`$. Assume moreover that $`\mathrm{\Phi }(\mathrm{})`$ is finite. Then for every increasing function $`u_0(x)`$ defined in an interval $`[a,b]`$ with $`u(a)=0`$, $`u(b)=1`$ and $`u_0^{}c>0`$ there exists a unique continuous function $`u(x,t)`$ which is defined in a set $`\mathrm{\Omega }`$ as above, is smooth and strictly monotone in $`x`$ for $`0<u<1`$, and there exist continuous curves $`l(t)`$ and $`r(t)`$, such that the triple ($`u,l,r`$) solves problem $`(\mathrm{P}_{\mathrm{II}})`$ in $`\mathrm{\Omega }_T`$. Besides, $`u_xc>0`$ whenever $`0<u<1`$.

We will also show in our construction that the boundary curves are monotone: $`l(t)`$ is nondecreasing, $`r(t)`$ nonincreasing. On the other hand, there is the problem of determining whether $`T`$ if finite or infinite. This depends on $`\mathrm{\Phi }`$ as we will see below.

### 6.2 Conjugate formulation

When dealing with smooth monotone solutions $`u_x>0`$ we can invert the variables $`x`$ and $`u`$ and write $`x=X(u,t)`$. Then $`u_xx_u=1`$, and after some computations we get the partial differential equation satisfied by $`x`$ a a function of $`u`$ and $`t`$:

(6.1)
$$x_t=(\mathrm{\Psi }(x_u))_u,$$

where $`\mathrm{\Psi }`$ is the conjugate flux function (conjugate to $`\mathrm{\Phi }`$), defined for $`s>0`$ as

(6.2)
$$\mathrm{\Psi }(s)=\mathrm{\Phi }(1/s).$$

Differentiation of equation (6.1) with respect to $`u`$ gives rise to the differentiated conjugate equation for $`w=X/u`$ as a function of $`u`$ and $`t`$:

(6.3)
$$w_t=(\mathrm{\Psi }(w))_{uu}.$$

We complete the list of related equations with the direct differentiated equation for $`v=u/x`$, already seen:

(6.4)
$$v_t=(\mathrm{\Phi }(v))_{xx},$$

and then $`v=1/w`$. It is important to point out that these relations are equivalent to the well-known Bäcklund transform \[BK\] between the main variables $`v`$ and $`w`$ of the second and fourth formulations. Indeed, we have

$$v(x,t)=\frac{1}{w(u,t)},u(x,t)=c+_\mathrm{\Gamma }(vdx+\mathrm{\Phi }(v)_xdt).$$

The reader is referred to \[RV95, V03c, V03d\] for other applications of this technique.

### 6.3 Posing and solving the conjugate problems

We now show how to use the conjugate formulations to solve the original problem. We assume that $`\mathrm{\Phi }`$ a flux function defined for all $`s>0`$ and such that $`\mathrm{\Phi }(\mathrm{})`$ is finite, say, $`\mathrm{\Phi }(\mathrm{})=0`$. Then $`\mathrm{\Psi }(s)`$ is continuous at $`s=0`$ and $`\mathrm{\Psi }(0)=0`$. This is the class of flux functions for which the conjugate problem looks simpler. Since we have assumed that $`\mathrm{\Phi }`$ is smooth, so is $`\mathrm{\Psi }`$ in its domain of definition.

(i) Since we are interested in solving the conjugate problem as an auxiliary step for Problem $`(\mathrm{P}_{\mathrm{II}})`$ we will relate the initial values for the function $`w(u,t)`$ to $`u_0`$ as follows. Assuming that $`u_0`$ is continuous and strictly monotone in the interval $`I=\{a<x<b\}`$, with $`u_0(a)=0`$, $`u_0(b)=1`$, and $`C^1`$ smooth inside $`I`$ with $`du_0/dx`$ bounded below away from zero, we define the inverse function $`x=h(u)=u_0^1(u):[0,1][a,b]`$.

(6.5)
$$w_0(u)=\frac{1}{u_{0,x}(h(u))},$$

which is defined for $`0u1`$ and is positive, bounded and smooth inside, i.e., for $`u(0,1)`$.

(ii) We then solve the conjugate problem

$$(\mathrm{P}_\mathrm{c})\{\begin{array}{cc}w_t=\mathrm{\Psi }(w)_{uu}\hfill & \text{ for }0<u<1,t>0\hfill \\ w(u,0)=w_0(u)\hfill & \text{ for }0u1\hfill \\ w(u,t)=0\hfill & \text{ for }u=0,1.\hfill \end{array}$$

As initial data we choose a nonnegative, bounded function $`w_0`$. Under these conditions Problem $`(\mathrm{P}_c)`$ has a unique solution by virtue of well-known nonlinear parabolic theories described in the previous section, cf. \[Be72, BC, BG\]; but note that now we are dealing with the homogeneous Dirichlet problem. The solution can be obtained as limit of the solutions $`w_\epsilon (u,t)\epsilon `$ of the nondegenerate problems with initial data $`w_{0,\epsilon }(u)=w_0(u)+\epsilon `$, $`\epsilon >0`$. In the monotone limit we get

$$\underset{\epsilon 0}{lim}w_\epsilon (u,t)=w(x,t)$$

which is nonnegative, continuous and bounded. Under the additional assumption that $`w_0`$ is locally bounded away from zero, it is easily proved that the solution $`w(u,t)`$ is positive, hence classical, in a strip

$$S_T=\{(u,t):0<u<1,0<t<T\}.$$

(iii) We also need to know something about the large time behaviour of the solutions to this problem. In the full generality ($`\mathrm{\Psi }`$ continuous at $`s=0`$), it can be proved rather easily that $`w(x,t)`$ goes to zero in uniform norm as $`t\mathrm{}`$. However, the rate depends on $`\mathrm{\Psi }`$, as we will see in detail in the next sections. Indeed, depending on $`\mathrm{\Psi }`$ it may happen that the solution vanishes identically after a finite time $`T>0`$ (so-called finite-time extinction). For power nonlinearities $`\mathrm{\Psi }^{}(s)=s^{q1}`$ this happens iff $`0<q<1`$, \[BH\]. On the other hand, it is well-known that for $`q1`$ the decay rate is $`O(t^{1/(q1)})`$, while for $`q=1`$ it is exponential.

(iv) Next, we pass to the integrated version using the formula

(6.6)
$$z(u,t)=_\mathrm{\Gamma }w𝑑u+\mathrm{\Psi }(w)_udt,$$

where $`\mathrm{\Gamma }`$ is any piece-wise smooth curve in $`(u,t)`$ space starting from a fixed point, say $`u=1/2,t=0`$ and arriving at a generic point $`(u,t)`$. In this way we obtain a solution of the integrated equation $`z_t=(\mathrm{\Phi }(z_u))_u`$, much as we did in the case of the original pair of formulations. Note that

(6.7)
$$z(1,t)z(0,t)=w𝑑u0$$

as $`tT`$.

(v) Let us examine the curves $`z_{(u)}(t)=z(u,t)`$ for fixed $`u[0,1]`$. It is clear from the smoothness of the solutions that these curves are $`C^{\mathrm{}}`$ smooth for every $`0<u<1`$. We are interested in the limit curves

(6.8)
$$l(t)=\underset{u0}{lim}z(u,t),r(t)=\underset{u1}{lim}z(u,t).$$

The limit is well-defined for every $`0<t<T`$ by monotonicity of $`z`$ as a function of $`u`$. These curves will show up in the next subsection as the interfaces of the original problem.

###### Lemma 6.2

The curves $`z=l(t)`$ and $`z=r(t)`$ are continuous and monotone for all $`0<t<T`$ with $`al(t)<r(t)b`$ ($`r(t)`$ is nonincreasing and $`l(t)`$ nonincreasing). As $`tT`$ we have $`r(t)l(t)0.`$

Proof. For $`0<s<t<T`$ and $`0<u_1<u_2<1`$ we have

$$z(u_1,t)z(u_1,s)=_{u_1}^{u_2}w(u,t)𝑑u+_{u_1}^{u_2}w(u,s)𝑑u+(z(u_2,t)z(u_2,s)).$$

Fix now $`s>0`$ and let $`t`$ be a bit larger than $`s`$. Since $`w`$ is bounded uniformly for $`ts`$ the two integrals are uniformly small as long as $`u_1u_2`$ is small. We will fix now $`u_20`$ and let $`u_1(0,u_2)`$ go to 0. It is clear then that $`z(u_1,t)z(u_1,s)`$ is uniformly small and goes to 0 as $`ts`$. Hence, $`l(t)`$ is continuous at $`t=s`$. The argument for $`r(t)`$ as $`u1`$ is the same.

For classical solutions the monotonicity of the limit curves is a consequence of the estimate

$$\frac{z}{t}=\frac{}{u}\mathrm{\Phi }(w)$$

which must be positive at $`u=0`$ and negative at $`u=1`$ because $`w0`$ and we have zero boundary conditions. For the general case we use the dependence of the solutions on $`\mathrm{\Phi }`$ as demonstrated in \[BC\]. Hence, $`l^{}(t)0`$, $`r^{}(t)0`$. We will see another proof below. Note that the rest of the curves $`z_{(u)}(t)`$ need not be monotone unless $`w_0`$ is a rearranged symmetrical function.

As for the last statement, it follows immediately from 6.7. It will mean for the original problem that the solution concentrates into a vertical profile, thus showing the formation of the desired front.

### 6.4 Inversion. Solutions of Type II

Thanks to the fact that $`z/u=w>0`$, we can invert the dependence between $`z`$ and $`u`$ in the previous construction to get a function $`u=u(z,t)`$ that is easily shown to satisfy the equation

$$u_t=(\mathrm{\Phi }(u_z))_z.$$

Besides, $`u`$ is a monotone function of $`z`$ and takes the values $`u=0`$ and $`u=1`$ respectively at the left and right endpoints of the domain of definition

$$\mathrm{\Omega }_z=\{l(t)<z<r(t)\},l(t)=z(0,t),r(t)=z(1,t).$$

where $`z(,t)`$ is the function defined in (6.6). Therefore, $`u(z,t)`$ is a candidate to solve our original problem if we identify the independent variable $`z`$ with $`xc`$, where $`c`$ is uniquely determined by the relation $`u_0(c)=1/2.`$

In order to check that we have solved the original problem we still have to check some particulars. It is clear that $`v=u_x`$ is related to the original $`w`$ by the formula

$$v(x,t)=\frac{1}{w(u,t)},$$

which simply states the derivative rule for the inverse function, and $`u`$ in the second member is given by $`u(z,t)`$, $`z=xc`$, as explained before. Here comes an important point: since $`w`$ takes on zero boundary values, $`v(x,t)`$ diverges at the endpoints of its domain of definition, $`\mathrm{\Omega }`$. In other words, the solutions of the original problem $`u=u(x,t)`$ enjoy the property of infinite gradients at the endpoints of the strip where they are defined. Since $`\mathrm{\Phi }(\mathrm{})=0`$ this also means zero flux at these points, a reasonable requirement, which explains why this condition has to be imposed on $`\mathrm{\Phi }`$.

As for the initial data, we have the mass formula

$$x=_{1/2}^uw_0(u)𝑑u+c,\text{with}u_0(c)=\frac{1}{2},$$

so that $`x`$ ranges over an interval $`[a,b]`$ when $`u`$ goes from 0 to 1, i.e., $`a=z(0,0)+c`$, $`b=z(1,0)+c`$. This rule is accompanied by the rule $`u=_{l(t)}^xv(x,t)𝑑x.`$

### 6.5 Uniqueness

Uniqueness of our class of monotone solutions works by translating any couple of solutions of Problem $`(\mathrm{P}_{\mathrm{II}})`$ with the same initial data into the conjugate formulation. They continue to have the same initial data. Uniqueness of weak solutions is well-known for that equation. The Theorem is proved.

### 6.6  Front formation

The asymptotic formation of a vertical front is a simple consequence of the fact that there exists the limit

$$\underset{tT}{lim}r(t)=\underset{tT}{lim}l(t)=x_{\mathrm{}}(a,b).$$

The existence of the common limit follows from Lemma 6.2.

## 7 Uniqueness and comparison for Type II solutions

The previous construction provides for existence and uniqueness of monotone solutions. We tackle next the property of comparison, which is stronger that uniqueness.

###### Theorem 7.1

Comparison applies to the solutions of Theorem 6.1 : let $`u_1(x,t)`$, $`i=1,2`$, be two solutions defined in a strip $`S_T=\times (0,T)`$ having initial data $`u_{0i}`$, where the solutions and the data have been extended by 0 the the left of the definition domain, by 1 to the right. If $`u_{01}(x)u_{02}(x)`$, then for all $`x`$ we have

(7.1)
$$u_1(x,t)u_2(x,t)\text{in }S_T.$$

Proof. The main idea is to use a contradiction argument at touching points of ordered solutions. We proceed in several stages and need a definition, taken from the theory of viscosity solutions in \[CV\]. We say that two solutions are strictly separated at a time $`t`$ if $`u_1>u_2`$ at all points intermediate to the bounding interfaces, and these are also separated.

(1) Given two solutions $`u_1`$ and $`u_2`$ with ordered data we must obtain approximations to which the contradiction argument can be applied. First, we separate the initial data of the solutions by displacing the first solution $`u_1`$ to the left (the problem is invariant under space displacements). Next, assume that we have modified $`u_1`$ so that it has a finite slope on the left-hand free boundary, and $`u_2`$ so that it has a finite slope on the hand-side free boundary. In this situation we argue on the first point and time $`(x_0,t_0)`$ where the graphs of both solutions touch and discover that: (i) it cannot happen with value $`u(0,1)`$ by the Strong Maximum Principle; (ii) not on the left-hand side where $`u_1=u_2=0`$, because at such a point we have the right derivative: $`u_{2,x}=+\mathrm{}`$ while $`u_{1,x}`$ is finite, which contradicts the fact that we still must have $`u_1(,t_0)u_2(,t_0)`$; (iii) same argument applies to the right-hand side where $`u_1=u_2=1`$. We conclude that such solutions $`u_1`$ and $`u_2`$ cannot touch, hence strict separation is preserved in time, and comparison $`u_1u_2`$ holds in the strong sense, i.e., with strict separation at all times.

(2) We must now prove that the original solutions can be approximated by solutions with finite derivatives on the lateral boundaries. This can be obtained by an easy modification of the previous construction. Thus, if $`u`$ is the constructed solution, it can be approximated at the level of the conjugate problem by putting the value $`w=\epsilon `$ on the left-hand side boundary $`u=0`$. It is rather standard monotonicity argument that $`w_\epsilon w`$ in the limit $`\epsilon 0.`$ Undoing the transformation this means that $`u_\epsilon `$ converges to the original constructed solution $`u`$. A similar construction applies to $`u_2`$ by modifying $`w_2`$ at the border $`u=1`$. After displacement of $`u_\epsilon `$ we get $`u_{1,\epsilon }(x+\delta ,t)u_{2,\epsilon }(x,t)`$. Let now $`\epsilon 0`$.

We derive next another proof of the monotonicity of the interfaces.

###### Corollary 7.2

The interfaces are monotone in all cases.

Proof. Given a solution $`u`$ with initial interfaces $`a`$ and $`b`$, we may place a very steep solution $`\widehat{u}`$ to the left, i.e., with $`\widehat{u}=0`$ for $`xa2\epsilon `$, $`\widehat{u}=1`$ for $`xa\epsilon `$. If $`\widehat{u}^{}`$ is symmetric then the interfaces are monotone, hence the left interface of $`\widehat{u}`$ lies to the right of $`a2\epsilon `$. Then $`\widehat{u}u`$ for all times, hence $`l(t)a2\epsilon `$, and in the limit $`l(t)a`$. The same argument applies taking origin of times at any $`y(0,T)`$, hence $`l(t)`$ is monotonically increasing in $`(0,T)`$. Likewise, $`r(t)`$ is proved to be monotone decreasing.

## 8 Self-similarity and asymptotics for power nonlinearities

In the power case, where $`\mathrm{\Phi }(s)=(1/m)s^m`$, the asymptotic condition “$`\mathrm{\Phi }(\mathrm{})`$ finite” means $`m<0`$. Then $`\mathrm{\Psi }(s)=(1/q)s^q`$ with $`q=m>0`$. Reminder: we consider only monotone solutions of (3.1) and nonnegative solutions of (4.1), and $`q=1+2\alpha `$ in the notation of the first section. Therefore, $`\alpha >1/2`$ means $`q>0`$.

There is a unique solution of problem $`(\mathrm{P}_c)`$ with initial data as in Section 6. The large time behaviour depends on the exponent. For $`q1`$ solutions are defined and nontrivial globally in time, while for $`0<q<1`$ they exist only in a finite interval $`0<t<T`$ and $`w(x,t)0`$ uniformly as $`tT`$.

As for the exact behaviour, if $`\alpha 0`$ the self-similar behaviour reduces to writing $`U=\mathrm{\Phi }(\xi )`$ with $`\xi =xt^{1/2\alpha }`$ and integrating the ODE

(8.1)
$$\mathrm{\Phi }^{\prime \prime }=\frac{1}{2\alpha }\xi (\mathrm{\Phi }^{})^{2+q}$$

to obtain in the monotone region

(8.2)
$$\mathrm{\Phi }^{}(\xi )=A(K^2\xi ^2)^{1/(2+2\alpha )},A=\left(2\alpha /(1+\alpha )\right)^{1/(2+2\alpha )}.$$

and a constant $`K>0`$ determined by the end values $`\mathrm{min}\mathrm{\Phi }(\xi )=0`$, $`\mathrm{max}\mathrm{\Phi }(\xi )=1`$. For the large-time behaviour of general solutions we have to make a study in different cases. We start with $`\alpha =0`$ which was not included before.

Case $`\alpha =0`$, hence $`q=1`$. The conjugate equation is the linear heat equation. It is well-known that the solutions behave for large times like

(8.3)
$$w(u,t)C_1W(u,t)C_1e^{\lambda t}f_1(u),$$

where $`\lambda =\pi ^2`$ is the first eigenvalue of the Laplacian in $`[0,1]`$, $`f_1(x)=\mathrm{sin}(\pi u)`$ is the first eigenfunction, and $`C_1>0`$ is a constant that depends on the initial data. Here and in the sequel we denote by capital letters the quantities corresponding to self-similar solutions. We get the intermediate step

(8.4)
$$xc=z(u,t)Z(u,t)Ce^{\lambda t}\mathrm{cos}(\pi u),$$

with $`C=C_1/\pi `$. This formula can be viewed as an implicit expression of $`u=u(x,t)`$ defined in the space between the interfaces. Note that $`u=1/2`$ for $`x=c`$ and all $`t>0`$. Putting $`c=0`$ we get for the self-similar solution

(8.5)
$$U(x,t)=F(xe^{\pi ^2t}/C),$$

where $`F^{}`$ behaves like $`d^{1/2}`$ near the end-values or interfaces, $`d`$ being the distance to these points. Moreover, the interfaces

(8.6)
$$R(t)=CR_1(t),L(t)=CR_1(t),\text{with}R_1(t)=e^{\pi ^2t}.$$

These estimates become the first order approximation when we consider general solutions:

(8.7)
$$u(x,t)F(xe^{\pi ^2t}/C),r(t),l(t)c\pm CR_1(t).$$

The gradients are given by

(8.8)
$$u_x(x,t)=\frac{1}{w}\frac{1}{C}e^{\lambda t}\frac{1}{\mathrm{sin}(\pi u)}.$$

They blow up exponentially as $`t\mathrm{}`$ and like $`O(d^{1/2})`$ at the interfaces, where $`d`$ is distance to the interface.

Case $`\alpha >0`$, hence $`q>1`$. This is the case treated in \[B01\] by the direct self-similar method. We recover the behaviour of general solutions from the conjugate problem. We have a self-similar asymptotic expression in the form of separation of variables:

(8.9)
$$w(u,t)W(u,t)t^{1/(q1)}f_q(u),$$

where $`f_q(u)>0`$ is the unique solution of the associated elliptic problem

$$(f^q)^{\prime \prime }+\mu f=0\text{in }[0,1],\mu =\frac{q}{q1},$$

and $`f=0`$ at the endpoints. This problem was studied in \[AP\]. In this case there is no free constant, i.e., the behaviour is universal. From it we get the intermediate asymptotic estimates

(8.10)
$$xc=z(u,t)t^{1/(q1)}g_q(u).$$

This formula defines implicitly $`u=u(x,t)`$ in the space between the interfaces as

(8.11)
$$uF_q((xc)t^{1/(q1)}).$$

Since $`f(u)`$ behaves exactly like $`O(u^{1/q})`$ near the end point $`u=0`$, and $`O((1u)^{1/q})`$ near $`u=1`$, it follows after integration and inversion that $`u=F`$ behaves like $`O(d^{(q/(q+1)})`$, where $`d`$ is the distance to the interfaces. These are given by

(8.12)
$$l(t)=cC_qt^{1/(q1)},r(t)=c+C_qt^{1/(q1)}.$$

The gradient is given by

(8.13)
$$u_x(x,t)t^{1/(q1)}\frac{1}{f_q(u)},$$

and blows up like $`O(t^{1/2\alpha })`$ as $`t\mathrm{}`$ and like $`O(d^\gamma )`$ at the interfaces, where $`\gamma =1/(q+1)=1/(2\alpha +2)`$.

Case $`1/2<\alpha <0`$, hence $`0<q<1`$. This is a new case, not included in the modelization of \[MS, B01\]. We get solutions from the conjugate problem which exist for a finite time $`T`$ and behave as $`tT`$ like the separation of variables formula:

(8.14)
$$w(u,t)(Tt)^{1/(1q)}f_q(u),$$

where $`f_q(x)>0`$ is the unique solution of the associated elliptic problem

$$(f^q)^{\prime \prime }+\mu f=0\text{in }[0,1],\mu =\frac{q}{q1},$$

and $`f=0`$ at the endpoints, with no free constant. We get by integration

(8.15)
$$xc=z(u,t)(Tt)^{1/(1q)}g_q(u).$$

This formula defines implicitly $`u=u(x,t)`$ in the space between the converging interfaces

(8.16)
$$l(t)=cC(Tt)^{1/(1q)},r(t)=c+C(Tt)^{1/(1q)}.$$

The gradient is given by

(8.17)
$$u_x(x,t)(Tt)^{1/(1q)}\frac{1}{f_q(u)},$$

and blows up like $`O((Tt)^{1/2\alpha })`$ as $`tT`$ and like $`O(d^\gamma )`$ at the interfaces, $`\gamma `$ as before.

## 9 Asymptotic estimates for more general $`\mathrm{\Phi }`$

Assumptions on $`\mathrm{\Phi }`$ to behave like a power at infinity guarantee that the asymptotic behaviour is as predicted in the power case. Thus, Bertsch and Peletier study in \[BP\] the asymptotic behaviour of equation (6.1) written in the form

$$\beta (v)_t=v_{xx}.$$

In our application the space variable $`x`$ becomes $`u`$ and $`\beta =\mathrm{\Psi }^1`$. They assume that $`\beta (0)=0`$, $`\beta ^{}(s)>0`$ for $`s>0`$ and

$$1\frac{s\beta ^{\prime \prime }(s)}{\beta ^{}(s)}\alpha $$

for all small $`0<s<s_0`$ and some $`\alpha (0,1)`$. They also need the more stringent condition that the limit

$$\sigma (s)=\underset{\epsilon 0}{lim}\frac{s\beta ^{\prime \prime }(\epsilon s)}{\beta ^{}(\epsilon )}$$

exists for every $`s>0`$. This limit is necessarily of the form $`\sigma (s)=s^m`$ for some $`m[0,1)`$. Those conditions cover in particular powers $`\mathrm{\Psi }(s)=s^q`$ with $`q=1/m>1`$, as well as the exponential $`\mathrm{\Psi }(s)=e^{1/s}`$, and any other function resembling such examples near $`s=0`$.

Under these restrictions, they establish an asymptotic separation-of-variables result that, when translated to our setting, means that the solution $`w(u,t)`$ of the Cauchy problem (6.1) behaves as $`t\mathrm{}`$ in the separate variable form $`wy(t)f(u)`$, more precisely,

(9.1)
$$\underset{t\mathrm{}}{lim}\frac{\beta (w)}{y(t)}=f(u),$$

where $`y(t)`$ is the solution $`y^{}=y^{1/q}`$, i.e. $`y(t)ct^{q/(q1)}`$ and $`f`$ is the profile obtained as in the previous section. This means that the analysis of the previous section is justified for almost power-like flux functions if they resemble a power with exponent $`q>1`$. The cases $`q=1`$ and $`q<1`$ should be justified in a similar manner. We note that a similar calculation for the Cauchy problem was done in \[PV\].

We point out that these restrictions on the nonlinearity are satisfied by the equation proposed in the image processing model, (1.4), hence the results apply to that model.

## 10  Numerical experiments

We have computed the solutions of the free boundary problem ($`\mathrm{P}_{\mathrm{II}})`$) by solving the conjugate problem, which is a homogeneous Neumann problem, in the cases $`\alpha =0`$, $`\alpha =1`$ and $`\alpha =1/2`$. The conjugate equations are $`p`$-Laplacian equations with $`p=2,4`$ and $`3/2`$. The computations have been done with symmetrical data $`x_0(u)=\mathrm{cos}(\pi u)`$ and asymmetrical $`x_0(u)=2u^4(1/2)u^61`$. The last two examples stabilize in finite time (since $`\alpha <0`$).

Figures 1 and 2. The case $`\alpha =0`$ with symmetrical and asymmetrical data

Figures 3 and 4. The case $`\alpha =1`$ with symmetrical and asymmetrical data

Figures 5 and 6. The case $`\alpha =1/2`$ with symmetrical and asymmetrical data

## 11  Discussion and conclusions

We have established the well-posedness of a free boundary problem that represents a one-dimensional version of the model for image contour enhancement. We have also established the asymptotic behaviour and its rates. The results apply to a large class of equations, which enlarges the scope of the results of \[B01, MS\]. We have performed a classification of the solutions and their properties according to the properties of the constitutive function $`\mathrm{\Phi }`$.

We have used as a framework the theory nonlinear parabolic equations of diffusive type, or nonlinear diffusion equations for short. These equations are typically used in describing processes of mass diffusion or thermal propagation (with or without additional effects, like convection or reaction). Here they appear in image processing. The types used here cover the model cases known as $`p`$-Laplacian equations and porous medium equations or their variants.

The analysis is technically performed by means of a series of remarkable transformations that lead to conjugate problems that are easier to analyze. These transformations are related to the Bäcklund transform.

There are a number of interesting consequences of the mathematical analysis that we sum up next.

$``$ Generality. The first one is the observation that solutions with gradient blowup can be obtained for the nonlinear diffusion equation with a large class of constitutive functions. Namely, $`\mathrm{\Phi }`$ can be any continuous nondecreasing function defined for $`s>0`$ such that $`\mathrm{\Phi }(\mathrm{})`$ is bounded. We have made for convenience the assumption of smoothness and strict monotonicity for $`0<s<\mathrm{}`$. Actually, the constitutive function can be more general, even a discontinuous maximal monotone graph as in \[Be72, Brz\]. But such a generality is not practical here, though the classification and results go through after some heavy work. A general presentation of that generality will be the object of the work \[Vbk\].

$``$ Focusing in finite time. In the case of powers $`\mathrm{\Phi }(s)s^m`$ the condition means $`m<0`$, and not only $`m1`$ as considered before. Moreover, the new range $`0>m>1`$ leads to a very fast evolution which arrives at a vertical front in finite time. This idea might have an applied interest.

$``$ Coexistence of types. We have also a so-called class of solutions of Type I that fall into the scope of the standard parabolic theory. They exist under the condition that the integral $`(\mathrm{\Phi }^{}(s)/s)𝑑s`$ converges at $`s=\mathrm{}`$, which in the case of powers means $`m>1`$. This fact leads to an interesting observation: there is a range of co-existence of solutions of Type I and II for $`1<m<0`$, i.e., $`1/2<\alpha <0`$. This may be of interest from the theoretical point of view: two different problems share the same equation and same initial data, but the solutions differ as a consequence of a further choice: of existence or not existence of a free boundary.

$``$ The pressure. There is a way of formally unifying the standard theory of free-boundary solutions for the porous medium equation with the present theory of free-boundaries in fast diffusion. This works by means of the new variable called pressure that for an equation of the form $`v_t=\mathrm{\Phi }(v)_{xx}`$ is defined as

(11.2)
$$p(v)=_a^v\frac{\mathrm{\Phi }^{}(s)}{s}𝑑s.$$

If the integral is convergent at $`v=0`$ then the choice $`a=0`$ is made. Otherwise, any value in the domain of $`\mathrm{\Phi }`$ is good. Indeed, the similarity in the behaviour of the pressure of Porous Medium case, $`\mathrm{\Phi }(u)=u^m`$, $`m>1`$, and our case (which lies in the power range $`m0`$) is seen when we write the partial differential equation for $`v`$ which reads

(11.3)
$$v_t=\sigma (v)\mathrm{\Delta }v+|v|^2,$$

with $`\sigma (v)=\mathrm{\Phi }^{}(u)`$. For $`\mathrm{\Phi }(u)=cu^m`$ we get $`v=cmu^{m1}/(m1)`$ and $`\sigma (v)=(m1)v`$. The paper \[V03c\] explores further the properties of free-boundary solutions for the pressure of fast diffusion equation. Note that, contrary to the porous medium case where the support of the pressure solutions expands in time, in our case it shrinks in time.

$``$ Grey levels as end values. Another variation of the main theme is to consider solutions losing the white or black level. This is investigated in the paper \[RV95\], where many solutions of the differentiated equation for $`v`$ are obtained in the range $`0m>1`$ with the same initial data by assigning fluxes at infinity:

(11.4)
$$\underset{x\mathrm{}}{lim}(v^{m1}v_x)=f(t),\underset{x\mathrm{}}{lim}(v^{m1}v_x)=g(t),$$

for bounded functions $`f,g0`$. In terms of $`u`$ it means in particular that we can obtain a free-boundary solution with a decreasing value of highest color level,

$$u(\mathrm{},t)=F(t)=1_0^tf(t)𝑑t<1.$$

Analogously, we can impose an increasing value for the lowest color level,

$$u(\mathrm{},t)=G(t)=_0^tg(t)𝑑t>0.$$

$``$ Non-monotone fronts. There is a gap in the theory we have developed, namely that we assumed that the front is monotone in the space direction, in other words, that $`u`$ is nondecreasing as a function of $`x`$. The limitations of fast diffusion equations to admit non-monotone solutions are not accidental. In the paper \[RV02\] it is proved that the model equation $`u_t=u_x^{m1}u_{xx}`$ does not admit non-monotone solutions if $`m0`$. In terms of the equation for $`v=u_x`$ it means that there are no solutions with both signs. However, there is hope when using functions $`\mathrm{\Phi }`$ like in (1.4) that are degenerate at infinity, i.e., $`\mathrm{\Phi }^{}(\mathrm{})=0`$, but regular at all other values, e.g., $`\mathrm{\Phi }C^1()`$. Indeed, the problem with the monotonicity happens because of the singularity at $`u_x=0`$ which forbids maxima or minima, while our main interest is in the free-boundary that is governed by the large values, $`u_x\mathrm{}`$. This question will be investigated separately.

$``$ Non-monotone nonlinearities. Note that in our problem setting the equation is forward parabolic, and the backward movement of the interface is due to the effect of the singular boundary condition, which happens to be compatible with the equation for the appropriate class of functions $`\mathrm{\Phi }`$ (e.g., for powers $`m<0`$). In the spirit of the Perona-Malik model, situations can be considered where $`\mathrm{\Phi }`$ is not monotone, hence the equation is not parabolic. A regularization is then needed so that one faces a regularized forward-backward diffusion problem. See in this respect \[BBDU\] and \[Wit\]. This subject is again under study.

Acknowledgment

The work of the first author was supported by Applied Mathematics subprogram of the US Department of Energy under Contract DE-AC03-76-SF00098. The second author has been partly supported by TMR Projects RTN2-2001-00349 and FMRX-CT98-0201, and Proyecto MCYT BFM2002-04572-C02-02. He is also grateful to the Lawrence Berkeley National Laboratory for supporting a research visit. The authors express their thanks to Prof. A. J. Chorin for his stimulating interest in our work. Computations were performed by Raúl Ferreira using Matlab.

Permanent addresses:

Grigori I. Barenblatt

Department of Mathematics and Lawrence National Laboratory

Berkeley, California 94720-3840

e-mail: gib@math.uberkeley.edu

Juan Luis Vázquez

Departamento de Matemáticas, Universidad Autónoma de Madrid

28046 Madrid, Spain.

e-mail: juanluis.vazquez@uam.es


text:
# Magnetic Phase Diagram of GdNi2B2C: Two-ion Magnetoelasticity and Anisotropic Exchange Couplings.

## I Introduction

The field-temperature ($`HT`$) magnetic phase diagrams of the intermetallic borocarbides $`R`$Ni2B2C ($`R`$=magnetic rare earth) manifest wide varieties of zero-field magnetic structures such as Neel-type, helimagnetic, amplitude-modulated, and squared-up states.1-8 Most of these states, in particular the modulated ones, are unstable against field and temperature variation, leading to a cascade of transformations Er-HT-diagram ; Ho-HT-Diagram ; Dy-HT-diagram  similar to the ones observed in the elemental rare-earths. It is the general opinion that the remarkable features observed in these $`HT`$ phase diagrams are governed by the combined influence of exchange couplings and crystalline electric field (CEF) interactions Theory-exchange-CEF-Phase-diagram : such an approach reproduced successfully the gross features of the magnetic $`HT`$ phase diagrams of HoNi2B2C. Nevertheless, various experimental observations suggest that two additional interactions, namely, anisotropic exchange interaction and magnetoelastic coupling, are necessary ingredients for the understanding of the behavior of the magnetic borocarbides. Magnetoelastic interactions are spectacularly manifested in the onset of a tetragonal-to-orthorhombic distortion at $`T_N`$ for $`R`$= Er, Ho, Dy, Tb compounds Er-magnetostriction ; Ho-distortion ; Dy-magnetostriction ; Tb-dichroism-HRMXRD , while the anisotropic exchange interactions were reported to be necessary for the description of the low-temperature magnetic properties of ErNi2B2C (Ref.Er-magnetostriction ) and GdNi2B2C (Ref.Gd-single-crystal-Tdynamics ). Obviously, the anisotropic exchange interactions and (two-ion) magnetoelastic couplings can be conveniently investigated in GdNi2B2C wherein CEF anisotropy is negligible and the de Gennes factor is the strongest.

This work reports a study of the magnetic and magnetoelastic properties of single crystal GdNi2B2C. The $`HT`$ magnetic phase diagrams for fields along the a- and c-axis were determined. Moreover, within the paramagnetic phase, the parastriction and paramagnetic measurements revealed the presence of weakly anisotropic exchange interaction as well as two-ion magnetoelastic couplings while within the ordered state these perturbations were observed to be strong enough to induce a noticeable modification in the nuclear and magnetic structures.

GdNi2B2C crystallizes in a body-centered tetragonal structure (space group I4/mmm, point symmetry of the Gd site is 4/mmm).R1221-structure  The zero-field magnetic state immediately below T$`{}_{N}{}^{}=19.5`$ K is a transverse, incommensurate, sine-modulated structure (moments along b-axis) with a wave vector $`\stackrel{}{k_a}`$ that decreases linearly from 0.551$`a^{}`$ at $`T_N`$ to 0.550$`a^{}`$ at $`T_R`$13.5 K. Gd-XRES ; Gd-single-crystal-Tdynamics ; Gd-poly-high-field ; Gd-chemical-composition ; Gd-MES  At $`T_R`$, a moment reorientation sets-in leading to an additional modulated component transversely polarized along the c-axis.Gd-XRES ; Gd-MES  On further temperature decrease, $`\stackrel{}{k_a}`$ increases monotonically reaching 0.553$`a^{}`$ at 3.5 K.

## II Theoretical Backgrounds

The $`{}_{}{}^{7/2}S`$-character of GdNi2B2C leads to a greater simplification of its magnetic Hamiltonian which is considered to consist mainly of an exchange, a Zeeman, a two-ion magnetoelastic, and an elastic term. Within the mean-field approximation (MFA) the tetragonal-invariant bilinear exchange interaction can be written as a sum of two terms:Gignoux-Schmitt-1990

$`_{ex}^{\alpha _1}`$ $`=(g__J\mu _B)^2n^{\alpha _1}<J>J,`$
$`_{ex}^{\alpha _2}`$ $`=(g__J\mu _B)^2n^{\alpha _2}\left(2<J_z>J_z<J_x>J_x<J_y>J_y\right)`$ (1)
$`_{ex}^{\alpha _1}`$ corresponds to the isotropic bilinear exchange coupling while $`_{ex}^{\alpha _2}`$ corresponds to an anisotropic bilinear coupling. The values of the isotropic ($`n^{\alpha _1}`$) and anisotropic ($`n^{\alpha _2}`$) effective exchange coefficients can be determined from the paramagnetic susceptibilities, measured with the field along ($`\chi _c`$) and perpendicular ($`\chi _c`$) to the c-axis:

$`\chi _c`$ $`={\displaystyle \frac{C}{T(n^{\alpha _1}n^{\alpha _2})C}}`$
$`\chi _c`$ $`={\displaystyle \frac{C}{T(n^{\alpha _1}+2n^{\alpha _2})C}}`$ (2)
$`C`$ is the Curie constant. In Gd-based compounds, the two-ion magnetoelasticity is related to the modification of the magnetic interactions by the strains. For strains that preserve the initial tetragonal symmetry, the two-ion magnetoelastic Hamiltonians can be written as:Gignoux-Schmitt-1990

$`_{ME}^{\alpha _1}`$ $`=\left(D_{\alpha _1}^{\alpha _1}ϵ^{\alpha _1}+D_{\alpha _1}^{\alpha _2}ϵ^{\alpha _2}\right)<J>J,`$
$`_{ME}^{\alpha _2}`$ $`=\left(D_{\alpha _2}^{\alpha _1}ϵ^{\alpha _1}+D_{\alpha _2}^{\alpha _2}ϵ^{\alpha _2}\right)\left(2<J_z>J_z<J_x>J_x<J_y>J_y\right)`$ (3)
where $`D_{\alpha j}^{\alpha _i}`$ are the two-ion magnetoelastic constants and $`ϵ^{\alpha _i}`$ are the normalized, symmetrized strains. $`ϵ^{\alpha _1}`$ is the volume strain while $`ϵ^{\alpha 2}`$ is the axial strain acting on the c/a ratio:

$`ϵ^{\alpha _1}`$ $`={\displaystyle \frac{1}{\sqrt{3}}}(ϵ_{xx}+ϵ_{yy}+ϵ_{zz})`$
$`ϵ^{\alpha 2}`$ $`={\displaystyle \frac{1}{\sqrt{6}}}(2ϵ_{zz}ϵ_{xx}ϵ_{yy})`$ (4)
The minimization of the free energy with respect to each strain yields the equilibrium values of, say, the $`\alpha `$-strains:

$`ϵ^{\alpha _1}`$ $`=M_{\alpha _1}^{\alpha _1}<J>^2+M_{\alpha _2}^{\alpha _1}\left(2<J_z>^2<J_x>^2<J_y>^2\right)`$
$`ϵ^{\alpha 2}`$ $`=M_{\alpha _1}^{\alpha _2}<J>^2+M_{\alpha _2}^{\alpha _2}\left(2<J_z>^2<J_x>^2<J_y>^2\right)`$ (5)
where $`M_{\alpha j}^{\alpha _i}`$ ($`i,j`$=1,2) is a combinations of the $`D_{\alpha j}^{\alpha _i}`$ and the symmetrized elastic constant $`C_0^{\alpha _k}`$($`k`$=1,2,12). The magnetoelastic couplings can be experimentally probed by measuring the macroscopic length change under an applied magnetic field. In a tetragonal symmetry, the relative length change, $`\lambda =\frac{l}{l}`$, measured along the ($`\beta _1,\beta _2,\beta _3`$) direction when an external field is along the (x,y,z) direction is given by:

$${}_{}{}^{\beta _1\beta _2\beta _3}\lambda _{xyz}^{}=\frac{1}{\sqrt{3}}ϵ^{\alpha _1}+\frac{1}{\sqrt{6}}ϵ^{\alpha _2}(2\beta _3^2\beta _1^2\beta _2^2)+\frac{1}{\sqrt{2}}ϵ^\gamma (\beta _1^2\beta _2^2)$$
(6)
where $`ϵ^\gamma `$ corresponds to the strain mode that leads to a symmetry lowering from tetragonal to orthorhombic. This $`ϵ^\gamma `$ mode can be determined from the appropriate magnetostriction measurements, namely:

$$ϵ^\gamma =\frac{1}{\sqrt{2}}(^{100}\lambda _{100}^{010}\lambda _{100})$$
(7)
The thermal and field variation of the measurable magnetostrictions arise from the $`T`$\- and $`H`$-dependence of $`<J_i>`$. The calculation of such a dependence is usually carried out either by proper diagonalization of the total Hamiltonian (mostly appropriate for the ordered state and higher applied fields) or by perturbation methods (appropriate within the paramagnetic phase and for small magnetic fields). Here we will be concerned only with the paramagnetic case (for more details see Ref.Gignoux-Schmitt-1990 ). Utilizing the susceptibility formalism Gignoux-Schmitt-1990 , one gets the following useful relations:

$`ϵ^{\alpha _i}`$ $`=\chi _{100}^{\alpha i}H^2`$
$`ϵ^{\alpha _i}`$ $`=\chi _{001}^{\alpha _i}H^2`$ (8)
where $`\chi _j^{\alpha _i}`$ ($`i`$=1,2) is the strain field susceptibility along the indicated direction. A complete description of the parastriction can be obtained by combining these equations and Eqs.5 together with the relation:

$$<J_i>=\frac{\chi _i}{N_Ag\mu _B}H;i=x,y,z;\text{ }N_A\text{ is the Avogadro number}$$
(9)
## III Experimental

Single crystals were grown by floating zone method.FZ-method  X-ray and magnetic measurements revealed a single phase of an excellent quality. The size and orientation of the rectangular spark-cut bars were chosen to suit the requirement of the magnetization or the magnetostriction measurements. The isofield and isothermal magnetization curves were measured by the extraction method for temperatures down to 1.5 K and fields up to 160 kOe. The sensitivity of the measurement is better than 5x10-5 emu.

The magnetostriction measurement were performed using a high-accuracy capacitance dilatometer which allows for measuring thermal expansion and forced magnetostriction in the temperature range of 3 - 250 K and under magnetic field up to 65 kOe. The typical resolution is better than 1 $`\stackrel{}{A}`$. In addition, the set-up allows for the rotation of the capacitance cell around the vertical axis of the cryostat while the field is maintained in a horizontal direction. In this way, the angle between the sample axis and the field direction can be varied across the range $`0180^{}`$. The relative length variation for temperature scan at fixed fields is defined as $`\lambda _H=\left[l(T)l(T_0)\right]/l(T_0)`$ while that for field scan at fixed temperatures is defined as $`\lambda _T=\left[l(H)l(H=0)\right]/l(H=0)`$. Diamagnetic contributions and demagnetizing fields are expected to be small and as such no correction for their influences on the magnetization or magnetostriction curves were attempted.

## IV Results and Analysis

Below we show representative magnetization and magnetostriction curves taken within the paramagnetic and the ordered phase. By applying the field along the (001) or (100) direction, it was possible to investigate separately the volume distortion $`ϵ^{\alpha 1}`$, the axial distortion $`ϵ^{\alpha 2}`$ or the orthorhombic distortion $`ϵ^\gamma `$. The paramagnetic and parastriction curves were analyzed according to susceptibility formalism given in Sec.II. Based on the overall analysis of the ordered state, $`HT`$ phase diagrams for $`H//\left(100\right)`$ and $`H//\left(001\right)`$ are constructed.

### IV.1  Paramagnetic Phase

#### IV.1.1 Magnetization

The paramagnetic susceptibilities were deduced at each temperature (25 K$`<T<`$300 K) from the Arrott plot ($`M^2`$ versus $`H/M`$). The thermal evolution of these susceptibilities along each of the (100) and (001) direction (shown in Fig.1) reveals a weak anisotropic character, similar to the one reported by Canfield et al Gd-single-crystal-Tdynamics . From the analysis of Fig.1 and according to Eq.2, we obtained $`n^{\alpha _1}`$=0.12($`\pm `$0.04) mole/emu and $`n^{\alpha _2}`$=-0.05($`\pm `$0.04) mole/emu. The experimental error involved in the determination of these $`n^{\alpha _1}`$ and $`n^{\alpha _2}`$ were dictated by the experimental resolution, nonetheless, the determination of their signs and their magnitudes are consistent with the results of Canfield et al Gd-single-crystal-Tdynamics : on the one hand, $`n^{\alpha _1}`$ is positive and weak indicative of a small effective $`\theta `$. On the other hand, the magnitude of $`n^{\alpha _2}`$ is twice smaller (indicative of a weak paramagnetic anisotropy) and its sign is negative (i.e. $`\chi _c<\chi _c`$) implying that it is energetically favorable for the spins to lie within the basal plane.

The effective moments deduced from the slope of the inverse susceptibilities were found to be $`\mu _{eff\left[100\right]}`$=8.33(4) $`\mu _B`$ and $`\mu _{eff\left(001\right)}`$=8.12(4)$`\mu _B`$; both are slightly higher than the theoretical value ($`\mu _{eff}`$=7.94$`\mu _B`$). Such a small discrepancy arises possibly due to matrix contribution: an exchange-enhancement of the Ni-sublattice susceptibility and/or a polarization of the conduction band.

#### IV.1.2 Parastriction

Representative curves of the parastrictions $`{}_{}{}^{100}\lambda _{100}^{}`$, $`{}_{}{}^{010}\lambda _{100}^{}`$, and $`{}_{}{}^{001}\lambda _{100}^{}`$ are shown in Fig.2. Similar curves (not shown) were obtained as well for $`{}_{}{}^{i}\lambda _{010}^{}`$and $`{}_{}{}^{i}\lambda _{001}^{}(i=(100),(010),(001))`$. In all cases, $`\lambda `$ is of the order of few 10-5, explaining the earlier report that magnetoelastic effects were apparently absent in the temperature-dependent XRD measurements.Gd-XRD-magnetoelastic

Fig.3 shows the symmetrized, normalized strains $`ϵ^{\alpha 1}`$ and $`ϵ^{\alpha 2}`$ which were constructed from parastrictions curves (as those in Fig.2) according to Eq.6. Both $`ϵ^{\alpha 1}`$ and $`ϵ^{\alpha 2}`$ are nonzero and follow faithfully the quadratic field dependence (see Eq.8). These results coupled with the tetragonal character of GdNi2B2C predict a magnetic contribution to the thermal expansion, in addition to the conventional contribution expected in the diamagnetic isomorph $`R`$Ni2B2C ($`R`$=Lu,Y, La). The symmetry lowering mode $`ϵ^\gamma `$ (not shown) is mostly field-independent ruling out any possibility for a field-induced orthorhombic distortion within the paramagnetic phase.

Fig.4 shows the thermal evolution of $`\chi _{100}^{\alpha _1}(T)`$ and $`\chi _{001}^{\alpha _1}(T).`$ These curves were obtained from the analysis of Fig.3 according to Eq.8. Following the same procedure, we obtained the $`\chi _{100}^{\alpha _2}(T)`$ and $`\chi _{001}^{\alpha _2}(T)`$ curves (shown in Fig.5). Within this paramagnetic range and according to Eqs.6 we get:-

$`ϵ_{H//100}^{\alpha _1}`$ $`=ϵ_{H//001}^{\alpha _1}`$
$`ϵ_{H//100}^{\alpha _2}`$ $`=ϵ_{H//001}^{\alpha _2}`$ (10)
then:

$`\chi _{100}^{\alpha _1}`$ $`=\chi _{001}^{\alpha _1}`$
$`\chi _{100}^{\alpha _2}`$ $`=\chi _{001}^{\alpha _2}`$ (11)
Within the experimental errors, these equalities were experimentally confirmed in Figs.4-5: the analysis of the $`\chi ^{\alpha 1}`$ curves gave:

$`1/\sqrt{\chi _{100}^{\alpha _1}(T)}`$ $`=361(30)+41(1)\text{x}T`$
$`1/\sqrt{\chi _{001}^{\alpha _1}(T)}`$ $`=324(73)+36(2)\text{x}T`$ (12)
while for the $`\chi ^{\alpha 2}`$ curves gave:

$`1/\sqrt{\chi _{100}^{\alpha _2}(T)}`$ $`=27(17)+20(1)\text{x}T`$
$`1/\sqrt{\chi _{001}^{\alpha _2}(T)}`$ $`=37(26)+21(1)\text{x}T`$ (13)
Thus, the thermal evolutions of the two $`\alpha _1`$\- curves (as well the $`\alpha _2`$-curves) are similar: the small numerical discrepancy is attributed to artefacts stemming from the very small values of $`\chi ^{\alpha _i}`$ ($`i=1,2`$).

From these expressions and the above-cited equations, $`M_{\alpha j}^{\alpha _i}`$ ($`i,j`$=1,2) can be straightforwardly calculated. However since the values of the symmetrized elastic constant $`C_0^{\alpha _k}`$($`k`$=1,2,12) are unavailable, these calculation are not sufficient for determining the extremely useful magnetoelastic constants $`D_{\alpha j}^{\alpha _i}`$.

### IV.2 The Ordered Phase

#### IV.2.1 Magnetization

Fig.6a shows representative magnetization isotherms along the a-axis, $`M_{100}`$. Each $`M_{100}(H)`$ isotherm increases monotonically with the field and, for a certain field range, such an increase is linear. Moreover, for low temperatures, $`M_{100}(H)`$ curves demonstrate a weak event at $`H_D`$ $``$12 kOe (defined as the field value at which $`\left(M/H\right)_T`$ is a maximum) and, on a further field increase, a saturation at a relatively high field, for instance $`H_S`$($`1.5`$ K$`)=128`$ kOe. It was observed that on increasing the temperature, the saturation process becomes less pronounced. The thermal evolution of each of $`H_D`$ and $`H_S`$ were plotted in Fig.11. The saturated moment (7.2$`\mu _B`$ at $`T=`$1.5K) is weakly enhanced in comparison with the theoretically expected value (7$`\mu _B`$). This enhancement is attributed to matrix contribution, just as the above-mentioned enhancement of the effective moment.

The thermal variation of the isofield $`M_{100}(T)`$ is shown in Fig.6b. Characteristic features signaling magnetic events at $`H_D`$ and $`H_S`$ were also observed and their values were found to agree with the values obtained from the magnetization isotherms of Fig.6a. This agreement was convincingly demonstrated by plotting in Fig.6b the thermal evolution of each of $`H_D`$ and $`H_S`$. Moreover, Fig.6b shows an extra curve, $`H_R`$, representing the spin reorientation process (see below).

Representative magnetization isotherms along the c-axis ($`M_{001}`$) are shown in Fig.7. These isotherms are very similar to the $`M_{100}`$ ones except that (i) there are only two events namely at $`H_R`$ and $`H_S`$, (ii) the field dependence of the magnetization at these events is more pronounced, and (iii) $`H_R`$ was observed to have a different temperature dependence than the corresponding one along the $`a`$-axis. However, for both orientation, the saturation is reached at the same $`H_S(T)`$ curve. The thermal evolution of $`H_S`$ and $`H_R`$ deduced from Fig.7 were plotted in Fig.12. It is worth noting that within the experimental conditions of the extraction technique, no hysteresis features or relaxation effects were observed.

#### IV.2.2 Magnetostriction

Fig.8a shows the forced magnetostriction curves, $`{}_{}{}^{100}\lambda _{100}^{}(H)`$, at various temperatures. For lower temperatures, a field increase is accompanied by a steep increase in $`{}_{}{}^{100}\lambda _{100}^{}`$, then a shallow maximum, and afterwards a steady decrease with a weak slope. The initial steep increase is very likely to be due to the domain-purifying influence of the applied field: the field derivative of $`\lambda `$ shows a maximum at $`H_D`$ (compare Fig.6a with Fig.8b). Surprisingly the $`{}_{}{}^{100}\lambda _{100}^{}`$ curves present a very small hysteresis suggesting that the domain distribution has a weak dependence on the sample history.

The thermal evolution of the isofield $`{}_{}{}^{100}\lambda _{100}^{}(T)`$ (shown in Fig.9) manifests magnetic events that are very much similar (though more pronounced) to the ones observed in the magnetization (Figs.6b) and forced magnetostriction curves (Fig8.a). The corresponding field and temperature values at these events were plotted in Fig.11.

It is interesting to note that for intermediate fields, e.g. $`H=15`$ kOe in Fig.9, on decreasing the temperature below $`T_N`$, $`{}_{}{}^{100}\lambda _{100}^{}(T)`$ shows at first a slow increase till 17.4 K where a decrease commences and continues till $`T_R.`$ Afterwards $`{}_{}{}^{100}\lambda _{100}^{}(T)`$ increases fast but later tends to saturate to 1.3x10-4 at $`T=0`$ K. In contrast, for lower (higher) fields, $`{}_{}{}^{100}\lambda _{100}^{}(T)`$ is monotonically decreasing (increasing) with decreasing temperature, tending to saturate at $`T`$=0 K. For all applied fields, the signature of the event at $`T_R(H)`$ is clearly observable but tends to fade away as higher fields are approached. It is also interesting to compare the magnitude of this two-ion effect with the one observed in typical single-ion compound such as ErNi2B2C (Ref.Er-magnetostriction ): for zero field, the saturation value of $`{}_{}{}^{100}\lambda _{100}^{}`$ in GdNi2B2C is $``$-0.5x10-4 while in ErNi2B2C is $``$3x10-4.

Fig.10a shows an impressive demonstration of the magnetic domain effects within the ordered state. In these curves, the angular dependence of various $`{}_{}{}^{100}\lambda _{\theta }^{}`$ isotherms was monitored by rotating the single crystal in a field of, say, 30 kOe. Furthermore, Fig.10b shows the orthorhombic strain $`ϵ^\gamma `$ for $`H=`$7.5, 30 kOe, constructed according to Eq.7. Under $`H`$=30 kOe and a decreasing temperature below $`T_N`$, $`ϵ^\gamma `$ commences to increase slowly till $`T_R(=`$15.5 K$`)`$ and afterwards with a faster rate till reaching saturation of -2.6x10$`^{\text{ 4}}`$. Across the borocarbides, this two-ion orthorhombic distortion should vary roughly with the de Gennes factor and as such for ErNi2B2C this partial contribution amounts to $`ϵ^\gamma `$0.4x10-4 which is one third of the total $`ϵ^\gamma `$ $`1.4`$x10-4 observed for ErNi2B2C at 2.5 K (Ref.Er-magnetostriction ): thus for borocarbides, the two-ion magnetoelastic contribution (in particular $`ϵ^\gamma `$) is not negligible as compared to that of the single-ion effect.

#### IV.2.3 The $`HT`$Phase Diagram

Based on the analysis of the magnetization and magnetostriction curves, we were able to construct the $`HT`$phase diagrams with the field either along the a-axis (shown in Fig.11) or along the c-axis (shown in Fig.12). In both diagrams, the boundary describing the saturation process (when the exchange field is completely counterbalanced) is well described by the relation $`H_S\left[1\left(T/T_N\right)^2\right]^{\frac{1}{2}}`$ kOe with $`H_S(T0)`$=128.5(5) kOe and $`T_N(H0)`$=19.5 K. The second boundary, denoted by $`H_R(T)`$, is attributed to the reorientation process: its evolution for $`H//c`$ (much simpler than that for $`H//a`$) is given by the relation $`H_R`$(1-$`T/T_R`$)$`^{\frac{1}{3}}`$ kOe where $`H_R(T0)`$=33.5(5) kOe and $`T_R(H0)`$=13.5 K. The extrapolated values of the characteristic points$`H_S,`$ $`T_N,`$ $`T_R`$ are in excellent agreement with the reported ones.Gd-XRES ; Gd-single-crystal-Tdynamics ; Gd-poly-high-field ; Gd-chemical-composition ; Gd-MES  The third boundary $`H_D(T),`$ appearing only in the a-axis phase diagram, (see Fig.11) is related to the domain effects as can be inferred from the angular dependence of Fig.10.a. This $`H_D(T)`$ has the following features: (i) the extrapolated $`H_D(T0)`$=12.0(5) kOe while $`T_D(H0)`$ $`=T_N(H0)`$ $`=`$19.5 K, (ii) it is almost constant for $`T<`$ 8 K, increases slowly up to a maximum around 16 K and afterwards decreases steadily till vanishes at $`T_N`$: such an increase of $`H_D`$ with temperature is not usual for domain effects.

## V Discussion

The modulated state in GdNi2B2C (and borocarbides in general) is usually related to the nesting features in their electronics structures.Generalized-X  Below, we discuss the stability of this state against field and temperature variation. First MFA arguments are applied to account for $`H_S`$, $`H_D`$, and the linearity of the magnetization isotherms. Next, anisotropic perturbation are introduced and and afterwards the anisotropic exchange interaction and magnetoelastic couplings are considered.

Among the borocarbides, GdNi2B2C has the strongest exchange interactions and the weakest magnetic anisotropy. The strength of the exchange interactions (denoted as $`H_{ex}`$, or $`J_{eff}`$) is reflected as a high value of $`T_N`$ and $`H_S`$ (see Figs.11-12) while the weakness of the magnetic anisotropy (denoted as $`H_{an}`$) is manifested as low $`H_D`$ and as a weak difference between the parallel and perpendicular paramagnetic susceptibility. MFA arguments can be used to calculate an effective value for each of $`H_{ex}`$, $`J_{eff}`$, and $`H_{an}`$ from the observed values of $`T_N`$, $`H_S,`$ and $`H_D`$. Let us, at first, ignore the anisotropic feature ($`H_{an}`$=0). Then, from the relation $`zJ_{eff}=\frac{3T_N}{2S(S+1)}`$ (z is an effective nearest neighbors), $`zJ_{eff}`$ was calculated $``$1.86 K. Substitution of this value into the relation $`H_S=2zJ_{eff}S(S+1)/g\mu _BS`$ gave $`H_S`$ $``$125 kOe, a reasonable value for the effective saturation field. Next, taking into account the case $`H_{an}0`$ and using Eqs.14, one is able to calculate an effective value for each of $`H_{ex}`$ and $`H_{an}`$ and furthermore to account for the observed linearity of $`M_{100}(H_D<H<H_S)`$ : Smart

$`H_S`$ $`=(H_{ex}H_{an})`$
$`H_D`$ $`=\sqrt{H_{an}H_{ex}H_{an}^2}`$
$`M_{100}/M_s`$ $`=H/H_S,\text{ for }H_D<H<H_S`$ (14)
Substituting the extrapolated $`H_S(T0`$ K) and $`H_D(T0`$ K) (see Fig.11 below) one obtains the satisfactory value $`H_{ex}(T0`$ K) = 129.6(5) kOe and $`H_{an}(T0`$ K)=1.1(3) kOe: as expected for a Gd-based compound, the former is much stronger than the latter.

In the above treatment, the zero-temperature magnetic structure is implicitly assumed as a collinear 3d AF state. Rather, for $`T_R(H=0)<T<T_N,`$ the magnetic structure is a modulated state with no c-component and the ordered moments are distributed among the a- or b-polarized domains. Then applying a field along the a-axis would force the a-polarized domains to align perpendicular to the field at $`H_D(T>T_R)`$ and later on to achieve saturation at $`H_S(T>T_R)`$. For the same temperature range, a field applied along the c-axis is already perpendicular to the domains and as such only the saturation transition would appear. On the other hand for $`T<T_R(H=0),`$ the spontaneous emergence of the c-component at $`H_R(T)`$ constitutes an extra event for both field orientations: thus for $`H//`$a-axis, there are three transitions (at $`H_D`$, $`H_R`$, and $`H_S)`$ while for the $`H//`$c-axis, there are only two transitions (at $`H_R`$ and $`H_S)`$.

Within the ordered phase, the manifestation of the magnetoelastic coupling and anisotropic exchange interactions is much more pronounced, leading to a substantial perturbation of the nuclear and the magnetic structures. To demonstrate the importance of these perturbations, let us divide the ordered phase into two regions depending on the magnitude of $`ϵ^\gamma `$. The first region is$`T_R<T<T_N`$ wherein the low-field $`ϵ^\gamma `$ is relatively weak (see Fig.10) and as such is hardly sufficient to modify (substantially) the nuclear structure or the modulated magnetic state. The second region is $`T<T_R`$ wherein the magnitude of the low-field $`ϵ^\gamma `$ is very large (see Fig.10) and increases whenever the temperature is decreased or the field is increased. As a consequence, the nuclear structure is distorted ($`\frac{a}{b}<1`$ and the four-fold symmetry of the basal plane is reduced to a two-fold symmetry). In addition, there is an accompanying dramatic modification of the magnetic structure: as the temperature is decreased, the anisotropic character of the exchange interactions (just as the magnitude of $`ϵ^\gamma `$) is enhanced and the magnetic structure is more and more perturbed towards a squared-up state wherein the magnetic moments approach equal amplitudes and progressively bunch towards an orientation within the b-c plane. While the square-up character can be understood based on entropy arguments, the bunching of the moments is most probably induced by the above mentioned perturbations.

By demonstrating the presence of magnetoelastic and anisotropic exchange couplings in the ordered phase of GdNi2B2C, this works confirms the far-sight analysis of Detlefs et al Gd-XRES  and Tomala et al Gd-MES . On the one hand, Detlefs et al attributed the strong asymmetry in the line shape of the resonant XRES below $`T_R(H=0)`$ as being due to a straining. Consequently, their XRES spectra were analyzed as being composed of two peaks (thus the asymmetric character) below $`T_R(H=0)`$ and as a single line above $`T_R(H=0)`$. It is assuring to note that the thermal evolution of the position and intensity of these XRES peaks correlate very well with $`ϵ^\gamma `$ (compare Fig.7 of Ref.Gd-XRES  with our Fig.10). Tomala et al, on the other hand, analyzed their 155Gd Mossbauer spectra assuming a transverse (moments along b-axis) sine-modulated state above $`T_R(H=0)`$ and a strongly squared-up and bunched state below $`T_R(H=0)`$.

## VI Conclusion

We presented experimental results on single crystal of GdNi2B2C. The parastriction and paramagnetic features, well accounted for within the susceptibility formalism, demonstrated the presence of weak anisotropic exchange couplings and two-ion $`ϵ^{\alpha _1}`$ and $`ϵ^{\alpha _2}`$ strain modes. The strength of the latter coefficients indicate that the strain derivative of the coupling constant at the Brillouin zone centre have significant amplitudes. The ordered state, on the other hand, is characterized by very rich field and temperature phase diagrams that manifest high values of $`T_N`$ and $`H_S`$ (indicative of strong exchange interactions), a spin reorientation process at $`H_R(T)`$, and a domain-purifying field $`H_D`$. Anisotropic exchange interactions and relatively strong magnetoelastic couplings (particularly the orthorhombic distortion mode) were observed, lending support to the reported modification of the magnetic structure below $`T_R`$.

It is worth recalling that the modulated state of GdNi2B2C within the temperature range $`T_R<T<T_N`$ is very much similar to that of ErNi2B2C within the range 2 K$`<T<T_N`$, even though the dominant aligning forces operating in each isomorphs are very much different. Since the exchange interaction and the magnetoelastic couplings are driven by the indirect RKKY couplings and that de Gennes scaling holds well among the heavy borocarbides,Canfield-physics-today  then both the anisotropy exchange interactions and magnetoelastic couplings must be present (though with varying strength depending on de Gennes factor) in other borocarbide magnets. As such, these perturbations must be explicitly considered when discussing the magnetic properties (in particular the $`HT`$phase diagrams) of these magnets.

###### Acknowledgement 1

M. ElM. acknowledge partial financial support from the Japan Society for the promotion of science, the University of Joseph Fourier of Grenoble, and the Brazilian agencies CNPq and FAPERJ.