text:
# 𝑌-Scaling Analysis of the Deuteron Within the Light-Front Dynamics Method

## I Introduction

High-energy electron scattering from nuclei can provide important information on the wave function of nucleons in the nucleus. In particular, using simple assumptions about the reaction mechanism, scaling functions can be deduced that, if shown to scale (i.e., if they are independent of the momentum transfer), can provide information about the momentum and energy distribution of the nucleons. Several theoretical studies  have indicated that such measurements may provide direct access to studies of short-range nucleon-nucleon (NN) correlation effects.

Since West’s pioneer work , there has been a growth of interest in $`y`$-scaling analysis, both in its experimental and theoretical aspects. This is motivated by the importance of extracting nucleon momentum distributions from the experimental data. West showed that in the impulse approximation, if quasielastic scattering from a nucleon in the nucleus is a dominant reaction mechanism, a scaling function $`F(y)`$ could be extracted from the measured cross section which is related to the momentum distribution of the nucleons. In the simplest approximation the corresponding scaling variable $`y`$ is the minimum momentum of the struck nucleon along the direction of the virtual photon. In principle, the scaling function $`F(y,Q^2)`$ depends on both $`y`$ and momentum transfer $`Q^2`$ (-$`Q^2`$ is the square of the four-momentum transfer), but at sufficiently high $`Q^2`$ values the dependence on $`Q^2`$ should vanish yielding scaling. However, any attempt to extract the momentum distribution from the $`y`$-scaling in electron-nucleus scattering faces the problem of estimating both effects from the final-state interactions (FSI) of the struck nucleon with the rest of the nucleus and from the nucleon binding. Previous calculations  suggest that the contribution from the final-state interactions should vanish at sufficiently high $`Q^2`$. The FSI lead to sizable scaling violation effects only at low values of the three-momentum transfer $`|𝐪|=\sqrt{Q^2+\nu ^2}`$ ($`\nu `$ is the energy loss) . The most important dynamical effects, such as binding corrections, which represent the fact that for complex nuclei the final spectator $`A1`$ system can be left in all possible excited states including the continuum, have been treated in  in terms of spectral functions. This problem has been solved in  by introducing a new scaling variable which gives direct, global and independent of $`A`$ link between the experimental data and the longitudinal momentum components.

Recently inclusive electron scattering has been studied at the Thomas Jefferson National Accelerator Facility (TJNAF) with 4.045 GeV incident beam energy from C, Fe and Au targets  to $`Q^2`$ $``$ 7 (GeV/c)2. Data were also taken using liquid targets of hydrogen and deuterium . The data presented in  represent a significant increase of the $`Q^2`$ range compared to previous SLAC measurement , in which an approach to the scaling limit for heavy nuclei is suggested but for low values of $`|y|<`$ 0.3 GeV at momentum transfers up to 3 (GeV/c)2, and, moreover, a scaling behaviour is observed for the first time at very large negative $`y`$ ($`y`$= -0.5 GeV/c). From theoretical point of view the extended region of $`y`$ measured at TJNAF is of significant importance since this is a regime where the nucleon momentum distribution is expected to be dominated by short-range NN correlations. On the other hand, it is interesting to note that contributions from short-range FSI may also result in a scaling-like behavior due to the small $`Q^2`$ dependence of these effects, and that these contributions are also dominated by short-range correlations. Obviously, a complete understanding of this electron-nucleus scattering requires a relativistic approach to the quantities related to the $`y`$-scaling analysis for a detailed comparison with the experimental data.

A relativistic $`y`$-scaling has been considered in  by generalizing the nonrelativistic scaling function to the relativistic case. Realistic solutions of the spinor-spinor Bethe-Salpeter (BS) equation for the deuteron with realistic interaction kernel were used for systematic investigation of the relativistic effects in inclusive quasielastic electron-deuteron scattering. The approach of $`y`$-scaling presented in  is fully covariant, with the deuteron being described by eight components, namely the $`{}_{}{}^{3}S_{1}^{++}`$, $`{}_{}{}^{3}D_{1}^{++}`$, $`{}_{}{}^{3}S_{1}^{}`$, $`{}_{}{}^{3}D_{1}^{}`$, $`{}_{}{}^{1}P_{1}^{+}`$, $`{}_{}{}^{1}P_{1}^{+}`$, $`{}_{}{}^{3}P_{1}^{+}`$, and $`{}_{}{}^{3}P_{1}^{+}`$ waves. The first two waves directly correspond to the $`S`$ and $`D`$ waves in the deuteron, with the waves with negative energy vanishing in the nonrelativistic limit. It has been demonstrated in  that, if the effects from the negative energy $`P`$-states are disregarded, the concept of covariant momentum distribution can be introduced.

Recently a successful relativistic description of the nucleon momentum distribution in deuteron has been done  within the light-front dynamics method . The most important result from the calculations in  is the possibility of the LFD method to describe simultaneously both deuteron charge form factors (that has been shown in ) and the momentum distribution. It is shown in  that after the projection of the Bethe-Salpeter amplitude on the light front, the six components of the LFD deuteron wave function are expressed through integrals over the eight components of the deuteron Bethe-Salpeter amplitude. Provided the nucleon-nucleon interaction is the same, these approaches incorporate by different methods the same relativistic dynamics. The wave functions in LFD are the direct relativistic generalization of the nonrelativistic ones in the sense that they are still the probability amplitudes. Therefore they can be used in the relativistic nuclear physics (e.g. ).

The aim of our work is using the nucleon momentum distribution $`n(k)`$ obtained with the LFD method to calculate the deuteron scaling function. The result for the asymptotic function is compared with the recent TJNAF data measured at six values of $`Q^2`$. In particular, the scaling behavior observed for very large negative $`y`$ providing momenta higher than those corresponding to existing experimental data for $`n(k)`$ may allow to distinguish the properties of the covariant LFD method from the potential approaches. The comparison with the BS result for the scaling function serves as a test for the consistency of both covariant approaches treating the deuteron relativistically in the case of $`y`$-scaling.

The paper is organized as follows. In Section II the definition and the physical meaning of the scaling variable and the scaling function are briefly reviewed together with some basic relations between the nucleon spectral function, the scaling function and the momentum distribution. The results for the nucleon momentum distribution in deuteron obtained within the LFD method are given in Section III. The calculated LFD asymptotic scaling function of the deuteron is presented in Section IV, where a comparison with the Bethe-Salpeter result and with some nonrelativistic calculations is also done. The summary of the present work is given in Section V.

## II Basic relations in the $`y`$-scaling method

The scaling function is defined as the ratio of the measured cross section to the off-shell electron-nucleon cross section multiplied by a kinematic factor:

$$F(q,y)=\frac{d^2\sigma }{d\mathrm{\Omega }d\nu }(Z\sigma _p+N\sigma _n)^1\frac{q}{[M^2+(y+q)^2]^{1/2}},$$
(1)
where $`Z`$ and $`N`$ are the number of protons and neutrons in the target nucleus, respectively, $`\sigma _p`$ and $`\sigma _n`$ are the off-shell cross sections, and $`M`$ is the mass of the proton. In analysing quasielastic scattering in terms of the $`y`$-scaling a new variable $`y=y(q,\nu )`$ is introduced. Then the nuclear structure function which is determined using the spectral function $`P(k,E)`$ as

$$F(q,\nu )=2\pi _{E_{min}}^{E_{max}(q,\nu )}𝑑E_{k_{min}(q,\nu ,E)}^{k_{max}(q,\nu ,E)}k𝑑kP(k,E),$$
(2)
can be expressed in terms of $`q`$ and $`y`$ rather than $`q`$ and $`\nu `$ (see Eq. (1)). In Eq. (2) $`E=E_{min}+E_{A1}^{}`$ is the nucleon removal energy with $`E_{A1}^{}`$ being the excitation energy of the final $`A1`$ nucleon system.

The most commonly used scaling variable $`y`$ is obtained  starting from relativistic energy conservation, setting $`k=y`$, $`\frac{\mathrm{𝐤𝐪}}{kq}=1`$ and the excitation energy $`E_{A1}^{}`$=0, and is defined through the equation

$$\nu +M_A=(M^2+q^2+y^2+2yq)^{1/2}+(M_{A1}^2+y^2)^{1/2},$$
(3)
where $`M_A`$ is the mass of the target nucleus and $`M_{A1}`$ is the mass of the $`A1`$ nucleus. Therefore, $`y`$ represents the longitudinal momentum of a nucleon having the minimum removal energy ($`E=E_{min}`$, i.e. $`E_{A1}^{}`$=0).

At high values of $`q`$ a pure scaling regime is achieved, where $`k_{min}|y(EE_{min})|`$ and Eq. (2) becomes

$$F(q,y)f(y)=2\pi _{E_{min}}^{\mathrm{}}𝑑E_{|y(EE_{min})|}^{\mathrm{}}k𝑑kP(k,E).$$
(4)
In Eq. (4) the particular behavior of $`P(k,E)`$ at large $`k`$ and $`E`$ is used in order to extend the upper limits of integration to infinity .

In the deuteron one always has $`E_{A1}^{}`$=0, so that the spectral function is entirely determined by the nucleon momentum distribution $`n(k)`$, i.e. $`P(k,E)=n(k)\delta (EE_{min})`$, and, consequently, $`k_{min}=|y|`$ for any value of $`q`$. The scaling function (4) reduces to the longitudinal momentum distribution

$$f(y)=2\pi _{|y|}^{\mathrm{}}k𝑑kn(k).$$
(5)
## III Nucleon momentum distribution in the light-front dynamics method

The relativistic deuteron wave function (WF) on the light-front $`\mathrm{\Psi }(\stackrel{}{k},\stackrel{}{n})`$ depends on two vector variables: i) the relative momentum $`\stackrel{}{k}`$ and ii) the unit vector $`\stackrel{}{n}`$ along $`\stackrel{}{\omega }`$ which determines the position of the light-front surface. Due to this, the WF is determined by six invariant functions $`f_i`$ ($`i`$=1,…,6) instead of two ($`S`$\- and $`D`$-waves) in the nonrelativistic case. Each one of these functions depends on two scalar variables $`k`$ and $`z=cos(\widehat{\stackrel{}{k},\stackrel{}{n}})`$. In LFD these six functions are calculated within the relativistic one-boson-exchange model. As shown in , in the nonrelativistic limit the functions $`f_{36}`$ become negligible, $`f_{1,2}`$ do not depend on $`z`$ and turn into $`S`$\- and $`D`$-waves ($`f_1u_S`$, $`f_2u_D`$) and the wave function $`\mathrm{\Psi }(\stackrel{}{k},\stackrel{}{n})`$ becomes the usual nonrelativistic wave function. One of the most important properties of the functions $`f_{16}`$ found in  is that for $`k`$ 2$`÷`$2.5 fm-1 the component $`f_5`$ (being related mainly to $`\pi `$-exchange) exceeds sufficiently the $`S`$\- and $`D`$-waves. This fact is very important in the calculations of $`n(k)`$ in deuteron as it will be shown below.

The LFD calculations have shown (for more details, see Ref. ) that, as expected, the most important contributions to the total $`n(k)`$ give terms related to the $`f_1`$, $`f_2`$ and $`f_5`$ functions

$$n(k)n_1(k)+n_2(k)+n_5(k).$$
(6)
The contributions of $`n_1`$, $`n_2`$, $`n_{12}=n_1+n_2`$ and $`n_5`$ are compared in Fig. 1. It can be seen that, while the functions $`f_1`$ and $`f_2`$ give a good description of the $`y`$-scaling data of $`n(k)`$ for $`k<2`$ fm-1 (like the $`S`$\- and $`D`$-wave functions in the nonrelativistic case), it is impossible to explain the high-momentum components of $`n(k)`$ at $`k>2`$ fm-1 without the contribution of the function $`f_5`$. We note that the deviation of the total $`n(k)`$ from the sum $`n_{12}=n_1+n_2`$ starts at $`k`$ around 1.8 fm-1. All this shows the important role of NN interactions which incorporate exchange of relativistic mesons in the case of the deuteron.

The nucleon momentum distribution in the deuteron using $`S`$\- and $`D`$-wave functions $`\mathrm{\Psi }_S(k)`$ and $`\mathrm{\Psi }_D(k)`$ corresponding to various NN potentials, such as the charge-dependent Bonn potential , the Argonne v18 , the Nijmegen - I ,- II and - Reid 93  and Paris 1980  is also computed in  by the expression:

$$n(k)=\frac{1}{4\pi }[\mathrm{\Psi }_S^2(k)+\mathrm{\Psi }_D^2(k)]n_S(k)+n_D(k)$$
(7)
with

$$n(k)d^3\stackrel{}{k}=1.$$
(8)
In Fig. 2 the result for $`n(k)`$ using the charge-dependent Bonn potential  is given and compared with the $`y`$-scaling data. As can be seen, the $`D`$-component of $`n(k)`$ is important but even its inclusion does not give a very good agreement with the data for $`k2`$ fm-1. In the next Fig. 3 we present the LFD result for $`n(k)`$ compared with the calculations using the WF’s corresponding to Nijmegen-I,-II, -Reid 93, Argonne v18 and Paris 1980 NN potentials and with the $`y`$-scaling data. We would like to note that: i) the results of the calculations using the NN potentials, such as Nijmegen-II, -Reid 93, Argonne v18 and Paris 1980 (shown in Fig. 3) are in better agreement with the $`y`$-scaling data than those using the charge-dependent Bonn potential (Fig. 2). This might be related to the fact that these potentials describe NN phase shifts up to larger energies (e.g. the Nijmegen-II potential gives reasonable $`pp`$ phase shifts up to 1.2 GeV, while the charge-dependent one-boson exchange Bonn potential fits the phase-shift data below 350 MeV); ii) It can be seen from Fig. 3 that there are small differences between the curves corresponding to different NN potentials for $`k`$ 3 fm-1 (which give a good description of the $`y`$-scaling data and almost coincide with the LFD result) and larger ones for $`k>3`$ fm-1. Large differences take place, however, between all of them and the LFD result for $`k>3`$ fm-1.

## IV Results for the asymptotic scaling function of deuteron

In this Section we present the results for the asymptotic scaling function of deuteron which is calculated by using of the LFD nucleon momentum distribution given in Section III.

The scaling function for deuteron calculated within the LFD method is shown in Fig. 4. It is compared with the TJNAF experimental data  for all measured angles. The $`Q^2`$ values are given for Bjorken $`x=Q^2/2M\nu =1`$ and correspond to elastic scattering from a free nucleon. It is seen from Fig. 4 that the relativistic LFD scaling function is in good agreement with the data in the whole region of negative $`y`$ available. As known, the scaling breaks down for values of $`y>0`$ due to the dominance of other inelastic processes beyond the quasielastic scattering. Our LFD deuteron scaling function is also compared in Fig. 4 with the scaling function obtained within the BS formalism . A small difference between the two results is observed for $`y<400`$ MeV/c but, at the same time, the theoretical LFD scaling function is closer to the experimental data in the same region of $`y`$. The fact that both LFD and BS functions reveal similar behavior is a strong indication in favor of the consistency of the two relativistic covariant approaches in case of the $`y`$-scaling.

In Fig. 5 the asymptotic relativistic LFD and BS scaling functions $`f(y)`$ are compared with the nonrelativistic ones, calculated with some realistic interactions. We have already shown (see Figs. 2 and 3 of the present work) that the results using these potentials explain very well (with the exception of the charge-dependent Bonn potential) all the available data for $`n(k)`$ up to $`k3`$ fm-1 exactly like the LFD method. However, it is concluded in  that for $`k>3`$ fm-1 the LFD results for $`n(k)`$ deviate strongly from those of the calculations using NN potentials. Here we would like to emphasize the existence of the same discrepancies between the scaling functions observed from the comparison in Fig. 5. It is shown that for $`|y|>`$ 400 MeV/c both LFD and BS curves start to deviate from the nonrelativistic scaling functions. The result for $`f(y)`$ calculated using the Nijmegen-II NN potential is in better agreement with the experimental data than those using other potentials. It is in accordance with the result for $`n(k)`$ shown in Fig. 3. For instance, by a thorough comparison between the relativistic Bethe-Salpeter and the nonrelativistic scaling functions of deuteron it has been found in  that the two functions start to sensibly differ also at $`|y|>`$ 400 MeV/c. Thus, the necessity to treat realistically the relativistic dynamics inside the deuteron in a way different from the potential approaches becomes apparent. In this sense, the results calculated for both momentum distribution and asymptotic scaling function confirm the abilities of the LFD method to describe with a good accuracy the experimental data measured at high momentum transfers.

## V Conclusions

In the present paper inclusive electron-deuteron scattering data have been analyzed in terms of the $`y`$-scaling function within the light-front dynamics method. For this purpose, the nucleon momentum distribution in deuteron has been used in order to calculate the asymptotic scaling function. For the trivial case of deuteron, for which the structure function (Eq. 4) coincides with the longitudinal momentum distribution (Eq. 5) we have found a good agreement of the calculated scaling function with the experimental data. Thus, the concept of relativistic $`y`$-scaling can be introduced in the LFD relativistic description of inclusive quasielastic $`eD`$ scattering, in the same way as it is done in the conventional nonrelativistic approach, i.e. by introducing a scaling function (which, in the scaling regime, is nothing but the nucleon longitudinal momentum distribution), and in terms of the same variable $`y`$. It has been pointed out that for $`|y|>`$ 400 MeV/c the differences between the LFD and the nonrelativistic scaling functions are very large.

Exploring the light-front dynamics, we continue in this paper our analysis of important deuteron characteristics. The effective inclusion of the relativistic nucleon dynamics and of short-range NN correlations can be better seen when analyzing electron scattering at high momentum transfer from complex nuclei, for which a proper theoretical $`y`$-scaling analysis is still lacking. Although scaling violation effects due to final-state interactions sharply decrease with increasing momentum transfer, a consistent treatment of both FSI and nucleon binding must be made in order to perform a precise comparison with the new TJNAF data. Such an investigation is in progress.

## ACKNOWLEDGMENTS

This work was partly supported by the Bulgarian National Science Foundation under the Contracts Nrs.$`\mathrm{\Phi }`$–809 and $`\mathrm{\Phi }`$–905.


format: text

text:
# Contact Line Instability and Pattern Selection in Thermally Driven Liquid Films

## I Introduction

Driven spreading of liquid films is a process which occurs in numerous industrial coating applications, such as spin-coating of hard drives, so understanding its dynamics and learning to control it is very important. For instance, instabilities which arise during the spreading of the liquid on the solid substrate can lead to nonuniform coverage, adversely affecting the quality of produced coating. Driven spreading of thin films and patterning also have important implications for microfluidics.

Driven spreading of liquid films under the action of gravity , centrifugal acceleration , thermocapillary effects , or combination thereof  has been extensively studied in the literature. Stability analysis of such flows has attracted the most attention and the mechanism of the linear instability is now well understood. Considerable progress has also been reached in feedback control of flat liquid layers , whose dynamics is governed by normal differential operators. The attempts to influence the stability of spreading films have so far been limited to non-feedback control achieved through either imposing an externally generated counterflow  or chemically patterning the substrate .

This study represents the first theoretical treatment of the feedback control problem for spreading films. The spatially nonuniform nature of spreading films and the presence of mean flow make the control problem much more difficult compared to the case of flat stationary films, because the dynamics of disturbances in spreading films is governed by a nonnormal evolution operator and thus requires a more careful analysis. For instance, liquid films flowing down an incline have been found to develop a contact line instability when the linear analysis predicted stable evolution . We derive the slip model of thermally driven spreading and use it to show that the contact line instability can be suppressed using adaptive thermal perturbations which depend on the distortion of the contact line. (This type of feedback is chosen because it is easiest to implement experimentally with sufficient spatial and temporal resolution via optical means .) Although the results of the following analysis should be applicable regardless of the driving force, we concentrate our attention on the case of thermal driving.

The layout of the paper is as follows. The slip model of thermally driven liquid films is derived in Section II and its linear stability analysis is conducted in Section III. The validity of linear stability analysis and transient effects are considered in Section IV. Section V presents the proposed algorithm for feedback control of the contact line instability and Section VI contains the conclusions.

## II Slip model for thermal spreading

We consider the spreading of a thin layer of partially wetting liquid on a horizontal substrate (see Fig. 1). The spreading process is conventionally described using the lubrication approximation , with the horizontal velocity governed by the Stokes equation

$$\mu _{zz}𝐯=\overline{p},$$
(1)
where $`\mu `$ is the dynamic viscosity, $`\overline{p}`$ is the modified pressure, $`=(_x,_y)`$ is the 2-dimensional gradient operator, and the vertical velocity is neglected.

It is well known  that the standard no-slip boundary condition at the liquid-solid interface results in a stress singularity at the contact line. The only approach explored in the literature for the thermally driven case was to relieve this singularity by introducing a thin precursor film . However, as we intend to use the position of the contact line in our control algorithm later, the precursor model becomes inconvenient. Instead we choose to specify a microscopic contact angle and relieve the stress singularity by employing a partial slip boundary condition

$$𝐯=\frac{\alpha }{3h}_z𝐯$$
(2)
at the bottom of the liquid layer, where $`h`$ is the local thickness of the film, and $`\alpha `$ is the phenomenological slip coefficient. The slip boundary condition (2) was originally introduced by Greenspan  for modeling the unforced spreading of liquid drops and later used to model the forced spreading of liquid films under the action of gravity . At the free surface the standard stress balance boundary condition

$$\mu _z𝐯=\sigma $$
(3)
applies, where $`\sigma `$ is the surface tension coefficient. Solving (1) subject to these boundary conditions we obtain the horizontal velocity

$$𝐯=\frac{1}{\mu }z\sigma \frac{1}{\mu }\left(\frac{\alpha }{3}+hz\frac{1}{2}z^2\right)\overline{p}.$$
(4)
In order to make the phenomenological boundary condition (2) consistent with the physics of the flow, in (4) we have dropped an unphysical term $`\alpha \sigma /3\mu h`$ which diverges for vanishing film thickness. This divergence is not necessarily a cause for alarm as the continuous approximation underlying the Stokes equation itself breaks down in this limit. The shear stress also becomes poorly defined for very thin films. Our choice of the functional form for the horizontal velocity (4), therefore, amounts to picking an appropriate phenomenological model in the region where the continuous description of the flow becomes invalid and fluctuations become important.

For a film which is sufficiently thin the hydrostatic pressure can be ignored, so that the modified pressure is given by the normal component of the surface tension $`\overline{p}=\sigma \kappa =\sigma ^2h`$. Substituting (4) into the mass conservation condition

$$_th=_0^h(𝐯)𝑑z$$
(5)
and integrating we obtain an evolution equation for the thickness:

$`_th=\left[{\displaystyle \frac{1}{2\mu }}h^2\sigma +{\displaystyle \frac{1}{3\mu }}(\alpha h+h^3)(\sigma ^2h)\right].`$ (6)
Now consider the situation which arises when the substrate covered by the liquid film is subjected to a linear temperature gradient in the $`x`$-direction. Assuming that the surface tension changes linearly with temperature $`\theta `$,

$$\sigma (x)=\sigma (\theta _0)+x_x\theta _\theta \sigma \sigma _0\tau x,$$
(7)
and neglecting the variation in $`\sigma `$ in the second term of (6), which produces subdominant contribution (see, e.g., the discussion in Ref. ), we obtain

$`_th={\displaystyle \frac{\tau }{2\mu }}_xh^2`$ $``$ $`{\displaystyle \frac{\sigma _0}{3\mu }}_x[(\alpha h+h^3)(_{xxx}h+_{xyy}h)]`$ (8)
$``$ $`{\displaystyle \frac{\sigma _0}{3\mu }}_y[(\alpha h+h^3)(_{xxy}h+_{yyy}h)].`$ (9)
We can absorb most parameters into the spatial and temporal scales by introducing the nondimensional variables $`t^{}=t/T`$, $`x^{}=x/X`$, $`y^{}=y/X`$, and $`h^{}=h/H`$. The vertical length scale $`H`$ is defined by the characteristic thickness of the film and sets both the horizontal length scale $`X=[2\sigma _0H^2/3\tau ]^{1/3}`$ and the time scale $`T=2\mu X/\tau H`$. After defining the dimensionless slip coefficient $`\alpha ^{}=\alpha /H^2`$ and dropping the primes we obtain

$`_th=_xh^2`$ $``$ $`_x[(\alpha h+h^3)(_{xxx}h+_{xyy}h)]`$ (10)
$``$ $`_y[(\alpha h+h^3)(_{xxy}h+_{yyy}h)].`$ (11)
The obtained equation has the same form as the one describing gravity driven rather than temperature driven films (see, e.g., equation (33) in Ref. ), with the exception that $`h^2`$ in the first term on the right-hand-side is replaced with $`\alpha h+h^3`$. This similarity of the structure of the governing equations suggests that the feedback control problem can be treated in the same way regardless of the nature of the driving force.

The liquid spreads in the direction opposite to the temperature gradient, so the motion of the contact line is most conveniently described in the reference frame moving with speed $`u`$ towards negative $`x`$. In this frame the evolution equation possesses a transversely uniform steady state solution, which gives the asymptotic film profile for constant flux boundary conditions. Substituting $`h(x,y,t)=h_0(x+ut)`$ into (10) and integrating once we obtain

$$uh_0h_0^2+(\alpha h_0+h_0^3)h_0^{\prime \prime \prime }=d,$$
(12)
where prime indicates the differentiation with respect to the $`x`$ coordinate. The constants $`u`$ and $`d`$ can be determined from the appropriate boundary conditions. Following Spaid and Homsy  we require that the film thickness vanish at the contact line, $`h_0(0)=0`$, and specify the slope $`h_0^{}(0)=c(X/H)\mathrm{tan}\varphi `$, where $`\varphi `$ is the microscopic dynamic contact angle in the dimensional variables. Furthermore, the constant flux boundary condition at the tail of the film far away from the contact line gives $`h_0(\mathrm{})=1`$ (this choice corresponds to using the dimensional tail thickness as the vertical lengthscale, $`H`$) and $`h_0^{\prime \prime \prime }(\mathrm{})=0`$. These boundary conditions require $`d=0`$ and $`u=1`$, and consequently

$$h_0^{\prime \prime \prime }=\frac{h_01}{h_0^2+\alpha }.$$
(13)
The solution to this equation describes the height profile of the spreading film once the distance from the contact line to the reservoir becomes sufficiently large.

At this point it is appropriate to make two comments regarding the structure of the lubrication equations produced by the slip model. First of all, from a mathematical perspective, introduction of partial slip at the liquid-solid interface is equivalent to a regularization procedure for a singular differential equation. Dropping the diverging term in (4) is equivalent to keeping the regularization parameter $`\alpha `$ only in the terms in (6) responsible for the singular behavior of the solution at the contact line. We could have kept all the terms in (4) and (6) just as well. Neither (10) nor (13) would have changed.

Second, (13) is very similar to the equation produced by the precursor model of Kataoka and Troian

$$h_0^{\prime \prime \prime }=\frac{(h_01)(h_0b)}{h_0^3},$$
(14)
where $`b`$ is the thickness of the precursor layer. In fact, the solutions of the two equations are virtually indistinguishable (see Fig. 2) for the proper choice of parameters (precursor thickness, slip coefficient, and contact angle). This is hardly a surprise, as the two equations become identical in the limit $`b0`$ and $`\alpha 0`$, and proves that the final result is essentially independent of the regularization procedure used to get rid of the stress singularity at the contact line.

The slip model has two independent parameters ($`c`$ and $`\alpha `$), while the precursor model has only one ($`b`$). A closer inspection shows that the gross features of the slip model are, in fact, determined by a single parameter (say, the height of the capillary ridge). As Fig. 3 shows, one can change the value of the slip coefficient by two orders of magnitude without causing noticeable changes in the film profile (the contact angle will also change in response to the changes in $`\alpha `$ to keep the height of the capillary ridge constant). The second parameter fine tunes the shape of the film in the immediate vicinity of the contact line. The slip model is, therefore, not only more convenient for the purpose of control (it explicitly defines the easy to measure position of the contact line), but it is to some extent more flexible than the precursor model.

Finally, let us explore the dependence of the film profile on the parameters. As we have just seen, the profile is rather insensitive to the changes in the slip coefficient, so let us fix it and vary the contact angle (see Fig. 4). The first important observation is that regardless of the magnitude of the contact angle, the capillary ridge never disappears (this turns out to be the case for any reasonable value of the slip coefficient, i.e., $`\alpha <1`$). This is in contrast with the results obtained for gravity driven films, where the capillary ridge completely disappears at small inclination angles . Second, the height of the ridge is a monotonically increasing function of $`c`$ and becomes essentially independent of the contact angle for $`c<0.3`$. The dependence on the slip coefficient is more subtle: the value of $`\alpha `$ controls the minimum height of the capillary ridge which grows with decreasing slip coefficient.

## III Contact line instability

Linear stability of the asymptotic solution $`h_0(x+ut)`$ can be determined in a standard way. Since this solution is uniform in the transverse direction, the linearized equation describing the dynamics of small disturbances can be partially diagonalized by Fourier transforming it in the $`y`$-direction. By substituting

$$h(x,y,t)=h_0(x+ut)+ϵg(x+ut,t)e^{iqy}$$
(15)
into (10) and retaining terms of order $`ϵ`$, we obtain

$$_tg=L(q)g,$$
(16)
where $`L(q)=L_0+q^2L_1+q^4L_2`$ is a $`4^{\mathrm{t}h}`$ order differential operator defined via

$`L_0g`$ $`=`$ $`[\{12h_0+(\alpha +3h_0^2)h_0^{\prime \prime \prime }\}g+(\alpha h_0+h_0^3)g^{\prime \prime \prime }]^{},`$ (17)
$`L_1g`$ $`=`$ $`[(\alpha +3h_0^2)h_0^{}g]^{}+2(\alpha h_0+h_0^3)g^{\prime \prime },`$ (18)
$`L_2g`$ $`=`$ $`(\alpha h_0+h_0^3)g.`$ (19)
The boundary condition on $`g`$ at the contact line

$$h_0^{\prime \prime }(0)g(0)=h_0^{}(0)g^{}(0)$$
(20)
can be obtained by requiring the perturbed solution to have the same contact angle as the unperturbed solution. The other three boundary conditions can all be imposed at the tail, say,

$$g(\mathrm{})=g^{}(\mathrm{})=g^{\prime \prime }(\mathrm{})=0.$$
(21)
Even though we cannot find the eigenvalues of $`L(q)`$ for arbitrary $`q`$ analytically, for long wavelength disturbances we can use perturbation theory to get the leading order (in $`q^2`$) terms. This requires finding the eigenfunctions of $`L_0`$ and its adjoint, $`L_0^{}`$. Taking the second derivative of (12) we obtain

$$L_0h_0^{}=0,$$
(22)
so that $`g_0=h_0^{}`$ is an eigenfunction of $`L_0`$ with eigenvalue $`\beta _0^0=0`$. The adjoint operator is defined via

$`L_0^{}f`$ $`=`$ $`[12h_0+(\alpha +3h_0^2)h_0^{\prime \prime \prime }]f^{}`$ (23)
$``$ $`[(\alpha h_0+h_0^3)f^{}]^{\prime \prime \prime },`$ (24)
so its respective eigenfunction is just a constant, say, $`f_0=1`$. In fact, these are generic results with deep physical meaning. Identical relations between the asymptotic state and the leading eigenfunctions were obtained, e.g., for gravity driven films using the precursor model . The relation for $`g_0`$ is due to the fact that equations for the asymptotic state are translationally invariant in the direction of the flow (this reflects an arbitrary choice in the position of the contact line), while the relation for $`f_0`$ is the consequence of the divergence form of (5), which reflects mass conservation.

According to the perturbation theory the leading eigenvalue has the following $`q`$-dependence:

$$\beta _0(q)=\beta _0^0+q^2\frac{_0^{\mathrm{}}f_0^{}L_1g_0𝑑x}{_0^{\mathrm{}}f_0^{}g_0𝑑x}+O(q^4).$$
(25)
Using (13) this can be reduced to

$$\beta _0(q)=q^2_0^{\mathrm{}}h_0(h_01)𝑑x+O(q^4).$$
(26)
As this is the largest eigenvalue, its sign determines the stability of the asymptotic state. It is easy to see that, if the asymptotic profile were monotonic, $`0<h_0<1`$, the integral would be strictly negative and the system would be linearly stable with respect to long wavelength disturbances. However, the presence of a capillary ridge makes the integral positive, showing that the increased mobility of the ridge provides the mechanism for the long wavelength instability in the thermally driven case. This mechanism has been originally conjectured by Kataoka and Troian based on the energy analysis of the precursor model . Equation (26) gives an explicit condition on the shape of the capillary ridge which determines the onset of instability, and echoes a similar result obtained for the case of gravity-driven flows .

Substituting $`g(x,t)=h_0^{}(x)\mathrm{exp}[\beta _0(q)t]`$ into (15) we notice that for small disturbances the right hand side represents the first two terms of the Taylor expansion of a distorted asymptotic state $`h_0(x+\xi +ut)`$, where

$$\xi (y,t)=ϵe^{iqy+\beta _0(q)t}.$$
(27)
is the deviation of the contact line from the mean. In fact, the marginal translational mode $`g_0`$ is not the only eigenfunction of $`L_0`$. There is an infinite discrete spectrum of eigenvalues $`\beta _n`$ and eigenfunctions $`g_n`$. Therefore, in the presence of an arbitrary disturbance (27) should be replaced with

$$\xi (y,t)=\frac{1}{c}\underset{n}{}g_n(0)_{\mathrm{}}^{\mathrm{}}ϵ_n(q)e^{iqy+\beta _n(q)t}𝑑q.$$
(28)
As the asymptotic state is unstable, the amplitude of the distortion will grow exponentially in time and eventually the contact line will form equally spaced “fingers”. In order to calculate the wavenumber of the pattern we numerically compute the eigenfunctions and eigenvalues of the evolution operator $`L(q)`$. This is most easily accomplished by discretizing the eigenvalue equation on a spatially nonuniform mesh to properly resolve the rapid change in the solutions near the contact line. As both the asymptotic state $`h_0`$ and the eigenfunctions $`g_n`$ exponentially flatten for large $`x`$, a truncated domain can be used, so that the boundary conditions (21) are imposed at finite distance $`l_x`$ away from the contact line (we used $`l_x=80`$ in most of the calculations). In order to compute the asymptotic state we used a shooting method on a nonuniform mesh with roughly $`10^4`$ points. As $`h=1`$ is an unstable fixed point of (13), numerical integration (we used the $`4^{\mathrm{t}h}`$ order Runge-Kutta method) cannot be performed beyond $`l_x30`$ using double precision arithmetics. To determine the solution for a longer interval we used a method  in which the numerical solution is extended using the analytical solution of (13) for $`h1`$:

$$h(x)=1+[a\mathrm{cos}(\sqrt{3}\kappa x)+b\mathrm{sin}(\sqrt{3}\kappa x)]e^{\kappa x}+\mathrm{},$$
(29)
where $`\kappa =(1+\alpha )^{1/3}/2`$ and $`a`$, $`b`$ are constants determined by a least squares fit. The eigenvalues and eigenfunctions where then computed using a built-in function in Matlab on a 1025-point mesh (finer resolutions did not change the eigenvalues by more than about 5%).

The results of our calculations for a typical choice of parameters are presented in Fig. 5. The fastest growing disturbance is found to have a transverse wavenumber $`q_{max}`$ lying between $`0.29`$ and $`0.34`$. This wavenumber decreases with the thickness of the capillary ridge $`h_{max}`$ and gives the characteristic wavelength of the fingering pattern $`\lambda _{max}=2\pi /q_{max}`$, which ranges between $`18.5`$ and $`21.7`$, in excellent agreement with the predictions of both the precursor model  and the experiments  which found the dimensionless wavelength to be between $`18`$ and $`22`$. The growth rate of the most unstable mode increases with $`h_{max}`$, varying between $`0.12`$ and $`0.17`$ for the range of contact angles considered here. Brzoska et al. have obtained an experimental value of about $`0.15`$, which is also consistent with the theory. As the thickness of the ridge has not been determined in experiments, it is impossible to make a more direct comparison, but nevertheless these results can be used to establish the ranges of parameters relevant for experimental conditions.

## IV Linear Stability and Transient Dynamics

We have determined earlier that the capillary ridge is present for any reasonable choice of parameters, so the asymptotic state of a thermally driven film is always linearly unstable and the contact line instability will inevitably set in as the asymptotic state is approached. However, as a quick comparison of (17) and (23) shows, the evolution operator $`L(q)`$ is nonnormal, and we have to consider the possibility that the transient dynamics associated with nonnormality could change the predictions of the linear stability theory concerning the growth rates of disturbances. In fact, transient effects could be quite significant. For instance, turbulence in channel flows arises at values of the Reynolds number well below the critical one predicted by the linear stability analysis . Both in the driven films and in channel flows nonnormality arises due to a significant mean flow, so it is natural to expect that transient behavior could be important for driven liquid films as well. Indeed, a disagreement between theoretical and experimental predictions of the most unstable wavelength in gravity driven films at low angles of inclination has been attributed to transient dynamics .

Linear stability analysis presented in the previous section is based on the assumption that the nonlinear terms are negligible at all times. If the disturbances were small and their dynamics were governed by a normal evolution operator, this assumption would have been well justified. For instance, when the system is stable, the eigenvalues predict both the short term and the long term dynamics. The situation can change dramatically when the evolution operator becomes nonnormal, as the eigenvalues become poor predictors of the short term dynamics. Inclusion of the nonlinear terms in (16) has two major consequences. First of all, nonlinear terms couple the dynamics of modes with different transverse wavenumbers . Second, nonlinear terms produce deviations from the asymptotic state which can be transiently amplified due to the nonnormality of the linearized evolution operator . A combination of these two effects can lead to a nonlinear instability which can compete with the linear instability.

There are two scenarios which could invalidate the results of linear stability analysis. In the simplest (purely linear) scenario , an initial disturbance with transverse wavenumber $`q_0`$ and magnitude $`g(x,0)=\xi `$ could be transiently amplified by a factor $`\gamma (q_0)`$ to produce a disturbance with magnitude $`\gamma \xi =O(1)`$. If this amplification occurs on a time scale shorter than $`1/\beta _0(q_{max})`$ and $`\xi \gamma ^1`$, the transient effects will dominate and a distortion of the contact line with the wavenumber $`q_0`$ rather than $`q_{max}`$ will result.

In a more complicated scenario , an initial disturbance is transiently amplified by the linear part of the evolution operator, while the nonlinear terms produce secondary disturbances which are further transiently amplified. This could lead to a positive feedback loop bootstrapping a nonlinear instability, provided the secondary disturbances contain the wavenumber of the initial disturbance and have the magnitude which is at least as large. It is easy to check that the nonlinear evolution operator will only contain terms quadratic, cubic, and quartic in $`g(x,t)`$ in addition to the linear terms kept in (16). Only the cubic terms will contain the original wavenumber, so the secondary disturbances produced by the quadratic and quartic terms will not be further transiently amplified. (Initial disturbances with $`q_0=0`$ represent the only exception, but they do not lead to distortion of the contact line.) An initial disturbance of magnitude $`\xi `$, which is transiently amplified by a factor $`\gamma `$, will produce a secondary disturbance of magnitude $`O((\gamma \xi )^3)`$ via the cubic nonlinearities. The secondary disturbance will exceed the primary disturbance only when $`(\gamma \xi )^3\xi `$, so that a self-sustaining nonlinear instability becomes possible for $`\xi \gamma ^{3/2}`$.

The maximal transient amplification in the stable band can be defined in a conventional manner :

$`\gamma (q,t)=\underset{g(x,0)}{sup}{\displaystyle \frac{g(x,t)}{g(x,0)}}=e^{Lt}.`$ (30)
In the unstable band the fast transient amplification will be followed by the slower exponential growth. Factoring out the exponential growth one obtains the following upper bound:

$`\gamma (q,t)=\underset{g(x,0)}{sup}{\displaystyle \frac{g(x,t)}{g(x,0)}}e^{\beta _0t}=e^{(L\beta _0)t}.`$ (31)
Numerical calculations show that in the unstable band $`\gamma (q,t)`$ is a monotonically increasing function of time, so the maximum is reached for $`t\mathrm{}`$. The time dependence for several values of the transverse wavenumber is shown in Fig. 6. As $`L\beta _0`$ and $`L`$ have the same eigenfunctions, we can easily calculate the matrix elements of the operator $`U(t)=\mathrm{exp}[(L\beta _0)t]`$:

$$U_{mn}(t)=_0^{\mathrm{}}f_m^{}e^{(L\beta _0)t}g_n𝑑x=e^{(\beta _n\beta _0)t}\delta _{mn}.$$
(32)
As only the element with $`m=n=0`$ survives for large times, the maximal transient amplification is achieved for the “optimal” initial disturbances equal to multiples of $`f_0`$. The evolution amplifies these disturbances and transforms them into multiples of the leading eigenfunction $`g_0`$. Exponential growth with rates predicted by the linear stability analysis sets in rather quickly as the time scales of the exponential and transient growth are of the same order of magnitude (roughly between $`2`$ and $`9`$ nondimensional units, depending on the wavenumber).

The ultimate test of the importance of nonlinear terms and transient dynamics is provided by a direct calculation of the transient amplification factor

$$\gamma (q)=\underset{t}{\mathrm{max}}\gamma (q,t).$$
(33)
The numerical results are presented in Fig. 7. For a fixed contact angle, $`\gamma (q)`$ is the largest for zero wavenumber disturbances (which do not lead to distortion of the contact line) and quickly decreases with increasing $`q`$. The maximum transient amplification increases with the contact angle and can become quite significant for typical experimental parameters ($`\gamma (0)70`$ for $`c=1`$, $`\gamma (0)100`$ for $`c=1.95`$). However, the effective value for a finite system will likely be in the range of a few tens. The strong transient amplification at $`q=0`$ can be easily traced to a very close alignment of a large group of eigenfunctions $`g_2`$ through $`g_{21}`$. As Fig. 8 shows, their shapes are extremely similar. It is, therefore, appropriate to associate the transient behavior with a whole group of stable eigenfunctions. The size of the group increases with $`c`$, increasing the degree of nonnormality. In contrast, for channel flows apparently only a couple of near-marginal eigenfunctions become closely aligned .

The transient amplification at small $`q`$ was found to depend rather sensitively on the length of the domain used in the numerical calculations. For instance, the values of $`\gamma `$ obtained for $`l_x=30`$ were generally about a half of those obtained for $`l_x=80`$. At about $`l_x=80`$ the dependence levelled off, and further increase in $`l_x`$ resulted in large fluctuations in the eigenvalues due to numerical inaccuracies resulting from strong nonnormality of the matrix produced by discretizing $`L(q)`$. This sensitivity can be explained by the shape of the adjoint eigenfunctions. While the eigenfunctions $`g_n(x)`$ decay exponentially for large $`x`$, their adjoints $`f_n(x)`$ grow exponentially (see Fig. 9). As a result, small inaccuracies in the boundary conditions at the tail of the film introduced by truncating the computational domain can have a significant effect on the adjoint eigenfunctions, thus affecting the transient amplification, which depends on both $`g_n(x)`$ and $`f_n(x)`$.

The minimal noise level required to trigger the nonlinear instability according to the second scenario is about 0.1% of the total film thickness for $`\gamma =100`$ (the first scenario requires an even higher level of noise). A more accurate calculation of the noise threshold is likely to raise this level much higher (to maybe a few percent) by taking into account the fact that (31) determines the maximal amplification achieved for a specially chosen initial condition, while the secondary disturbances will generically be amplified less strongly. In typical experimental conditions the noise is likely to be substantially smaller than one percent of the film thickness. The existing experimental data  agrees with the predictions of the linear theory rather well, supporting the conclusion that the transient effects are relatively weak and, therefore, the modal linear stability analysis of the previous section accurately describes the dynamics.

## V Feedback control of the contact line instability

Now that we understand the limits of the linear stability analysis let us consider the control problem. Can the contact line instability be suppressed, or alternatively, can a pattern with a desired wavelength be imposed by applying feedback? In principle, the answer seems to be clear, as feasibility of feedback control of several other types of instability (buoyancy , thermocapillary , evaporative ) in liquid layers has been demonstrated. Although one might hope that the developed control methods could be adapted for suppressing the contact line instability, in reality the spreading films turn out to be dramatically different.

The existing control methods have been developed for stabilizing flat films with no mean flow, i.e., steady, uniform target states. The evolution operators describing the dynamics of small disturbances about such states have eigenfunctions whose horizontal dependence is given by Fourier modes, which are normal to each other for any realistic boundary conditions (in principle, one can design the boundary conditions such that the eigenfunctions will not be normal even in this case ). Not only does it mean that all horizontal modes become uncoupled, so any mode can be controlled independently of the others, there is no nonnormality, so there are no transients and the linear stability analysis is unconditionally valid.

The target state $`h_0(x+ut)`$ in the present problem is nonuniform in the direction of the flow. As a result, the differential operator $`L`$ does not fully diagonalize and the control problem becomes vastly more complicated. Not only are all the modes in the system coupled, feedback applied to suppress one mode generally affects all other modes, so an infinite-dimensional problem has to be considered from the outset. However, even if these problems are resolved and a feedback making the dynamics asymptotically stable is found, there is no guarantee that the transient effects will not invalidate the whole analysis.

Let us repeat the linear stability analysis of the spreading film, but now in the presence of feedback. First we make use of the simplification afforded by the uniformity of the target state in the transverse direction, which allows us to partially diagonalize the evolution operator. For the moment we restrict our attention to monochromatic disturbances $`ϵg(x+ut,t)\mathrm{exp}(iqy)`$. Since the flow is driven by the gradient in the temperature (and hence surface tension), the stability of the flow is most easily altered by varying the temperature field behind the contact line. Suppose we modify the temperature profile by adding a perturbation

$$\mathrm{\Delta }\theta (x,y,t)=ϵ\tau (_\theta \sigma )^1s(t)w(x+ut)e^{iqy},$$
(34)
where the transverse wavelength $`q`$ is the same as that of the disturbance and $`s(t)`$ and $`w(x)`$ are some functions determining the temporal and spatial profile of the perturbation, which will be determined later. Consequently, (7) and hence (16) will be modified to account for the variation in the surface tension transversely to, as well as along, the direction of the flow. At order $`ϵ`$ instead of (16) we obtain

$`_tg`$ $`=`$ $`L_0g+s[N_1(h_0)w+N_2(h_0)w^{}]^{}`$ (35)
$`+`$ $`q^2[L_1gsN_2(h_0)w]+q^4L_2g,`$ (36)
where the influence functions

$`N_1(h_0)`$ $`=`$ $`{\displaystyle \frac{2}{3}}(h_0^2h_0),`$ (37)
$`N_2(h_0)`$ $`=`$ $`h_0^2+{\displaystyle \frac{2}{3}}(\alpha h_0+h_0^3)h_0^{\prime \prime }`$ (38)
determine the effect of the imposed thermal perturbation.

To get a sense of the dynamics of different modes in the presence of feedback, we expand the disturbance in the basis formed by the eigenfunctions of $`L_0`$,

$$g(x,t)=\underset{m}{}G_m(t)g_m(x),$$
(39)
and make the strength of the applied perturbation proportional to the magnitude of the distortion of the contact line (with a proportionality constant $`k`$, called the gain, to be determined later),

$$s(t)=k\frac{g(0,t)}{c}=\frac{k}{c}\underset{m}{}G_m(t)g_m(0).$$
(40)
Multiplying (35) by $`f_n^{}`$ and integrating from $`0`$ to $`\mathrm{}`$ we obtain an infinite system of ODEs describing the dynamics of individual modes:

$$\dot{G}_n=\beta _n^0G_n+\underset{m}{}(A_{nm}+q^2B_{nm}+q^4C_{nm})G_m,$$
(41)
where (assuming that all adjoint eigenfunctions are normalized such that $`_0^{\mathrm{}}f_n^{}g_m𝑑x=\delta _{nm}`$)

$`A_{nm}`$ $`=`$ $`k{\displaystyle \frac{g_m(0)}{c}}{\displaystyle _0^{\mathrm{}}}f_n^{}[N_1(h_0)w+N_2(h_0)w^{}]^{}𝑑x,`$ (42)
$`B_{nm}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}f_n^{}L_1g_m𝑑xk{\displaystyle \frac{g_m(0)}{c}}{\displaystyle _0^{\mathrm{}}}f_n^{}N_2(h_0)w𝑑x,`$ (43)
$`C_{nm}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}f_n^{}L_2g_m𝑑x.`$ (44)
As Fig. 5 shows, the uncontrolled system possesses a single unstable eigenvalue well separated from the rest of the spectrum. One can, therefore, expect that the dynamics of small disturbances should be well described by a single mode truncation of (41) in the absence of feedback. The same is not generally true when the feedback is applied. As all modes are coupled, the feedback designed to suppress the leading mode will always affect, and can potentially destabilize, the subleading modes. As our numerical calculations show, such destabilization does indeed occur, unless the spatial profile $`w(x)`$ of the thermal perturbation is chosen carefully to avoid this. Had the evolution operator been normal, we could have always chosen the feedback in such a way that different modes became uncoupled, so only the stability of a few independent unstable modes had to be considered. The problem of controlling the contact line instability turns out to be quite delicate in comparison.

Assuming that $`w(x)`$ is chosen such that all subleading modes remain stable, we can truncate the system (41) by discarding all modes except the leading one. The stability of the single-mode truncation

$$\dot{G}_0=(A_{00}+q^2B_{00}+q^4C_{00})G_0$$
(45)
is a necessary condition for the stability of the full system (41). The truncated system is stable when the following three conditions are satisfied:

$`A_{00}`$ $`=`$ $`kw^{}(\mathrm{})0,`$ (46)
$`B_{00}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}h_0(h_01)𝑑xk{\displaystyle _0^{\mathrm{}}}N_2(h_0)w𝑑x0,`$ (47)
$`C_{00}`$ $`=`$ $`{\displaystyle \frac{1+2\alpha }{4}}0.`$ (48)
The first condition is satisfied for any feedback localized near the contact line, because in this case $`A_{00}=0`$, while the second condition can always be satisfied with the proper choice of the gain $`k`$. The third condition is clearly satisfied as well. Since (41) is valid for all values of $`q`$ for which the governing equation (10) is valid and the choice of $`k`$ in (46) is independent of $`q`$, we can immediately generalize to non-monochromatic disturbances by integrating over all $`q`$, such that the feedback will be given by

$$\mathrm{\Delta }\theta (x,y,t)=k\tau (_\theta \sigma )^1w(x+ut)\xi (y,t),$$
(49)
where $`\xi (y,t)`$ is the instantaneous deviation of the contact line from its mean position.

How can we choose the heating profile in the flow direction which will stabilize the unstable mode without destabilizing any of the initially stable modes? As Fig. 8 shows, the leading mode is localized under the capillary ridge. Therefore, to suppress the instability one has to apply a perturbation which will be similarly localized to within the region occupied by the ridge. The inspection of matrix elements (42) shows that such a localized perturbation will affect all modes whose adjoint eigenfunctions are not small in that region. According to Fig. 9 both $`f_0`$ and $`f_1`$ are of the same order of magnitude there, while the other subleading modes are several orders of magnitude smaller. As a result, only the stability of the leading and the first subleading mode may change in response to feedback. The numerically computed spectra of the controlled system support this conclusion.

The stability of the two leading modes is determined by the signs and magnitudes of the matrix elements $`A_{nm}`$ and $`B_{nm}`$ with $`n,m=0,1`$. A few general comments about these matrix elements can be made based on the structure of the eigenfunctions. As Fig. 10 shows, $`f_0`$ and $`f_1`$ are both nearly constant under the ridge and have opposite signs, while $`g_0(0)`$ and $`g_1(0)`$ have the same sign. As a result, a decrease in $`B_{00}`$ is necessarily accompanied by a commensurate increase in $`B_{11}`$. Furthermore, since $`f_0`$ is constant, we have $`A_{0n}=0`$ for any $`n`$, so at $`q=0`$ the eigenvalues of the controlled system are $`\beta _0^{}(0)=0`$ and $`\beta _1^{}(0)\beta _1(0)+A_{11}`$. The matrix element $`A_{11}`$ is generally nonzero and changes in response to the strength of the applied feedback. An inappropriate choice of the profile $`w(x)`$ or the gain constant $`k`$ can make $`A_{11}`$ large enough to cause destabilization at long wavelengths. However, even when $`w(x)`$ is chosen such that $`A_{11}`$ is negative (but small), destabilization of the $`n=1`$ mode at short wavelenghts can occur due to the increase in $`B_{11}`$, if the feedback gain is too large.

Additional insights can be gained by considering the effect of feedback from the physical point of view. The action of feedback in the direction transverse to the flow is described by the influence function $`N_2(h_0)`$. The first and second term of this function describe the motion of the liquid under the action of, respectively, surface forces and pressure produced by the local gradients in the surface tension. For a convex region of the film, where $`h_0^{\prime \prime }<0`$, these two effects will compete with each other. For instance, a local maximum of surface tension will induce the flow along the surface toward that location and the flow in the interior of the liquid layer away from that location. The first effect will dominate for low capillary ridges (small contact angles), the second one for high capillary ridges (large contact angles).

As the instability is caused by the increased mobility of the capillary ridge, one might envision enforcing control by changing the thickness of the film. The local thickness of the capillary ridge could, in principle, be modified by locally heating or cooling it to redistribute the liquid in such a way as to decrease the thickness, and hence the mobility, where we need to slow down the motion of the contact line and increase the thickness and mobility, where we need to speed it up to compensate for the deviation. This can be achieved by choosing $`w(x)`$ which is localized under the ridge and does not change sign. For instance, one could pick a Gaussian profile representing the effect of thermal spreading in the solid substrate

$$w(x)=\mathrm{exp}\left[\frac{(xx_0)^2}{2\mathrm{\Delta }x^2}\right],$$
(50)
where $`x_0`$ and $`\mathrm{\Delta }x`$ are chosen such that $`w(x)`$ is centered under the capillary ridge and has a comparable width. A sample profile is shown in Fig. 11(a) for a special choice of parameters. The numerically computed spectrum, Fig. 11(b), demonstrates that the stabilization can indeed be achieved by this method for driven films with a relatively low capillary ridge.

This simple approach, however, does not always achieve the desired result. In fact, it only succeeds when the largest growth rate $`\beta _0(q_{max})`$ in the uncontrolled system is small compared to $`|\beta _1(0)|`$, i.e., when the feedback required to stabilize the $`n=0`$ mode is too weak to destabilize the $`n=1`$ mode. One can already notice the sign of approaching trouble by looking at Fig. 11(b): the eigenvalue of the $`n=1`$ mode at $`q=0`$ starts to creep upward due to the increase in the matrix element $`A_{11}`$. If this approach is used to stabilize a flow with $`\beta _0(q_{max})`$ comparable to $`|\beta _1(0)|`$, the feedback required to suppress the $`n=0`$ mode becomes strong enough to destabilize the $`n=1`$ mode at low wavenumbers. The numerically computed spectra for higher capillary ridges ($`h_{max}>1.7`$) show that the low wavenumbers are destabilized before the high wavenumbers are stabilized for any choice of the width $`\mathrm{\Delta }x`$ of the thermal perturbation.

We are thus forced to look for an alternative solution. Changing the overall thickness of the capillary ridge, even locally, is a rather ineffective procedure, especially for small $`q`$, as the liquid has to be redistributed over large distances in the transverse direction. One could instead apply a local force to the ridge, redistributing the liquid between its front and back. For instance, by heating the front of the ridge and cooling its back one creates the pressure gradient enhancing the flow. Reversing the sign of the applied perturbation impedes the flow, directly affecting the propagation velocity. The corresponding thermal perturbation should have a profile $`w(x)`$ which changes the sign near the highest point of the ridge. For instance, if one chooses the profile to be anti-symmetric

$$w(x)=(xx_0)\mathrm{exp}\left[\frac{(xx_0)^2}{2\mathrm{\Delta }x^2}\right],$$
(51)
the matrix element $`A_{11}`$ can be made large and negative, so the $`n=0`$ mode can be stabilized without destabilizing the $`n=1`$ mode. Fig. 12(a) shows that such a thermal perturbation with the position $`x_0`$ and width $`\mathrm{\Delta }x`$ tuned to be roughly the same as those of the capillary ridge, should have the largest effect on $`B_{00}`$. Indeed, the influence function $`N_2(h_0)`$, dominated for high capillary ridges by its second term, also changes sign near the highest point of the capillary ridge. The numerically computed spectrum, Fig. 12(b) shows that one can again successfully suppress the instability.

Another alternative is to exploit the narrow concave region of the film near the contact line. One can again use the Gaussian thermal profile (50) centered at the contact line (see Fig. 13(a)). Heating this region and thus lowering the surface tension one produces gradients in both the pressure and surface tension, which induce the secondary flow away from the contact line. Respectively, cooling this region draws the liquid towards it, providing a direct way to locally control the propagation speed of the film. Making the amount of heating or cooling proportional to the displacement of the contact line again allows one to suppress the contact line instability. The numerically computed spectrum of the system with and without feedback is shown in Fig. 13(b).

Finally, let us look at the transient amplification of disturbances in the presence of feedback. It is unclear a priori what effect the control would have on the transient dynamics. Numerical calculations show that, depending on the choice of the thermal profile $`w(x)`$, feedback can either increase or decrease the degree of nonnormality in the system. This observation is consistent with the theory. As the matrix elements (42) show, feedback affects different modes in a different way, and small changes in the eigenvalues can have a large effect on the transient amplification. For instance, thermal perturbations with an anti-symmetric profile, such as (51), can decrease the transient amplification at small wavenumbers by almost an order of magnitude (see Fig. 14). Generally, direct control of the propagation velocity via longitudinal surface tension gradients leads to a decrease in the transient amplification, while indirect control via transverse surface tension gradients affecting the mobility increases the transient amplification for long wavelength disturbances, which is also consistent with naive expectations.

We can thus conclude that an appropriately chosen feedback can make the dynamics asymptotically stable without increasing the transient amplification of disturbances, such that the contact line instability is suppressed in the presence of noise characteristic of typical experimental conditions. Moreover, the feedback is capable of reducing the transient amplification as well, so we can expect that feedback control can be effective in suppressing the contact line instability even when it is caused by the nonlinear effects. Finally, once the instability is suppressed, selective patterning can be achieved by removing feedback and/or introducing additional forcing at a wavenumber corresponding to a desired pattern.

The proposed control algorithm has been verified experimentally. Although no systematic investigation of different thermal perturbation profiles has been attempted so far, the proof-of-principle experiments have shown that by heating the advanced regions of the film and cooling the retarded regions one can completely suppress the contact line instability. It has also been shown that transversely modulated thermal perturbations applied near the contact line can be used to achieve selective patterning, producing perfectly periodic patterns of rivulets. The details of the experiments will be presented in a separate publication (N. Garnier, R. O. Grigoriev, and M. F. Schatz, “Optical manipulation of microscale fluid flow”, in preparation).

## VI Discussion

To summarize our results, we have determined that the evolution operator governing the dynamics of spontaneous disturbances for thermally driven films, both with and without feedback, is significantly nonnormal and can transiently amplify those disturbances. The strongest transient amplification occurs for the zero wavenumber which does not lead to contact line instability. However, even for nonzero wavenumbers transient amplification is unlikely to produce an instability for levels of noise characteristic of typical experimental conditions. Therefore, linear stability analysis accurately describes both the controlled and the uncontrolled dynamics. In contrast, for gravity driven films at small inclination angles the transient amplification could be much stronger, with the maximum achieved at a nonzero wavenumber , providing an alternative mechanism for instability.

We have also shown that the contact line instability in thermally driven films can be effectively suppressed by locally heating or cooling the liquid behind the contact line. Such thermal perturbation can be easily imposed experimentally with sufficient spatial and temporal resolution by radiatively heating the substrate . This approach offers significant advantages in controlling the dynamics of microflows compared to the one based on chemical patterning of the substrate . First of all, no preparation of the substrate is needed, while the patterns can be dynamically reconfigured, offering potential for a significant increase in flexibility. Second, feedback control can be used to achieve extremely small feature size, if high intensity radiation is used to drive the flow on a thin substrate with moderate thermal conductivity , opening up new prospects for microfluidics and microfabrication applications. Finally, feedback control provides a unique opportunity for studying the dynamics of subdominant modes and even unstable states of the system. For instance, it can be used to experimentally measure the growth (or decay) rates of monochromatic disturbances with wavenumbers $`q_0q_{max}`$ by using a wavenumber dependent gain to suppress the disturbances which would otherwise obscure the dynamics of the mode of interest.

We also expect that feedback control can be equally effective in suppressing nonlinear instabilities such as those occurring in gravity driven spreading at small inclination angles. Indeed, we have seen that the profile of the thermal perturbations can be tuned to decrease the degree of nonnormality. Therefore, by suppressing the transient growth feedback can also quench the bootstrapping mechanism leading to a nonlinear instability. However, because small disturbances at the contact line could be transiently amplified to produce $`O(1)`$ changes in the thickness of the capillary ridge , it is possible that the control algorithm will have to use direct measurements of the thickness rather than the much easier to monitor position of the contact line.


text:
# Tracing and Explaining Execution of CLP(FD) Programs

## 1 Introduction

In this paper, we present new ideas in the area of tracing CLP(FD) programs. New trends such as tracing global constraints and explaining solver actions are investigated. Since these issues are still in their development stages and not systematically used, we have not tried to create a universal trace model. Instead, we focus on CLP(FD) and in particular the CLP(FD) solver of SICStus Prolog.

The idea of tracing constraint programs originates from Meier . Also, the search-tree visualization tool in CHIP  uses a trace log to be able to display propagation events. More current work is carried out by Langevine et al. who presented a CLP(FD) trace model in . Their work is naturally influenced by the work of Ducassé  on Opium.

Previous work such as the above mainly focuses on providing information about control of execution and domain modification. Of course, this information is needed, but if nothing else is provided, important information about the execution is lost. We have looked at the execution of a CLP(FD) program from the following perspectives:

1. Control of Execution; information about posting new constraints, waking up delayed constraints, constraint entailment/failure.
2. Domain Modification; information about variable domain narrowing (pruning).
3. Declarative Aspect; information about why certain actions are taken by the solver and in which context they occur.
4. Procedural Aspect; information about filtering algorithms and their entities.

We build the trace of events as a nested block structure. This means that most of the trace events appear inside a block surrounded by a starting event and an ending event. One example of this is the posting of a new constraint, which gives rise to a block surrounded by the trace events `begin_new_ctr` and `end_new_ctr`. Inside this block, trace events describing the posting of demon(s) for the new constraint and initial prunings are generated. Building the trace like this has several advantages. First of all, it makes it easier to get a structured overview of the trace. Second, similar actions may be grouped together inside the same block. Furthermore, extra information about the action(s) taking place may be added to the opening and closing events as well as through added events inside the block.

We identify variables and other entities by their context, i.e. their positions in the constraints, and not only by some identifier. This is useful when one needs to distinguish a specific occurrence among several of the same variable. In addition to that, it makes it possible to also clearly identify entities such as integers, atoms and compound terms.

Also, corresponding to point 3 above, we provide explanations  for why certain actions are taken by the solver. These explanations are incorporated in the trace by specific explanation trace events. Actions for which explanations are provided include domain narrowing, failure, demon enqueuing and internal constraint posting. We extend previous work by also providing compact explanations for global constraints. By doing this, we obtain a stronger link between the consistency rules of the global constraints and the trace.

Moreover, corresponding to point 4 above, we provide trace events containing detailed information about the filtering algorithms of global constraints. This concerns identifying and presenting the different entities inside the algorithms. This too increases the link between the global constraints and the trace. However, this is more related to the procedural aspect.

An example of a simple CLP(FD) program in SICStus Prolog is shown in Ex. 1. It will be used throughout the paper. Line `(2)` initializes the domains of the included variables. Lines `(3)` and `(4)` post the well-known constraints `all_different/1`222This is the naive version of all\_different, achieving the same pruning as a set of pairwise inequality constraints. and `element/3`. Line `(5)` calls SICStus’s builtin search predicate with a `leftmost`333Leftmost variable, smallest value first. search strategy.

In the following, Sect. 2 presents a CLP(FD) execution model. This is needed in order to be able to present the trace structure in its relevant context. Sect. 3 presents the different trace events. First, some design choices concerning how to represent different entities such as constraints and variables are introduced. Following that, the trace events are presented in chunks; trace events corresponding to the same perspective, as noted above, are presented together. After that, a short description of the implementation made in SICStus Prolog is given. Sect. 4 discusses possible areas of usage. Finally, Sect. 5 concludes the paper.

## 2 Execution Model

This section presents an execution model of a CLP(FD) kernel; the different data structures it contains, and how these interact. By doing this, we can later present the trace structure in its relevant context. The model is influenced by the constraint solvers of SICStus Prolog  and Choco .

#### 2.0.1 Data Structures.

The data structures contain lists and queues for storing information regarding constraints, variables, demons and propagation events.

Each constraint is associated with a list of domain variables, `VList`, and a list of demons, `DList`. `VList` contains the variables that occur in the constraint. `DList` contains the demons that are responsible for removing inconsistent values from the variables in `VList`.

Each variable is associated with a domain, `Domain`, containing the values the variable may take. These domains are narrowed by propagation events specified below. Furthermore, each variable is associated with a list of constraints, `CList`, which it is involved in. As soon as the domain of a variable changes, one or more of these constraints may wake up and activate their demons.

Each demon contains a wake-up condition. When the wake-up condition is fulfilled, the demon’s filtering algorithm is activated. This filtering algorithm may exclude some values from the domains of some variables. A filtering algorithm consists of a set of methods, $`\{𝙼_\mathrm{𝟷},\mathrm{},𝙼_𝚗\}`$, where each $`𝙼_𝚒`$ corresponds to some consistency rule of the constraint that the demon belongs to. If some method, $`𝙼_𝚒`$, notices that the domains of some variables are not consistent with the constraint, $`𝙼_𝚒`$ generates propagation events which specify domain narrowings on these variables.

The propagation events have the form `X in_set S`, which constrains the variable `X` to be a member of the set `S` of integers.

Two global lists are used for storing current constraints and variables in the system; `ConstraintList` and `VariableList`. Furthermore, two global queues are needed to manage the actions taken by the kernel. `ReadyQueue` contains the demons that are about to wake up. `PropagationQueue` contains the propagation events that have been created by active demons.

Fig. 1 shows the different data structures and how they are connected to each other. The state is taken from Ex. 1, just before the execution of line `(4)`. `C1` and `C2` are the `all_different/1` and `element/3` constraints respectively. Detailed information about variables `X` and `Y` as well as about constraint `C2` is shown in three boxes.

#### 2.0.2 Services.

Our CLP(FD) kernel provides the following services:

* connectCtr(C), adds the constraint `C` to the kernel. This includes adding `C` to `ConstraintList`, adding the variables in `C` to `VariableList`, `adding` `C` to the `CList` of all its variables, creating the demons444Our trace structure supports the use of more than one demon for each constraint. However, since the examples in this paper are created using SICStus Prolog, and it currently only associates one demon with each constraint, more than one demon per constraint will never be shown. associated with `C` and running each demon’s filtering algorithm for initial propagation.
* disconnectCtr(C), removes the constraint `C` from the kernel. This is done when the constraint is entailed, i.e. when `C` is true no matter what values its variables take.
* enqueueCtr(C,E), puts last in `ReadyQueue` the demons associated with the constraint `C` that are to wake up due to the propagation event `E`. If any of these demons are already in `ReadyQueue`, it will not be enqueued again.
* dispatchCtr(), dequeues the demon that is first in `ReadyQueue` and activates its filtering algorithm.
* enqueueEvent(E), puts `E` last in `PropagationQueue`.
* dispatchEvent(), dequeues the event `E` that is first in `PropagationQueue`, performs the narrowing on some domain variable `X` it specifies, and calls `enqueueCtr(C,E)` for all constraints in `X`’s constraint list.

If the functions `dispatchEvent()` and `dispatchCtr()` can both be executed, the former will have higher priority.

#### 2.0.3 The Active System: An Example.

The easiest way to see how the above described data structures interact is by looking at an example. Assume that the state of the kernel is as in Fig. 1. There is no propagation event to be performed, but the demon `D2` of constraint `C2`, the `element/3` constraint, should be woken up. The kernel calls the function `dispatchCtr()`, which removes $`\mathrm{𝙳𝟸}`$ from `ReadyQueue` and activates its filtering algorithm. A first method, `M1`, generates a propagation event, `E1`, which restricts `X` to take a value inside the range `1..4`. Equally, a second method, `M2`, generates a propagation event, `E2`, which restricts `Y` to take a value inside the set `{2,3,4}`. Also, a third method, `M3`, generates a propagation event, `E3`, which restricts `X` from taking the value `4`555Since the value 8 on position 4 in the list is not in the domain of Y.. When the filtering algorithm of `D2` has finished, the kernel notices that `PropagationQueue` has changed to contain the propagation events `E1`, `E2` and `E3`. The kernel calls the function `dispatchEvent()` three times which does the following each time:

* The first propagation event in `PropagationQueue` is dequeued.
* The narrowing of the variable specified by the propagation event is performed by the kernel.
* The demons of the constraints `C1` and `C2`, that are supposed to wake up due to the narrowing performed, are put in `ReadyQueue` by calling the function `enqueueCtr()` twice.

When `PropagationQueue` is empty, the kernel starts over by calling the function `dispatchCtr()` again. This is repeated until both queues are empty. When this is the case, the control of execution is given back to the Prolog level. In our case, this means that labeling takes place which will fix a variable and therefore trigger more propagation.

## 3 Trace Structure

In this section, the trace structure and the information it contains is presented. We start by introducing our ideas about how to identify the different entities in a CLP(FD) program such as constraints and variables. Following that, the trace events are presented together with some examples. Finally, our implementation in SICStus Prolog is described.

### 3.1 Representing and Identifying Entities in CLP(FD) Programs

#### 3.1.1 Identifiers.

The id of a constraint is a Prolog atom created by prepending the atom `ctr_` to the source code functor and appending a unique666Unique for the specific type of constraint. number to this. The id of a variable is a Prolog atom of the form `fdvar_N` where `N` is a unique number. For example, the following identifiers are produced by Ex. 1 where `C0`, `C1` and `C2` are the constraints posted on lines `(2)`, `(3)` and `(4)` respectively:

> | { | $`𝚇\mathrm{𝚏𝚍𝚟𝚊𝚛}\mathrm{\_}\mathrm{𝟷},𝚈\mathrm{𝚏𝚍𝚟𝚊𝚛}\mathrm{\_}\mathrm{𝟸},\mathrm{𝚅𝟷}\mathrm{𝚏𝚍𝚟𝚊𝚛}\mathrm{\_}\mathrm{𝟹},\mathrm{𝚅𝟸}\mathrm{𝚏𝚍𝚟𝚊𝚛}\mathrm{\_}\mathrm{𝟺},`$ |
> | --- | --- |
> |  | $`\mathrm{𝙲𝟶}\mathrm{𝚌𝚝𝚛}\mathrm{\_}\mathrm{𝚍𝚘𝚖𝚊𝚒𝚗}\mathrm{\_}\mathrm{𝟷},\mathrm{𝙲𝟷}\mathrm{𝚌𝚝𝚛}\mathrm{\_}\mathrm{𝚊𝚕𝚕}\mathrm{\_}\mathrm{𝚍𝚒𝚏𝚏𝚎𝚛𝚎𝚗𝚝}\mathrm{\_}\mathrm{𝟷},\mathrm{𝙲𝟸}\mathrm{𝚌𝚝𝚛}\mathrm{\_}\mathrm{𝚎𝚕𝚎𝚖𝚎𝚗𝚝}\mathrm{\_}\mathrm{𝟷}`$} |

#### 3.1.2 Constraint Representation.

A constraint is represented similarly to its source code representation. It has the same functor atom and structure, with all variables replaced by their given identifiers. In Ex. 1 for instance, the constraint

> `all_different([X,Y,3,V1,8,V2])`

is represented as

> `all_different([fdvar_1,fdvar_2,3,fdvar_3,8,fdvar_4])`.

#### 3.1.3 Location of Variables.

It is not enough to be able to identify variables by their identifiers. While this allows to distinguish different variables, it will not make it possible to, within one constraint, uniquely determine a specific occurrence among several of the same variable.

We have solved this with what we call a path. A path uniquely determines a position or a list of positions in a constraint. When referencing arguments of a constraint, for example when displaying a variable’s domain, a path referring to that variable is given instead of its id. This path must always have a context, a corresponding constraint, to be interpreted in.

A path is represented as a Prolog list containing (lists of) fixed values or other paths. The first element of a path identifies the topmost position in the corresponding constraint. For example, any path referring to some position(s) in the `element/3` constraint in Ex. 1 has one of the integers 1, 2, 3 or one of the lists `[1,2]`, `[1,3]`, `[2,3]`, `[1,2,3]` as first element. Lines a and b in Table 1 show two of these in paths with only one element and the position(s) they refer to.

The i+1:st element of a path refers to some position(s) inside whatever the i:th element refers to. This means that the position(s) the i:th element refers to must be a compound term or a list of compound terms. Line c in Table 1 illustrates this. This example also introduces one more symbol in the path term; the `#` sign. It is used whenever the term(s) we refer to is inside a list (a non-list compound term would not have a preceding `#` sign).

Some more syntax is needed in order for the path to be more expressive and compact. First of all, we will use `[*]` to express “all positions” of a compound term. It may also be useful to identify almost all positions in a compound term. We will use the `\` sign for this purpose, denoting set subtraction777In our case, the sets are lists.. Line d in Table 1 illustrates this.

In addition to the above, we may apply functions to the entities referred to by a path. `P/F` has the meaning that the function `F` is applied to each entity referred to by the path `P`. Possible values for `F` are, among others, `min`, `max` and `length`. Table 2 shows some examples of this. If `P` refers to a list of positions of size $`n`$, `P/F` returns a list of values of size $`n`$. However, one element lists are simplified as shown on lines a and b in Table 2.

### 3.2 Trace Events

Now, it is time to present the trace events. They are organized in chunks, each chunk containing trace events for a specific purpose. The chunks of trace events are presented in tables followed by examples generated by Ex. 1.

All trace events have an `EventId` attribute, an increasing number. We will often identify them by the value of this attribute and simply call it id. The `EventId` is always the first attribute of the Prolog fact. Furthermore, when we talk about a whole block of trace events, we mean the opening and closing trace events of the block and all trace events in between. For example, a `prune` block includes the trace events `begin_prune` and `end_prune`, and all trace events between these.

Within our illustrations of Ex. 1, since the constraint representations occur very frequently, we will use the following short notations:

* `ALLDIFF = all_different([fdvar_1,fdvar_2,3,fdvar_3,8,fdvar_4])`
* `ELEMENT = element(fdvar_1,[2,4,3,8],fdvar_2)`

Also, when mentioning variables, we will use the variable name from the source code of Ex. 1, even if the trace events contain the given identifier.

#### 3.2.1 Control of Execution.

In this section, trace events describing the control of execution are presented. This concerns things like new constraint posting, demon enqueuing and demon wakening. The trace events that describe these actions are presented in Table 3.

Lines `(3)` and `(4)` in Ex. 1 post two constraints to the system. This gives rise to the following trace events to be created999There are several gaps, identified by ‘…’, since we only present the trace events relevant to the control of execution at this stage.:

```
begin_new_ctr(9,ctr_all_different_1,ALLDIFF)
  new_demon(10,ctr_all_different_1,all_different_1,[[1,#[*]]-val],
            ALLDIFF)
  begin_wake_demon(11,ctr_all_different_1,all_different_1)
  ...
  end_wake_demon(21,delay)
end_new_ctr(22,delay)
begin_new_ctr(23,ctr_element_1,ELEMENT)
  ...
  begin_new_ctr(25,ctr_in_1,fdvar_1 in 1..4)
  ...
  end_new_ctr(32,entail)
  new_demon(33,ctr_element_1,element_3,
            [[1]-dom,[2,#[*]]-minmax,[3]-minmax],ELEMENT)
  begin_wake_demon(34,ctr_element_1,element_3)
  ...
  end_wake_demon(57,delay)
end_new_ctr(58,delay)

```

The first five trace events are created due to the posting of the constraint `all_different/1`. This gives rise to the creation of a new demon (id `10`) and the wakening of that demon (id:s `11` and `21`). The last one of these trace events (id `22`) contains the result `delay`, meaning that the constraint will not wake up again until the state of at least one of its variables changes. Moving on, the `element/3` constraint is introduced in much the same way. The difference is the creation of an internal `in/1` constraint (id:s `25` and `32`) to ensure that `X` takes a value in the range `1..4`.

#### 3.2.2 Domain Modification.

Now, we introduce trace events describing domain modification, i.e. pruning of variables. Such trace events are presented in Table 4.

For Ex. 1, the following trace events are added due to the posting of the `all_different/1` constraint on line `(3)`:

```
begin_prune(14,remove_value(3),ALLDIFF)
  ...
  before_prune(16,[1,#[1,2,4,6]],[[[1|6]],[[1|6]],[[1|6]],[[1|6]]],
               ALLDIFF)
  prune(17,ctr_all_different_1,[1,#[1,2,4,6]],remove_value(3),ALLDIFF,
        succeed)
  after_prune(18,[1,#[1,2,4,6]],[[[1|2],[4|6]],[[1|2],[4|6]],
                                 [[1|2],[4|6]],[[1|2],[4|6]]],ALLDIFF)
end_prune(19,succeed)

```

These will occur between the trace events with id:s `11` and `21`, the wakening of the `all_different/1` demon. The first trace event contains information about the intended domain narrowing. In this case the value `3` is to be removed from some variables. The path to these variables is shown in the following three trace events. The first and last of these (id:s `16` and `18`) contain domain information about the variables before and after pruning. The middle trace event (id `17`) contains information about the actual pruning. Finally, the trace event with id `19` contains the result of the pruning. Since it did not produce any empty domains, the result is `succeed`.

The posting of the `element/3` constraint generates the following trace events containing information about pruning the variable `X`:

```
begin_prune(50,remove_values([[4|4]]),ELEMENT)
  ...
  before_prune(52,[1],[[[1|2],[4|4]]],ELEMENT)
  prune(53,ctr_element_1,[1],remove_value(4),ELEMENT,succeed)
  after_prune(54,[1],[[[1|2]]],ELEMENT)
end_prune(55,succeed)

```

#### 3.2.3 Explanations.

The trace events presented in this section provide explanations  for other trace events, i.e. reasons for why certain actions are taken by the solver. These trace events are presented in Table 5.

First of all, we will make clear what we mean with explanations: An explanation $`e`$ for some kernel action $`a`$ is a logical formula $`e`$ such that $`ea`$, i.e. if $`e`$ holds, $`a`$ will occur. In Prolog, we represent an explanation with a pair, `N-CondList`, where `N` is an integer and `CondList` is a list of `cond(M,Pa,PL)` terms. An explanation $`e=`$ `N-CondList` holds iff at least `N` of the `cond/3` terms in `CondList` hold. Each `cond/3` term contains an integer `M`, a path `Pa` and a list of properties `PL`. The property on index $`i`$ in `PL` is associated with the entity on index $`i`$ in the list of positions referred to by `Pa`. For example, assume that `Pa` $`=`$ `[[1,3]]`, `PL` $`=`$ `[eq(2),neq(3)]` and that `Pa` is associated with the `element/3` constraint in Ex. 1. Then `X` is associated with `eq(2)` and `Y` is associated with `neq(3)`. When `PL` contains a single property, it is associated with each entity referred to by `Pa`. Also, the list structure is skipped and `PL` will denote the value of the single property. Possible properties in `PL` are, among others, `eq(I)`, `neq(I)`, `inset(S)` and `notinset(S)` where `I` is an integer and `S` is a finite set of integers. A `cond/3` term holds iff at least `M` of the entities referred to by `Pa` fulfill their respective properties.

Now, let us give an illustrative example. Assume that

> $`e=`$ `1-[cond(1,[1,#3],eq(3))]`

explains the solver action $`a`$. Also, assume that $`e`$ is associated with a pruning generated due to the `all_different/1` constraint in Ex. 1. This gives us the information that $`a`$ occurs since at least one of the `cond/3` terms holds. For the only one, it holds since at least one of the entities referred to by `[1,#3]` can take the value `3`, `eq(3)`. In our case, this is the position marked with gray in all\_different(\[X,Y,3,V1,8,V2\]). In order for $`a`$ not to occur, the least thing one must do is to change the value on that position to some value distinct from `3`. Actually, this explanation occurs in the trace event

```
prune_because(15,ctr_all_different_1,1-[cond(1,[1,#3],eq(3))],ALLDIFF).

```

Hence, $`e`$ is an explanation for removing the value `3` from all variables in the `all_different/1` constraint.

The following are some more explanation trace events from Ex. 1:

```
new_ctr_because(24,ctr_element_1,1-[cond(1,[2]/length,eq(4))],ELEMENT)
...
prune_because(51,ctr_element_1,1-[cond(1,[3],notinset([[8|8]]))],ELEMENT)

```

The first one is generated before the `begin_new_ctr` trace event with id `25`. It contains an explanation for posting the `X in 1..4` constraint inside the creation of the `element/3` constraint. The explanation is that the length of the list on the second position in the constraint is `4`.

Similarly, the id of the last one tells us that it is generated before the `before_prune` trace event with id `52`. It contains an explanation for pruning the variable `X`. This too occurs inside the creation of the `element/3` constraint. The explanation is that the third argument of the constraint, the variable `Y`, does not intersect with the set consisting of the single value `8`. Due to this, `X` cannot take the value `4`.

For some solver actions, there are several different reasons for why they occur. In order to show this, we need to introduce one more example.

The `all_distinct/1` constraint on line `(4)` in Ex. 2 uses the complete filtering algorithm of Régin . When, during the posting of that constraint, its demon is woken up, it will produce a failure since the variables `X2`, `X3`, `X4` and `X5` are all in the range `1..2`. For this failure, we generate the following explanation:

> `1-[cond(all,[1,#[2,3,4,5]],inset([[1|2]]))]`

where `all` is a keyword for identifying all entities referred to by the associated path. By looking at this explanation, we conclude that, in order for the failure not to occur, we need to enlarge the domains of any two of the variables `X2`, `X3`, `X4` and `X5` with at least two distinct values. We see this since there are $`4`$ variables taking values in a set of integers of size $`2`$ and $`42=2`$.

#### 3.2.4 Filtering Algorithms.

The trace events presented in this section describe the execution of the filtering algorithms associated with the demons. Since each filtering algorithm is built up by a set of methods, information about each active method is needed. Table 6 presents the relevant trace events.

Each `prune` block is surrounded by a `method` block, and possibly preceded by an `info_method` trace event. For the first example concerning domain modification, the following trace events are added:

```
begin_method(12,propagate_ground_variable)
  info_method(13,ground_variable,[1,#3]-3,ALLDIFF)
  ...
end_method(20,succeed)

```

By looking at the id:s, we see that they are generated before and after the `prune` block describing the removal of the value `3` from the variables `X`, `Y`, `V1` and `V2`. The `begin_method` trace event contains a name of the active method. In this case, `propagate_ground_variable` denotes the method that removes a ground value from all non-ground domain variables in the `all_different/1` constraint. Following that, the `info_method` trace event contains more detailed information about the actual entities involved in the method. More specifically, it tells us that the ground variable is on the position referred to by `[1,#3]` and its value is `3`. Finally, the `end_method` trace event contains the result, `succeed`, of the whole `method` block. The result will be `fail` iff at least one failure inside the block is detected, otherwise `succeed`.

The second example concerning domain modification, the one associated with the `element/3` constraint, is surrounded by the following `method` block:

```
begin_method(49,prune_x)
...
end_method(56,succeed)

```

### 3.3 Implementation

The described trace structure has been implemented in a research branch of the CLP(FD) solver in SICStus Prolog version 3.9.

Each trace event is represented by an abstract datatype (ADT) in C. Complex information, such as compound terms and lists, is represented using C strings. Each ADT has creation and deletion primitives as well as a primitive for converting the information to a Prolog fact.

On the Prolog side, it is possible to create and access the trace events through some predicates. For the creation, each trace event (with two exceptions131313The trace events push\_demon and push\_demon\_because are created at a central place within the C part of the kernel.) have a corresponding constructor predicate. For retrieving, a predicate `get_events(+Selector, -TraceEvents)` is available. `Selector` is an integer, an atom or a compound term specifying a set of trace events. `TraceEvents` is a list of trace events corresponding to `Selector`.

## 4 Discussion

Tracing CLP(FD) programs is still an open area. In our case, more testing is needed in order to evaluate the trace model. This includes building debugging tools, using the different kinds of information the trace provides.

One interesting tool to build would be a search-tree visualizer such as the Oz Explorer  by Schulte. Such a tool would explore the information in the trace concerning the control of execution aspect. Also, as was done for the Oz Explorer, one could connect to it other tools using a plug-in scheme, displaying various information concerning other aspects at the nodes of the search-tree. An example of such a tool for the Oz Explorer is the Constraint Investigator  by Müller.

Another tool, exploring the information concerning the domain modification aspect, would be a variable domain visualizer. Domains could for instance be visualized in a way similar to as in .

The information concerning the declarative aspect could be used by a tool for solving over-constrained problems. Explanations would help to point out the set of inconsistent constraints. Jussien and Ouis presented a tool like this in .

For educational purposes, a tool that provides information about filtering algorithms would be useful. Such a tool could connect the consistency rules of the constraints with the different methods in the corresponding filtering algorithms.

As a first use of the trace structure, we have implemented an extension to the Prolog debugger in SICStus. At the usual breakpoints of the debugger, it is possible to ask queries about the trace and the entities in it. The user may for example display certain trace events, ask queries about domain variables and generate explanations for specific domain narrowings.

## 5 Conclusion

In this paper, we presented new ideas in the area of tracing CLP(FD) programs. In addition to providing information about control of execution and domain modification, we also provide information concerning the declarative and procedural aspects. We did this by incorporating the idea of explanations in the trace, as well as coming up with new ideas about tracing filtering algorithms. By doing this, we have obtained a link between the consistency rules and the filtering algorithms in constraint solving.

Furthermore, we presented some new ideas about how to present the trace, both regarding the structure and the information contained. This includes building the trace as a nested block structure, as well as using a path for identifying the different entities in a CLP(FD) program.

Implementing the trace structure in a high-end system like SICStus is not an easy task. This is especially true for the task of adapting filtering algorithms to generate explanations. Sometimes, the added overhead it meant to provide the information needed for sharp explanations was too large. In this case, some approximation explanation141414An explanation that is true but not necessarily the one that expresses the least dependency. has been used. In the long run, providing explanations probably means that filtering algorithms need to be changed in order to, more naturally, provide the necessary information.