text:
# Partially fluidized shear granular flows: Continuum theory and MD simulations

## I Introduction

In the last few year there have been many experimental behringer99 ; gollub98 ; gollub00 ; chicago ; pouliquen99 ; daerr ; Tsai02  and theoretical grest91 ; wolf96 ; grest02 ; sparks99 ; sparks02 ; lemaitre ; ertas  studies that explored a broad range of granular flow conditions from rapid dilute flows to slow dense flows, as well as the details of the shear-driven fluidization transition. While dilute granular flows can be well described by the kinetic theory of dissipative granular gases Jenkins85 , dense granular flows still present significant difficulty in formulation of a continuous theory. In Ref. savage98 , Savage proposed a continuum theory for slow dense granular flows based on the so-called associated flow rule that relates the strain rate and the shear stress in plastic frictional systems. Averaging strain rate fluctuations yields a Bingham-like constitutive relation in which the shear stress has a viscous as well as a strain-rate independent parts. According to this theory, the stress and strain rate tensors are always co-axial and, furthermore, it also postulates that the viscosity diverges as the density approaches close packing limit. Losert et al. gollub00  (see also bocquet02 ) proposed a similar hydrodynamic model based on a Newtonian stress-strain constitutive relation with density dependent viscosity without strain-rate independent component. As observed in Ref. gollub00 , the ratio of the full shear stress to the strain rate diverges at the fluidization threshold. This was also interpreted in Ref. gollub00  as a divergence of the viscosity coefficient when the volume fraction approaches the randomly packed limit. This description works only in a fluidized state and can not properly account for hysteretic phenomena in which static and fluidized states co-exist under the same external load, such as stick-slip oscillations gollub98 , avalanching daerr , or shear band formation.

In many granular flows of interest static and dynamics regions co-exist under the same external load conditions. Examples of such hysteretic phenomena include stick-slip oscillations gollub98 , avalanching daerr , or shear band formation. This calls for a unified theory which would be applicable both in the flowing regime and in the static regime.

In our recent papers AT1 ; AT2  we proposed a different approach based on the order parameter description of the granular matter. The value of the order parameter specifies the ratio between static and fluid parts of the stress tensor. The order parameter was assumed to obey dissipative dynamics governed by a free energy functional with two local minima. This description based on the separation of static and fluid components of the shear stress, calls for an alternative definition of viscosity as a ratio of the fluid part of the shear stress to the strain rate. Since the fluid shear stress vanishes together with the strain rate, the viscosity coefficient in our theory is expected to remain finite at the fluidization threshold. We assumed the simplest Newtonian friction law, so the viscosity coefficient is a constant. This model yielded a good qualitative description of many phenomena occurring in granular flows, such as hysteretic transition to chute flow, stick-slip regime of a driven near-surface flow, structure of avalanches in shallow chute flows, etc.

However, several important issues have not been addressed: mesoscopic definition of the order parameter, quantitative specification of the order parameter dynamics and the constitutive relation. In this paper we set out to perform two dimensional (2D) molecular dynamics simulations which should provide us with the way to achieve these goals. We classify all contacts as either “fluid-like” or “solid-like” and define the order parameter as a mesoscopic space-time average fraction of solid-like contacts. Using this order parameter we obtain the constitutive relation and the relaxational dynamics of the order parameter directly from simulations of granular flow in a thin Couette geometry at zero gravity. Preliminary account of our results is presented in Ref. volfson .

The paper is organized as follows. Section II outlines the standard granular hydrodynamics theory based on the kinetic theory of dissipative granular gases. Section III introduces the continuum description of partially fluidized flows based on the relaxational dynamics of the order parameter. In Section IV we describe our 2D molecular dynamics simulations and define measurement protocol for the order parameter and fluid and solid components of the stress tensor. In Section V we study Couette flow in a thin granular layer to obtain the free energy controlling the relaxational order parameter dynamics and to extract the constitutive relations. In Section VI the obtained set of equations is used to calculate the stress and velocity distributions in a different system, a thick granular layer under non-zero gravity driven by a moving heavy upper plate.

## II Granular hydrodynamics

In this Section we outline the standard continuum description of granular flows based on continuity equations for mass, momentum and fluctuation kinetic energy (or “granular temperature”). This description is usually applied to dilute granular gases where it can be rigorously derived from the kinetic theory Jenkins85 , although slightly modified hydrodynamics based on kinetic theory often works reasonably well for relatively dense flows, even though the kinetic theory itself is not applicable to these conditions.

The mass, momentum and energy conservation equations have the usual form

$`{\displaystyle \frac{D\nu }{Dt}}`$ $`=`$ $`\nu 𝐮,`$ (1)
$`\nu {\displaystyle \frac{D𝐮}{Dt}}`$ $`=`$ $`𝝈+\nu 𝐠,`$ (2)
$`\nu {\displaystyle \frac{DT}{Dt}}`$ $`=`$ $`𝝈:\dot{𝜸}𝐪\epsilon ,`$ (3)
where $`\nu `$ is the density, $`𝐮`$ is the velocity field, $`T=\stackrel{~}{u}_i\stackrel{~}{u}_i/2`$ is the granular temperature, $`D/Dt=_t+(𝐮)`$ is the material derivative, $`𝐠`$ is the gravity acceleration, $`\sigma _{\alpha \beta }`$ is the stress tensor, $`𝐪`$ is the energy flux vector, $`\dot{\gamma }_{\alpha \beta }=_\alpha u_\beta +_\beta u_\alpha `$ is the strain rate tensor, and $`\epsilon `$ is the energy dissipation rate.

These three equations have to be supplemented by the constitutive relations for the stress tensor $`𝝈`$, energy flux $`𝐪`$, and the energy dissipation rate $`\epsilon `$. For dilute systems, a linear relations between stress $`𝝈`$ and strain rate $`\dot{𝜸}`$ is obtained,

$`\sigma _{\alpha \beta }`$ $`=`$ $`[p+(\mu \lambda )\text{Tr}\dot{𝜸}]\delta _{\alpha \beta }\mu \dot{\gamma }_{\alpha \beta },`$ (4)
$`𝐪`$ $`=`$ $`\kappa T,`$ (5)
In the kinetic theory of granular gases Jenkins85 , these equations are closed with the following equation of state

$$p=\frac{4\nu T}{\pi d^2}[1+(1+e)G(\nu )],$$
(6)
and the expressions for the shear and bulk viscosities

$`\mu `$ $`=`$ $`{\displaystyle \frac{\nu T^{1/2}}{2\pi ^{1/2}dG(\nu )}}\left[1+2G(\nu )+\left(1+{\displaystyle \frac{8}{\pi }}\right)G(\nu )^2\right],`$ (7)
$`\lambda `$ $`=`$ $`{\displaystyle \frac{8\nu G(\nu )T^{1/2}}{\pi ^{3/2}d}},`$ (8)
the thermal conductivity

$$\kappa =\frac{2\nu T^{1/2}}{\pi ^{1/2}dG(\nu )}\left[1+3G(\nu )+\left(\frac{9}{4}+\frac{4}{\pi }\right)G(\nu )^2\right],$$
(9)
and the energy dissipation rate

$$\epsilon =\frac{16\nu G(\nu )T^{3/2}}{\pi ^{3/2}d^3}(1e^2).$$
(10)
Here $`0<e<1`$ is restitution coefficient and $`d`$ is particle diameter. The function $`G(\nu )`$ which enters these relations is the spatial particle-particle correlation function, and for a dilute 2D gas of elastic hard disks was derived by Carnahan and Starling carnahan ,

$$G_{CS}(\nu )=\frac{\nu (17\nu /16)}{(1\nu )^2}$$
(11)
This formula is expected to work for densities roughly below 0.7. For high density granular gases, this function has been calculated using free volume theory kirkwood51 ,

$$G_{FV}=\frac{1}{(1+e)\left[(\nu _c/\nu )^{1/2}1\right]}$$
(12)
where $`\nu _c0.82`$ is the density of the random close packing limit. Luding Luding01  proposed a global fit

$$G_L=G_{CS}+(1+\mathrm{exp}((\nu \nu _0)/m_0))^1)(G_{FV}G_{CS})$$
(13)
with empirically fitted parameters $`\nu _00.7`$ and $`m_010^2`$. However, even with this extension, the continuum theory comprised of Eqs.(1)-(10) cannot describe the force chains which transmit stress via persistent contacts remaining in the dense granular flows, as well as the transition from solid to static regimes and coexisting solid and fluid phases.

## III Order parameter description of partially fluidized granular flows

In this Section we review briefly our continuum theory AT1 ; AT2  which provides an alternative approach to the formulation of the constitutive relations in partially fluidized granular flows. In the dense flow regime, the granular matter can be considered incompressible, so Eq. (1) can be replaced by $`𝐮=0`$ and the density $`\nu =\nu _c`$. This also allows us to drop the energy equation (3) and the equation of state (6), as for $`\nu \nu _c`$, $`G(\nu )\mathrm{}`$ and $`T0`$, so that $`G(\nu )Tconst`$. This of course leads to the familiar divergence of the viscosity coefficient $`\mu G(\nu )T^{1/2}`$.

Next, we separate the stress tensor $`𝝈`$ into two parts, a static (contact) part $`𝝈^s`$, and a fluid part $`𝝈^f`$. The latter is assumed to take a purely Newtonian form

$$\sigma _{\alpha \beta }^f=p_f\delta _{\alpha \beta }\mu _f\dot{\gamma }_{\alpha \beta }$$
(14)
where $`p_f`$ is the “partial” fluid pressure, $`\mu _f`$ is the viscosity coefficient associated with the fluid stress tensor which is different from $`\mu _0`$ introduced for the full stress tensor. As we shall see in the following, unlike $`\mu _0`$, $`\mu _f`$ does not diverge as $`\nu \nu _c`$. In our original model AT1 ; AT2  we simply assumed $`\mu _f=const`$.

We postulated note1  that the fluid part of the off-diagonal components of the stress tensor is proportional to the off-diagonal components of the full stress tensor with the proportionality coefficient being a function of the order parameter $`\rho `$,

$$\sigma _{yx}^f=q(\rho )\sigma _{yx}.$$
(15)
This assumption stipulates the equation for the static part of the off-diagonal stress components,

$$\sigma _{yx}^s=(1q(\rho ))\sigma _{yx}.$$
(16)
Both fluid and solid parts of the stress tensor are assumed symmetric, $`\sigma _{yx}^{f,s}=\sigma _{xy}^{f,s}`$. This assumption is confirmed by our numerical simulations (see below). We choose a fixed range for the order parameter such that it is zero in a completely fluidized state and one in a completely static regime. Thus the function $`q(\rho )`$ has the property $`q(0)=1,q(1)=0`$. In Refs.AT1 ; AT2  for simplicity we took $`q(\rho )=1\rho `$.

A similar relationship can be postulated for the diagonal terms of the stress tensor,

$`\sigma _{xx}^f`$ $`=`$ $`q_x(\rho )\sigma _{xx},\sigma _{yy}^f=q_y(\rho )\sigma _{yy},`$ (17)
$`\sigma _{xx}^s`$ $`=`$ $`(1q_x(\rho ))\sigma _{xx},\sigma _{yy}^s=(1q_y(\rho ))\sigma _{yy},`$ (18)
where the scaling functions $`q_{x,y}(\rho )`$ can differ from $`q(\rho )`$.

Combining Eqs.(14)-(18), we obtain the constitutive relation in the closed form,

$$\sigma _{\alpha \beta }=p_f\delta _{\alpha \beta }/q_\alpha (\rho )\mu _f\dot{\gamma }_{\alpha \beta }/q(\rho )$$
(19)
where $`\alpha ,\beta \{x,y\}`$.

The order parameter itself was not related to any microscopic properties of granular assemblies in Refs. AT1 ; AT2 . We simply assumed that because of strong dissipation in dense granular flows it has purely relaxational dynamics controlled by the Ginzburg-Landau equation,

$$\frac{D\rho }{Dt}=D^2\rho \frac{(\rho )}{\rho }$$
(20)
Here $`D`$ is the diffusion coefficient and $`(\rho )`$ is the free energy density which was assumed to have a quartic polynomial form to account for the bistability near the solid-fluid transition,

$$(\rho )=^\rho \rho (\rho 1)(\rho \delta )𝑑\rho .$$
(21)
The control parameter $`\delta `$ is determined by the stress tensor which in Ref. AT1 ; AT2  was taken to be a linear function of $`\varphi =\mathrm{max}|\sigma _{mn}/\sigma _{nn}|`$, where the maximum is sought over all possible orthogonal directions $`m`$ and $`n`$. It is easy to see that in the interval $`0<\delta <1`$ Eq. (20) has two stable uniform solutions $`\rho =0,1`$ corresponding to fluid and solid states and one unstable solution $`\rho =\delta `$.

Momentum conservation equation (2) together with Eqs. (19)-(21) represent a closed set of continuous equations which after being augmented by appropriate boundary conditions, can describe a variety of interesting granular flows such as avalanches in thin chute flows, drum flows, stick-slip oscillation in surface-driven flows, etc. AT1 ; AT2 . Since a rigorous derivation of the continuum model of dense partially fluidized flows from first principles does not seem feasible, the assumptions made ad hoc in the formulation of the model (2),(19)-(21) have to be checked against available experimental and/or numerical data. Unfortunately, it appears to be extremely difficult to directly measure the order parameter in the bulk of granular material, as it is determined by subtle changes in the contact fabric. In this paper we attempt to extract the properties of the order parameter from molecular dynamics simulations and on this basis make a quantitative fit of the order parameter model.

## IV Molecular Dynamics Simulations

To model the interaction of individual grains we use the so-called soft-contact approach. The grains are assumed to be non-cohesive, dry, inelastic disk-like particles. Two grains interact via normal and shear forces whenever they overlap. For the normal impact we employ *spring-dashpot* model wolf96 . This model accounts for repulsion and dissipation; the repulsive component is proportional to the degree of the overlap, and the velocity dependent damping component simulates the dissipation. The model for shear force is based upon the technique developed by Cundall and Strack cundall79 . It incorporates tangential elasticity and Coulomb laws of friction. The elastic restoring force is proportional to the integrated tangential displacement during the contact and limited by the product of the friction coefficient and the instantaneous normal force. The grains possess two translational and one rotational degrees of freedom. The motion of a grain is obtained by integrating the Newton’s equations with the forces and torques produced by its interactions with all the neighboring grains and walls of the container.

Consider a grain $`i`$ of radius $`R^i`$ located at $`𝐫^i`$ moving with translational velocity $`𝐯^i`$ and angular velocity $`\omega ^i`$. This grain is in contact with grain $`j`$ whenever overlap $`\delta _n=R^i+R^j|𝐫^i𝐫^i|>0`$. The relative velocity at the contact point and its normal and tangent components are given by

$`𝐯^{ij}`$ $`=`$ $`𝐯^i𝐯^j+(R^i\omega ^i+R^j\omega ^j)𝐭^{ij},`$ (22)
$`v_n`$ $`=`$ $`(𝐯^i𝐯^j)𝐧^{ij},v_t=(𝐯^i𝐯^j)𝐧^{ij}+(R^i\omega ^i+R^j\omega ^j),`$ (23)
where $`𝐧^{ij}=(𝐫^i𝐫^i)/r^{ij}=(n_x,n_y)`$ is the inward normal to the surface of $`i`$ at the contact point with $`j`$, and the direction of the tangent $`𝐭^{ij}=(n_y,n_x)`$ is chosen so that $`𝐭^{ij}\times 𝐧^{ij}`$ is co-linear with the angular velocity. Then the grain $`i`$ is subjected to the contact force due to the interaction with $`j`$ with normal and tangential components

$`F_n`$ $`=`$ $`k_n\delta _n2\gamma _nm_ev_n,`$ (24)
$`F_t`$ $`=`$ $`\text{sign}(\delta _t)\mathrm{min}(|k_t\delta _t|,|\mu F_n|),\delta _t(t)={\displaystyle _{t_0}^t}v_t(\tau )𝑑\tau ,`$ (25)
where $`k_{n,t}`$ are the corresponding spring constants, $`\gamma _n`$ is the normal damping coefficient, $`m_e=m^im^j/(m^i+m^j)`$ is the reduced mass, $`\mu `$ is coefficient of friction, and $`\delta _t`$ is tangential displacement since the moment $`t_0`$ of the initial contact between $`i`$ and $`j`$. When the static yield criterion Eq. (25) is satisfied, the magnitude of $`\delta _t`$ is adjusted to an instantaneous equilibrium value providing $`F_t=\mu F_n`$. According to the analytical solution of the linear spring-dashpot model the coefficient of restitution and the duration of heads-on collision are $`e=\mathrm{exp}(\gamma _nt_n)`$, $`t_n=\pi /(k_n/m_e\gamma _n^2)^{1/2}`$.

The advantages and limitations of the employed contact force model were thoroughly studied by a number of authors wolf96 ; grest91 ; grest02 . In fact this is the simplest model which allows us to account for both static and dynamic friction.

When all forces acting on grain $`i`$ from other grains, boundaries and perhaps external fields are computed, the problem is reduced to the integration of the Newton’s equations for translational and rotational degrees of freedom,

$`m^i{\displaystyle \frac{d^2𝐫^i}{dt^2}}`$ $`=`$ $`m^i𝐠+{\displaystyle \underset{c}{}}𝐅^{ic},`$ (26)
$`I^i{\displaystyle \frac{d\omega ^i}{dt}}`$ $`=`$ $`R^i{\displaystyle \underset{c}{}}F_t^{ic},`$ (27)
where the mass of particle $`i`$ is denoted by $`m^i`$ and its moment of inertia is $`I^i=1/2m^iR_{}^{i}{}_{}{}^{2}`$, $`m^i𝐠`$ stands for an external gravity field, and the sums in Eqs. (26),(27) run over all contacts of particle $`i`$.

The results of the simulations reported here are presented in dimensionless form. All quantities are normalized by an appropriate combination of the average particle diameter $`d`$, mass $`m`$, and gravity $`g`$.

Eqs. (26),(27) were integrated using the fifth-order predictor-corrector allen87  with a constant time step $`\delta t`$. The spring constant $`k_n`$ and damping coefficient $`\gamma _n`$ were chosen to provide the desired value of the restitution coefficient $`e`$ and guarantee an accurate resolution of an individual collision. Typically, we used $`e=0.92`$, $`\delta t=10^4`$, $`t_c=50\delta t`$, $`k_n=\mathrm{2.0\; 10}^5`$, $`\gamma _n=16.7`$, $`k_s=1/3k_n`$.

The computational domain spans $`L_x\times L_y`$ area, and is periodic in horizontal direction $`x`$. Unless otherwise mentioned the granulate is slightly polydisperse to avoid crystallization effects grest91 . We assume that the grain diameters are uniformly distributed around mean with relative width $`\mathrm{\Delta }_r`$. To provide a link between micro-mechanical quantities obtained through simulations and continuous fields, we define the following coarse-graining procedure. Since all experiments described below deal with steady quantities, the procedure consists of two steps: space and then time averaging. The space averaging of a field x is performed over horizontal bins (along the flow direction) of size $`V_bL_x\times d`$ and is denoted as $`\text{x}`$ in the following. Contributions of those particles which only partially belong to a certain bin are weighted by the fraction of their area. After a simulation has reached a steady state, instantaneous profiles are averaged over a suitable number of time-snapshots. We shall denote the time-averaging with $`\overline{\text{x}}`$.

For example, steady solid fraction profile is given by

$$\nu (y)=\overline{\nu (y,t)},\nu (y,t)=V_b^1\underset{iV_b}{}w^i(y)V^i,$$
(28)
where the summation runs over grains at least partially in a bin $`V_b`$ centered at $`y`$, $`V^i`$ is the area of grain $`i`$, and $`w^i(y)`$ is a corresponding fraction of a grain’s area within $`V_b`$. The coarse-grained velocity field is

$$𝐮(y)=\overline{𝐮(y,t)},𝐮(y,t)=\nu (y,t)^1𝐯^i(y,t)=\nu (y,t)^1V_b^1\underset{iV_b}{}w^i(y)V^i𝐯^i.$$
(29)
For simplicity we extend the space-time averaging technique described above for other quantities, such as the stress tensor

$$\sigma _{\alpha \beta }(y)=\overline{\sigma _{\alpha \beta }(y,t)},\sigma _{\alpha \beta }(y,t)=\frac{1}{2}\underset{ci}{}r_\alpha ^{ic}F_\beta ^{ic}+m^i\stackrel{~}{v}_\alpha ^i\stackrel{~}{v}_\beta ^i,$$
(30)
where $`\alpha ,\beta =\{x,y\}`$, $`r_\alpha ^{ic}=𝐫^{ic}𝐞_\alpha `$, $`F_\beta ^{ic}=𝐅^{ic}𝐞_\beta `$, $`\stackrel{~}{v}_\alpha ^i=v_\alpha ^i\overline{v_\alpha ^i}`$. The stress tensor in Eq. (30) has two distinct components. The first one – virial or contact – describes pairwise interactions of grains. The second one – kinetic or Reynolds – is due to velocity fluctuations.

### IV.1 Order parameter for granular fluidization: static contacts vs. fluid contacts

At any moment of time all contacts are classified as either “fluid” or “solid”. A contact is considered “solid” if it is in a stuck state ($`F_t<\mu F_n`$) and its duration is longer than a typical time of collision $`t_{}`$. The first requirements eliminates long-lasting sliding contacts, and the second requirement excludes short-term collisions pertinent for completely fluidized regimes. We choose a typical collision to last $`t^{}=1.1t_c`$. When either of the requirements is not fulfilled, the contact is assumed “fluid”.

We define the order parameter as the ratio between space-time averaged numbers of “solid” contacts $`\overline{Z_s}`$ and all contacts $`\overline{Z}`$ within a sampling area silbert02 ,

$$\rho (y)=\overline{Z_s^i}/\overline{Z^i}.$$
(31)
This definition endures at least two limiting cases: when a granulate in a static state and when it is strongly agitated, i.e. completely fluidized. In fact, in the former quiescent state all contacts are stuck and $`\rho =1`$. In the fluidized case $`\overline{Z_s}`$ is zero and $`\overline{Z}`$ is small but finite, therefore $`\rho =0`$.

Let us note that the order parameter $`\rho `$ just introduced is expected to be as sensitive to the degree of fluidization as the stress tensor. A small rearrangement of the force network may result in strong fluctuations of either field, while such quantities as the solid fraction or granular temperature will remain virtually constant. It is known that granular aggregates exhibit rigidity phase transition near the critical volume fraction $`\nu _c`$ sparks99 . When the volume fraction is near or above $`\nu _c`$, the granular material has elasto-plastic rheology and below this critical value it behaves as a fluid. Near and above the critical volume fraction for the same boundary conditions the granular aggregate may exhibit either slow creeping flow, solid-like behavior, or both. The latter case, known as stick-slip phenomenon, is observed both experimentally and numerically behringer96 ; gollub98 ; sparks99 . In an experiment under constant volume boundary conditions a series of stick-slip events occur without a significant change in the solid fraction, which is insensitive to global rearrangements relieving the accumulated stress. We expect the order parameter to be able to reflect such rearrangements and ultimately describe the corresponding phase transition.

### IV.2 Stress tensor

The full stress tensor (30) consists of a contact (virial) part and the Reynolds part. In turn, the contact stress tensor can be split into the “solid contact” component, $`\sigma ^s`$, and the “fluid contact” component, ($`\sigma ^f`$), in the same fashion as was done with contacts themselves. Combining the “fluid contact” component with the Reynolds stress, we obtain the full stress tensor as a sum of two parts

$`\sigma =\sigma ^f+\sigma ^s,`$ (32)
$`\sigma _{\alpha \beta }^s={\displaystyle \frac{1}{2}}{\displaystyle \underset{c}{\overset{}{}}}r_\alpha ^cF_\beta ^c,\sigma _{\alpha \beta }^f={\displaystyle \frac{1}{2}}{\displaystyle \underset{c}{\overset{\prime \prime }{}}}r_\alpha ^cF_\beta ^c+m^i\stackrel{~}{v}_\alpha \stackrel{~}{v}_\beta `$ (33)
where summation in $`^{}`$ and $`^{\prime \prime }`$ is restricted to “solid” and “fluid” contacts respectively.

The “fluid” part of the stress tensor is due to short-term collisional stresses and the Reynolds stresses, whereas the solid part accounts for persistent force chains. The Reynolds contribution to the stress is negligibly small in the vicinity of the phase transition, but comes into play when the granular aggregate is highly fluidized. In the system which is neither completely rigid nor completely fluidized we expect the coexistence (in time and space) of both phases. A particular grain may have both types of contacts at the same time thus contributing to both $`\sigma ^f`$ and $`\sigma ^s`$. This picture is reminiscent of the concept of bimodal character of stress transmission in static contact network introduced by Radjai et al moreau96 ; moreau98 .

## V Testbed system: Couette flow in a thin granular layer

### V.1 Bifurcation diagram

We studied the fluidization transition in a thin ($`50\times 10`$) granular layer between two “rough plates” under fixed pressure $`P`$ and zero gravity conditions (Figure1). The parameters of simulations were $`\mu _0=0.5,e=0.92,\mathrm{\Delta }_r=0.2,k_n=210^5,k_s/k_n=1/3,\gamma _n=16.7`$. We used the same material parameters of grains throughout this Section. The rough plates were simulated by two straight chains of large grains (twice as large as an average particle diameter) glued together. The layer was chosen narrow enough, so the properties of the granular layer (shear rate, shear stress components, order parameter, etc) were constant across the layer. Two opposite forces $`𝐅_1=𝐅_2`$ were applied to the plates along the horizontal $`x`$ axis to induce shear stress in the bulk. We started with weak forces not sufficient to initiate a shear flow and slowly ramped them up in small increments well above the critical yield force at which the granular flow started. After that we ramped the shear forces down until the granular layer was jammed again. At every “stop” we measured all stress components, strain rate, and the order parameter, and averaged the data over the whole layer and over time of each step. Figure 2 shows the strain rate $`\dot{\gamma }`$ (we drop subscripts $`yx`$ at the strain rate) and the order parameter $`\rho `$ as functions of the shear stress $`\sigma _{yx}`$ which is approximated by the applied force $`F`$ normalized by the layer length $`L_x`$, for $`P=40`$. As it is to be expected, the strain rate remains zero, and the order parameter is one until the shear stress reaches a certain critical value $`\sigma _112.6`$. This value differs slightly for different runs because of the finite system size and absence of self-averaging in the static regime. Above the yield stress, the strain rate abruptly jumps to a finite value $`0.35`$, and the order parameter drops to $`0.15`$. At larger $`|\sigma _{yx}|`$, the strain rate increases faster than linearly, and the order parameter rapidly approaches zero. The return curve corresponding to the diminishing of the shear stress follows roughly the same path, and then continues to another (smaller) value of the shear stress ($`\sigma _29.4`$). At this value the layer jams, the strain rate returns to zero, and the order parameter jumps back to one.

The most striking feature of this figure is the hysteretic behavior of both the strain rate and the order parameter as a function of the shear stress. This hysteresis was anticipated in our order-parameter model AT1 ; AT2 , however now we are in a position to fit the model equations quantitatively. We repeated these simulations at several different values of the compressing pressure $`P`$. Data for different pressure values in the flow regime fall onto the same universal curve if one normalizes the shear stress by the pressure (see Figure 3). Assuming that there is an (unobserved) unstable branch of the bifurcation diagram which merges with the stable branch at $`\delta 0.26`$, we can make a simple analytic fit of this curve as

$$F(\rho ,\delta )=(1\rho )\left(\rho ^22\rho _{}\rho +\rho _{}^2\mathrm{exp}[A(\delta ^2\delta _{}^2)]\right)=0$$
(34)
with $`\rho _{}=0.6,A=25,\delta _{}=0.26`$ (see Figure 3, line) and use it in the polynomial expansion of the free energy density which enters the order parameter equation (20):

$$(\rho )=^\rho F(\rho ,\delta )𝑑\rho .$$
(35)
We also measured the density and the granular temperature of the grains as we decreased the shear force. These measurements could only be performed in the range $`\rho <0.55`$ since for larger $`\rho `$ the partially fluidized state is unstable. Note that that density of grains stays almost constant in a wide range of the order parameter $`\rho >0.1`$ (see Figure 4a). The granular temperature (which is defined as $`T=\stackrel{~}{u}_i\stackrel{~}{u_i}/2`$) normalized by the applied pressure $`P`$ appears to be a unique function of the order parameter (Figure 4b).

### V.2 Relaxation dynamics of the order parameter

We probed the relaxation dynamics of the order parameter by performing the following numerical experiment. The granular layer was prepared as in the previous Section. Lateral shear forces were increased adiabatically until the granular system reached a metastable solid (jammed) state within a hysteretic region. Then the layer was perturbed by applying random forces to a small randomly selected fraction of particles. The dynamics of the order parameter vary depending on the magnitude of the perturbation. Figure 5 shows an example of the evolution from the same jammed state for two different magnitudes of initial perturbation. Interestingly, the relaxation back to the jammed state is very fast, whereas the relaxation toward stable shear flow is much slower. We believe that it has to do with inertia of grains and the upper plate, so the intrinsic time scale of the order parameter relaxation is rather small, $`O(1)`$.

Unfortunately, our thin Couette flow system does not allow us to probe the local coupling of the order parameter since the order parameter is uniformly distributed throughout the system. In the absence of this data, this coupling was modeled by the linear diffusion term in (20) with the constant diffusion coefficient. As we will see in the following Section, this approximation indeed provides a good description for spatially non-uniform near-surface flow, however the value of the diffusion coefficient appears to be a function of local stress.

### V.3 Fitting the constitutive relation

The next step is to fit the constitutive relation from MD simulations. To this end, we use the same Couette flow simulations, but now we analyze the “fluid stress” $`\sigma _{\alpha \beta }^f`$ and the “static stress” $`\sigma _{\alpha \beta }^s`$ separately during our ramp-down simulations at three different values of $`P`$. Figure 6 shows a sample of these data for $`P=30`$. At large $`\delta =F/(L_xP)`$, when the order parameter $`\rho `$ is low, the total stress is dominated by the fluid component, but as the flow stops and $`\rho `$ approaches unity, the fluid stress turns to zero, and the total stress is equal to the static stress. Plotting $`\sigma _{yx}^f/\sigma _{yx}`$ as a function of the order parameter $`\rho `$ for different $`P`$ (Figure 7a), we observe that all data collapse onto a single curve which is well fitted by $`q(\rho )=(1\rho )^{2.5}`$. The lines in Figure 6 show the fit of the fluid and static stress tensors using (15),(16) with $`q(\rho )=(1\rho )^{2.5}`$.

The fluid as well as solid parts of the stress tensor are nearly symmetric, $`\sigma _{yx}^{f,s}=\sigma _{xy}^{f,s}`$, so the ratio $`\sigma _{xy}^f/\sigma _{xy}`$ is described by the same scaling function $`q(\rho )`$. On the other hand, the same procedure for the diagonal elements of the stress tensor yields a noticeably different scaling (see Figure7,b). Furthermore, a small but noticeable difference is evident between $`\sigma _{xx}^f/\sigma _{xx}`$ and $`\sigma _{yy}^f/\sigma _{yy}`$. More detailed analysis shows that in fact fluid parts of the diagonal components of the stress tensor $`\sigma _{xx}^f`$ and $`\sigma _{yy}^f`$ are nearly identical, and the difference is due to the solid part of the normal stresses (see Figure 8). This observation is consistent with the fact that the diagonal terms of the static stress tensor are determined by the details of the external loading. On the other hand, in a completely fluidized state the diagonal terms are all equal to the hydrodynamic pressure $`p_f`$ walton86 . In a partially fluidized regime, the diagonal terms of the shear stress can be expressed as

$$\sigma _{xx}=p_f/q_x(\rho ),\sigma _{yy}=p_f/q_y(\rho ).$$
(36)
Both functions $`q_{x,y}(\rho )`$ should approach 1 as $`\rho 0`$, but they may have different functional form to reflect the anisotropy of the static stress tensor. If the fluid pressure $`p_f`$ approaches zero in the solid state as $`p_f(1\rho )^\alpha `$, then $`q_{x,y}(1\rho ^{\beta _{x,y}})^\alpha `$, so $`\sigma _{xx,yy}\beta _{x,y}^\alpha `$. In our simple Couette flow, the diagonal stress tensor components can be well fitted by $`q_x(\rho )(1\rho )^{1.9}`$ and $`q_y(\rho )(1\rho ^{1.2})^{1.9}`$ (see Figure 7,b). So we observe that even in a partially fluidized regime, the “fluid phase” component indeed behaves as a real fluid with a well-behaved “partial” pressure $`p_f`$ which is zero in a solid state at $`\rho =1`$ and is becoming the full pressure in a completely fluidized state $`\rho =0`$.

Plotting the fluid shear stress versus the strain rate, we can test the validity of the Newtonian model for the stress-strain relation (14). Figure 9 shows $`\sigma _{yx}^f`$ vs $`\dot{\gamma }`$ for three different pressures $`P=20,30,40`$. At small $`\dot{\gamma }`$ all three lines are close to the same straight line $`\sigma _{yx}^f=12\dot{\gamma }`$, which indicates that the Newtonian scaling for fluid shear stress holds reasonably well. The deviations at large $`\dot{\gamma }`$ are evidently caused by variations of temperature and density in the dilute regime. Note that in contrast, the full shear stress does not goes to zero as $`\dot{\gamma }0`$ (see Figure 9b), so a viscosity coefficient conventionally defined as the ratio of the full shear stress to the strain rate diverges at the fluidization threshold as observed in Ref. gollub00 .

Combining the Newtonian law for the fluid stress-strain dependence with the order parameter scaling of the fluid stress tensor, we arrive at the relationship between the full stress tensor and the strain rate tensor (19) with $`\mu =12,q_x(\rho )=(1\rho )^{1.9},q_y(\rho )=(1\rho ^{1.2})^{1.9},q(\rho )=(1\rho )^{2.5}`$.

### V.4 Toward dilute granular flows - granular temperature revisited

While the above fittings have been made for the regime of a slow dense flow with $`\nu \nu _c`$, it is tempting to generalize the theoretical model so it smoothly crosses over to the standard kinetic continuum theory (1)-(10) for $`\rho 0`$. This goal can be achieved by including the equation of state (6) and the equation for the granular temperature (3) back into the theory. The important difference with respect to the standard kinetic theory is that the pressure which we calculate with Eq.(6) is not the total pressure, but the partial pressure associated with the fluid part of the stress tensor. Of course, as $`\rho 0`$, the static part of the stress tensor disappears, and the partial pressure becomes the total pressure.

We can test the relevance of this combined approach by calculating the spatial correlation function $`G(\nu )=(1+e)^1[\pi d^2p_f/4\nu T1]`$ with the values of fluid pressure $`p_f`$, temperature $`T`$ and density $`\nu `$ calculated in our testbed Couette flow at different external pressures $`P`$ and comparing it with the theoretical functions $`G_{CS}(\nu ),G_{CS}(\nu ),G_L(\nu )`$ (11)-(13). Figure 10 shows that the Carnahan-Sterling formula works very well in the dilute range $`\nu <0.67`$ as expected. In the high density regime $`G(\nu )`$ approaches the free volume result (12), and the overall dependence is in agreement with the global interpolation $`G_L(\nu )`$ by Luding Luding01 .

We can carry this analysis one step further and test the kinetic theory prediction for the shear viscosity. If we scale $`\sigma _{yx}^f`$ by the shear viscosity calculated using the kinetic formula (7) with globally fitted correlation function $`G_L(\nu )`$ and actual temperature and density values from corresponding runs, all three lines in Figure 9a collapse onto the same straight line dependence $`\mu _f^1\sigma _{yx}^f=\dot{\gamma }`$ (see Figure 9a, inset).

One more test of the kinetic theory predictions can be performed by analyzing the granular temperature as a function of the shear strain rate. According to Ref. savage98 , $`T^{1/2}=A\dot{\gamma }`$ for a plane parallel shear flow where $`A`$ is a material constant weakly dependent on volume fraction. Our numerical data for three different external pressures shown in Figure 11, are consistent with this scaling law although small deviations can be observed at small $`\dot{\gamma }`$.

We have not done a similar comparison for the bulk viscosity, thermal conductivity, and the energy loss, however the presented data strongly suggest that the constitutive relations for the fluid part of the stress tensor are well described by the standard granular hydrodynamics.

### V.5 Order parameter description of partially fluidized granular flows - Take 2

Let us summarize the equations of the continuum theory as specified on the basis of the 2D molecular dynamics simulation of the thin Couette flow.

The mass, momentum, and energy conservation conditions are expressed by Eqs. (1)-(3). The order parameter equation now has the form

$$\frac{D\rho }{Dt}=D^2\rho (\rho 1)\left(\rho ^22\rho _{}\rho +\rho _{}^2\mathrm{exp}[A(\delta ^2\delta _{}^2)]\right)$$
(37)
with $`\rho _{}=0.6,\delta _{}=0.26,A=25`$. As mentioned before, for the lack of simulation data we assume linear diffusion coupling of the order parameter with a constant non-dimensional diffusion coefficient $`D`$. The constitutive relation now reads

$$\sigma _{\alpha \beta }=[p_f/q_\alpha (\rho )+(\mu _f\lambda _f)\text{Tr}\dot{𝜸}/q(\rho )]\delta _{\alpha \beta }\mu _f\dot{\gamma }_{\alpha \beta }/q(\rho )$$
(38)
with $`q_x(\rho )=(1\rho )^{1.9},q_y(\rho )=(1\rho ^{1.2})^{1.9},q(\rho )=(1\rho )^{2.5}`$.

The equation of state and expressions for viscosity, thermoconductivity and the energy dissipation have the same functional form as (6)-(10), but they are now written for the “fluid” parameters $`p_f,\mu _f,\lambda _f`$, and $`\epsilon `$.

## VI Surface-driven shear granular flow under gravity

In this Section we apply the theoretical description which was formulated in the previous Section on the basis of numerical simulations of a thin Couette flow with no gravity to another model problem. We consider shear granular flow in a thick granular layer under gravity driven by the upper plate which is pulled in a horizontal direction, see Figure 12. A similar system has been studied experimentally by Nasuno et al. gollub98 , as well as by Tsai et al. Tsai02 .

We simulated up to 20,000 particles in a rectangular box under a heavy plate which was moved either with a constant speed $`V`$ or a constant force $`F`$. Periodic boundary conditions were assumed in a horizontal direction. After a transient, a quasi-stationary fluidization and shear flow established in the near-surface layer, while near the bottom grains remained in a nearly static jammed regime. Here, we show the results of several different runs with small/large pressure and small/large shear. The details of these runs are given in Table 1. The vertical profiles of the density, flow velocity, and the order parameter are shown in Figure 13. The density of grains remain nearly constant very close to the maximum random packing density value, except for a narrow near-surface boundary layer. The most significant density variations were observed for the case $`P10V50`$ corresponding to a light fast moving upper plate. The horizontal velocity decays roughly exponentially off the plate in agreement with experimental evidence gollub98 ; Tsai02 . The vertical profiles of the order parameter demonstrate a well-defined transition from fluid state near the upper plate to solid state below. The width of the “interface” grows with the applied pressure $`P`$, which indicates the stress dependence of the diffusion coefficient for the order parameter. Surprisingly, we found that the order parameter does not approach 1 at large depths, but instead seems to saturate at some value slightly below 1. We believe that this behavior is inherent to our 2D geometry with periodic boundary conditions on side walls. The moving upper plate oscillates vertically and produces vibrations in the bulk of granular layer. These slowly decaying with depth vibrations break weak contacts between particles which are not strongly pressed against each other (e.g. lying under arches). We believe that in full 3D simulations with more realistic boundary conditions this effect may be less pronounced. In principle, it may be included in the theoretical description by proper averaging of fluctuations of the strain tensor in the spirit of Savage savage98 .

By averaging velocity fluctuations and forces acting on individual particles, we calculated vertical profiles of the fluid and solid parts of the shear and normal stress components $`\sigma _{\alpha \beta }`$ (see Section IV). These profiles for Run P10V50 are shown in Figure 14. Strong fluctuations of the horizontal component of the stress tensor $`\sigma _{xx}`$ are mostly related to the static component of the tensor. According to Eqs. (15),(16), fluid and solid parts of the shear stress components should be related to the shear component of the full stress tensor (which in this geometry is roughly independent of $`y`$) via the function $`q(\rho )`$. Figure 15,a depicts this function as a parametric plot of $`\sigma _{yx}^f(y)/\sigma _{yx}`$ vs. $`\rho (y)`$ made using vertical profiles of stresses and the order parameter. As seen from this Figure, the same fit $`q(\rho )=(1\rho )^{2.5}`$ approximates the data quite well. However, unlike the zero-gravity case of the thin Couette flow, the normal stress components seem to be isotropic $`\sigma _{xx}=\sigma _{yy}`$ and they both are well described by $`q_y(\rho )`$ (see Figure 15,b).

The dependence of the fluid part of the shear stress component on the strain rate (Figure 16) shows the same behavior as for the thin Couette system: at small $`\dot{\gamma }`$ the viscosity is nearly constant $`\mu _f12`$, and at larger $`\dot{\gamma }`$ shear thinning is observed. Interestingly, the dependence of the local Reynolds shear stress on the local strain rate is well described by the Bagnold scaling $`\sigma _{yx}^R\dot{\gamma }_{yx}^2`$. Eventually, at large $`\dot{\gamma }`$, this scaling should dominate the full stress-strain rate relationship. Figure 17 compares the behavior of the viscosity coefficient $`\mu =\sigma _{yx}/\dot{\gamma }`$ and $`\mu _f=\sigma _{yx}^f/\dot{\gamma }`$ calculated along the vertical profiles of $`\sigma `$ and $`\dot{\gamma }`$, as a function of density $`\nu `$ and order parameter $`\rho `$. While the former diverges as $`\nu \nu _c`$ and $`\rho 1`$, the latter approaches the constant value $`\mu _f12`$, in agreement with results of Section V.

As in the previous Section, we can extract the particle-particle correlation function by calculating $`G(\nu )=(1+e)^1[\pi d^2p_f/4\nu T1]`$ using the vertical profiles of $`\sigma _{yx}^f,T`$, and $`\nu `$. Again, we obtain a good agreement with theoretical predictions based on the kinetic theory for fluid component of the stress tensor (Figure 18).

Finally, we can compare the stationary vertical profiles of the order parameter and the horizontal velocity with theoretical predictions. In most of our numerical simulations we specified the velocity of the upper plate rather than the applied force. That allowed us to study the regimes of slow dense flows which would be unstable had we applied a constant shear force. The shear stress tensor component $`\sigma _{yx}`$ in the stationary regime was indeed constant across the layer (see for example Figure14,c). However, due to slippage near the moving plate the relation between the plate speed and the shear stress is complicated. We do not address the issue of boundary conditions here as it is a subject of a separate study (see for example hui84 ). Here we simply use the values which are obtained in numerical simulations (the last column in Table 1), as parameters in our theoretical model.

In the stationary regime, the relevant stress tensor components are specified as follows:

$$\sigma _{yy}=P+(Hy),$$
(39)
$$\sigma _{yx}=\sigma _{xy}=const$$
(40)
where $`P`$ is the external pressure applied to the upper wall.

We also need to specify the boundary conditions for the order parameter at the top and bottom plates. This is a serious issue in its own right which we will address elsewhere. Here we simply impose no-flux boundary conditions both at the top and the bottom plate for the order parameter, $`_y\rho (0)=_y\rho (H)=0`$.

We limit ourselves with the case of slow dense flow regime, when the granular temperature plays a minor role, and the granular flow can be considered incompressible. This allows us to use the reduced set of equations (2),(37),(38) with the fixed viscosity $`\mu _f=12`$. The stationary shear flow solution of the continuum equations can be found numerically as follows. Since the components of the full stress tensor are assumed known, we solve the time-dependent order parameter equation (37) using the pseudo-spectral method until the solution reaches a stationary state. The resulting solution for the order parameter is then used to obtain the velocity profile by integrating the constitutive relation (38) from the bottom ($`y=0`$) up. Since the grains are strongly compressed near the rough bottom plate due to gravity, we assume the no-slip boundary condition for the horizontal velocity at $`y=0`$. The momentum conservation equation (2) is satisfied automatically. Thus obtained profiles of velocity and the order parameter were compared with our 2D molecular dynamics simulations.

The results of the comparison between the velocity and order parameter profiles obtained in simulation and by using the continuum theory are shown in Figure 19 for four runs P10V5 (a), P10V50 (b), P50V5 (c), P50V50 (d). The only fitting parameter used was the diffusion constant $`D`$ in the order parameter equation, which has not been determined in our testbed analysis. We used $`D=1`$ for runs P10V5 and P10V50, $`D=5`$ for P50V5 and $`D=10`$ for P50V50. From this we can conclude that the diffusion coefficient depends on the the local stress tensor, however more elaborate numerical experiments are needed to pinpoint this dependence more quantitatively. All other parameters were identical for all four cases as specified in Section V.5. The vertical profiles of the order parameter and the horizontal velocities are reasonably well described by the theory. However, for low pressure runs P10V5 and P10V50, the horizontal velocity profiles deviate from the numerical data presumably because the viscosity coefficient is no longer a constant in a dilute region near the top plate.

## VII Conclusions

In this paper we performed a series of numerical simulations of 2D wall-driven granular flows. These simulations were designed with a specific goal, to quantify the continuum theory of partially fluidized granular flows which was introduced ad hoc earlier AT1 ; AT2 . We defined the order parameter as a ratio of the number of static contacts to the total coordination number averaged over a small mesoscopic volume. Using simulations of a thin Couette flow between two rough plates, we determined the free energy density for the order parameter. Simulations confirmed that the ratio of the shear to the normal stress in the bulk of the granular flow can parametrize the stationary states of the order parameter equation. The same simulations allowed us to determine the detailed structure of the constitutive relation. We split the total stress tensor into the fluid and solid components, in which the former is comprised of the Reynolds stresses and the stresses transmitted through short-term collisions, while the latter is formed by the force chains through persistent contacts. The ratio of fluid and solid stress components is indeed determined by the order parameter through scaling functions $`q(\rho ),q_{x,y}(\rho )`$. Remarkably, the fluid component of the stress tensor is a linear function of the strain rate $`\dot{\gamma }`$ in the slow dense flow regime. This justifies the Newtonian scaling of the stress-strain relationship adopted in the theory.

Using the calibrated theory, we studied the flow structure of a thick surface driven Couette granular flow under gravity. We found the theoretical predictions to be in a good quantitative agreement with simulations.

The evidence presented here suggests an intriguing interpretation for the order parameter description of dense and slow granular flows. The granular material under shear stress appears to be similar to a multi-phase system with the fluid phase “immersed” into the solid phase. The fluid phase behaves as a simple Newtonian fluid for small shear rates when the density is almost constant, but exhibit shear thinning at larger shear rates when the density begins to drop. We observed that the Reynolds part of the fluid shear stress obeys the Bagnold scaling $`\sigma _{yx}^R\dot{\gamma }^2`$. We anticipate that for very large shear rates when the Reynolds stress becomes dominant the overall stress tensor should exhibit Bagnold scaling locally.

While this theory is primarily intended for dense and slow granular flows, we have shown that it can be combined with existing models of rapid granular flows based on the kinetic theory of granular gases. This requires to drop the assumption of incompressibility and include the equation for the granular temperature. Our simulations showed that the kinetic theory works well for the fluid part of the stress tensor in the whole range of densities from dilute regime to the critical random close packing density.

Many issues still remain open. The spatially non-uniform dynamics of the order parameter requires a more detailed study. We found that the diffusion constant postulated in Eq.(37) appears to be a function of the normal shear stress as well as the local strain rate, however, we do not have sufficient numerical data to provide a quantitative description of this dependence. It would be of interest to analyze the propagation of a fluidization front in a granular layer prepared in a meta-stable static regime. Such simulations could provide an insight into the mechanisms of the local coupling of the order parameter.

The molecular dynamics algorithm employed is based on a number of approximations. These approximations, however well tested and widely accepted wolf96 ; grest91 ; grest02 , directly affect the results of our fitting the continuum model. For example, if one replaces the Hookian model of particle interaction with a Hertzian one, an appreciable difference in the structure of the order parameter may be observed. More numerical work is needed to quantify the relationships between the microscopic parameters of the system (nature of collisions, restitution coefficient, friction, etc) and the parameters of the continuum model.

Finally, our simulations were limited by 2D systems, and of course the resulting continuum theory can only be directly applicable to 2D systems. While we anticipate that the structure of the model to remain in 3D systems, the specific form of the fitting functions should change. This future work will allow us to perform a comparison of the 3D model not only with numerical simulations, but also with experimental data.

The authors are indebted to B. Beringer, P. Cvitanovic, J. Gollub, T. Halsey, and J. Viñals for useful discussions, and to J. C. Tsai for sharing his unpublished experimental data. This work was supported by the Office of the Basic Energy Sciences at the US Department of Energy, grants W-31-109-ENG-38, and DE-FG03-95ER14516. Simulation were performed at the National Energy Research Scientific Computing Center.


format: text

text:
# The diffusive pair contact process and non-equilibrium wetting. 11footnote 1Unpublished notes intended as a basis for further research, disseminated exclusively on the cond-mat preprint server.

## 1 Introduction

One of the major challenges in non-equilibrium statistical physics is the classification of phase transitions from fluctuating phases into absorbing states . It is believed that the critical behavior of absorbing phase transitions can be associated with a finite number of universality classes. So far only few universality classes are firmly established, the most important ones being directed percolation (DP) , the parity-conserving (PC) class , voter-type transitions , and the general epidemic process . Searching for further universality classes the pair contact process with diffusion (PCPD)

$`2A3A`$ $`\text{with rate }\sigma `$
$`2A\mathrm{}`$ $`\text{with rate }\nu `$ (1)
diffusion of individual particles $`\text{with rate }D,`$
also called annihilation-fission process, is currently one of the most promising candidates as it exhibits a continuous phase transition with an unusual type of critical behavior which has not been seen before. These exceptional properties may be related to the fact that the PCPD is a binary spreading process, i.e., two particles have meet in order to generate offspring or annihilate.

The unusual critical behavior of binary spreading processes was first observed by Grassberger in 1982 . The problem was then rediscovered 15 years later by Howard and Täuber , who proposed a bosonic field theory for the 1+1-dimensional PCPD which turned out to be unrenormalizable. More recently Carlon et al  investigated a ‘fermionic’ lattice model of the PCPD model, in which the occupancy per site is restricted by an exclusion principle. Their paper trigerred a series of numerical and analytical studies  and released a debate concerning the asymptotic critical behavior at the transition. Currently several viewpoints are being discussed, stating that the PCPD

1. represents a new universality class with a unique set of critical exponents ,
2. represents two different universality classes depending on the diffusion rate ,
3. can be interpreted as a cyclically coupled DP and annihilation process ,
4. may be regarded as a marginally perturbed DP process with continuously varying critical exponents ,
5. may cross over to DP after very long time .

Each of these explanations has been supported to a different extent by physical arguments, mean field approaches, DMGR methods, and state-of-the-art simulations. The surprising variety of viewpoints demonstrates that the PCPD is a highly non-trivial process and that the resolution of these open questions is an exciting challange of non-equilibrium statistical physics.

The purpose of these notes is to point out that the Langevin equation of the PCPD is related to the problem of non-equilibrium wetting, leading to conclusions in favor of a slow crossover to DP. However, I would like to emphasize that the arguments presented here are partly speculative and need to be substantiated. Therefore these notes do not present fully validated results, rather they are intended as a basis for further research and discussions.

## 2 Why DP?

Currently most authors believe that the PCPD represents a new universality class. Depending on the model under consideration, it is observed that the asymptotic scaling regime is only reached after a long time of $`10^4\mathrm{}10^6`$ Monte Carlo steps. The estimates for $`\delta =\beta /\nu _{}`$ seem to be close to $`0.21`$, while the dynamic exponenent $`z1.7`$ is clearly smaller than $`2`$, indicating superdiffusive spreading at criticality.

This conjecture, however, poses a fundamental problem. As shown in Refs. , binary spreading processes are characterized by two different modes (or sectors) of spreading, namely, a high-density mode dominated by self-reproducing and annihilating pairs of particles, and a low-density mode of solitary diffusing particles. The interplay of the two modes in a critical binary spreading process is illustrated in Fig. 1, where pairs and solitary particles are represented as red and blue pixels, respectively. Plotting $`x/L^{1/2}`$ versus $`\mathrm{log}_{10}t`$ the figure covers four decades in time. As can be seen, patches of high activity (red) are connected by lines of diffusing solitary particles (blue). Obviously this interplay is present on all scales up to $`10^6`$ time steps.

The problem arises precisely at this point: Even after $`10^6`$ time steps the solitary particles perform simple random walks over large distances. However, such a random walk is always characterized by the dynamic exponent $`z=2`$, while the process as a whole spreads superdiffusively with $`z<2`$. Therefore, the effective diffusion constant for solitary particles has to vary slightly under rescaling, meaning that a cluster such as in Fig. 1 cannot be scaling-invariant. Therefore it seems that the process is still far away from the asymptotic scaling regime, even after $`10^6`$ time steps.

Another hint can be found in the paper by Noh and Park , who measured the life time distribution $`F(\tau )`$ of solitary particles in a critical binary spreading process, finding an approximate power-law behavior $`F(\tau )\tau ^\theta `$ with an exponent $`\theta =2.25(5)`$. Since this distribution decays faster than $`\tau ^2`$, their result would imply that the mean life time

$$\overline{\tau }=\frac{F(\tau )\tau 𝑑\tau }{F(\tau )𝑑\tau }$$
(2)
is finite, introducing a non-trivial time scale in the critical PCPD. The existence of such a typical life time indicates that the true asymptotic critical behavior may only be seen on extremely large scales, where $`\overline{\tau }`$ is virtually invisible. This scaling regime may be far beyond the accessible range of today’s numerical simulations.

Observing that the numerical estimates for the critical exponents seem to move in the direction of DP values with increasing numerical effort, I suggested that a very slow crossover to DP should not be ruled out . Very recently Carlon and Barkema  supported this point of view by a quantitative Monte Carlo and density matrix renormalization group study. In the present notes the DP hypothesis is supported in a completely different way by relating the PCPD to a non-equilibrium wetting process. However, as mentioned before, some of the arguments presented in the following are still speculative so that the hypothesis of an asymptotic DP behavior should be regarded as one out of many possible scenarios.

## 3 Langevin equation for the PCPD

In Ref.  the Langevin equation for the ($`d`$+1)-dimensional PCPD was derived rigorously by introducing a bosonic operator formalism and performing the continuum limit. Using a simplified notation this Langevin equation reads

$$\frac{}{t}\rho (\stackrel{}{x},t)=b\rho ^2(\stackrel{}{x},t)c\rho ^3(\stackrel{}{x},t)+D^2\rho (\stackrel{}{x},t)+\rho (\stackrel{}{x},t)\xi (\stackrel{}{x},t),$$
(3)
where $`\rho (\stackrel{}{x},t)`$ is a coarse-grained particle density and $`\xi (\stackrel{}{x},t)`$ denotes a white Gaussian noise with the correlations

$$\xi (\stackrel{}{x},t)\xi (\stackrel{}{x}^{},t^{})=\mathrm{\hspace{0.33em}2}\mathrm{\Gamma }\delta ^d(\stackrel{}{x}^{}\stackrel{}{x})\delta (t^{}t).$$
(4)
The four terms on the r.h.s. of Eq. (3) can be interpreted as follows. Dividing the discrete lattice of the PCPD into boxes which are much larger than the lattice spacing but much smaller than the system size, $`\rho (\stackrel{}{x},t)`$ may be understood as a coarse-grained average density of particles in a box at position $`\stackrel{}{x}`$. Assuming the particles in each box to be uncorrelated, the interplay of the binary reactions $`2A3A`$ and $`2A\mathrm{}`$ leads to a quadratic term $`b\rho ^2(\stackrel{}{x},t)`$, where $`b`$ is essentially determined by the difference $`\sigma 2\nu `$ of the two reaction rates. For so-called ‘fermionic’ models with an exclusion principle we added a cubic term $`c\rho ^3(\stackrel{}{x},t)`$ by hand which prevents the particle density in the active phase from diverging. Moreover, there is a diffusion term and a noise field accounting for density fluctuations.

Note that the amplitude of the noise in Eq. (3) is proportional to the density $`\rho (\stackrel{}{x},t)`$. This type of noise, which is known as multiplicative noise in the literature (see e.g. ), can be motivated as follows. Since the noise accounts for fluctuations of the particle density in each box, it is primarily generated by the binary reactions $`2A3A`$ and $`2A\mathrm{}`$ so that number of noise-generating sites in each box will be proportional to $`\rho ^2(\stackrel{}{x},t)`$. Thus, according to the central limit theorem, the total noise generated in the box is Gaussian and its intensity is expected to be proportional to $`\rho (\stackrel{}{x},t)`$.

Analyzing the Langevin equation by simple power-counting one can compute the mean field critical exponents and the upper critical dimension (see e.g. ). Neglecting diffusion and noise, the homogeneous stationary solution is $`\rho =b/c`$, hence the mean field critical point is $`b_c=0`$. According to the standard scaling theory of absorbing phase transitions, invariance under rescaling yields the mean-field critical exponents

$$\beta ^{MF}=1,\nu _{}^{}{}_{}{}^{MF}=1,\nu _{}^{}{}_{}{}^{MF}=2$$
(5)
and the upper critical dimension

$$d_c=2.$$
(6)
For $`d>d_c`$ the coefficient $`\mathrm{\Gamma }`$ scales to zero, meaning that the noise will be irrelevant on large scales so that the critical exponents are given by their mean field values (5). In fact, recent high-precision simulations in two spatial dimensions  confirm this prediction for various values of the diffusion rate. In $`d<d_c`$ dimensions, however, fluctuation effects lead to a non-trivial critical behavior. For this reason the present study is restricted to the (1+1)-dimensional case.

Let us first recall the main results of Ref. . For the unrestricted PCPD, where the cubic term is absent, the bare coefficient $`b`$ and the noise amplitude $`\mathrm{\Gamma }`$ are related to the reaction rates in Eq. (1) by

$$b=\sigma 2\mu ,\mathrm{\Gamma }=2\sigma \mu .$$
(7)
Remarkably, in the unrestricted PCPD the critical point is always $`b_c=0`$, even in the presence of fluctuation effects below the upper critical dimension. This implies that the noise amplitude $`\mathrm{\Gamma }=3\sigma /4`$ is positive at criticality so that the transition is characterized by ‘real’ noise in the sense of Ref. . Moreover, the average particle density at the critical point was found to be constant. Regarding numerical simulations we note that this observation ubiquitously requires the Langevin equation to be iterated in the Ito sense, i.e., the density-dependent prefactor of the noise has to be evaluated before the update is carried out. Using the Stratonovich scheme one would have to introduce an additional linear drift term $`\mathrm{\Gamma }\rho (\stackrel{}{x},t)`$ in Eq. (3).

In the inactive phase $`b<0`$ the annihilation process $`2A\mathrm{}`$ dominates so that $`\rho (t)`$ decays algebraically. In this case the noise amplitude is expected to become negative in the renormalization group sense, i.e., the system crosses over to ‘imaginary’ noise after some time. For $`b>0`$, however, the particle density grows without limit and diverges exponentially.

## 4 The PCPD as a non-equilibrium wetting process

Following Ref.  we perform a Cole-Hopf transformation

$$h(\stackrel{}{x},t)=\mathrm{ln}\rho (\stackrel{}{x},t)$$
(8)
which maps Eq. (3) to

$$\frac{}{t}h(\stackrel{}{x},t)=\mathrm{\Gamma }be^{h(\stackrel{}{x},t)}+ce^{2h(\stackrel{}{x},t)}+D^2h(\stackrel{}{x},t)D[h(\stackrel{}{x},t)]^2+\xi (\stackrel{}{x},t)$$
(9)
where $`\xi (\stackrel{}{x},t)`$ is a non-multiplicative real Gaussian noise with the same correlations as in Eq. (4). Note that the constant drift term $`\mathrm{\Gamma }`$ is a consequence of the Ito interpretation in Eqs. (3)-(4). Introducing the notations $`\lambda =2D`$, $`a=\frac{\lambda }{2}(h)^2`$, and $`v_0=\mathrm{\Gamma }a`$ this equation may be rewritten as a KPZ equation

$$\frac{}{t}h(\stackrel{}{x},t)=a\frac{\delta V[h(\stackrel{}{x},t)]}{\delta h(\stackrel{}{x},t)}+D^2h(\stackrel{}{x},t)+\frac{\lambda }{2}[h(\stackrel{}{x},t)]^2+\xi (\stackrel{}{x},t)$$
(10)
in a potential

$$V[h]=\frac{ce^{2h}}{2}be^hv_0h,$$
(11)
which has been studied recently in the context of non-equilibrium wetting . Note that the constant $`a`$ has been defined in such a way that it compensates the average drift caused by the KPZ nonlinearity so that $`v_0`$ can be regarded as the average velocity of a freely evolving interface.

Let us now turn to the question how the order parameters are related in both cases. The order parameter of the PCPD is the density of particles $`\rho (\stackrel{}{x},t)`$. According to Eq. (8) the corresponding order parameter in the wetting process is

$$\rho (\stackrel{}{x},t)=\mathrm{exp}[h(\stackrel{}{x},t)].$$
(12)
Since it is known from numerical simulations of the PCPD that higher moments of the density $`\rho ^n(t)`$ scale in the same way as $`\rho (t)`$, we may approximate the exponential function by a step function

$$\mathrm{exp}[h]\{\begin{array}{c}1\text{ if }0<h1\hfill \\ 0\text{ if }h>1\hfill \end{array}$$
(13)
which – in a model with discrete heights – is essentially the density of sites at zero height

$$n_0(\stackrel{}{x},t)=\delta _{h(\stackrel{}{x},t),0}.$$
(14)
Comparing non-equilibrium wetting and the PCPD we therefore expect both quantities $`\rho (\stackrel{}{x},t)`$ and $`n_0(\stackrel{}{x},t)`$ to exhibit essentially the same type of asymptotic scaling behavior. Roughly speaking, the sites where the interface touches the substrate can be regarded as the active sites of the PCPD.

## 5 Interpretation of the phase transition in the unrestricted PCPD

In the case of the unrestricted PCPD, where multiple occupancy per site is allowed, the cubic term in Eq. (3) vanishes so that the field theory of Ref.  applies. As mentioned before it was shown that the transition takes place at $`b=0`$ even below the upper critical dimension. Moreover, the density of particles at criticality was found to be constant.

Interpreting the PCPD as a wetting process these results are easy to understand. In the inactive phase $`b<0`$ the potential $`V[h]=be^h`$ may be regarded as a lower wall representing a hard-core substrate on which the wetting layer is deposited (see Fig. 3). As shown in  the presence of a lower wall leads to a continuous wetting transition with a critical point where the propagation velocity $`v_0`$ of a freely evolving interface is zero. Apparently the mapping ensures that after renormalization this velocity vanishes automatically, i.e., the unrestricted PCPD is mapped onto the phase transition line of the corresponding wetting problem. Therefore, starting with a flat interface at $`h=0`$ (corresponding to a fully occupied lattice in the PCPD) the interface is neither pinned nor does it propagate uniformly, rather it roughens close to the wall.

In the active phase $`a>0`$ the potential is simply turned upside down so that an exponentially increasing force pulls the interface downwards, corresponding to a rapidly increasing particle density in the PCPD (see Fig. 3). Therefore, in the unrestricted PCPD the transition results from a changing sign in the potential, turning the repulsive force into an attractive one, hence the transition takes place exactly at $`b_c=0`$. Obviously this mechanism works in any dimension and does not depend on fluctuation effects.

## 6 Interpretation of the phase transition in the restricted PCPD

Adding a cubic term with $`c>0`$ in the Langevin equation (3) the particle density in the active phase does no longer diverge. Such a cubic term emerges, e.g., in ‘fermionic’ lattice models with an exclusion principle, where multiple occupancy per site is forbidden. The cubic term can also be implemented in models with unrestricted occupancy per site by choosing the update rule in such a way that the effective fission rate decreases with increasing particle density .

Although the restricted PCPD still exhibits a phase transition, its physical properties are very different:

* In contrast to the unrestricted case the particle density at criticality is no longer constant, instead it decays slowly, probably as a power law with strong corrections.
* In the active phase two different stationary states coexist, namely, the absorbing state (with $`0`$ or $`1`$ particles) and a fluctuating state with a finite density of particles.

Mapping the restricted PCPD to a wetting process we expect that the presence of a cubic term does not change the sign of the renormalized noise amplitude, i.e., right at the transition we are still dealing with ‘real’ noise so that the Cole-Hopf transformation remains valid. As shown in Fig. 3 the cubic term gives rise to an additional potential well at zero height. This potential well may be interpreted as an attractive short-range force between substrate and wetting layer . As a main result it was observed that such a force, if strong enough, may turn the continuous wetting transition into a discontinuous one. Moreover, in those parts of the phase diagram, where the coefficient $`\lambda `$ of the KPZ nonlinearity is negative, an extended region emerges, where the bound and the moving phase coexist. The main motivation of the present notes is to relate this type of phase coexistence in the wetting process with the aforementioned coexistence of fluctuating and absorbing states in the corresponding PCPD.

The phase coexistence observed in non-equilibrium wetting works as follows. Keeping $`c>0`$ fixed the parameter $`b`$ controls essentially the depth of the potential well. If $`b`$ is sufficiently small the transition is not affected, i.e., it is still continuous and takes place at $`v_0=0`$ (see left panel of Fig. 4). However, if the potential well is deep enough, i.e., if $`b`$ exceeds a certain critical threshold $`b^{}`$, the transition becomes first order and a phase coexistence region emerges in those parts of the phase diagram where $`\lambda <0`$, as shown in the right panel of Fig. 4. Within this region the short-range force is strong enough to bind the interface to the substrate although a freely evolving interface far away from the wall would already advance with the velocity $`v_0>0`$.

The mechanism, which ensures the stability of both phases in the thermodynamic limit, can be understood as follows (see Ref. ): If the interface detaches partly from the substrate due to a large fluctuation, it first advances because of $`v_0>0`$. The island continues to grow until its edges reach a certain critical slope, from where on the negative KPZ-nonlinearity suppresses further growth. The resulting pyramidial island then shrinks laterally at constant pace and eventually disapears (see Fig. 5). The velocity at which the island shrinks is maximal at the wetting transition line $`v_0=0`$ and tends to zero at the upper boundary of the phase coexistence region. Obviously, this mechanism requires $`\lambda `$ to be negative. We note that in Eq. (10) this is indeed the case.

Regarding the Cole-Hopf transformation we now postulate that for $`c>0`$ the PCPD corresponds to a wetting process for which the velocity of a freely evolving interface $`v_0`$ is positive. Only then the absorbing state of the PCPD, which corresponds to a completely detached interface far away from the wall, is thermodynamically stable. The postulate $`v_0>0`$ should emerge from a renormalization group calculation, probably due to a shift of the noise amplitude when the cubic term is introduced. It implies that the restricted PCPD corresponds to a point above the horizontal wetting transition line in Fig. 4. Depending on the value of $`b`$ – the critical parameter of the PCPD – this point is located either inside the phase coexistence region or above. We conjecture that the coexistence region corresponds to the active phase of the PCPD and that the phase transition takes place at its upper boundary (the green line in Fig. 4).

## 7 Critical properties at the borderline of phase coexistence

Approaching the upper boundary of the phase coexistence region the attractive short-range force becomes so weak that the velocity at which the islands shrink tends to zero. Consequently the average size of the islands in the stationary state increases, whereas the typical slope of their edges remains almost the same.

A schematic illustration of a typical interface configuration in the coexistence region close to the upper boundary is shown in Fig. 6. As can be seen, the particle density $`\rho =e^h`$ in the corresponding PCPD is indeed proportional to the density of interface sites at zero height. To understand the critical properties of the PCPD, it is therefore essential to analyze the dynamics of bottom layer sites in the corresponding wetting problem at the upper boundary of the coexistence region. Fig. 7 shows a spatio-temporal plot of the wetting model described in Appendix A, where the bottom layer sites are represented as black pixels. If the previous assumptions are correct, these pixels should display essentially the same critical behavior as the active sites of the PCPD. Therefore, the question arises to what extent the critical dynamics of the black pixels is universal.

A possible answer is given in a recent preprint by Muñoz and Pastor-Satorras . Considering the problem of synchronization transitions in extended coupled maps they are led to exactly the same Langevin equation as in Eq. (10). Discretizing space-time and analyzing the critical behavior at the upper boundary of the coexistence region they find numerical evidence of a directed percolation transition. This observation, together with the postulates of the present work, would imply that for $`b>b^{}`$ the transition of the one-dimensional PCPD belongs to the DP class.

Using the language of non-equilibrium wetting the conclusion by Muñoz and Pastor-Satorras seems to be reasonable. The evolution of the interface is in fact dominated by the dynamics of sites at zero height. As in DP, they can spontaneously generate offspring and disappear. The large islands between those sites do not mediate effective long-range interactions, instead they seem to follow passively the dynamics at the bottom layer, adjusting their size quickly whenever their base grows or shrinks. Once the interface detaches completely, it advances with constant velocity $`v_0>0`$, meaning that the PCPD has entered the absorbing state.

In addition the authors of Ref.  find a regime of first-order transitions. However, such a regime can only exist in those parts of the phase diagram where $`\lambda `$ is positive. In the present case, where $`\lambda `$ is negative, the transition belongs either to the class of multiplicative noise $`(b<b^{})`$ or to DP $`(b>b^{})`$, separated by a tricritical point at $`b=b^{}`$ (see Ref. ). Nevertheless one may observe a transient first-order behavior for $`\lambda <0`$, which crosses over to DP after very long time.

## 8 Conclusions

In these notes I have propsed a relation between the diffusive pair contact process and non-equilibrium wetting, arriving at the following conclusions:

1. Roughly speaking the particles in the PCPD correspond to interface sites and zero height in the corresponding wetting problem.
2. The transitions in the unrestricted (bosonic) and the restricted (fermionic) PCPD rely on fundamentally different mechanisms: In the unrestricted case the transition in the corresponding wetting process is caused by a simple change of sign in the potential, whereas in the restricted case the transition emerges at the borderline of a region of phase-coexistence.
3. Referring to recent simulation results by Muñoz and Pastor-Satorras one can conclude that for $`b>b^{}`$ the phase transition in the one-dimensional PCPD belongs to DP universality class, whereas for $`b<b^{}`$ it belongs to the usual continuous wetting (multiplicative noise) universality class. Thus there may be an open door for the scenario of two universality classes depending on the models parameters.

Various question arise:

1. Where is the counterpart of the diffusive background of solitary particles?
This is probably the weakest point of the continuum approach discussed here. Obviously the continuum description does not account for single particles. Here it would be useful to study the difference between discrete and continuum models in more detail, as it has been done in the context of synchronization transitions.
2. Does the wetting model reproduce the algebraic decay inactive phase?
No, instead one obtains an exponential decay. This failure may be related to the fact that the PCPD crosses over from ‘real’ to ‘imaginary’ noise when it enters the annihilation-dominated regime. By contrast the noise in the wetting model is always real.
3. Suppose that the PCPD belongs to the DP class, why is $`d_c=2`$ and not $`4`$ ?
For $`d>2`$ the nonlinear term in the KPZ equation is irrelevant (unless it is very large) and thus the phase coexistence region in the wetting problem no longer exists. Hence the DP regime cannot be accessed in $`d>2`$ and the model falls into the universality class of multiplicative noise.

Clearly the ideas presented in these notes are still speculative. Firstly, it is assumed that the Langevin description of the PCPD (including the cubic term) is valid and that the noise at criticality is Gaussian and real. Moreover, it is assumed that the Cole-Hopf transformation can be applied as usual. Finally, we postulate that the parameters renormalize in such a way that for $`c>0`$ the velocity of a freely evolving interface is positive. All these assumptions have to be verified. However, in my opinion the main problem of the continuum description is the missing notion of ‘solitary particles’. As single particles play an important role in the PCPD in restarting avalanches of high activity, it may well happen that one of the main features of the model, namely, the discrete nature of the particle density, is lost by introducing the Langevin equation. Nevertheless it is interesting to see that both problems, the PCPD and non-equilibrium wetting, are closely related and I hope that these notes may stimulate further research in this direction.

Acknowledgment: I would like to thank D. Mukamel for fruitful discussions.

## Appendix A A minimal model for non-equilibrium wetting

The probably simplest model for non-equilibrium wetting, which may be regarded as a realization of Eqs. (10)-(11), has been introduced some time ago in Ref. . The model is defined as a restricted solid-on-solid deposition-evaporation process of a growing interface in which the substrate is implemented as a hard-core wall at zero height. The dynamic rules involve three different elementary processes (see Fig. 8), namely,

> * deposition of atoms on the substrate at rate $`q_0`$,
> * deposition of atoms on top of islands at rate $`q`$,
> * evaporation from the edges of islands at rate $`1`$, and
> * evaporation from the middle of plateaus at rate $`p`$.

The reduced growth rate $`q_0`$ at the bottom layer accounts for the attractive short-range force between substrate and wetting layer and determines the depth of the potential well.

The phase diagram of the model without attractive force $`q_0=q`$ is shown in the left panel of Fig. 9. The moving phase and the bound phase are separated by a second-order transition line, where the velocity of a freely evolving interface vanishes. The line has two special points. For $`p=0`$ the model exhibits a special critical behavior since the transition is driven by a DP process at the bottom layer, as discussed in Refs. . Another special transition point is located at $`p=q=1`$, where the dynamic rules and are symmetric under reflection $`hh`$ so that the nonlinear term in the KPZ equation vanishes, corresponding to an Edwards-Wilkinson (EW) equation in a potential. Comparing the velocities of a horizontal and an artificially tilted interface it is possible to determine a line where the effective coefficient $`\lambda `$ of the nonlinear term vanishes (shown as a dotted line in Fig. 9). As expected this line intersects the phase transition line at $`p=q=1`$.

Introducing an attractive short-range force between substrate and surface layer by lowering $`q_0`$ the critical point $`q_c(p)`$ is not changed, i.e., the transition line remains the same. However, if $`q_0`$ is smaller than a certain threshold $`q_0^{}(p)`$ the transition may become first order (bold line in the right panel of Fig. 9). The first order line ends in a tricritical point (denoted as TCP) which moves along the phase transition line as $`q_0`$ is varied. In those parts of the phase diagram, where $`p<1`$ (i.e., $`\lambda <0`$), a phase coexistence region emerges above the transition line.

The parameter $`q`$ controls the growth rate and may be associated with $`v_0`$, while $`q_0`$ determines the strength of the short range force which is related to the parameter $`b`$ in the Langevin equation. The parameter $`p`$ can be used to control the effective value of $`\lambda `$, the coefficient of the nonlinear term in the KPZ equation.

References: