text:
# 1 Introduction

## 1 Introduction

One of the strong arguments in favour of supersymmetry is that R-parity conserving supersymmetric models have a cold dark matter candidate(CDM), the lightest supersymmetric particle(LSP). The preferred candidate is a neutralino. The contribution of the LSP to the relic density is however very model dependent and varies by several orders of magnitude over the whole allowed parameter space of the mininal supersymmetric standard model (MSSM). The measurements of the relic density then imposes stringent constraints on the parameters of the MSSM often favouring solutions with light supersymmetric particles. There are basically two mechanisms that can significantly reduce the estimate of the relic density leading to acceptable values while having a not so light supersymmetric spectrum: annihilation via a s-channel resonance and coannihilations where the LSP interacts with slightly heavier sparticles. Special care must be taken to treat these two cases carefully if one wants to make predictions for the relic density of CDM at the few percent level. Although, at present, the generally assumed range for the relic density of CDM is $`.1<\mathrm{\Omega }h^2<.3`$ , the recent data from BOOMERANG coupled with some constraints already indicates a more restrictive range at 2$`\sigma `$, $`.11<\mathrm{\Omega }h^2<.15`$. As new measurements will be performed in the near future, improvements over the present limits are expected.

There exist many calculations of the relic density in the MSSM using various approximations both for the evaluation of the thermally averaged cross-section and for solving the Boltzmann equation for the density of dark matter particles \[3-8\]. Among these, Neutdriver and DarkSusy are publicly available. Our purpose was to provide a tool that evaluates with high accuracy the annihilation cross-sections even in regions near poles and thresholds, that is both flexible and upgradable and that goes beyond DarkSusy as far as the calculation of the relic density is concerned.

The first calculations of the relic density used an expansion of the annihilation cross-section in power series of neutralino velocities, $`\sigma =a+bv^2`$. While this approximation works well in large regions of parameter space, it fails when annihilation through s-channel resonance is important. When the masses are such that the neutralinos can annihilate through a s-channel Higgs resonance  or a s-channel Z resonance, the annihilation rate increases significantly often bringing the relic density in an acceptable range as one gets close to the s-channel pole. In fact as one gets very close to the pole, the annihilation rate becomes so fast that the neutralinos cannot constitute the only source of dark matter. These effect are especially important in models with large $`\mathrm{tan}\beta `$ where the heavy Higgs resonances are very wide. The proper relativistic formalism for the evaluation of the thermally averaged cross-section was then introduced by . The relativistic formalism, generalized to the case of coannihilations , was implemented in the codes of  and , for the case of gaugino coannihilations. We follow basically this formalism, although contrary to DarkSusy, we still rely on approximations for the solution of the relic density equations and the determination of the freeze-out temperature. This allows to significantly increase the speed of the program and proves to be very useful when scanning over a large parameter space.

Coannihilation processes where the LSP interacts with slightly heavier sparticles can occur in principle with any supersymmetric particle, although in SUGRA models, the most common coannihilations are with gauginos or right-handed sleptons. The importance of the coannihilation channels were emphasized before both for gauginos , sleptons or stops . In micrOMEGAs we include ALL coannihilation channels, in all more than 2800 processes not counting charged conjugate processes. The tree-level cross-sections are calculated exactly including the full set of diagrams contributing to each process. The calculations of cross-sections are based on CompHEP, an automatic program for the evaluation of tree-level Feynman diagrams. Furthermore we include also some higher order effects, namely the two-loop corrections to the Higgs masses and the one-loop corrections to the Higgs widths. The latter turns out to be important in the large $`tan\beta `$ region. Although we will generally assume that the neutralino is the LSP, micrOMEGAs can be used to compute the relic density with any supersymmetric particle as the LSP, in particular the sneutrino. This is because all (co-)annihilation of any pairs of supersymmetric particles into any pairs of standard model or Higgs particles are included.

The main characteristics of micrOMEGAs, are

* Complete tree-level matrix elements for all subprocesses
* Includes all coannihilation channels with gauginos, sleptons and squarks.
* Loop-corrected Higgs masses and widths
* Speed of calculation

After the important equations for the calculation of the relic density are summarized, we give a short description of the parameters of the supersymmetric model and of the package. Finally we present some results and comparisons with another program in the public domain, DarkSusy.

## 2 Calculation of the relic density

The relic density at present is calculated from

$`\mathrm{\Omega }_{\stackrel{~}{\chi }_1^0}h^2=2.755\times 10^8{\displaystyle \frac{m_{\stackrel{~}{\chi }_1^0}}{GeV}}Y_0`$ (2.1)
where $`Y_0`$ is the abundance of the LSP today. To find $`Y_0=Y(T=T_0)`$, one needs to solve the evolution equation for $`Y`$

$$\frac{dY}{dT}=\sqrt{\frac{\pi g_{}(T)}{45G}}<\sigma v>(Y^2Y_{eq}^2)$$
(2.2)
$`g_{}(T)`$ is a degrees of freedom parameter derived from the thermodynamics describing the state of the universe  and $`Y_{eq}=Y_{eq}(T)`$ represents the thermal equilibrium abundance

$$Y_{eq}(T)=\frac{45}{4\pi ^4h_{eff}(T)}\underset{i}{}g_i\frac{m_i^2}{T^2}K_2(\frac{m_i}{T})$$
(2.3)
where we sum over all supersymmetric particles $`i`$ with mass $`m_i`$ and $`g_i`$ degrees of freedom. $`K_n`$ is the modified Bessel function of the second kind of order $`n`$ . Note that $`Y_{eq}`$ falls rather rapidly as the temperature decreases. $`<\sigma v>`$ is the relativistic thermally averaged annihilation cross-section

$$<\sigma v>=\frac{\underset{i,j}{}g_ig_j\underset{(m_i+m_j)^2}{}𝑑s\sqrt{s}K_1(\sqrt{s}/T)p_{ij}^2\sigma _{ij}(s)}{2T\left(\underset{i}{}g_im_i^2K_2(m_i/T)\right)^2},$$
(2.4)
where $`\sigma _{ij}`$ is the total cross section for annihilation of a pair of supersymmetric particles into some Standard Model particles, and $`p_{ij}`$ is the momentum of the incoming particles in their center-of-mass frame. The summation is over all supersymmetric particles. Integrating Eq. 2.2 from $`T=\mathrm{}`$ to $`T=T_0`$ would lead $`Y_0`$.

Although one can solve for $`Y`$ numerically, the procedure is extremely time consuming especially when scanning over a large parameter space and when we include a great number of processes. It is therefore important to seek as good an approximation as possible to speed up the code. We follow the usual procedure of defining a freeze-out temperature $`T_f`$. which can be extracted by solving iteratively

$$\frac{dln(Y_{eq})}{dT}=\sqrt{\frac{\pi g_{}(T)}{45G}}<\sigma v>Y_{eq}\delta (\delta +2)$$
(2.5)
where $`\delta `$ is some small constant number. The freeze-out temperature is defined from $`Y_f=Y(T_f)=(1+\delta )Y_{eq}(T_f)`$. Starting from a typical value $`T_f=m_{\stackrel{~}{\chi }_1^0}/25`$ only a few iterations are necessary to find a solution to this equation with $`\delta =1.5\pm .2`$.

In the second regime, where $`YY_{eq}`$, one can neglect $`Y_{eq}^2`$ completely. One finds

$`{\displaystyle \frac{1}{Y(0)}}={\displaystyle \frac{1}{Y_f}}+\sqrt{{\displaystyle \frac{\pi }{45G}}}{\displaystyle _{T_0}^{T_f}}g_{}^{1/2}(T)<\sigma v>𝑑T,`$ (2.6)
Furthermore we find that the solution (2.6) does not depend significantly on $`\delta `$.

Typical freeze-out temperatures vary between 1GeV and 10GeV, in this temperature range, the $`h_{eff}`$ and $`g_{}`$ functions can vary by about 20% . To achieve an accuracy better than 10%, one cannot use a constant value for these functions. We use the numerical tables of the DarkSUSY package  for a precise evaluation of $`h_{eff}`$ and $`g_{}`$.

### 2.1 Numerical integration and summation

In order to find $`T_f`$ by solving Eq. 2.5, we have to evaluate several integrals and perform a summation over different annihilation channels. In the evaluation of the thermally averaged cross-section we have included all two-body subprocesses involving two LSP’s, the LSP and a co-annihilating SUSY particle (all the particles of the MSSM) as well as all subprocesses involving two coannihilating SUSY particles. The final states include all possible standard model and Higgs particles that contribute to a given process at tree-level. The total number of processes exceeds 2800 not including charged conjugate processes. In practice processes involving the heavier SUSY particles contribute only when there is a near mass degeneracy with the LSP since there is a strong Boltzmann suppression factor $`B_f`$ in Eq. 2.4.

$$B_f=\frac{K_1((m_i+m_j)/T_f)}{K_1(2m_{\stackrel{~}{\chi }_1^0}/T_f)}e^{\frac{(m_i+m_j2m_{\stackrel{~}{\chi }_1^0})}{T_f}}$$
(2.7)
where $`m_i,m_j`$ are the masses of the incoming particles. To speed up the program a given subprocess is removed from the sum (2.4) if the total mass of the incoming particles is such that $`B_f`$ is below some limit, $`B_ϵ`$ defined by the user. Alhtough it would be sufficient to use $`B_ϵ=10^2`$ to give a precision of $`1\%`$ when $`\sigma _{coan}\sigma _{\stackrel{~}{\chi }_1^0\stackrel{~}{\chi }_1^0}`$, we use a more restrictive value, $`B_ϵ=10^6`$, to allow for cases where $`\sigma _{coan}\sigma _{\stackrel{~}{\chi }_1^0\stackrel{~}{\chi }_1^0}`$. This can occur for example for coannihilation processes with squarks which depend on $`\alpha _s`$, for processes with poles or in some regions of parameter space where $`\sigma _{\stackrel{~}{\chi }_1^0\stackrel{~}{\chi }_1^0}`$ is suppressed.

In our program we provide two options to do the integrations, the fast one and the accurate one. The fast mode already gives a precision of about 1% which is good enough for all practical purposes. In the accurate mode the program evaluates all integrals by means of an adaptative Simpson program. It automatically detects all singularities of the integrands and checks the precision. In the case of the fast mode the accuracy is not checked. We integrate the squared matrix elements over $`\mathrm{cos}\theta `$, the scattering angle, by means of a 5 point Gauss formula. For integration over $`s`$, Eq. 2.4 we use a restricted set of points which depends whether we are in the vicinity of a s-channel Higgs(Z,W) resonance or not. We increase the number of poins if the Boltzmann factor corresponding to $`m_{pole}`$ is larger than $`0.01B_ϵ`$.

## 3 MSSM parameters, Higgs mass and widths

The input parameters are the ones of the soft SUSY Lagrangian defined at the weak scale, using the same notation as in CompHEP/SUSY models. In the model used, the masses of fermions of the first generation are set to zero. The masses of the quarks of the second generation are also set to zero, so that there is no mixing of the squarks of the first two generations. However we have kept the mass for the muon as well as the trilinear coupling $`A_\mu `$ as this is relevant for the calculation of the muon anomalous magnetic moment. Although first generation sleptons are pure left/right states, they are properly ordered according to their masses so that the correct coannihilating particle corresponding to the lightest slepton is taken into account. The sign convention for the parameters $`\mu `$ and $`A`$, is given explicitly in .

The calculation of the Higgs masses are done with FeynHiggsFast. For the Higgs widths, it is necessary to take into account the QCD corrections to the partial widths, $`h(H,A)b\overline{b}`$. This is particularly important at large $`\mathrm{tan}\beta `$ where the $`b\overline{b}`$ mode is the dominant one. The QCD corrections can be very large for heavy Higgses, easily a factor of two above the tree-level value for a Higgs of $`1`$TeV, due mostly to the running of the quark mass at high scale. To take these corrections into account we have redefined the vertices $`hq\overline{q},Hq\overline{q}`$ and $`Aq\overline{q}`$ using an effective mass that reproduces the radiatively corrected width. We used HDECAY to produce a table of mass-dependent QCD-corrected Higgs partial widths. From this, the effective quark masses $`m_b(m_H)`$ are extracted and simple interpolation can reproduce the one-loop corrected width for any value of $`m_H`$. In the region of physical interest, the precision on the width is at the per-mil level safe for the neutral Higgses partial widths near the $`t\overline{t}`$ threshold. However, this region does not contribute significantly to the neutralino cross-section. For the charged Higgs, one can extract from HDECAY both an effective $`m_t`$ as well as an effective $`m_b`$ using high and low $`tan\beta `$ values. This way we could reproduce at better than 1% level the partial width $`H^+t\overline{b}`$. The effective quark mass $`m_b(Q)`$ evaluated at $`Q=m_1+m_2`$, the sum of the initial masses, is used by default in the $`(h,H,A)b\overline{b}`$ vertices. This will lead to the correct result when the contribution of the Higgs resonance is very important ($`M_H2m_{\stackrel{~}{\chi }_1^0}`$).

## 4 Description of micrOMEGAs

micrOMEGAs is a C program that also calls some external FORTRAN functions. micrOMEGAs
relies on CompHEP for the definition of the parameters and the evaluation of all cross-sections. The program is contained in a package that lets the user choose between weak scale parameters or parameters of SUGRA models as input parameters. The latter is achieved through a link with ISASUGRA/Isajet. The package can be obtained at
http://wwwlapp.in2p3.fr/lapth/micromegas

As we have already mentionned, due to a strong Boltzmann suppression factor, only a small fraction of the available processes are needed, those with a sparticle close in mass to the LSP. In principle, compilation of the full set of subprocesses is possible, but such a program would be huge and could not be distributed easily. To avoid this problem, we include in our package the program CompHEP which generates, while running, the subprocesses needed for a given set of MSSM parameters. The generated code is linked during the run to the main program and executed. The corresponding “shared” library is stored on the user disk space and is accessible for all subsequent calls, thus each process is generated and compiled only once. Such approach can be realized only on Unix platforms which support dynamic linking.

The complete list of input soft SUSY parameters at the weak scale can be found in . When using the SUGRA option, the 7 input parameters are the usual 5 parameters of a SUGRA model, defined at the GUT scale, $`m_0,m_{1/2},A_0,\mathrm{tan}\beta ,sign(\mu )`$. In addition one must specify the value of the top mass as well as the the Isajet option (model $`=1,2`$ ) for a SUGRA model with or without gauge coupling unification at the GUT scale. When this option is used, the value of the weak scale soft supersymmetric Lagrangian are then extracted from ISASUGRA. Whether the input parameters are defined at the weak scale or at the GUT scale, one can always choose to redefine additional parameters, such as the standard model parameters or the widths of Higgses or SUSY particles. In processes with t-channel poles it is sometimes necessary to specify a width for some particles such as gauginos. By default these widths have been set to 1GeV. Even though the LSP is assumed to be stable, to avoid any spurious pole it is necessary to introduce also a small width. Numerically, the results do not depend on the exact value chosen for these widths.

After reading the input parameters and calculating the physical parameters needed for the evaluation of the cross-sections using the functions of CompHEP, the calculation of the relic density is performed. The value of $`\mathrm{\Omega }h^2`$ as well as the list of channels that give the most significant contribution to it are given. For all specifications on the functions available and the options that can be set see . In addition we provide subroutines that calculate various constraints on the MSSM parameters: direct limits from colliders, $`\mathrm{\Delta }\rho `$, $`bs\gamma `$ and $`(g2)_\mu `$. All these constraints can be updated or replaced easily.

## 5 Results and Comparisons

The micrOMEGAs code was extensively tested against another public package for calculating the relic density, DarkSUSY. As discussed previously, the two codes differ somewhat in the numerical method used for solving the density equations. micrOMEGAs includes more subprocesses(e.g. all coannihilations with sfermions), loop-corrected Higgs widths, and complete tree-level matrix elements for all processes. DarkSUSY includes on the other hand some loop induced processes such as $`\chi \chi gg,\gamma \gamma `$ which are generally small. Whenever the coannihilation channels with sfermions and the Higgs pole are not important we expect good agreement with DarkSUSY. We have first compared micrOMEGAs with a version of DarkSUSY where we have replaced the matrix elements by CompHEP matrix elements. As expected, we found excellent numerical agreement between the two programs (at the $`2\%`$ level) for all points tested. We have then made comparisons with the original version of DarkSUSY. The results of these comparisons are displayed for a few test points in Table 1. We found in general good agreement between the two codes. However we have observed some discrepancies that could reach up to ($`30\%`$) in particular in the process $`\chi _1^0\chi _1^+t\overline{b}`$(model B in Table 1)111The CompHEP result for this matrix element agrees with the result of GraceSUSY.. As displayed in the line $`\mathrm{\Omega }_{tree}^\chi `$, when removing non-gaugino coannihilation channels and reverting to the tree-level treatment of the Higgs width we recover results similar to DarkSUSY. The impact of these extra channels, model C for sleptons and model D for squarks can be as large as an order of magnitude and depends critically on the mass difference with the lightest neutralino 222Extensive comparisons of micrOMEGAs with an improved version of DarkSusy including slepton coannihilation channels were also performed recently in , complete agreement was found..

In Fig. 1 the variation of the relic density as function of the mass difference between the neutralino $`\stackrel{~}{\chi }_1^0`$ and the Next-to-Lightest Supersymmetric particle(NLSP) is displayed. In Fig. 1a, the SUSY parameters correspond to model C, with either the right-handed slepton masses Mr1=Mr2=Mr3 or the left-handed slepton masses Ml1=Ml2=Ml3 as free parameters (in this case we have also fixed Mr3=289 GeV). For the former choice of parameters, the $`\stackrel{~}{\tau }`$ is the NLSP while for the latter the NLSP is a sneutrino. Model C corresponds to $`\mathrm{\Delta }M_{\stackrel{~}{\tau }\stackrel{~}{\chi }}=10.4`$GeV. In a model where the neutralino is the only source of darkmatter, the lower bound on the relic density then implies a minimum mass difference between the LSP and the NLSP($`>510`$GeV), avoiding the difficulties of detecting a nearly degenerate NLSP at colliders. In Fig.1b, the dependence of the relic density on the mass difference between a coanihilating squark and the LSP is displayed. The SUSY parameters are those of model D with either the right-handed u-squark or d-squark masses kept as free parameters. These two cases lead to a $`\stackrel{~}{t}`$ or a $`\stackrel{~}{b}`$ NLSP respectively. Model D corresponds to $`\mathrm{\Delta }M_{\stackrel{~}{t}\stackrel{~}{\chi }}=26.9`$GeV. Due to the large cross-sections in coannihilation channels involving strongly interacting particles, One witnesses a sharp drop in the relic density as soon as $`\mathrm{\Delta }M_{NLSPLSP}`$ drops below 50GeV.

As stressed above, the effect of the Higgs width is particularly important at large $`\mathrm{tan}\beta `$ with the enhanced contribution of the $`b`$-quark coupling to the heavy scalar Higgs. However, the one-loop QCD corrections to the widths amount to a reduced effective b-quark mass and a much smaller width especially at large values of $`m_H`$. If it was not for the strong Boltzmann suppression factor singling out the contribution at $`\sqrt{s}2m_\chi `$ there will be little difference after integrating over the peak for the one-loop or tree-level result. However the effect observed can be as much as a factor 2. For $`m_\chi M_A/2`$ the narrower resonance (1-loop) suffers less from the Boltzmann reduction factor leading to $`<\sigma ^{1loop}>/<\sigma _{tree}>>1`$ and $`\mathrm{\Omega }_{1loop}<\mathrm{\Omega }_{tree}`$. Further away from the pole however one catches the contribution from the wider resonance without excessive damping from the Boltzmann factor, as expected $`\mathrm{\Omega }_{1loop}>\mathrm{\Omega }_{tree}`$. This is illustrated both for the parameters of models E and F while varying $`M_A`$ in a wide range around the values $`M_A=2m_\chi `$.

Our numerical results were also compared with Ref.. Qualitative agreement is found in the case of a SUGRA model although we use a different RGE code. Precise comparisons necessitates a careful tuning of parameters to make sure we have the same parameters at the weak scale. A random scan over $`m_0m_{1/2}`$ for $`\mu >0`$ shows the typical shape of the allowed region ($`.1<\mathrm{\Omega }h^2<.3`$) in SUGRA models for moderate values of $`\mathrm{tan}\beta `$. For large $`\mathrm{tan}\beta `$, one clearly see the dramatic effect of the heavy Higgs pole in the central band that goes up to very large values of $`m_0`$, $`m_{1/2}`$. Anywhere in this band or in the coannihilation tails a generally heavy supersymmetric particles spectrum can be expected.

## 6 Conclusion

The package micrOMEGAs allows to calculate the relic density of the LSP in the most general MSSM with $`R_P`$ conservation. This is the first program that includes all possible coannihilation channels333Another code that also uses CompHEP for calculating all processes appeared very recently.. The package is self-contained safe for the ISASUGRA package that is required when using the SUGRA option. All possible channels for coannihilations are included and all matrix elements are calculated exactly at tree level with the help of CompHEP. Loop corrections for the masses of Higgs particles (two-loop) and the width of the Higgs (QCD one-loop) are implemented. Good agreement with existing calculations is found when identical set of channels are included. Future versions will include interfaces to other codes that use the renormalization group equations to calculate the weak scale parameters. Including loop corrections to neutralino masses is also planned. Even though these corrections are only a few GeV’s they can alter significantly the calculation of the relic density when there is a near mass degeneracy with the next to lightest supersymmetric particle that contributes to a coannihilation channel . Although the loop processes are in general small, we plan to include $`\chi \chi \gamma \gamma ,\gamma Z,gg`$ in an update of micrOMEGAs.

Acknowledgements

This work was supported in part by the PICS-397, Calcul en physique des particules.


format: text

text:
# A Lattice Boltzmann method for simulations of liquid-vapor thermal flows

## I Introduction

After years of research, the Lattice Boltzmann Methods (LBM) has become an established numerical approach in computational fluid dynamics (CFD). Many models and extensions have been formulated that cover a wide range of complex fluids and flows. Furthermore, LBM has been extended to include turbulence models that have already had a direct and substantial impact on engineering applications.

Among many desirable LBM features such as simplicity, parallelizability, and robustness in dealing with complex boundary conditions, one recognized advantage is its capability of simulating fluid flows with multiple phases. The core mechanism in LBM modeling of multiphase flows is its microscopic level realization of non-ideal gas equations of state. As a result, at sufficiently low temperature and proper pressure, liquid-vapor like first order phase transitions are spontaneously generated. There is no need to explicitly tracking the interfaces between immiscible phases. Furthermore, unlike static statistical physical models, LBM also contains the momentum conservation, so that bubbles and liquid droplets are formed along with fluid hydrodynamic processes. The success and simplicity of LBM for multiphase flows has led to various applications that include simulations of oil-water mixtures through porous media, Rayleigh-Taylor problems, and many more. On the other hand, there is an crucial missing piece. That is, so far all the existing multiphase LBM models are limited to regimes in which the temperature dynamics is either negligible or its effect on flow is unimportant. This limitation, along with the overall unavailability in CFD, has prevented us from dealing with an important class of flows, namely multiphase flows involving strong couplings with thermodynamics. Specific examples of such type of flows range from the common water boiling processes to thermal nuclear reactor applications. Thus from both fundamental and practical point of views, extensions of the existing CFD and LBM methods to simulation of thermal multiphase flows is extremely important.

The LBM is originally evolved from lattice gas models obeying fundamental conservation laws and symmetries. Now it has also been shown to be systematically derivable from the continuum Boltzmann equation. The most commonly known lattice Boltzmann equation (LBE) has the following form (adopting the lattice units convention in which $`\mathrm{\Delta }t=\mathrm{\Delta }x=1`$),

$$f_i(𝐱+\widehat{c}_i,t+1)f_i(𝐱,t)=𝒞_i,$$
(1)
where time $`t`$ takes on only positive integer values, and the particle velocity takes on a finite set of discrete vector values (speeds), $`\{\widehat{c}_i;i=0,\mathrm{},b\}`$. These speeds form links among nodes on a given lattice. The collision term on the right hand side of eqn.(1) now often uses the so called Bhatnagar-Gross-Krook (BGK) approximation,

$$𝒞_i=\frac{f_if_i^{eq}}{\tau },$$
(2)
having a single relaxation time parameter, $`\tau `$. Here, $`f_i^{eq}`$ is the local equilibrium distribution function that has an appropriately prescribed functional dependence on the local hydrodynamic properties. The basic hydrodynamic quantities, such as fluid density $`\rho `$ and velocity $`𝐮`$, are obtained through simple moment summations,

$`\rho (𝐱,t)`$ $`=`$ $`{\displaystyle \underset{i}{}}f_i(𝐱,t),`$ (3)
$`\rho 𝐮(𝐱,t)`$ $`=`$ $`{\displaystyle \underset{i}{}}\widehat{c}_if_i(𝐱,t).`$ (4)
In addition, one can also define a fluid temperature $`T`$ from,

$$\rho \frac{D}{2}T(𝐱,t)=\underset{i}{}\frac{1}{2}(\widehat{c}_i𝐮(𝐱,t))^2f_i(𝐱,t),$$
(5)
where $`D`$ is the dimension of the momentum space of the discrete lattice velocities. It has been theoretically shown that the hydrodynamic behavior produced from LBE obeys the Navier-Stokes fluid dynamics at a long wave-length and low frequency limit. The resulting equation of state is that of an ideal gas fluid, namely the pressure $`p`$ obeys a linear relation with density and temperature,

$$p=\rho T.$$
(6)
The kinematic viscosity of the fluid is related to the relaxation parameter by

$$\nu =(\tau \frac{1}{2})T.$$
(7)
LBM has been extended to simulations of multiphase flows. The key step is to introduce an additional term, $`\mathrm{\Delta }f_i(𝐱,t)`$ on the right hand side of eqn.(1), to represent a body-force. This force term is self-consistently generated by the neighboring distribution functions around each lattice site, and it does not either violate the local mass conservation nor the global momentum conservation. However, the local momentum is altered by an amount,

$$𝐅(𝐱,t)=\underset{i}{}\widehat{c}_i\mathrm{\Delta }f_i(𝐱,t).$$
(8)
The appearance of the body force term can be physically, in a mean-field sense, attributed to a non-local interacting potential $`U`$ among the particles. The existence of such an interacting potential is the essential mechanism in the non-ideal gas type of fluids. Hence with a suitable choice of $`U`$, spontaneous phase separations can be produced, and one can use it conveniently to numerically study multiphase flow phenomena. Through the years, there have been many progresses in LBM models for multiphase flows. On the other hand, as pointed out at the beginning, all the existing attempts are limited to isothermal (or “a-thermal”) situations in which the dynamics of temperature in the fluid is suppressed. That is, $`T`$ is either assumed a constant or, at best, a prescribed function of space (or time).

## II Multiphase flows with incorporation of thermodynamics

In this paper, we present an extension of multiphase LBM to include the full thermodynamic process.

The most natural extension in LBM for thermodynamics has been to introduce a conserved energy degree of freedom. This is relatively straightforward for the ideal gas type of models in which only pointwise collisions are involved and only kinetic energy is considered. When a sufficient number of particle speeds is used, one can theoretically show that LBM leads to the correct full set of thermohydrodynamic equations of an ideal-gas fluid. Unfortunately, besides being considerably more expensive computationally than the isothermal LB models, such an approach cannot be easily generalized to multiphase thermodynamic flows. The most obvious obstacle is the difficulty in tracking the energy evolution while maintaining a total energy conservation: For a non-ideal gas system, the total energy also contains an interacting energy part that is a function of the relative positions among the particles. Without a total energy conservation, a temperature variable cannot be defined fully self-consistently at the microscopic level. In addition, it has been shown that, unlike the isothermal models, an LBE with an energy degree of freedom does not guarantee a global $`H`$-theorem. As a consequence, the system can exhibit significantly less stability. Other undesirable features in this direct approach include 1) difficulty in changing Prandtl number value from unity, unless a substantial generalization to the BGK collision term is made; and 2) a rather limited temperature range (with the maximal allowable value only about twice the minimal value), unless significantly more speeds are added. Because of these reasons, the progress in LBM for thermal multiphase flows has been rather slow.

Here we present a new LBM approach that can essentially avoid all of the above mentioned drawbacks. The fundamental idea can be briefly summarized: First of all, the fluid dynamics part (i.e., the density and momentum evolution) is represented by a modified isothermal LBE, while the energy evolution part is determined by an additional scalar energy transport equation. The latter can be solved either via a finite difference scheme or an auxiliary LBE. Secondly, the coupling of the two parts is through a properly defined body force term in the LBE (and the compression and dissipation terms in the energy equation). As we shall realize below, although conceptually rather simple, the new model produces the correct full thermohydrodynamic equations together with a non-ideal gas equation of state.

We choose a common isothermal LBE (e.g., D3Q19,) as a starting basis. As discussed earlier, an isothermal LBE model for fluid density and velocity evolution is considerably simpler compared to its energy conserving counterpart. This is certainly desirable for doing efficient fluid flow simulations. Furthermore, the equilibrium distribution in an isothermal LB model is only a function of fluid density and velocity. The lack of temperature dynamics in the equilibrium distribution is the key for achieving a higher stability in LBE. Having these in mind, it is much desirable to introduce a macroscopic mechanism to recover thermodynamics. Specifically, instead of letting temperature to influence the equilibrium distributions in an LB system, the thermodynamic effect is obtained via a temperature-dependent body-force. Because of the “external” nature of the coupling, the LB system and its equilibrium property remain to be microscopically isothermal. Nonetheless, as explained below, this alternative way of coupling achieves the desired thermodynamics at the macroscopic level.

Ignoring the higher order contributions, the body-force term can be simply expressed as,

$$\mathrm{\Delta }f_i(𝐱,t)=\frac{w_i}{T_0}\widehat{c}_i𝐅,(𝐱,t)$$
(9)
where the constant weights $`w_i`$ and $`T_0`$ are directly determined by the LBE model (e.g., D3Q19;, in which $`T_0=1/3`$). One can easily verify that this gives rise to eqn.(8). The global momentum conservation is preserved as long as $`𝐅(𝐱,t)`$ is expressed as a spatial gradient of a scalar function,

$$𝐅(𝐱,t)=U(𝐱,t).$$
(10)
It is straightforward to implement this condition in a discrete space by proper finite-difference procedures. Based on the consideration for higher order isometry in surface tension, we choose (for D3Q19) the following specific form,

$$U(𝐱,t)\underset{i}{}\frac{D}{bc_i^2}\widehat{c}_iU(𝐱+\widehat{c}_i,t)$$
(11)
With the additional body-force term, one can easily recognize that the overall effective pressure in the resulting fluid momentum equation has become,

$$p=\rho T_0+U$$
(12)
where the first term is a result of the isothermal LBM. From (12), one can obtain any form of equation of state, $`p=p(\rho ,T)`$, simply by making a corresponding choice for $`U`$,

$$U(𝐱,t)=p(\rho (𝐱,t),T(𝐱,t))\rho (𝐱,t)T_0.$$
(13)
The above quantity is determined once the local values of $`\rho (𝐱,t)`$ and $`T(𝐱,t)`$ are provided. Obviously the resulting fluid is no longer isothermal if the temperature $`T(𝐱,t)`$ varies. In other words, because of the macroscopic way of coupling, the resulting fluid dynamics is no longer isothermal. In addition, with the high flexibility of choosing the equation of state, this approach can be applied to simulation of non-ideal gas fluids and multiphase flows. Indeed, we confirmed this basic feature through a set of spinodal decomposition tests based on a Van de Waals gas model (Carnahan-Starling equation of state). Similar to the other multiphase LBM, a spontaneous phase separation process is well observed at sufficiently low temperature values.

The evolution of the temperature $`T(𝐱,t)`$ in this new approach is obtained from solving a supplemental scalar energy transport equation,

$$\rho (_t+𝐮)e=p𝐮+\kappa T+\mathrm{\Psi }$$
(14)
where $`e=c_vT`$ is the internal energy, and $`c_v`$ is the specific heat at constant volume of the fluid. The overall pressure $`p`$ is defined by the equation of state (12), and $`\kappa `$ is the heat conductivity that can be specified flexibly. The term $`\mathrm{\Psi }`$ represents the viscous dissipation of flow and the contribution of surface tension. The energy evolution equation (14) is a standard macroscopic description for thermal fluids. The computation of an isothermal LB model along with a scalar energy equation is considerably less expensive than any microscopic attempts: for it neither requires many particle speeds nor complicated tracking of the energy evolution. Moreover, the difficulties in stability and Prandtl number associated with the original thermal LBM are not issues in this approach. Solving a scalar transport equation is rather straightforward. There are many finite-difference schemes for accurately and efficiently solving the scalar transport equation. In our particular simulations, we have used an extended Lax-Wendroff scheme. The combination of eqns. (1)-(4), (9)-(11), and (13)-(14) forms the new LBM approach for modeling multiphase thermodynamic fluid flows. The thermal boundary condition can be realized via standard numerical procedures so that,

$$\kappa \widehat{n}T|_w=q,$$
(15)
with a prescribed heat flux $`q`$, that can either be fixed or a function of local properties in order to achieve a fixed wall temperature. The unit vector $`\widehat{n}`$ denotes the surface normal direction.

There is one more important feature in this model worth pointing out. That is, the new approach avoids the fundamental limitation on the temperature range that has constrained the other thermal LB models. Notice the temperature only appears in the body-force term in a gradient function form. Hence, unlike that for equilibrium distribution functions, there is no absolute upper or lower bound on the temperature values except that it should not change too rapidly across a given resolution scale. Moreover, there is obviously no absolute bound on temperature in the energy equation.

Based on the above discussions, one can realize that the new LB model generates a fully macroscopically consistent description for thermodynamic flows involving generalized equations of state. Therefore, this approach offers a convenient and efficient numerical tool for studying thermal multiphase flow problems.

## III Simulation of the liquid-vapor boiling process

In this section we present computational results of a typical multiphase thermodynamic flow simulation with the new LB approach. In particular, a liquid-vapor boiling process involving Rayleigh-Benard like convections, phase changes, together with a complex temperature dynamics is simulated successfully, albeit qualitative. Although representative of a wide range of important applications, boiling flow problems have received little success from CFD in general. Consequently, it is very important that the new approach can demonstrate such a fundamental capability.

The Rayleigh-Benard convection process has been widely used as a benchmark for many fluid computations. It is the simplest representation of a boiling phenomenon in which a complex buoyancy-driven convection process occurs at various values of the Rayleigh number. On the other hand, most of the boiling processes occur in nature also involve evolutions of multiple thermodynamic phases. That is, besides thermal convection, the fluid under goes a phase transition process in which liquid droplets and vapor bubbles are generated. The most obvious practical examples include the common water boiling in a pot.

We choose a standard Rayleigh-Benard setup, in which both the upper and the lower solid plates obey no-slip boundary conditions, while the horizontal boundary condition is periodic. To achieve more stable and second-order accurate numerical results, we have also applied the new scheme to the modified LB discretization formulation of He et al. As discussed above, the Carnahan-Starling equation of state is used here for convenience. The mean density value $`\rho =1.36\rho _c`$. The temperature on the upper wall is fixed at $`T_u=0.795T_c`$ while on the lower wall $`T_l=0.954T_c`$. Here $`T_c`$ ($`=0.55`$ in lattice units for the choice of the model) is the critical temperature for formation of two phases. The initial temperature is set to be linearly distributed between the two plates and is consistent with their temperature boundary conditions. The simulation volume is $`L\times H=256\times 128`$ grid points. The gravity value $`g=5\times 10^6`$ (lattice units) is used. In order to avoid unnecessary complications, a weak surface tension effect is also applied in the LBM flow simulation, so that the surface tension contribution to the energy evolution can thus be neglected. Specifically, we choose the surface tension coefficient $`\sigma `$ to be $`0.01`$ for these simulations. The kinematic viscosity and the Prandtl number are set at $`\nu =0.02`$ and $`Pr=10`$, respectively. For simplicity, the heat capacity $`c_v`$ in our multiphase flow is chosen as a constant ($`=1`$). All the other fluid parameters are the same as in He et al. Based on the choice of these parameter values, the resulting Rayleigh number is $`R_a3.0\times 10^5`$, which is much higher than the first threshold ($`Ra_c=1708`$) for an onset of convection in the conventional single-phase Rayleigh-Benard system.

The simulation starts from a uniform density distribution with one percent random fluctuations. To enhance bubble formation, small temperature fluctuations are added to the equation of state in the first grid point layer near the bottom solid plate. Due to the higher temperature, a lower density fluid is been produced that subsequently leads to small vapor bubbles at the bottom plate, and then these merge with each other to form larger size ones. The phase formation process in the simulation is combined with the convection process in which the hotter and lighter vapor phase rises while the colder and heavier liquid phase descends due to gravity. With the above choice of parameter values, such a full thermohydrodynamic cycle is seen to be able to sustain itself indefinitely. Figure 1 shows the density distributions at some representative nondimensionalized (by $`\sqrt{H/g}`$) times. One can observe the formation of two streams of bubbles near the bottom plate. One can also observe that two pairs of counter rotating convection rolls are locked between the bubble streams. In addition, the convection rolls are seen to pinch off the small bubbles from the bottom plate. As bubbles rise, their sizes are seen to increase slightly. Since small bubbles move faster than large bubbles, collision and coalescence often occur among them. Figure 2 shows streamlines of the velocity field at time 39.5, from which one can clearly see the two pairs of counter rotating convection rolls. Figure 3 depicts the corresponding temperature deviation from the linear distribution at the same time. Interestingly, the temperature exhibits a non-trivial behavior: Its value is seen to be relatively lower in the vapor phase domains near their interfaces. As expected physically, this phenomenon is due to the $`p𝐮`$ term in the energy equation (14) associated with the volume expansion from water to vapor phases. All the above simulated phenomena are qualitatively correct for a realistic two phase thermodynamic flow.

We also ran another simulation with the exact same setup as the above, except that the surface tension is made ten times stronger. As seen in figure 4, only one pair of counter rotating convection rolls is produced at its asymptotic state. Different from the standard single phase Bernard convection, we see that other intrinsic properties such as the surface tension can also alter the thermal convection characteristics in a multiphase flow.

## IV Discussions

In this paper, we present a novel approach that combines a multiphase lattice Boltzmann method with a scalar temperature equation. The coupling is realized macroscopically via a self-consistent body-force. The basic formulation is applicable to both 2D and 3D flow situations. It is directly verifiable theoretically that new LB model obeys the correct full thermohydrodynamic equations with a non-ideal gas equation of state. The incorporation of thermohydrodynamics in the new approach allows for simulations of complex multiphase flows coupled with temperature dynamics. Other important features of the new scheme include its simplicity, efficiency and robustness for simulating thermal multiphase flow processes. The latter has shown to be extremely difficult with the other LBM based schemes or the more conventional numerical methods. Furthermore, like other LBMs, the new approach can handle complex physical boundary conditions . All these are highly desirable for a practical computational model.

We have demonstrated the capability of our new approach for doing boiling flow simulations. This type of flows has so far not been successfully handled via other methods. Thus the new approach has opened a promising opportunity for numerically studying thermal multiphase flow problems. On the other hand, further improvements in the direction of the new approach are necessary in order for it to become a more useful and quantitative computational tool. Without going into detail, the major remaining issues include: 1) incorporating more realistic equations of state instead of the van der Waals type of gas model; 2) a more physical treatment of the heat capacity and latent heat; 3) further understanding and modeling of surface tension and near interface physics; 4) further enhancement of numerical stability to achieve significantly higher density ratio and lower viscosity; and 5) generalization of boundary conditions. Among all the above tasks, we consider (4) to be the most challenging one for LBM in general.

## V Acknowledgments

The authors wish to thank Drs. Xiaoyi He, YueHong Qian, Ilya Staroselsky, Adrian Tentner, Steven Orszag, Sauro Succi, Rick Shock and Nick Martys for their useful discussions. This work is supported in part by DOE SBIR grant No. 65016B01-I.


text:
# Corner transfer matrix renormalisation group method for two-dimensional self-avoiding walks and other 𝑂(𝑛) models

## Abstract

We present an extension of the corner transfer matrix renormalisation group (CTMRG) method to $`O(n)`$ invariant models, with particular interest in the self-avoiding walk class of models ($`O(n=0)`$). The method is illustrated using an interacting self-avoiding walk model. Based on the efficiency and versatility when compared to other available numerical methods, we present CTMRG as the method of choice for two-dimensional self-avoiding walk problems.

The self-avoiding walk class of models on the two-dimensional square lattice, with a variety of possible interactions, has mobilised the scientific community for about half a centuryclois ; vander . The number of exact results for such models is limited, and numerical studies are hard. A clear illustration of the numerical difficulty is the disagreement which existed over the numerical determination of critical temperature and exponents for the standard $`\theta `$-point model, see for example referencesthetanum .

To date the numerical methods available for the study of interacting self-avoiding walks in two dimensions are series expansions of walks of lengths of a few tens of stepsenumer , transfer matrices for lattice widths up to about 12derr  and increasingly complicated Monte-Carlo simulation methodsmcslow , limited in practice to only a portion of the phase diagram.

Motivated by these numerical difficulties, we decided to extend the Corner Transfer Matrix Renormalisation Group (CTMRG) methodnishino . The CTMRG method is based on White’s Density Matrix Renormalisation Group method (DMRG)white  and Baxter’s corner matrix formalismbaxter . To date the CTMRG method has only been applied to discrete spin models, where it is shown to be computationally efficientnishino .

Our extension to interacting self-avoiding walk models exploits the connection betweeen these models and the $`O(n)`$ invariant spin modelsnienon , which contains as special cases the Ising model ($`n=1`$), the XY model ($`n=2`$) and the Heisenberg model ($`n=3`$). The method therefore has applications well beyond the self-avoiding walk type models ($`n=0`$).

The $`O(n)`$ spin model is defined through the partition functionon

$$𝒵_{O(n)}=\underset{\{\stackrel{}{s}_i\}}{}\mathrm{exp}\left(\frac{1}{2}\beta J\underset{i,j}{}\stackrel{}{s}_i\stackrel{}{s}_j\right),$$
(1)
where $`i,j`$ refers to a sum over nearest neighbour spins. The spin $`\stackrel{}{s}_i`$ has $`n`$ components, and is normed such that $`s_i^2=1`$. Another formulation of $`O(n)`$ invariant models, with the same critical behaviour, is:

$$𝒵_{O(n)}=\underset{\{\stackrel{}{s}_i\}}{}\underset{i,j}{}\left(1+K\stackrel{}{s}_i\stackrel{}{s}_j\right).$$
(2)
where the spins are now placed on the lattice bondsnienon . A diagrammatic expansion of Equation 2 follows if we identify the $`1`$ as the weight of an empty bond between the sites $`i`$ and $`j`$ and the $`K`$ as the weight of an occupied bond. This expansion may be expressed in terms of graphs $`𝒢`$ of non-intersecting loops (collisions at sites are however allowed)nienon . The partition function may then be written:

$$𝒵_{O(n)}=\underset{𝒢}{}n^{l(𝒢)}K^{b(𝒢)},$$
(3)
where $`l`$ is the number of loops and $`b`$ is the number of occupied bonds.

The parameter $`n`$ is now a fugacity controlling the number of loops, and need no longer be taken as an integer. This fugacity corresponds to a long ranged interaction, since the loops may be of any size. This non-locality is undesirable for our purposes. We would like to express $`n`$ as a product over local weights. This may be achieved as follows. Each loop may be followed clockwise or anticlockwise. The loops in Equation 3 are not oriented, but may be oriented by associating $`2^{l(𝒢)}`$ oriented graphs with each non-oriented graph. By associating a loop fugacity $`n_+`$ ($`n_{}`$) with the clockwise (anticlockwise) oriented loops, the partition function may be rewrittennien2 :

$$𝒵_{O(n)}=\underset{𝒢}{}(n_++n_{})^lK^b=\underset{𝒢^{}}{}n_+^{l_+}n_{}^l_{}K^b$$
(4)
where $`𝒢^{}`$ is the set of oriented loop graphs, and $`l_+`$ ($`l_{}`$) is the number of clockwise (anticlockwise) oriented loops. Setting $`n_+=\mathrm{exp}(i\theta )`$ and $`n_{}=\mathrm{exp}(i\theta )`$ gives $`n=2\mathrm{cos}(\theta )`$. The oriented loop factor is now broken up into local weights by associating a corner weight $`w_i=\mathrm{exp}(i\theta /4)`$ with every clockwise corner and $`w_i=\mathrm{exp}(i\theta /4)`$ with every anticlockwise corner. On the square lattice, there must be four more corners with one orientation, compared with the other orientation, in order to close a loop. The product of the local weights will then give the correct weight for the oriented loops.

The partition function may now be rewritten in terms of a vertex modelnienon :

$$𝒵_{O(n)}=\underset{𝒢^{}}{}\underset{i}{}v_i$$
(5)
where $`v_i`$ is the weight of the vertex at site $`i`$. The derivation given here is only for the simplest case, but we may freely change the weights of the vertex configurations in order to generate different interactions in the original model (see figure 1).

The limit $`n0`$ corresponds to the self-avoiding walk modelpgdg . In the generalised form presented here this corresponds to an interacting self-avoiding walk model due to Blöte and Nienhuisnienon ; blbatnien . A step fugacity $`K`$ and an attractive short ranged interaction $`\epsilon <0`$ are introduced. In the standard $`\theta `$ point model the interactions are between non-consecutively visited nearest-neighbour sites and a given site may only be visited onceclois ; vander . In our current models, this last constraint is relaxed; the walk may collide at a site, but not cross, and remains self-avoiding for the bonds. The interaction is now assigned to doubly visited sites. An additional weight $`p`$ is added for sites which are visited by a straight section of walk (i.e. do not sit on a corner). The partition function may then be written:

$$Z=\underset{\mathrm{walks}}{}(Kp)^L\tau ^{N_I}p^{N_c},$$
(6)
where $`\tau =\mathrm{exp}(\beta \epsilon )`$, $`N_I`$ is the number of site collisions and $`N_c`$ are the number of corners in the walk. This model gives rise to the vertex weights shown in figure 1. The standard self-avoiding walk model is found setting $`\tau =0`$ and $`p=1`$ and when $`p=0`$ the model has the same critical behaviour as the standard $`\theta `$ point modelblbatnien . During the remainder of this letter we shall illustrate the CTMRG method for $`O(n)`$ models using this Blöte-Nienhuis walk model.

Following Baxterbaxter , the partition function of a two dimensional lattice model may be written in terms of the product of four matrices representing the four quarters of the lattice. The inputs and outputs of the matrices are the configurations at the seams of the four quarters. These matrices are known as corner transfer matrices. In general the four matrices are different, but may often be related by lattice symmetries. For our model the four matrices are the same up to a complex conjugation operationpreprint .

It is usually not possible to calculate explicitly these matrices for systems with a large number of sites. This is where the CTMRG method comes in; the matrices for larger lattices are calculated from smaller lattices iterativelynishino . This is done as follows. An initial system, consisting of a small number of sites, is mapped exactly onto a prototype system made up of four $`m`$-state spins. At each iteration the system is enlarged by adding sites, this enlarged system is then projected back onto the prototype system in some optimal way, so as to minimise the loss of information. The value of $`m`$ determines the amount of information which may be carried forward at each iteration, the larger the value of $`m`$ the better the approximation. For details see nishino ; preprint .

As with the DMRG method, the innermost sites of the lattice are treated exactly. This means that the CTMRG method is most appropriate for the calculations of one-point functions such as the site free energy, the site density of monomers, specific heats etc. In figure 2 we show the density, $`\rho `$, as a function of the step-fugacity $`K`$ for fixed values of $`\tau `$. We chose to fix $`p=0`$ since for this value the Nienhuis-Blöte model has the same critical behaviour as the standard $`\theta `$ point modelblbatnien . When $`\tau <2`$ the model is expected to be in the self-avoiding walk universality class and present a critical transition where the density changes continuously from $`\rho =0`$ for $`K<K_c`$ to $`\rho >0`$ for $`K>K_c`$. When $`\tau =2`$ the transition is expected to be tricritical, in the same universality class as the $`\theta `$ point, with $`K_c=1/2`$ exactly. For $`\tau >2`$ the transition is expected to be first order. These three types of behaviour may be clearly seen in figure 2.

Whilst the raw finite-size data presented in figure 2 already gives a fairly precise idea of the behaviour of the system (order of the transition and first estimate of $`K_c`$), practical calculations of critical temperature and critical exponent estimates rely on finite-size analysis. In transfer matrix calculations one usually uses a phenomenological renormalisation group based on the correlation lengths due to Nightingalempn76 . In the CTMRG method an effective transfer matrix may be determined directly. However, as already mentioned, the highest precision is obtained for the one point correlation functions. Since the number of sizes available is an order of magnitude larger than that available for transfer matricies it is advantageous to exploit the finite size scaling laws for the one point functions. In particular it is expected that the singular part of the density scales as

$$\rho _s(K,L)=L^{1/\nu 2}\stackrel{~}{\rho }(|KK_c|L^{1/\nu }).$$
(7)
This scaling behaviour implies that the function

$$\phi _{L,L^{}}(K)=\frac{\mathrm{log}\left(\rho _s(K,L)/\rho _s(K,L^{})\right)}{\mathrm{log}\left(L/L^{}\right)}$$
(8)
takes the value $`\phi _{L,L^{}}=1/\nu 2`$ when $`K=K_c`$, independantly of $`L`$ and $`L^{}`$. Naturally there are additional finite size corrections which should be taken into account, but the conclusion is that if the function $`\phi _{L,L^{}}(K)`$ is plotted for various values of $`L`$ and $`L^{}`$ then it will converge to a fixed point given by $`\phi (K_c)=1/\nu 2`$. In what follows we have set $`L^{}=L2`$ and looked for solutions of the equation

$$\phi _{L,L2}(K_c^L)=\phi _{L2,L4}(K_c^L),$$
(9)
where $`K_C^L`$ is the estimated critical temperature from sizes $`L`$, $`L2`$ and $`L4`$. If such a solution does exist then $`K_c=lim_L\mathrm{}K_c^L`$ and $`\nu =lim_L\mathrm{}1/(\phi _{L,L2}(K_c^L)2)`$.

Figure 3 shows $`\phi `$ plotted for several values of $`L`$ for the pure self-avoiding walk model ($`p=1`$, $`\tau =0`$). The different curves cross at a point defining $`K_c`$ and $`\nu `$. The corresponding estimates of $`K_c`$ and $`\nu `$ are shown in Figure 4 plotted as a function of $`1/L`$. The extrapolations of $`K_c^L`$ and $`\nu ^L`$ to $`L\mathrm{}`$ are given in Table 1 along with prelimenary estimates for $`p=0`$ and different values of $`\tau `$. The accuracy with which the critical points are determined is in general an order of magnitude better than found with transfer matricies. Full details and definitive estimates will be given elsewhere preprint .

In (numerically) exact methods, such as the transfer matrix method or the exact enumeration method, the main source of uncertainty in the results is due to the extrapolation of a small number of points. In the CTMRG method the extrapolation problems are largely removed since lattice sizes an order of magnitude larger have been reached. The uncertainty now lies in the precision related to the calculation of each density point.

In this article results for the Blöte-Nienhuis interacting self-avoiding walk are presented since good quality results already exist, providing a good test of the efficiency of the CTMRG method. In particular we have focused on the self-avoiding walk model ($`p=1`$, $`\tau =0`$) for which the asymptotic limit is well described by transfer matrix calculations and yet we still find a substantial increase in precision (see Table 1). It is clear that in circumstances where larger system sizes are required to extract the scaling behaviour, CTMRG should far exceed the numerically exact methods in performance.

The quality of results presented in this article is virtually unattainable with such ease by any other numerical method we know of (there is a very small number of exceptional cases where better accuracy was obtainedexceptional ) and so we present CTMRG as the method of choice for two-dimensional self-avoiding walk models. CTMRG may easily be extended to more complicated interacting self-avoiding walk models, such as the hydrogen-bonding self-avoiding walkfs .


text:
# Partially Ordered generalized patterns and 𝑘-ary words

## Abstract

Recently, Kitaev \[Ki2\] introduced partially ordered generalized patterns (POGPs) in the symmetric group, which further generalize the generalized permutation patterns introduced by Babson and Steingrímsson \[BS\]. A POGP $`p`$ is a GP some of whose letters are incomparable. In this paper, we study the generating functions (g.f.) for the number of $`k`$-ary words avoiding some POGPs. We give analogues, extend and generalize several known results, as well as get some new results. In particular, we give the g.f. for the entire distribution of the maximum number of non-overlapping occurrences of a pattern $`p`$ with no hyphens (that allowed to have repetition of letters), provided we know the g.f. for the number of $`k`$-ary words that avoid $`p`$.

## 1. Introduction

Let $`[k]^n`$ denote the set of all the words of length $`n`$ over the (totally ordered) alphabet $`[k]=\{1,2,\mathrm{},k\}`$. We call these words by $`n`$-long $`k`$-ary words. A *generalized pattern* $`\tau `$ is a word in $`[\mathrm{}]^m`$ (possibly with hyphens between some letters) that contains each letter from $`[\mathrm{}]`$ (possibly with repetitions). We say that the word $`\sigma [k]^n`$ *contains* a generalized pattern $`\tau `$, if $`\sigma `$ contains a subsequence isomorphic to $`\tau `$ in which the entries corresponding to consecutive entries of $`\tau `$, which are not separated by a hyphen, must be adjacent. Otherwise, we say that $`\sigma `$ *avoids* $`\tau `$ and write $`\sigma [k]^n(\tau )`$. Thus, $`[k]^n(\tau )`$ denotes the set of all the words in $`[k]^n`$ that avoid $`\tau `$. Moreover, if $`P`$ is a set of generalized patterns then $`|[k]^n(P)|`$ denotes the set all the words in $`[k]^n`$ that avoid each pattern from $`P`$ simultaneously.

###### Example 1.1.

A word $`\pi =a_1a_2\mathrm{}a_n`$ avoids the pattern $`13\text{-}2`$ if $`\pi `$ has no subsequence $`a_ia_{i+1}a_j`$ with $`j>i+1`$ and $`a_i<a_j<a_{i+1}`$. Also, $`\pi `$ avoids the pattern $`121`$ if it has no subword $`a_ia_{i+1}a_{i+2}`$ such that $`a_i=a_{i+2}<a_{i+1}`$.

Classical patterns are generalized patterns with all possible hyphens (say, $`2\text{-}1\text{-}3\text{-}4`$), in other words, those that place no adjacency requirements on $`\sigma `$. The first case of classical patterns studied was that of permutations avoiding a pattern of length 3 in $`\mathrm{SS}_3`$. Knuth \[Knuth\] found that, for any $`\tau \mathrm{SS}_3`$, $`|\mathrm{SS}_n(\tau )|=C_n`$, the $`n`$th Catalan number. Later, Simion and Schmidt \[SS\] determined the number $`|\mathrm{SS}_n(P)|`$ of permutations in $`\mathrm{SS}_n`$ simultaneously avoiding any given set of patterns $`P\mathrm{SS}_3`$. Burstein \[Bu\] extended this to $`|[k]^n(P)|`$ with $`P\mathrm{SS}_3`$. Burstein and Mansour \[BM1\] considered forbidden patterns with repeated letters. Also, Burstein and Mansour \[BM2, BM3\] considered forbidden generalized patterns with repeated letters.

Generalized permutation patterns were introduced by Babson and Steingrímsson \[BS\] with the purpose of the study of Mahonian statistics. Claesson \[C\] and Claesson and Mansour \[CM\] considered the number of permutations avoiding one or two generalized patterns with one hyphen. Kitaev \[Ki1\] examined the number of $`|\mathrm{SS}_n(P)|`$ of permutations in $`\mathrm{SS}_n`$ simultaneously avoiding any set of generalized patterns with no hyphens. Besides, Kitaev \[Ki2\] introduced a further generalization of the generalized permutation patterns namely partially ordered generalized patterns.

In this paper we introduce a further generalization of the generalized patterns namely partially ordered generalized patterns in words (POGPs), which is an analogue of POGPs in permutations \[Ki2\]. A POGP is a generalized pattern some of whose letters are incomparable. For example, if we write $`\tau =1\text{-}1^{}2^{}`$, then we mean that in occurrence of $`\tau `$ in a word $`\sigma [k]^n`$ the letter corresponding to the $`1`$ in $`\tau `$ can be either larger, smaller, or equal to the letters corresponding to $`1^{}2^{}`$. Thus, the word $`113425[5]^6`$ contains seven occurrence of $`\tau `$, namely $`113`$, $`134`$ twice, $`125`$ twice, $`325`$, and $`425`$.

Following \[Ki2\], we consider two particular classes of POGPs – shuffle patterns and multi-patterns, which allows us to give an analogue for all the main results of \[Ki2\] for $`k`$-ary words. A multi-pattern is of the form $`\tau =\tau ^0\text{-}\tau ^1\text{-}\mathrm{}\text{-}\tau ^s`$ and a shuffle pattern of the form $`\tau =\tau ^0\text{-}a_1\text{-}\tau ^1\text{-}a_2\text{-}\mathrm{}\text{-}\tau ^{s1}\text{-}a_s\text{-}\tau ^s`$, where for any $`i`$ and $`j`$, the letter $`a_i`$ is greater than any letter of $`\tau ^j`$ and for any $`ij`$ each letter of $`\tau ^i`$ is incomparable with any letter of $`\tau ^j`$. These patterns are investigated in Sections 3 and 4.

Let $`\tau =\tau ^0\text{-}\tau ^1\text{-}\mathrm{}\text{-}\tau ^s`$ be an arbitrary multi-pattern and let $`A_{\tau ^i}(x;k)`$ be the generating function (g.f.) for the number of words in $`k`$-letter alphabet that avoid $`\tau ^i`$ for each $`i`$. In Theorem 4.7 we find the g.f., in terms of the $`A_{\tau ^i}(x;k)`$, for the number of $`k`$-ary words that avoid $`\tau `$. In particular, this allows us to find the g.f. for the entire distribution of the maximum number of non-overlapping occurrences of a pattern $`\tau `$ with no hyphens, if we only know the g.f. for the number of $`k`$-ary words that avoid $`\tau `$. Thus, in order to apply our results in what follows we need to know how many $`k`$-ary words avoid a given ordinary generalized pattern with no hyphens. This question was examined, for instance, in \[BM1, Sections 2 and 3\]\[BM2, Section 3\] and \[BM3, Section 3.3\].

## 2. Definitions and Preliminaries

A partially ordered generalized pattern (POGP) is a generalized pattern where some of the letters can be incomparable.

###### Example 2.1.

The simplest non-trivial example of a POGP that differs from the ordinary generalized patterns is $`\tau =1^{}\text{-}2\text{-}1^{\prime \prime }`$, where the second letters is the greatest one and the first and the last letters are incomparable to each other. The word $`\sigma =31421`$ has five occurrences of $`\tau `$, namely $`342`$, $`341`$, $`142`$, $`141`$, and $`121`$.

Let $`A_\tau (x;k)=_{n0}a_\tau (n;k)x^n`$ denote the generating function (g.f.) for the numbers $`a_\tau (n;k)`$ of words in $`[k]^n`$ avoiding the pattern $`\tau `$. For $`\tau =1^{}\text{-}2\text{-}1^{\prime \prime }`$, we have

(2.1)
$$A_{1^{}\text{-}2\text{-}1^{\prime \prime }}(x;k)=\frac{1}{(1x)^{2k1}}\underset{j=1}{\overset{k1}{}}\frac{x}{(1x)^{2j}}.$$

Indeed, if $`\sigma [k]^n`$ avoids $`\tau `$, and $`\sigma `$ contains $`s>0`$ copies of the letter $`k`$, then the letters $`k`$ appear as leftmost or rightmost letters of $`\sigma `$. If $`\sigma `$ contains no $`k`$ then $`\sigma [k1]^n`$. So, for all $`n0`$, we have

$$a_\tau (n;k)=a_\tau (n;k1)+2a_\tau (n1;k1)+3a_\tau (n2;k1)+\mathrm{}+(n+1)a_\tau (0;k1),$$

since there are $`(i+1)a_\tau (ni;k1)`$ possibilities to place $`i`$ letters $`k`$ into $`\sigma `$, for $`0in`$. Hence, for all $`n2`$,

$$a_\tau (n;k)2a_\tau (n1;k)+a_\tau (n2;k)=a_\tau (n;k1),$$

together with $`a_\tau (0,k)=1`$ and $`a_\tau (1,k)=k`$. Multiplying both sides of the recurrence above with $`x^n`$ and summing over all $`n2`$, we get Equation 2.1.

###### Definition 2.2.

If the number of words in $`[k]^n`$, for each $`n`$, that avoid a POGP $`\tau `$ is equal to the number of words that avoid a POGP $`\varphi `$, then $`\tau `$ and $`\varphi `$ are said to be equivalent and we write $`\tau \varphi `$.

The reverse $`R(\sigma )`$ of a word $`\sigma =\sigma _1\sigma _2\mathrm{}\sigma _n`$ is the word $`\sigma _n\mathrm{}\sigma _2\sigma _1`$. The complement $`C(\sigma )`$ is the word $`\theta =\theta _1\theta _2\mathrm{}\theta _n`$ where $`\theta _i=k+1\sigma _i`$ for all $`i=1,2,\mathrm{},n`$. For example, if $`\sigma =123331[3]^6`$, then $`R(\sigma )=133321`$, $`C(\sigma )=321113`$, and $`R(C(\sigma ))=311123`$. We call these bijections of $`[k]^n`$ to itself trivial. For example, the number of words that avoid the pattern $`12\text{-}2`$ is the same as the number of words that avoid the patterns $`2\text{-}21`$, $`1\text{-}12`$, and $`21\text{-}1`$, respectively.

Following \[Ki2\], it is convenient to introduce the following definition.

###### Definition 2.3.

Let $`\tau `$ be a generalized pattern without hyphens. A word $`\sigma `$ quasi-avoids $`\tau `$ if $`\sigma `$ has exactly one occurrence of $`\tau `$ and this occurrence consists of the $`|\tau |`$ rightmost letters of $`\sigma `$, where $`|\tau |`$ denotes the number of letters in $`\tau `$.

For example, the word $`5112234`$ quasi-avoids the pattern $`1123`$, whereas the words $`5223411`$ and $`1123345`$ do not.

###### Proposition 2.4.

Let $`\tau `$ be a non-empty generalized pattern with no hyphens. Let $`A_\tau ^{}(x;k)`$ denote the g.f. for the number of words in $`[k]^n`$ that quasi-avoid $`\tau `$. Then

(2.2)
$$A_\tau ^{}(x;k)=(kx1)A_\tau (x;k)+1.$$

###### Proof.

Using the similar arguments as those in the proof of \[Ki2, Proposition 4\], we get that, for $`n1`$,

$$a_\tau ^{}(n;k)=ka_\tau (n1;k)a_\tau (n;k),$$

where $`a_\tau ^{}(n;k)`$ denotes the number of words in $`[k]^n`$ that quasi-avoid $`\tau `$. Multiplying both sides of the last equality by $`x^n`$ and summing over all natural numbers $`n`$, we get the desired result. ∎

###### Definition 2.5.

Suppose $`\{\tau ^0,\tau ^1,\mathrm{},\tau ^s\}`$ is a set of generalized patterns with no hyphens and

$$\tau =\tau ^0\text{-}\tau ^1\text{-}\mathrm{}\text{-}\tau ^s,$$

where each letter of $`\tau ^i`$ is incomparable with any letter of $`\tau ^j`$ whenever $`ij`$. We call such POGPs multi-patterns.

###### Definition 2.6.

Suppose $`\{\tau ^0,\tau ^1,\mathrm{},\tau ^s\}`$ is a set of generalized patterns with no hyphens and $`a_1a_2\mathrm{}a_s`$ is a word of $`s`$ letters. We define a shuffle pattern to be a pattern of the form

$$\tau =\tau ^0\text{-}a_1\text{-}\tau ^1\text{-}a_2\text{-}\mathrm{}\text{-}\tau ^{s1}\text{-}a_s\text{-}\tau ^s,$$

where each letter of $`\tau ^i`$ is incomparable with any letter of $`\tau ^j`$ whenever $`ij`$, and the letter $`a_i`$ is greater than any letter of $`\tau ^j`$ for any $`i`$ and $`j`$.

For example, $`1^{}\text{-}2\text{-}1^{\prime \prime }`$ is a shuffle pattern, and $`1^{}\text{-}1^{\prime \prime }`$ is a multi-patterns. From definitions, we obtain that we can get a multi-pattern from a shuffle pattern by removing all the letters $`a_i`$.

There is a connection between multi-avoidance of the generalized patterns and the POGPs. In particular, to avoid $`1^{}\text{-}2\text{-}1^{\prime \prime }`$ is the same as to avoid simultaneously the patterns $`1\text{-}2\text{-}1`$, $`1\text{-}3\text{-}2`$, and $`2\text{-}3\text{-}1`$. A straightforward argument leads to the following proposition.

###### Proposition 2.7.

For any POGP $`\tau `$ there exists a set $`T`$ of generalized patterns such that a word $`\sigma `$ avoids $`\tau `$ if and only if $`\sigma `$ avoids all the patterns in $`T`$.

For example, if $`\tau =1^{}2^{}\text{-}3\text{-}1^{\prime \prime }`$, then to avoid $`\tau `$ is the same to avoid $`5`$ patterns, $`12\text{-}3\text{-}1`$, $`12\text{-}3\text{-}2`$, $`12\text{-}4\text{-}3`$, $`13\text{-}4\text{-}2`$, and $`23\text{-}4\text{-}1`$. Moreover, the following proposition holds:

###### Proposition 2.8.

Suppose $`\tau =\tau _1\text{-}a\text{-}\tau _2`$ (resp. $`\varphi =\varphi _1\text{-}\varphi _2`$) is a shuffle pattern (resp. a multi-pattern) such that $`\tau _1,\varphi _1[r_1]^\mathrm{}_1`$, $`\tau _2,\varphi _2[r_2]^\mathrm{}_2`$ and each letter of $`[r_1]`$ is incomparable with any letter of $`[r_2]`$. Also, without lose the generality, suppose $`r_1r_2`$. Then to avoid $`\tau `$ (resp. $`\varphi `$) is the same as to avoid $`{\displaystyle \underset{i=0}{\overset{r_2}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{r_1}{i}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{r_2}{i}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{r_1+r_2i}{r_1}}\right)`$ generalazed patterns. In particular, the number of generalized patterns does not depend on the lengths $`\mathrm{}_1`$ and $`\mathrm{}_2`$.

###### Proof.

Obviously, to prove the statement, we need to find the number of ways to make a total order on $`[r_1][r_2]`$ (the letter $`a`$ does not play any roll, since it is always the greatest letter). Any total order on $`[r_1][r_2]`$ is an alphabet that can consist of $`r_1+r_2i`$ letters, where $`i`$ is the number of letters in $`[r_2]`$ that supposed to coincide with some letters in $`[r_1]`$. Clearly, $`0ir_2`$ and we can choose coinciding letters in $`\left(\genfrac{}{}{0pt}{}{r_1}{i}\right)\left(\genfrac{}{}{0pt}{}{r_2}{i}\right)`$ ways. Now, after choosing the coinciding letters, we can make a total order in $`\left(\genfrac{}{}{0pt}{}{r_1+r_2i}{r_1}\right)`$ ways, which is given by \[Ki2, Theorem 8\]. ∎

## 3. The shuffle pattern

We recall that according to Definition 2.6, a shuffle pattern is a pattern of the form

$$\tau =\tau ^0\text{-}a_1\text{-}\tau ^1\text{-}a_2\text{-}\mathrm{}\text{-}\tau ^{s1}\text{-}a_s\text{-}\tau ^s,$$

where $`\{\tau ^0,\tau ^1,\mathrm{},\tau ^s\}`$ is a set of generalized patterns with no hyphens, $`a_1a_2\mathrm{}a_s`$ is a word of $`s`$ letters, for any $`i`$ and $`j`$ the letter $`a_i`$ is greater than any letter of $`\tau ^j`$ and for any $`ij`$ each letter of $`\tau ^i`$ is incomparable with any letter of $`\tau ^j`$.

Let us consider the shuffle pattern $`\varphi =\tau \text{-}\mathrm{}\text{-}\tau `$, where $`\mathrm{}`$ is the greatest letter in $`\varphi `$ and letters each letter in the left $`\tau `$ is incomparable with any letter in the right $`\tau `$.

###### Theorem 3.1.

Let $`\varphi `$ be the shuffle pattern $`\tau \text{-}\mathrm{}\text{-}\tau `$ described above. Then for all $`k\mathrm{}`$,

$$A_\varphi (x;k)=\frac{1}{(1xA_\tau (x;k1))^2}\left(A_\varphi (x;k1)xA_\tau ^2(x;k1)\right).$$

###### Proof.

We show how to get a recurrence relation on $`k`$ for $`A_\varphi (x;k)`$, which is the g.f. for the number of words in $`[k]^n(\varphi )`$. Suppose $`\sigma [k]^n(\varphi )`$ is such that it contains exactly $`d`$ copies of the letter $`k`$. If $`d=0`$ then the g.f. for the number of such words is $`A_\varphi (x;k1)`$. Assume that $`d1`$. Clearly, $`\sigma `$ can be written in the following form:

$$\sigma =\sigma ^0k\sigma ^1k\mathrm{}k\sigma ^d,$$

where $`\sigma ^j`$ is a $`\varphi `$-avoiding word on $`k1`$ letters, for $`j=0,1,\mathrm{},d`$. There are two possibilities: either $`\sigma ^j`$ avoids $`\tau `$ for all $`j`$, or there exists $`j_0`$ such that $`\sigma ^{j_0}`$ contains $`\tau `$ and for any $`jj_0`$, the word $`\sigma ^j`$ avoids $`\tau `$. In the first case, the number of such words is given by the g.f. $`x^dA_\tau ^{d+1}(x;k1)`$, whereas in the second case, by $`(d+1)x^dA_\tau ^d(x;k1)(A_\varphi (x;k1)A_\tau (x;k1))`$. In the last expression, the multiple $`(d+1)`$ is the number of ways to choose $`j`$, such that $`\sigma ^j`$ has an occurrence of $`\tau `$, and $`A_\varphi (x;k1)A_\tau (x;k1)`$ is the g.f. for the number of words avoiding $`\varphi `$ and containing $`\tau `$.

Therefore,

$$A_\varphi (x;k)=A_\varphi (x;k1)+\underset{d1}{}(d+1)x^dA_\tau ^d(x;k1)A_\varphi (x;k1)\underset{d1}{}dx^dA_\tau ^{d+1}(x;k1),$$

equivalently,

$$A_\varphi (x;k)=A_\varphi (x;k1)+A_\varphi (x;k1)\frac{2xA_\tau (x;k1)x^2A_\tau ^2(x;k1)}{(1xA_\tau (x;k1))^2}\frac{xA_\tau ^2(x;k1)}{(1xA_\tau (x;k1))^2}.$$

The rest is easy to check. ∎

###### Example 3.2.

Let $`\varphi =1^{}\text{-}2\text{-}1^{\prime \prime }`$. Here $`\tau =1`$, so $`A_\tau (x;k)=1`$ for all $`k1`$, since only the empty word avoids $`\tau `$. Hence, according to Theorem 3.1, we have

$$A_\varphi (x;k)=\frac{A_\varphi (x;k1)x}{(1x)^2},$$

which together with $`A_\varphi (x;1)=\frac{1}{1x}`$ (for any $`n`$ only the word $`\underset{ntimes}{\underset{}{11\mathrm{}1}}`$ avoids $`\varphi `$) gives Equation 2.1.

More generally, we consider a shuffle pattern of the form $`\tau ^0\text{-}\mathrm{}\text{-}\tau ^1`$, where $`\mathrm{}`$ is the greatest element of the pattern.

###### Theorem 3.3.

Let $`\varphi `$ be the shuffle pattern $`\tau \text{-}\mathrm{}\text{-}\nu `$. Then for all $`k\mathrm{}`$,

$$A_\varphi (x;k)=\frac{1}{(1xA_\tau (x;k1))(1xA_\nu (x;k1))}\left(A_\varphi (x;k1)xA_\tau (x;k1)A_\nu (x;k1)\right).$$

###### Proof.

We proceed as in the proof of Theorem 3.1. Suppose $`\sigma [k]^n(\varphi )`$ is such that it contains exactly $`d`$ copies of the letter $`k`$. If $`d=0`$ then the g.f. for the number of such words is $`A_\varphi (x;k1)`$. Assume that $`d1`$. Clearly, $`\sigma `$ can be written in the following form:

$$\sigma =\sigma ^0k\sigma ^1k\mathrm{}k\sigma ^d,$$

where $`\sigma ^j`$ is a $`\varphi `$-avoiding word on $`k1`$ letters, for $`j=0,1,\mathrm{},d`$. There are two possibilities: either $`\sigma ^j`$ avoids $`\tau `$ for all $`j`$, or there exists $`j_0`$ such that $`\sigma ^{j_0}`$ contains $`\tau `$, $`\sigma ^j`$ avoids $`\tau `$ for all $`j=0,1,\mathrm{},j_01`$ and $`\sigma ^j`$ avoids $`\nu `$ for any $`j=j_0+1,\mathrm{},d`$. In the first case, the number of such words is given by the g.f. $`x^dA_\tau ^{d+1}(x;k1)`$. In the second case, we have

$$x^d\underset{j=0}{\overset{d}{}}A_\tau ^j(x;k1)A_\nu ^{dj}(x;k1)(A_\varphi (x;k1)A_\tau (x;k1)).$$

Therefore, we get

$$\begin{array}{c}A_\varphi (x;k)=A_\varphi (x;k1)+A_\varphi (x;k1)\underset{d1}{}x^d\underset{j=0}{\overset{d}{}}A_\tau ^j(x;k1)A_\nu ^{dj}(x;k1)\hfill \\ \underset{d1}{}x^d\underset{j=1}{\overset{d}{}}A_\tau ^j(x;k1)A_\nu ^{d+1j}(x;k1),\hfill \end{array}$$

equivalently,

$$A_\varphi (x;k)=(A_\varphi (x;k1)xA_\tau (x;k1)A_\nu (x;k1))\underset{d0}{}x^d\underset{j=0}{\overset{d}{}}A_\tau ^j(x;k1)A_\nu ^{dj}(x;k1).$$

Hence, using the identity $`{\displaystyle \underset{n0}{}}x^n{\displaystyle \underset{j=0}{\overset{n}{}}}p^jq^{nj}={\displaystyle \frac{1}{(1xp)(1xq)}}`$ we get the desired result. ∎

We now give two corollaries to Theorem 3.3.

###### Corollary 3.4.

Let $`\varphi =\tau ^0\text{-}\mathrm{}\text{-}\tau ^1`$ be a shuffle pattern, and let $`f(\varphi )=f_1(\tau ^0)\text{-}\mathrm{}\text{-}f_2(\tau ^1)`$, where $`f_1`$ and $`f_2`$ are any trivial bijections. Then $`\varphi f(\varphi )`$.

###### Proof.

Using Theorem 3.3, and the fact that the number of words in $`[k]^n`$ avoiding $`\tau `$ (resp. $`\nu `$) and $`f_1(\tau )`$ (resp. $`f_2(\nu )`$) have the same generating functions, we get the desired result. ∎

###### Corollary 3.5.

For any shuffle pattern $`\tau \text{-}\mathrm{}\text{-}\nu `$, we have

$$\tau \text{-}\mathrm{}\text{-}\nu \nu \text{-}\mathrm{}\text{-}\tau .$$

###### Proof.

Corollary 3.4 yields that the shuffle pattern $`\tau \text{-}\mathrm{}\text{-}\nu `$ is equivalent to the pattern $`\tau \text{-}\mathrm{}\text{-}R(\nu )`$, which is equivalent to the pattern $`R(\tau \text{-}\mathrm{}\text{-}R(\nu ))=\nu \text{-}\mathrm{}\text{-}R(\tau )`$. Finally, we use Corollary 3.4 one more time to get the desired result. ∎

## 4. The multi-patterns

We recall that according to Definition 2.5, a multi-pattern is a pattern of the form $`\tau =\tau ^0\text{-}\tau ^1\text{-}\mathrm{}\text{-}\tau ^s`$, where $`\{\tau ^0,\tau ^1,\mathrm{},\tau ^s\}`$ is a set of generalized patterns with no hyphens and each letter of $`\tau ^i`$ is incomparable with any letter of $`\tau ^j`$ whenever $`ij`$.

The simplest non-trivial example of a multi-pattern is the multi-pattern $`\varphi =1\text{-}1^{}2^{}`$. To avoid $`\varphi `$ is the same as to avoid the patterns $`1\text{-}12`$, $`1\text{-}23`$, $`2\text{-}12`$, $`2\text{-}13`$, and $`3\text{-}12`$ simultaneously. To count the number of words in $`[k]^n(1\text{-}1^{}2^{})`$, we choose the leftmost letter of $`\sigma [k]^n(1\text{-}1^{}2^{})`$ in $`k`$ ways, and observe that all the other letters of $`\sigma `$ must be in a non-increasing order. Using \[BM1\], for all $`n1`$, we have

$$|[k]^n(1\text{-}1^{}2^{})|=k\left(\genfrac{}{}{0pt}{}{n+k2}{n1}\right).$$

The following theorem is an analogue to \[Ki2, Theorem 21\].

###### Theorem 4.1.

Let $`\tau =\tau ^0\text{-}\tau ^1`$ and $`\varphi =f_1(\tau ^0)\text{-}f_2(\tau ^1)`$, where $`f_1`$ and $`f_2`$ are any of the trivial bijections. Then $`\tau \varphi `$.

###### Proof.

First, let us prove that the pattern $`\tau =\tau ^0\text{-}\tau ^1`$ is equivalent to the pattern $`\varphi =\tau ^0\text{-}f(\tau ^1)`$, where $`f`$ is a trivial bijection. Suppose that $`\sigma =\sigma ^1\sigma ^2\sigma ^3[k]^n`$ avoids $`\tau `$ and $`\sigma ^1\sigma ^2`$ has exactly one occurrence of $`\tau ^0`$, namely $`\sigma ^2`$. Then $`\sigma ^3`$ must avoid $`\tau ^1`$, so $`f(\sigma ^3)`$ avoids $`f(\tau ^3)`$ and $`\sigma _f=\sigma ^1\sigma ^2f(\sigma ^3)`$ avoids $`\varphi `$. The converse is also true, if $`\sigma _f`$ avoids $`\varphi `$ then $`\sigma `$ avoids $`\tau `$. Since any word either avoids $`\tau ^0`$ or can be factored as above, we have a bijection between the class of words avoiding $`\tau `$ and the class of words avoiding $`\varphi `$. Thus $`\tau \varphi `$.

Now, we use the considerations above as well as the properties of trivial bijections to get

$$\begin{array}{c}\tau \tau ^0\text{-}f_2(\tau ^1)R(\tau ^0\text{-}f_2(\tau ^1))R(f_2(\tau ^1))\text{-}R(\tau ^0)\hfill \\ R(f_2(\tau ^1))\text{-}f_1(R(\tau ^0))R(f_2(\tau ^1))\text{-}R(f_1(\tau ^0))f_1(\tau ^0)\text{-}f_2(\tau ^1).\hfill \end{array}$$

Using Theorem 4.1, we get the following corollary, which is an analogue to \[Ki2, Corollary 22\].

###### Corollary 4.2.

The multi-pattern $`\tau ^0\text{-}\tau ^1`$ is equivalent to the multi-pattern $`\tau ^1\text{-}\tau ^0`$.

###### Proof.

From Theorem 4.1, using the properties of the trivial bijection $`R`$, we get

$$\tau ^0\text{-}\tau ^1\tau ^0\text{-}R(\tau ^1)R(R(\tau ^1))\text{-}R(\tau ^0)\tau ^1\text{-}R(R(\tau ^0))\tau ^1\text{-}\tau ^0.$$

Using induction on $`s`$, Corollary 4.2, and proceeding in the way proposed in \[Ki2, Theorem 23\], we get

###### Theorem 4.3.

Suppose we have multi-patterns $`\tau =\tau ^0\text{-}\tau ^1\text{-}\mathrm{}\text{-}\tau ^s`$ and $`\varphi =\varphi ^0\text{-}\varphi ^1\text{-}\mathrm{}\text{-}\varphi ^s`$, where $`\tau ^1\tau ^2\mathrm{}\tau ^s`$ is a permutation of $`\varphi ^1\varphi ^2\mathrm{}\varphi ^s`$. Then $`\tau \varphi `$.

The last theorem is an analogue to \[Ki2, Theorem 23\]. As a corollary to Theorem 4.3, using Theorem 4.1 and the idea of the proof of \[Ki2, Corollary 24\], we get the following corollary which is an analogue to \[Ki2, Corollary 24\].

###### Corollary 4.4.

Suppose we have multi-patterns $`\tau =\tau ^0\text{-}\tau ^1\text{-}\mathrm{}\text{-}\tau ^s`$ and $`\varphi =f_0(\tau ^0)\text{-}f_1(\tau ^1)\text{-}\mathrm{}\text{-}f_s(\tau ^s)`$, where $`f_i`$ is an arbitrary trivial bijection. Then $`\tau \varphi `$.

The following theorem is a good auxiliary tool for calculating the g.f. for the number of words that avoid a given POGP. For particular POGPs, it allows to reduce the problem to calculating the g.f. for the number of words that avoid another POGP which is shorter. We recall that $`A_\tau ^{}(x;k)`$ is the generating function for the number of words in $`[k]^n`$ that quasi-avoid the pattern $`\tau `$.

###### Theorem 4.5.

Suppose $`\tau =\tau ^0\text{-}\varphi `$, where $`\varphi `$ is an arbitrary POGP, and the letters of $`\tau ^0`$ are incomparable to the letters of $`\varphi `$. Then for all $`k1`$, we have

$$A_\tau (x;k)=A_{\tau ^0}(x;k)+A_\varphi (x;k)A_{\tau ^0}^{}(x;k).$$

###### Proof.

Suppose $`\sigma =\sigma ^1\sigma ^2\sigma ^3[k]^n`$ avoids the pattern $`\tau `$, where $`\sigma ^1\sigma ^2`$ quasi-avoids the pattern $`\tau ^0`$, and $`\sigma ^2`$ is the occurrence of $`\tau ^0`$. Clearly, $`\sigma ^3`$ must avoid $`\varphi `$. To find $`A_\tau (x;k)`$, we observe that there are two possibilities: either $`\sigma `$ avoids $`\tau ^0`$, or $`\sigma `$ does not avoid $`\tau ^0`$. In these cases, the g.f. for the number of such words is equal to $`A_{\tau ^0}(x;k)`$ and $`A_\varphi (x;k)A_{\tau ^0}^{}(x;k)`$ respectively (the second term came from the factorization above). Thus, the statement is true. ∎

###### Corollary 4.6.

Let $`\tau =\tau ^1\text{-}\tau ^2\text{-}\mathrm{}\text{-}\tau ^s`$ be a multi-pattern such that $`\tau ^j`$ is equal to either $`12`$ or $`21`$, for $`j=1,2,\mathrm{},s`$. Then

$$A_\tau (x;k)=\frac{1\left(1+\frac{kx1}{(1x)^k}\right)^s}{1kx}.$$

###### Proof.

According to \[BM2\], $`A_{12}(x;k)=A_{21}(x;k)=\frac{1}{(1x)^k}`$. Using Theorem 4.5, Proposition 2.4 and induction on $`s`$, we get the desired result. ∎

More generally, using Theorem 4.5 and Proposition 2.4, we get the following theorem that is the basis for calculating the number of words that avoid a multi-pattern, and therefore is the main result for multi-patterns in this paper.

###### Theorem 4.7.

Let $`\tau =\tau ^1\text{-}\tau ^2\text{-}\mathrm{}\text{-}\tau ^s`$ be a multi-pattern. Then

$$A_\tau (x;k)=\underset{j=1}{\overset{s}{}}A_{\tau ^j}(x;k)\underset{i=1}{\overset{j1}{}}((kx1)A_{\tau ^i}(x;k)+1).$$

## 5. The distribution of non-overlapping generalized patterns

A descent in a word $`\sigma [k]^n`$ is an $`i`$ such that $`\sigma _i>\sigma _{i+1}`$. Two descents $`i`$ and $`j`$ overlap if $`j=i+1`$. We define a new statistics, namely the maximum number of non-overlapping descents, or $`MND`$, in a word. For example, $`MND(33211)=1`$ whereas $`MND(13211143211)=3`$. One can find the distribution of this new statistic by using Corollary 4.6. This distribution is given in Example 5.2. However, we prove a more general theorem:

###### Theorem 5.1.

Let $`\tau `$ be a generalized pattern with no hyphens. Then for all $`k1`$,

$$\underset{n0}{}\underset{\sigma [k]^n}{}y^{N_\tau (\sigma )}x^n=\frac{A_\tau (x;k)}{1y((kx1)A_\tau (x;k)+1)},$$

where $`N_\tau (\sigma )`$ is the maximum number of non-overlapping occurrences of $`\tau `$ in $`\sigma `$.

###### Proof.

We fix the natural number $`s`$ and consider the multi-pattern $`\mathrm{\Phi }_s=\tau \text{-}\tau \text{-}\mathrm{}\text{-}\tau `$ with $`s`$ copies of $`\tau `$. If a word avoids $`\mathrm{\Phi }_s`$ then it has at most $`s1`$ non-overlapping occurrences of $`\tau `$. Theorem 4.7 yields

$$A_{\mathrm{\Phi }_s}(x;k)=\underset{j=1}{\overset{s}{}}A_\tau (x;k)\underset{i=1}{\overset{j1}{}}((kx1)A_\tau (x;k)+1).$$

So, the g.f. for the number of words that has exactly $`s`$ non-overlapping occurrences of the pattern $`\tau `$ is given by

$$A_{\mathrm{\Phi }_{s+1}}(x;k)A_{\mathrm{\Phi }_s}(x;k)=A_\tau (x;k)((kx1)A_\tau (x;k)+1)^s.$$

Hence,

$$\underset{n0}{}\underset{\sigma [k]^n}{}y^{N_\tau (\sigma )}x^n=\underset{s0}{}A_\tau (x;k)((kx1)A_\tau (x;k)+1)^s=\frac{A_\tau (x;k)}{1y((kx1)A_\tau (x;k)+1)}.$$

All of the following examples are corollaries to Theorem 5.1.

###### Example 5.2.

If we consider descents (the pattern $`12`$) then $`A_{12}(x;k)=\frac{1}{(1x)^k}`$ (see \[BM2\]), hence the distribution of $`MND`$ is given by the formula:

$$\underset{n0}{}\underset{\sigma [k]^n}{}y^{N_{12}(\sigma )}x^n=\frac{1}{(1x)^k+y(1kx(1x)^k)}.$$

###### Example 5.3.

The distribution of the maximum number of non-overlapping occurrences of the pattern $`122`$ is given by the formula:

$$\underset{n0}{}\underset{\sigma [k]^n}{}y^{N_{122}(\sigma )}x^n=\frac{x}{(1x^2)^k+x1+y(1kx^2(1x^2)^k)},$$

since according to \[BM3, Theorem 3.10\], $`A_{122}(x;k)=\frac{x}{(1x^2)^k(1x)}`$.

###### Example 5.4.

If we consider the pattern $`212`$ then $`A_{212}(x;k)=\left(1x{\displaystyle \underset{j=0}{\overset{k1}{}}}{\displaystyle \frac{1}{1+jx^2}}\right)^1`$ (see \[BM3, Theorem 3.12\]), hence the distribution of the maximum number of non-overlapping occurrences of the pattern $`212`$ is given by the formula:

$$\underset{n0}{}\underset{\sigma [k]^n}{}y^{N_{212}(\sigma )}x^n=\frac{1}{1x{\displaystyle \underset{j=0}{\overset{k1}{}}}{\displaystyle \frac{1}{1+jx^2}}+xy\left({\displaystyle \underset{j=0}{\overset{k1}{}}}{\displaystyle \frac{1}{1+jx^2}}k\right)}.$$

###### Example 5.5.

Using \[BM3, Theorem 3.13\], the distribution of the maximum number of non-overlapping occurrences of the pattern $`123`$ is given by the formula:

$$\underset{n0}{}\underset{\sigma [k]^n}{}y^{N_{123}(\sigma )}x^n=\frac{1}{{\displaystyle \underset{j=0}{\overset{k}{}}}a_j\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right)x^j+y\left(1kx{\displaystyle \underset{j=0}{\overset{k}{}}}a_j\left({\displaystyle \genfrac{}{}{0pt}{}{k}{j}}\right)x^j\right)},$$

where $`a_{3m}=1`$, $`a_{3m+1}=1`$, and $`a_{3m+2}=0`$, for all $`m0`$.