text:
# Recent BCP progress in Taiwan

## 1 Introduction

The collaboration on $`B`$ physics and CP violation (BCP) is one of the most active groups in Taiwan. In this talk I will briefly review theoretical progresses on BCP, which were made in Taiwan recently. For experimental review, please refer to Dr. Wang’s talk in this workshop. I will concentrate on five approaches to exclusive $`B`$ meson decays based on factorization assumption (FA), $`SU(3)`$ symmetry, perturbative QCD (PQCD) factorization theorem, QCD factorization (QCDF), and light-front QCD (LFQCD) formalism. Abundant results of exclusive $`B`$ meson decays have been produced and important dynamics has been explored.

The effective Hamiltonian for $`b`$ quark decays, for example, the $`bs`$ transition, is given by

$`H_{\mathrm{eff}}={\displaystyle \frac{G_F}{\sqrt{2}}}{\displaystyle \underset{q=u,c}{}}V_q\left[C_1(\mu )O_1^{(q)}(\mu )+C_2(\mu )O_2^{(q)}(\mu )+{\displaystyle \underset{i=3}{\overset{10}{}}}C_i(\mu )O_i(\mu )\right],`$ (1)
with the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements $`V_q=V_{qs}^{}V_{qb}`$. Using the unitarity condition, the CKM matrix elements for the penguin operators $`O_3`$-$`O_{10}`$ can also be expressed as $`V_u+V_c=V_t`$. The definitions of the operators $`O_i`$ are referred to.

According to the Wolfenstein parametrization, the CKM matrix upto $`𝒪(\lambda ^3)`$ is written as

$`\left(\begin{array}{ccc}V_{ud}& V_{us}& V_{ub}\\ V_{cd}& V_{cs}& V_{cb}\\ V_{td}& V_{ts}& V_{tb}\end{array}\right)=\left(\begin{array}{ccc}1\frac{\lambda ^2}{2}& \lambda & A\lambda ^3(\rho i\eta )\\ \lambda & 1\frac{\lambda ^2}{2}& A\lambda ^2\\ A\lambda ^3(1\rho i\eta )& A\lambda ^2& 1\end{array}\right),`$ (2)
with the parameters $`\lambda =0.2196\pm 0.0023`$, $`A=0.819\pm 0.035`$, and $`R_b\sqrt{\rho ^2+\eta ^2}=0.41\pm 0.07`$. The unitarity angle $`\varphi _3`$ is defined via

$$V_{ub}=|V_{ub}|\mathrm{exp}(i\varphi _3).$$
(3)
One of the important missions of $`B`$ fatories is to determine the angles $`\varphi _1`$, $`\varphi _2`$ and $`\varphi _3`$. The angle $`\varphi _1`$ can be extracted from the CP asymmetry in the $`BJ/\psi K_S`$ decays, which arises from the $`B`$-$`\overline{B}`$ mixing. Due to similar mechanism of CP asymmetry, the decays $`B^0\pi ^+\pi ^{}`$ are appropriate for the extraction of the angle $`\varphi _2`$. However, these modes contain penguin contributions such that the extraction suffers large uncertainty. It has been proposed that the angle $`\varphi _3`$ can be determined from the decays $`BK\pi `$, $`\pi \pi `$. Contributions to these modes involve interference between penguin and tree amplitudes, and their analyses rely heavily on theoretical formulations.

## 2 Factorization Assumption

The conventional approach to exclusive nonleptonic $`B`$ meson decays is based on FA, in which nonfactorizable and annihilation contributions are neglected and final-state-interaction (FSI) effects are assumed to be absent. Factorizable contributions are expressed as products of Wilson coefficients and hadronic form factors, which are then parametrized by models. With these approximations, the FA approach is simple and provides qualitative estimation of various decay branching ratios.

To extract the angle $`\varphi _3`$, we consider the ratios $`R`$ and $`R_\pi `$ defined by

$`R={\displaystyle \frac{\mathrm{Br}(B_d^0K^\pm \pi ^{})}{\mathrm{Br}(B^\pm K^0\pi ^\pm )}},`$ $`R_\pi ={\displaystyle \frac{\mathrm{Br}(B_d^0K^\pm \pi ^{})}{\mathrm{Br}(B_d^0\pi ^\pm \pi ^{})}},`$ (4)
where $`\mathrm{Br}(B_d^0K^\pm \pi ^{})`$ represents the CP average of the branching ratios $`\mathrm{Br}(B_d^0K^+\pi ^{})`$ and $`\mathrm{Br}(\overline{B}_d^0K^{}\pi ^+)`$, and the definition of $`\mathrm{Br}(B^\pm K^0\pi ^\pm )`$ is similar. It has been shown that the data $`R1`$ imply $`\varphi _390^o`$.

To explain the data of $`R_\pi 4`$, a large angle $`\varphi _3130^o`$ must be postulated. It is easy to observe from Eqs. (2) and (LABEL:ckp) that the products of the CKM matrix elements $`V_{us}^{}V_{ub}`$ and $`V_{ud}^{}V_{ub}`$ have the same weak phase, and that the real parts of $`V_{ts}^{}V_{tb}`$ and $`V_{td}^{}V_{ub}`$ are opposite in sign. That is, the tree-penguin interference in the decays $`BK\pi `$ and $`B\pi \pi `$ is anti-correlated. A $`\varphi _3>90^o`$ then leads to constructive interference between the tree and penguin contributions in $`BK\pi `$, and a large $`R_\pi `$. The determination $`\varphi _3114^o`$ from the global fit to charmless $`B`$ meson decays, located between the two extreme cases $`90^o`$ and $`130^o`$, is then understood. On the other hand, a large $`B\rho `$ form factor $`A_0^{B\rho }0.48`$ has been extracted, which accounts for the large $`B\rho \pi `$ branching ratios. In the modified FA approach with an effective number of colors $`N_c^{\mathrm{eff}}`$, a large unitarity angle $`\varphi _3105^o`$ is also concluded.

The above $`\varphi _3`$, located in the second quadrant, contradict the extraction from other measurements, such as $`B_s`$-mixing. The best fit to experimental data of semileptonic $`B`$ meson decays, $`B`$-$`\overline{B}`$ mixing, and $`ϵ_K`$ indicates that $`\varphi _3`$ is located in the first quadrant. An improvement of FA has been considered, in which possible strong phases produced via FSI are introduced as arbitrary parameters. Performing the best fit to data, a large $`\varphi _3`$ is still required and strong phases are found to be large, which generate significant CP asymmetries in $`K\pi `$ and $`\pi \pi `$ modes.

## 3 $`SU(3)`$ Symmetry

The model-dependent determination of the angle $`\varphi _3`$ from FA seems not to be satisfactory. A more model-independent approach based on $`SU(3)`$ symmetry has been proposed, in which the light quarks $`u`$, $`d`$ and $`s`$ form a $`SU(3)`$ triplet, while the heavy quarks $`c`$, $`b`$, and $`t`$ form $`SU(3)`$ singlets. According to the above assumption, the $`B`$ mesons $`B_u`$, $`B_d`$, and $`B_s`$ form a $`SU(3)`$ triplet at the hadronic level. Pseudo-scalar mesons $`P`$ and vector mesons $`V`$ also possess definite $`SU(3)`$ structures,

$`M_j^i`$ $`=`$ $`\left(\begin{array}{ccc}\frac{\pi ^0}{\sqrt{2}}\frac{\eta _8}{\sqrt{6}}+\frac{\eta _1}{\sqrt{3}}& \pi ^{}& K^{}\\ \pi ^+& \frac{\pi ^0}{\sqrt{2}}\frac{\eta _8}{\sqrt{6}}+\frac{\eta _1}{\sqrt{3}}& \overline{K}^0\\ K^+& K^0& \frac{2\eta _8}{\sqrt{6}}+\frac{\eta _1}{\sqrt{3}}\end{array}\right),`$ (5)
$`V_j^i`$ $`=`$ $`\left(\begin{array}{ccc}\frac{\rho ^0}{\sqrt{2}}+\frac{\omega }{\sqrt{2}}& \rho ^{}& K^{}\\ \rho ^+& \frac{\rho ^0}{\sqrt{2}}+\frac{\omega }{\sqrt{2}}& \overline{K}^0\\ K^+& K^0& \varphi \end{array}\right).`$ (6)
Similarly, the four-fermion operators in the weak Hamiltonian can be decomposed into operators with definite $`SU(3)`$ structures. For example, the penguin operators $`O_{36}`$ are labelled as $`\overline{3}`$ states,

$`\overline{q}b(\overline{u}u+\overline{d}d+\overline{s}s)\overline{3},`$ (7)
since $`\overline{q}b`$ forms a triplet $`\overline{3}`$ and $`\overline{u}u+\overline{d}d+\overline{s}s`$ froms a singlet. The operators $`O_{1,2}`$ are written as

$`\overline{q}u\overline{u}b`$ $``$ $`\overline{3}\times 3\times \overline{3}=\overline{3}+\overline{3}+6+\overline{15}.`$ (8)
Following Eqs. (7) and (8), the effective Hamiltonian in Eq. (1) is decomposed into operators carrying different $`SU(3)`$ structures, such as $`H(\overline{3})`$, $`H(6)`$ and $`H(\overline{15})`$, whose coefficients are the linear combinations of the Wilson coefficients.

Employing the above results, we formulate various decay amlitudes of the tree and annihilation topologies for $`BPP`$ modes. For example, the parameter $`C_{\overline{3}}`$ associated with the contraction $`B_iM_k^iM_j^kH(\overline{3})^j`$ represents a tree amplitude. The parameter $`A_{\overline{3}}`$ associated with the contraction $`B_iH(\overline{3})^iM_l^kM_k^l`$ represents an annihilation amplitude. Collecting all expressions for the branching ratios and the CP asymmetries of the $`B\pi \pi `$, $`K\pi `$ and $`KK`$ modes, there are totally 13 free parameters. This number is too big for a global analysis of currently available data. As an approximation, annihilation contributions are neglected. 8 parameters, the absolute values of $`C_{\overline{3}}^T`$, $`C_{\overline{3}}^P`$, $`C_6^T`$ and $`C_{\overline{1}5}^T`$, the phases $`\delta _{\overline{3}}^P`$, $`\delta _6^T`$, and $`\delta _{\overline{1}5}^T`$, and the CKM phase $`\varphi _3`$, are then left, where $`T`$ ($`P`$) denotes the tree operators $`O_{1,2}`$ (the penguin operators $`O_{310}`$).

The best fit to data gives

$`\varphi _3=70^o,\rho =0.17,\eta =0.37,`$
$`C_{\overline{3}}^T=0.28,C_{\overline{3}}^P=0.14,C_6^T=0.33,C_{\overline{1}5}^T=0.14,`$
$`\delta _{\overline{3}}^P=12^o,\delta _6^T=6^o,\delta _{\overline{1}5}^T=74^o.`$ (9)
If the $`SU(3)`$ symmetry breaking effect from $`f_K/f_\pi 1`$ is taken into account, $`f_K`$ ($`f_\pi `$) being the kaon (pion) decay constant, the results are shifted only a bit. Hence, there is no indication that $`\varphi _3`$ should be located in the second quadrant. Certainly, the allowed range of the above parameters is still large.

## 4 Perturbative QCD

It has been shown that the decay amplitudes and strong phases discussed in the previous sections can be evaluated in the PQCD framework, and that it is possible to extract $`\varphi _3`$ from the $`BK\pi `$ data. According to PQCD factorization theorem, a $`B`$ meson decay amplitude is expressed as convolution of a hard $`b`$ quark decay amplitude with meson wave functions. A meson wave function, absorbing nonperturbative dynamics of a QCD process, is not calculable, while a hard amplitude is.

In perturbation theory nonperturbative dynamics is reflected by infrared divergences in radiative corrections. It has been proved to all orders that these infrared divergences can be separated and absorbed into meson wave functions. A formal definition of wave functions as matrix elements of nonlocal operators has been constructed, which, if evaluated perturbatively, reproduces the infrared divergences. The gauge invariance of the above factorization has been proved in. A meson wave function must be determined by nonperturbative means, such as lattice gauge theory and QCD sum rules, or extracted from experimental data. A salient feature of PQCD factorization theorem is the universality of nonperturbative wave functions. Because of universality, a $`B`$ meson wave function extracted from some decay modes can be employed to make predictions for other modes. This is the reason PQCD factorization theorem possesses a predictive power.

In the practical calculation small parton transverse momenta $`k_T`$ are included, which are essential for smearing the end-point singularities from small momentum fractions. Because of the inclusion of $`k_T`$, double logarithms $`\mathrm{ln}^2(Pb)`$ are generated from the overlap of collinear and soft enhancements in radiative corrections to meson wave functions, where $`P`$ denotes the dominant light-cone component of a meson momentum and $`b`$ is a variable conjugate to $`k_T`$. The resummation of these double logarithms leads to a Sudakov form factor $`\mathrm{exp}[s(P,b)]`$, which suppresses the long-distance contributions in the large $`b`$ region, and vanishes as $`b=1/\mathrm{\Lambda }`$, $`\mathrm{\Lambda }\mathrm{\Lambda }_{\mathrm{QCD}}`$ being the QCD scale. This suppression guarantees the applicability of PQCD to exclusive decays around the energy scale of the $`b`$ quark mass.

The hard amplitude contains all possible Feynman diagrams, such as factorizable diagrams, where hard gluons attach the valence quarks in the same meson, and nonfactorizable diagrams, where hard gluons attach the valence quarks in different mesons. The annihilation topology is also included, and classified into factorizable or nonfactorizable one. Therefore, FA for two-body $`B`$ meson decays is not necessary. It has been shown that factorizable annihilation contributions are in fact important, and give large strong phases in PQCD.

We emphasize that the hard amplitude is characterized by the virtuality of internal particles, $`t\sqrt{\overline{\mathrm{\Lambda }}M_B}1.5`$ GeV, $`\overline{\mathrm{\Lambda }}=M_Bm_b`$. The RG evolution of the Wilson coefficients $`C_{4,6}(t)`$ dramatically increase as $`t<M_B/2`$, such that penguin contributions are enhanced. With this penguin enhancement, the observed branching ratios of the $`BK\pi `$ and $`B\pi \pi `$ decays can be explained in PQCD using a smaller angle $`\varphi _3=90^o`$. That is, the data of $`R_\pi `$ do not demand large $`\varphi _3`$. Such a dynamical enhancement of penguin contributions does not exist in the FA approach.

Our predictions for the branching ratio of each $`K\pi `$ mode corresponding to $`\varphi _3=90^o`$,

$`\mathrm{Br}(B^+K^0\pi ^+)=21.72\times 10^6,\mathrm{Br}(B^{}\overline{K}^0\pi ^{})=21.25\times 10^6,`$
$`\mathrm{Br}(B_d^0K^+\pi ^{})=24.19\times 10^6,\mathrm{Br}(\overline{B}_d^0K^{}\pi ^+)=16.84\times 10^6,`$
$`\mathrm{Br}(B^+K^+\pi ^0)=14.44\times 10^6,\mathrm{Br}(B^{}K^{}\pi ^0)=10.65\times 10^6,`$
$`\mathrm{Br}(B_d^0K^0\pi ^0)=11.23\times 10^6,\mathrm{Br}(\overline{B}_d^0\overline{K}^0\pi ^0)=11.84\times 10^6,`$ (10)
are consistent with the CLEO data,

$`\mathrm{Br}(B^\pm K^0\pi ^\pm )=(18.2_{4.0}^{+4.6}\pm 1.6)\times 10^6,`$
$`\mathrm{Br}(B_d^0K^\pm \pi ^{})=(17.2_{2.4}^{+2.5}\pm 1.2)\times 10^6,`$
$`\mathrm{Br}(B^\pm K^\pm \pi ^0)=(11.6_{2.71.3}^{+3.0+1.4})\times 10^6,`$
$`\mathrm{Br}(B_d^0K^0\pi ^0)=(14.6_{5.13.3}^{+5.9+2.4})\times 10^6,`$
$`A_{CP}(B_d^0K^\pm \pi ^{})=0.04\pm 0.16,`$
$`A_{CP}(B^\pm K^0\pi ^\pm )=0.17\pm 0.24.`$ (11)
In the above expressions $`B(B_d^0K^\pm \pi ^{})`$ represents the CP average of the branching ratios $`B(B_d^0K^+\pi ^{})`$ and $`B(\overline{B}_d^0K^{}\pi ^+)`$.

## 5 QCD Factorization

Recently, Beneke, Buchalla, Neubert, and Sachrajda proposed the QCDF formalism for two-body nonleptonic $`B`$ meson decays. They claimed that factorizable contributions, for example, the form factor $`F^{B\pi }`$ in the $`B\pi \pi `$ decays, are not calculable in PQCD, but nonfactorzable contributions are in the heavy quark limit. Hence, the former are treated in the same way as FA, and expressed as products of Wilson coefficients and $`F^{B\pi }`$. The latter, calculated as in the PQCD approach, are written as the convolutions of hard amplitudes with three $`(B,\pi ,\pi )`$ meson wave functions. Annihilation diagrams are neglected as in FA, but can be included as $`1/M_B`$ correction. Values of form factors at maximal recoil and nonperturbative meson wave functions are treated as free parameters.

Here I mention some essential differences between the QCDF and PQCD approaches. For more detailed comparisions, refer to. Because of the neglect of annihilation diagrams in QCDF, strong phases and CP asymmetries are much smaller than those predicted in PQCD. In QCDF the leading-order diagrams are those that contain vertex corrections to the four-fermion operators. These diagrams, however, appear at the next-to-leading order in PQCD. This difference implies different characteristic scales in the two approaches: the former is characterized by the $`b`$ quark mass $`m_b`$, while the latter is characterized by the virtuality $`t`$ of internal particles, which leads to the penguin enhancement emphasized above. Without penguin enhancement, a large $`\varphi _3`$ is still necessary to account for the large ratio $`R_\pi `$.

The $`B\varphi K`$ decays have been analyzed in the QCDF formalism, and branching ratios much smaller than experimental data have been obtained. Since these modes are dominated by penguin contributions, the penguin enhancement may be crucial for explaining the data.

## 6 Light-front QCD

The evaluation of a form factor is simple in the LFQCD formalism, which is written as an overlap integral of initial- and final-state meson wave functions. Various form factors have been computed, such as the $`B\pi `$, $`\rho `$, $`K`$ and $`K^{}`$ form factors  and the $`\mathrm{\Lambda }_b\mathrm{\Lambda }`$ form factors. The results have been employed to predict the decay spectra of the $`BK\mu \mu (\tau \tau )`$ and $`\mathrm{\Lambda }_b\mathrm{\Lambda }\mu \mu (\tau \tau )`$ modes. This formalism has been also applied to the radiative leptonic $`B`$ meson decays $`Bl^+l^{}\gamma `$. These predictions can be compared with data in the future.

## 7 Conclusion

In this talk I have briefly summarized the theoretical progresses on exclusive $`B`$ meson decays, which were made by Taiwan BCP community recently. With the active collaboration, more progresses are expected in the near future.

## Acknowledgments

This work was supported in part by the National Science Council of R.O.C. under the Grant No. NSC-89-2112-M-006-033, and in part by Grant-in Aid for Special Project Research (Physics of CP Violation) and by Grant-in Aid for Scientific Exchange from Ministry of Education, Science and Culture of Japan.


format: text

text:
# Complex diffusion Monte-Carlo method: test by the simulation of the 2𝐷 fermions

## I Introduction

There are several objects of the quantum mechanics where the investigated system’s wave function is essentially complex(having imaginary part). The examples of these objects are a two-dimensional system of electrons in the external uniform magnetic field when vector potential has a central-symmetric form and system of anyons. The complexity of wave function does not give a possibility for the simulations of such kind of systems by using well-known Green Function Monte-Carlo method (see review ). This method essentially demand the reality of the system’s wave function which considered as probability weight in during stochastic process.

There had been undertaken the several attempts  for the construction a method for the simulation of quantum systems with complex wave function. Authors of the works  suggested to take as a probability weight the module of the complex distribution function. All quantities are calculated in  by averaging over this complex distribution function. The phase of this function $`\gamma `$ is taken accordingly $`\gamma =\gamma _i`$, where $`\gamma _i`$ is the phase at the $`i`$-time step of the stochastic process.

Another quantum Monte Carlo method using algorithm without branching for the simulating of complex problems was developed in . The main trouble of this approach was an increasing of the statistical error as a function of the whole time of the simulation.

In the previous paper  it was suggested new Complex Diffusion Monte-Carlo (CDMC) method allowing to simulate the quantum systems with complex wave function. Our method is rather close to the fixed - phase Diffusion Monte-Carlo method developed for the simulation of the two-dimensional electrons in the magnetic field . In this work there had been an exact simulation only of the module of system’s wave function. The phase of the wave function was considered as fixed and equal to phase of the Laughlin’s wave function . In contrast, our method include also the simulation of the phase factor of the wave function.

The basic difficulty of the numerical simulation of the fermions is essentially the same as for the systems with the complex wave function. Their wave function can change the sign and therefore can not be used as the probability weight in the simulation process.

For the simulation of the continuum (not on the lattice) fermionic systems it were developed widely used fixed node Monte-Carlo method ( see also reviews ) and recently proposed constrained path Monte-Carlo method . In these both methods it was assumed the restriction on the random walks connected with the uncertainties in the space localization of the wave function node surfaces. The comparison of these methods was done in .

In the $`2D`$ space we have unique tool – so-called anyons, particles (bosons or fermions) with additional gauge interaction which provide needed statistical properties . So, by tuning the coupling constant we may arrive to the fermions starting from the bosons. It looks very attractive to apply the anyons for the simulation of the $`2D`$ fermions, because statistical vector potential that gives anyonic property is relatively smooth and continuous . The remaining main problem is that the wave function of anyonic system is essentially complex, i.e. it contains an imaginary part.

The CDMC method was successfully tested by the simulation of ground state of one electron in magnetic field . In the present work we briefly outline the idea of CDMC (the details of the algorithm one can find in ) and provide a further test it by simulating the ground state of fermions in $`2D`$ parabolic well and comparing with well known analytical answer for the ground state energy as a function of the number of fermions. The simulation is done for the fermionic systems with the number of particles from two to ten. It was found a good agreement of the numerical results with analytical ones for the numbers of the particles great than four. This observation of the fermionic simulation allows us to hope that the CDMC is a good tool for the simulation of $`2D`$ fermion systems.

## II Model and CDMC method for the simulation

Hamiltonian of anyons in the $`2D`$ parabolic well has a form:

$$\widehat{H}=\frac{1}{2M}\underset{i=1}{\overset{N}{}}(\stackrel{}{p}_i+\stackrel{}{A}(\stackrel{}{r}_i))^2+\frac{Mw_o^2}{2}\underset{i=1}{\overset{N}{}}\stackrel{}{r}_i^2.$$
(1)
Here $`M`$ is the mass of particle, $`\stackrel{}{p}=i\mathrm{}\stackrel{}{}`$, where $`\stackrel{}{}=\stackrel{}{i}{\displaystyle \frac{}{x}}+\stackrel{}{j}{\displaystyle \frac{}{y}}`$, $`w_o`$ stands for characteristic frequency of free particles in $`2D`$ parabolic well and $`\stackrel{}{r}_i`$ is radius vector of $`i`$-th particle. The number of particles is $`N`$.

Vector potential for anyons $`\stackrel{}{A}(\stackrel{}{r}_i)`$  in (1) is

$$\stackrel{}{A}(\stackrel{}{r}_i)=\mathrm{}\nu \underset{j>i}{\overset{N}{}}\frac{\stackrel{}{z}\times \stackrel{}{r}_{ij}}{|\stackrel{}{r}_{ij}|^2}.$$
(2)
Here $`\stackrel{}{z}`$ is unit vector perpendicular to $`2D`$ plane and $`\nu `$ characterizes the form of fractional statistics ( spin of the anyon ). In this work we want to employ a bosonic representation of anyons, therefore $`\nu =0`$ corresponds to noninteracting bosons and $`\nu =1`$ is the case of our interest - free fermions.

In CDMC a Monte Carlo simulation is carried out with the complex distribution function

$$f(\stackrel{}{R},t)=\mathrm{\Psi }_T^{}(\stackrel{}{R})\mathrm{\Psi }(\stackrel{}{R},t),$$
(3)
where $`\stackrel{}{R}`$ stands for the coordinates of all particles, $`\mathrm{\Psi }_T(\stackrel{}{R})`$ is a trial wave function and the wave function $`\mathrm{\Psi }(\stackrel{}{R},t)`$ satisfies a Schrödinger equation with imaginary time (expressed in $`\mathrm{}`$ units)

$$\frac{\mathrm{\Psi }(\stackrel{}{R},t)}{t}=(\widehat{H}E_T)\mathrm{\Psi }(\stackrel{}{R},t).$$
(4)
With a suitable choice of a trial energy $`E_T`$ it is possible to arrange the convergence $`\mathrm{\Psi }(\stackrel{}{R},t)\mathrm{\Psi }_0(\stackrel{}{R})`$ \- an exact stationary wave function of ground state of Hamiltonian (1) at $`t\mathrm{}`$ .

For a bosonic representation of anyons we take (conjugated) trial wave function in the form:

$$\mathrm{\Psi }_T^{}(\stackrel{}{R})=\underset{i=1}{\overset{N}{}}\mathrm{\Psi }_T^{}(\stackrel{}{r}_i).$$
(5)
Good choice of single particle trial wave function $`\mathrm{\Psi }_T^{}(\stackrel{}{r}_i)`$ ( see below ) as energy $`E_T`$ must provide a convergence of simulation process to exact ground state.

The distribution function $`ff(\stackrel{}{R},t)`$ satisfies a diffusion equation:

$`{\displaystyle \frac{f}{t}}`$ $`=`$ $`D{\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{\Delta }_if+D{\displaystyle \underset{i=1}{\overset{N}{}}}\stackrel{}{}_i(f\mathrm{Re}\stackrel{}{F}_Q(\stackrel{}{r}_i))+`$ (7)
$`+i{\displaystyle \underset{i=1}{\overset{N}{}}}[\stackrel{}{}_i(Df\mathrm{Im}\stackrel{}{F}_Q(\stackrel{}{r}_i)){\displaystyle \frac{\mathrm{}}{M}}\stackrel{}{A}(\stackrel{}{r}_i)\stackrel{}{}_if]+(E_L(\stackrel{}{R})E_T)f.`$
Here $`D={\displaystyle \frac{\mathrm{}^2}{2M}}`$, $`\mathrm{\Delta }_i=\stackrel{}{}_i^2`$. When a time step integration of diffusion equation (7), $`\tau `$, goes to zero, then Green function for this equation has a form:

$`G(\stackrel{}{R},\stackrel{}{R}^{};\tau )`$ $`=`$ $`{\displaystyle \frac{\mathrm{exp}[D\tau {\displaystyle \underset{i=1}{\overset{N}{}}}\stackrel{}{A}_Q^2(\stackrel{}{r}_i,\stackrel{}{r}_i^{})]}{(4\pi D\tau )^N}}\times `$ (10)
$`\times \mathrm{exp}[{\displaystyle \frac{{\displaystyle \underset{i=1}{\overset{N}{}}}(\stackrel{}{r}_i\stackrel{}{r}_i^{}D\tau \mathrm{Re}\stackrel{}{F}_Q(\stackrel{}{r}_i^{}))^2}{4D\tau }}]\mathrm{exp}[\tau (E_L(\stackrel{}{R})E_T)]\times `$
$`\times \mathrm{exp}\left[i{\displaystyle \underset{i=1}{\overset{N}{}}}\stackrel{}{A}_Q(\stackrel{}{r}_i,\stackrel{}{r}_i^{})(\stackrel{}{r}_i\stackrel{}{r}_i^{}D\tau \mathrm{Re}\stackrel{}{F}_Q(\stackrel{}{r}_i^{}))\right].`$
In the expressions (7) and (10)

$$\stackrel{}{F}_Q(\stackrel{}{r}_i)=2\mathrm{\Psi }_T^1(\stackrel{}{r}_i)\stackrel{}{}_i\mathrm{\Psi }_T^{}(\stackrel{}{r}_i),$$
(11)
$$E_L(\stackrel{}{R})=\mathrm{\Psi }_T^1(\stackrel{}{R})\widehat{H}^{}\mathrm{\Psi }_T^{}(\stackrel{}{R}),$$
(12)
$$\widehat{H}^{}=\frac{1}{2M}\underset{i=1}{\overset{N}{}}(\stackrel{}{p}_i\stackrel{}{A}(\stackrel{}{r}_i))^2+\frac{Mw_o^2}{2}\underset{i=1}{\overset{N}{}}\stackrel{}{r}_i^2.$$
(13)
In general case $`\mathrm{\Psi }_T^{}(\stackrel{}{r}_i)`$ is a complex wave function and

$$\stackrel{}{F}_Q(\stackrel{}{r}_i)=\mathrm{Re}\stackrel{}{F}_Q(\stackrel{}{r}_i)+i\mathrm{Im}\stackrel{}{F}_Q(\stackrel{}{r}_i).$$
(14)
In (10) (see also ) we have introduced a new quantity

$$\stackrel{}{A}_Q(\stackrel{}{r}_i,\stackrel{}{r}_i^{})=\frac{1}{2}\mathrm{Im}\stackrel{}{F}_Q(\stackrel{}{r}_i^{})\frac{1}{\mathrm{}}\stackrel{}{A}(\stackrel{}{r}_i).$$
(15)
In the expressions (10) and (15) vectors $`\stackrel{}{R}`$ and $`\stackrel{}{r}_i`$ correspond to time point $`t+\tau `$ and vectors $`\stackrel{}{R}^{}`$, $`\stackrel{}{r}^{}`$ to time point $`t`$.

One can see that the Green function $`G(\stackrel{}{R},\stackrel{}{R}^{};\tau )`$ (10) as and distribution function $`f`$ in (3) is a complex. These both quantities are satisfying to the usual integral relation

$$f(\stackrel{}{R},t+\tau )=𝑑\stackrel{}{R}^{}G(\stackrel{}{R},\stackrel{}{R}^{};\tau )f(\stackrel{}{R}^{},t).$$
(16)
From relation (16) it is following that the module and phase of the distribution function at consequent time point are determining by the module and the phase of the Green function and of ones of the distribution function at previous time point of integration of diffusion equation (7).

As the energy $`E_L(\stackrel{}{R})`$ is a complex, real and imaginary part of $`E_L(\stackrel{}{R})`$ contribute separately to ones of the Green function. Therefore the expression of the last two exponents in (10) has a form:

$`\begin{array}{c}\mathrm{exp}[\tau (\mathrm{Re}E_L(\stackrel{}{R})E_T)]\times \\ \times \mathrm{exp}\left[i{\displaystyle \underset{i=1}{\overset{N}{}}}\stackrel{}{A}_Q(\stackrel{}{r}_i,\stackrel{}{r}_i^{})(\stackrel{}{r}_i\stackrel{}{r}_i^{}D\tau \mathrm{Re}\stackrel{}{F}_Q(\stackrel{}{r}_i^{}))i\tau \mathrm{Im}E_L(\stackrel{}{R})\right].\end{array}`$ (19)
The mean value of the some quantity $`B(\stackrel{}{R}(t))`$ on the complex distribution function $`f(\stackrel{}{R},t)`$ at the time $`t`$ of the simulating process is calculated by the formula (see )

$$<B(t)>=\frac{{\displaystyle \underset{i=1}{\overset{M}{}}}\mathrm{exp}[i\gamma (\stackrel{}{R}_i(t)]B(\stackrel{}{R}_i(t))}{{\displaystyle \underset{i=1}{\overset{M}{}}}\mathrm{exp}[i\gamma (\stackrel{}{R}_i(t)]}.$$
(20)
Here $`\stackrel{}{R}_i(t)`$ is the coordinates of all particles from the $`i`$-th configuration of the $`M`$ \- configurations at time $`t`$. The coordinates $`\stackrel{}{R}_i(t)`$ of particles in (20) are weighted with probability $`|f(\stackrel{}{R}_i,t)|`$ and the quantity $`\gamma (\stackrel{}{R}_i(t))`$ is a phase of the distribution function $`f(\stackrel{}{R}_i,t)`$.

From integral expression (16) it is seen that phase of the distribution function $`\gamma (\stackrel{}{R}_i(t+\tau ))`$ connected with the phase of the Green function $`\gamma _G(\stackrel{}{R}_i(t+\tau ),\stackrel{}{R}_j(t))`$ and the phase $`\gamma (\stackrel{}{R}_j(t))`$ of the distribution function by the relation:

$$\gamma (\stackrel{}{R}_i(t+\tau ))=\gamma _G(\stackrel{}{R}_i(t+\tau ),\stackrel{}{R}_j(t))+\gamma (\stackrel{}{R}_j(t)).$$
(21)
Here $`\stackrel{}{R}_j(t)`$ is $`j`$\- th configuration at time $`t`$. So, we can write

$`\begin{array}{c}<B(t+\tau )>={\displaystyle \frac{{\displaystyle \underset{i=1}{\overset{M}{}}}e^{i\gamma _G(\stackrel{}{R}_i(t+\tau ),\stackrel{}{R}_j(t))+i\gamma (\stackrel{}{R}_j(t))}B(\stackrel{}{R}_i(t+\tau ))}{{\displaystyle \underset{i=1}{\overset{M}{}}}e^{i\gamma _G(\stackrel{}{R}_i(t+\tau ),\stackrel{}{R}_j(t))+i\gamma (\stackrel{}{R}_j(t))}}}=\\ ={\displaystyle \frac{{\displaystyle \underset{i=1}{\overset{M}{}}}e^{i\gamma _G(\stackrel{}{R}_i(t+\tau ),\stackrel{}{R}_j(t))}B(\stackrel{}{R}_i(t+\tau ))}{{\displaystyle \underset{i=1}{\overset{M}{}}}e^{i\gamma _G(\stackrel{}{R}_i(t+\tau ),\stackrel{}{R}_j(t))}}}.\end{array}`$ (24)
The expression (24) is hold, because all quantities in numerator and denominator are weighted with probability $`|f(\stackrel{}{R}_i(t+\tau ),t+\tau )|`$, i.e. at time moment $`t+\tau `$ with new configurations $`\stackrel{}{R}_i(t+\tau )`$. Hence, a mean quantity $`<B(t+\tau )>`$ is determined only by phase $`\gamma _G`$ of the Green function. It is a key conclusion, providing the calculation of mean quantities in our CDMC.

Next, by taking into account the expression for the simulation of the random displacement for the particle $`\stackrel{}{r}_i=\stackrel{}{r}_i^{}+D\tau \mathrm{Re}\stackrel{}{F}_Q(\stackrel{}{r}_i^{})+2(D\tau )^{1/2}\chi `$ (see  and ), where $`\chi `$ is a gaussian random number, the expressions for $`\stackrel{}{F}_Q(\stackrel{}{r}_i)`$ (11), for $`E_L(\stackrel{}{R})`$ (12) and for $`\stackrel{}{A}_Q(\stackrel{}{r}_i,\stackrel{}{r}_i^{})`$ (15) and also the expression (19), one can show that at $`\tau 0`$ the phase of the Green function $`\gamma _G`$ tends to zero as $`\tau ^{1/2}`$. This property of $`\gamma _G`$ together with the formula (24) allows CDMC successfully simulate all mean quantities.

The simulation process starts with mean field approximation of Fetter, Hanna and Laughlin . In this approximation the mean field $`\overline{\stackrel{}{A}}`$ generated by the average density $`\rho `$ of the particles:

$$\overline{\stackrel{}{A}}(\stackrel{}{r})=\rho \pi \mathrm{}\nu (\stackrel{}{z}\times \stackrel{}{r})=\frac{1}{2}\stackrel{}{B}\times \stackrel{}{r}.$$
(25)
Here $`\stackrel{}{B}=2\pi \rho \mathrm{}\nu \stackrel{}{z}`$ is fictitious uniform mean magnetic field that defines corresponding magnetic length $`a_H=\sqrt{\mathrm{}/B}`$ and cyclotron frequency $`\omega _c=B/M`$. Particles move in this mean magnetic field. Energy spectrum of Hamiltonian (1) with $`\overline{\stackrel{}{A}}`$ and without $`2D`$ parabolic well is Landau levels. We have to simulate anyonic ground state. Therefore we take a wave function $`\mathrm{\Psi }_T^{}(\stackrel{}{r}_i)`$ in the following form:

$$\mathrm{\Psi }_T^{}(\stackrel{}{r}_i)=C\mathrm{exp}\left(\alpha \frac{(x_i^2+y_i^2)}{2R_o^2}\right)\mathrm{exp}\left(\frac{(x_i^2+y_i^2)}{4a_H^2}\right).$$
(26)
Here C is a normalization constant, $`\alpha `$ is numerical parameter and $`R_o=(\mathrm{}/Mw_o)^{1/2}`$. In this expression first exponential function with $`\alpha =1`$ corresponds to the ground state wave function of particle in $`2D`$ parabolic well and second one to the ground state of Landau levels.

Average density of particles $`\rho `$ is given by $`\rho =1/\pi r_o^2`$ where $`r_o`$ is the mean distance between particles in system. We assume that $`r_o`$ equals to $`R_o`$. The reason of this is that the interaction between two anyons at any non-zero $`\nu `$ has repulsive character and the area $`S_n`$ that occupied by free particle at $`n`$-th $`2D`$ oscillator’s level is $`S_n=n\pi R_o^2`$. Last relation for $`S_n`$ follows by using virial theorem for the quantum energetic states for the free particle in $`2D`$ parabolic well.

By substituting expression $`B`$ into expression $`a_H`$ and taking account that $`\rho =1/\pi r_o^2`$, we find the following relation between $`a_H`$ and $`r_o`$: $`a_H=\sqrt{1/2\nu }r_o`$. By similar way one can show that $`\mathrm{}\omega _c/2=\nu \mathrm{}\omega _o`$.

It is natural to use $`\mathrm{}\omega _c/2`$ as energy unit, $`a_H`$ as length unit and $`2/\mathrm{}\omega _c`$ as time unit $`\tau `$. In these units the last term of the Hamiltonian (1) has a coefficient $`1/4\nu ^2`$ and the argument of the first exponential function in (26) has a coefficient $`\alpha /4\nu `$. By using above relations between $`\mathrm{}\omega _c/2`$ and $`\mathrm{}\omega _0`$ and between $`a_H`$ and $`r_0`$ it is easily to express the energies in terms of $`2D`$ oscillator quanta $`\mathrm{}\omega _0`$ and the lengths in terms of $`r_0`$.

## III Discussion of the results

The simulation of the ground state of the systems of several fermions in $`2D`$ parabolic well was performed. By setting the coupling $`\nu =1`$ in (2) we get fermions from anyons. This value of $`\nu `$ is fixed below. As in , the total simulation time $`t`$ is divided into several (in our case we took 40) time blocks $`\mathrm{\Delta }t`$. In every time block $`\mathrm{\Delta }t`$ initial number of configurations $`N_c`$ is chosen equal to 1000 and number of time steps $`\tau `$ equal to 30.

We simulate the energy $`E_{SIM}`$, mean radius $`\overline{r}`$ and mean square of radius of $`\overline{r^2}`$ of the fermionic systems. The calculation of these quantities are performed by the formula (24). Here we take into account that the energy $`E_L(\stackrel{}{R})`$ has as a real part $`\mathrm{Re}E_L(\stackrel{}{R})`$, and also an imaginary part $`\mathrm{Im}E_L(\stackrel{}{R})`$. The number of fermions is varied from two to ten.

The $`\tau `$ dependence of the energy $`E_{SIM}`$ has a general structure $`A_0+A_1\tau ^{1/2}`$  where $`A_0`$ and $`A_1`$ are the numerical constants.

We have obtained that at a very small $`\tau `$ the mean deviation of $`E_{SIM}`$ larger than the fluctuations of $`E_{SIM}`$. At large $`\tau `$ further increasing of it leads to the increasing of the fluctuation of $`E_{SIM}`$ and thus to the growing of the population number $`N_p`$ (the population number $`N_p`$ is connected with first exponential function in (19) and was also introduced in ). The very large values of $`\tau `$ give a dramatical increasing of $`N_p`$ that leads to the instability of $`E_{SIM}`$ and whole stochastic process. So, we have to choose optimal time step $`\tau _0`$.

CDMC  has an artifact that the averaged real quantities become complex, but real parts of them are essentially bigger than imaginary parts because the imaginary parts are controlled by the phase of the Green function $`\gamma _G`$ depending as $`\tau _0^{0.5}`$. Then the imaginary part of $`E_{SIM}`$ for almost all number of fermions $`N`$ is less than 0.15 percent of its real part. For $`\overline{r}`$ and $`\overline{r^2}`$ they are less than 0.07 percent of their real parts.

Table 1 present the exact analytical result for the ground state energy of the fermionic systems in the $`2D`$ parabolic well $`E_{EXACT}`$, $`E_{SIM}`$ (in $`\mathrm{}\omega _o`$ units) and the simulated values of $`\overline{r}`$ and $`\overline{r^2}`$ ( in $`r_o`$ and $`r_o^2`$ units respectively). The values of $`E_{EXACT}`$ was calculated by occupying the available lowest states. The energy of $`K`$-th quantum number state of the $`2D`$ oscillator is $`E_K=\mathrm{}\omega _o(K+1)`$ and degeneracy number of this state is equal $`(K+1)`$ .

$`E_{SIM}`$ is depended on value of numerical parameter $`\alpha `$ in the trial wave function $`\mathrm{\Psi }_T^{}(\stackrel{}{r}_i)`$ (26). We take $`\alpha `$ providing minimal value of the energy $`E_{SIM}`$. The numerical data of $`\alpha `$ that give a minimal energy $`E_{SIM}`$ are outlined also in the Table 1. We have compared the CDMC dependence of the parameter $`\alpha `$ versus number of particles $`N`$ with the variational calculation of $`\alpha `$ . If at the changing of $`N`$ from two to ten CDMC $`\alpha `$ changes from -0.593 to -0.084, so variational $`\alpha `$ changes from -0.293 to -0.683 (see ). For the big number fermions (as one can see from table) CDMC $`\alpha `$ is almost the same for all systems. Thus a mean field approximation of Fetter, Hanna and Laghlin  with numerical parameter in a wave function can be considered as a good start approximation for the simulation of a big number anyons. Table 1 also shows the numerical values of the population number $`N_p`$ and optimal time steps $`\tau _0`$. Good agreement of the simulation results for the $`E_{SIM}`$ with the exact ground state energies $`E_{EXACT}`$ for the numbers of the particles bigger than four allows us to hope that CDMC  gives a possibility for the exact simulation of the big number $`2D`$ fermions.


text:
# Exact solution of a class of one-dimensional nonequilibrium stochastic models*footnote **footnote * to appear in Physical Review E

## I Introduction

Stochastic reaction-diffusion models play an important role in the description of classical interacting many-particle nonequilibrium systems in physics and in interdisciplinary areas (chemistry, biology, economics, …). Usually physical systems are much too complex to be treatable analytically or even numerically. However, in the context of critical phenomena simple toy models have been shown to be useful in determining universal properties and understanding possible relationships between microscopically different processes (see e.g.  and references therein).

In this work we study, on a periodic chain, models which are prototypes of one-dimensional diffusion-limited reactions. Our main purpose is to present a novel approach (based on a duality transformation  combined with generating function techniques ) to analyze the DPAC and related models. In this paper we illustrate this approach by recovering known results and deriving, in a simple and systematic manner, some new ones and show that the method employed here is complementary to the traditional (see e.g.  and references therein) ones. We also solve the dynamics of the dual of the DPAC model which is a biased generalization of Glauber’s model. We believe that the method used here is particularly suited to site-depend and/or disordered  systems. A more detailed presentation of the method, as well as other applications, will be given elsewhere .

Diffusion-annihilation and coagulation models (in their free-fermionic version) properly describe the kinetics of excitons in several materials: the dynamics of photoexcited solitons and polarons in, e.g., chains of $`[Pt(en)_2][Pt(en)_2Cl_2](BF_4)_4`$, called $`(MX)`$, where $`(en)`$ denotes ethylenedyamine, or the fusion of photogenerated excitons in chains of tetramethylammonium manganese trichloride (TMMC) . Models of pair-annihilation and creation are useful to decribe problems of dimer adsorption and desorption .

The motivation for studying such one-dimensional systems is not only their experimental relevance, but also their theoretical importance in the understanding of the fluctuations of low-dimensional systems. In both one- and two- dimensional systems, diffusive mixing is inefficient and leads to the building of large-scale correlations. The mean-field is not adequate and in this sense exact results are desirable.

The paper is organized as follows: In the sequel of this section we recall the stochastic Hamiltonian formulation of Markovian processes obeying a master equation. In section II we introduce two models: the diffusion-limited pair-annihilation-creation (DPAC) model and the diffusion-limited pair-annihilation-creation model connected with a source (DPACI model). In section III we map through a duality transformation the DPAC and DPACI models on other stochastic models. The section IV and V are devoted to the detailed study of the DPAC model. In section IV, we evaluate the exact generating function of the dual DPAC model. Correlation functions of the DPAC model are studied in section V. In section VI, we study the DPACI model via the duality transformation. In section VII we show how to extend the solution obtained for the diffusion-limited pair-annihilation (DPA) model to the diffusion-limited (DC) model. The section VIII is devoted to the conclusion.

In this work we study some one-dimensional two-states nonequilibrium systems for which the dynamics takes place on a periodic chain with $`L`$ (even) sites. The dynamics of particles (of a single species) is governed by a Master equation. Each site of the lattice can be either empty or occupied by a particle at most (the hard-core constraint), say, of species $`A`$ . It is known that such nonequilibrium problems can be reformulated as field-theoretical many-body problems. Below we briefly recall the basis of the field-theoretical formulation of the master equation.

The state of the system is represented by the ket $`|P(t)=_{\{n\}}P(\{n\},t)|n`$, where the sum runs over the $`2^L`$ configurations. At site $`i`$ the local state is specified by the ket $`|n_i=(\mathrm{1\hspace{0.33em}0})^T`$ if the site $`i`$ is empty, $`|n_i=(\mathrm{0\hspace{0.33em}1})^T`$ if the site $`i`$ is occupied by a particle of type $`A`$. We define the left vacuum $`\stackrel{~}{\chi }|`$, where $`\stackrel{~}{\chi }|_{\{n\}}n|`$. The master equation can be rewritten formally as an imaginary-time Schrödinger equation: $`\frac{}{t}|P(t)=H|P(t)`$, where $`H`$ is the Stochastic Hamiltonian which governs the dynamics of the system . In general, it is neither hermitian nor normal. Its construction from the master equation is a well established procedure (see e.g.  and references therein). Here we specifically focus on two-state systems so that the stochastic Hamiltonian can be rewritten in term of Pauli’s operator as a (pseudo-)spin chain (see below). Probability conservation (stochasticity of $`H`$)$`yields\stackrel{~}{\chi }|H=0`$.

## II The DPAC and DPACI models

In this section we present two models which we will consider in the sequel of the work. We begin with the diffusion-limited pair-annihilation-creation (DPAC) model.

We consider, without loss of generality, a periodic chain of $`L`$ (even) sites on which an even number of particles with hard-core repulsion move according to a Master equation Because the DPAC model preserves the parity of the number of particles, we should split the dynamics in a even and an odd sector. Here we focus specifically, and without loss of generality, on the even sector. The treatment of the odd sector is similar excepted for the boundary conditions (see e.g. ). The parity of the number of particles does not play an important role for the dynamics of the model. The restriction to the even sector is performed projecting out the odd sector from the initial state, say $`|P(0)`$, with help of the projector $`\frac{1}{2}\left(𝟙+\right)`$, with $`Q=_m\sigma _m^z`$. In the sequel, for simplicity, we will often call $`|P(0)`$ its even-sector projection, i.e. $`|P(0)^{even}\frac{1}{2}(𝟙+)|(\mathrm{𝟘})`$..

In the DPAC model, at each time-step four events can occur:

A particle can jump from site $`r`$ to $`r+1`$ (provided the latter was previously empty) with a probability rate $`h^{}`$.

A particle can jump from site $`r`$ to $`r1`$ (provided the latter was previously empty) with a probability rate $`h`$.

A pair of particles can be created at sites $`r`$ and $`r+1`$ (provided the latter were previously empty) with a probability rate $`ϵ^{}`$.

A pair of particles can be annihilated at sites $`r`$ and $`r+1`$ (provided the latter were previously occupied) with a probability rate $`ϵ`$.

The Master equation associated to these processes can be formulated as an imaginary time Schrödinger equation. Identifying a vacancy with a (pseudo-) spin up and a particle with a (pseudo-) spin down the stochastic Hamiltonian can be written in term of Pauli (pseudo-)spin$`1/2`$ operators, as follows:

$`H^{DPAC}={\displaystyle \underset{r=1}{\overset{L}{}}}H_{r,r+1}^{DPAC}`$ (1)
$`H_{r,r+1}^{DPAC}`$ $`=`$ $`\left({\displaystyle \frac{ϵ^{}ϵ+(hh^{})}{4}}\right)\sigma _r^z+\left({\displaystyle \frac{ϵ^{}ϵ(hh^{})}{4}}\right)\sigma _{r+1}^z+\left({\displaystyle \frac{ϵ^{}+ϵ+h+h^{}}{4}}\right)(1\sigma _r^x\sigma _{r+1}^x)`$ (2)
$`+`$ $`\left({\displaystyle \frac{ϵ^{}+ϵhh^{}}{4}}\right)(\sigma _r^z\sigma _{r+1}^z+\sigma _r^y\sigma _{r+1}^y)+i\left({\displaystyle \frac{ϵ^{}ϵ(hh^{})}{4}}\right)\sigma _r^x\sigma _{r+1}^y+i\left({\displaystyle \frac{ϵ^{}ϵ+(hh^{})}{4}}\right)\sigma _r^y\sigma _{r+1}^x,`$ (3)
where $`\sigma _j^x,\sigma _j^y,\sigma _j^z`$ are the usual Pauli’s matrices acting on site $`j`$ and $`2\sigma _j^\pm \sigma _j^x\pm i\sigma _j^y`$. The operator $`n_j\frac{1}{2}\left(𝟙\sigma _𝕛^𝕫\right)`$ is the local (at site $`j`$) density operator ($`𝟙`$ denotes the identity operator). This model (with a particle-hole transformation) was introduced by Grynberg et al. , and corresponds to a biased generalization of Lushnikov’s model

The constraint

$`ϵ+ϵ^{}=h+h^{},`$ (4)
is the free-femionic condition, for which the model becomes soluble Recently we have proposed  an approach to study analytically the DPAC model beyond the free-fermion case (4). Although the approach  provides the correct long-time behaviour of the density and correlation fuctions (including the subdominant terms), it turns out to be less general than expected . It is not the purpose of this work to study the limits of the approach devised in .. In the sequel, we assume that the condition (4) is fulfilled. When, in (1), no pairs are created, i.e. $`ϵ^{}=0`$, the model under consideration is called $`DPA`$ (diffusion-limited pair-annihilation) model.

Let us now consider the following stochastic Hamiltonian (with periodic boundary conditions):

$`H^{input}=\gamma {\displaystyle \underset{r=1}{\overset{L}{}}}\left[(\sigma _r^++\sigma _r^{})(\sigma _L^++\sigma _L^{})1\right]=\gamma {\displaystyle \underset{r=1}{\overset{L}{}}}\left[\sigma _r^x\sigma _L^x1\right]`$ (5)
This stochatic Hamiltonian term corresponds to a single particles “source” which injects in the system, with rate $`\gamma `$, particles at site $`j<L`$ and $`L`$, whenever both sites were previously vacant. If both sites $`j<L`$ and $`L`$ were previously occupied, the “source” annihilates (with rate $`\gamma `$) both particles. When one of the sites $`j<L`$ or $`L`$ is occupied and the other empty, the effect of $`H^{input}`$ (with rate $`\gamma `$) is to fill in the previously vacant site and to empty the previously occupied one. In (5), the term in $`\sigma _L^x`$ is an artefact due to the periodicity of the problem (because of the duality transformation, see next section) and preserves the parity of the system described by the stochastic Hamiltonian $`H^{DPACI}H^{DPAC}+H^{input}`$. The steady-state as well as the density relaxation time of such a system (in its translationally invariant version) have been studied  for the case where there was input of particles (i.e. $`\gamma >0`$), but no pair-creation (i.e. $`ϵ^{}=0`$), for an infinite system (where there is no problem of boundary conditions and one can consider simply $`H^{input}=\gamma _r\left[\sigma _r^x1\right]`$). Here we obtain the exact density of the DPACI model and study the interplay between the pair-creation and the source (see (47)). We also obtain the two-point non-instantaneous correlation functions of the (unbiased) DPACI model in the absence as well in the presence of pair-creation.

## III The duality transformation

Following , we introduce a set of operators which forms a periodic Temperley-Lieb algebra (a quotient of a Hecke algebra). Further, we define an unitary transformation which allows to map the DPAC model on another stochastic model for which we are able to calculate the generating function.

Let us define the following operators in term of Pauli spin-matrices:

$`e_{2j1}{\displaystyle \frac{1}{2}}(1+\sigma _j^x),\mathrm{\hspace{0.33em}1}jL;e_{2j}{\displaystyle \frac{1}{2}}(1+\sigma _j^z\sigma _{j+1}^z),\mathrm{\hspace{0.33em}1}jL1;e_0e_{2L}={\displaystyle \frac{1}{2}}(1+\widehat{C}\sigma _L^z\sigma _1^z),`$ (6)
where $`\widehat{C}_{j=1\mathrm{}L}\sigma _j^x`$.

We define the (unitary) duality transformation $`V`$ by

$`V\mathrm{exp}\left({\displaystyle \frac{i\pi }{4}}{\displaystyle \underset{j=1}{\overset{L}{}}}\sigma _j^y\right)\left\{{\displaystyle \underset{k=1}{\overset{2L1}{}}}[(1+i)e_k1]\right\}`$ (7)
The dual (unitary) transformation $`V`$ is also called domain-wall transformation  because of its physical interpretation. In  this duality transformation has been introduced to study the zero-temperature Ising model with help of its free-fermionic dual counterpart (the DPA model, where $`ϵ^{}=0`$). The transformation (7) connects a stochastic model with periodic boundary conditions to another stochastic models with the same boundary conditions .

Let us define the dual version of the stochastic Hamiltonian (1):

$`\widehat{H}_r^{DPAC}V^1H_{r,r+1}^{DPAC}(ϵ,ϵ^{},h,h^{})V;\widehat{H}^{DPAC}{\displaystyle \underset{r=1}{\overset{L}{}}}\widehat{H}_r^{DPAC}`$ (8)
with periodic boundary conditions: $`\sigma _{L+r}^\mathrm{\#}=\sigma _r^\mathrm{\#},(r<L)`$

From now on, we work with the dual model (8) and obtain information on the original model (1) (with the constraint (4)). Notice that the dual model (8) is still a stochastic one because, by construction, $`\stackrel{~}{\chi }|H^{DPAC}=\stackrel{~}{\chi }|H_r^{DPAC}=\stackrel{~}{\chi }|\widehat{H}^{DPAC}=\stackrel{~}{\chi }|\widehat{H_r}^{DPAC}=0`$, since $`\stackrel{~}{\chi }|V\widehat{\stackrel{~}{\chi }}|=\alpha \stackrel{~}{\chi }|`$, where $`\alpha =\frac{(1)^{L1}}{\sqrt{2}}e^{i\frac{\pi }{4}(L1)}`$  is a constant which plays no role in the following and therefore will be omitted.

Here we consider the domain-wall duality transformation (7) and, formally, obtain (when (4) is fulfilled) a generating functions for the dual model. The latter will allow to solve the DPAC model completely. In the third part of this work, we translate our results into the language of the diffusion-limited coagulation model.

The duality transformation (7) maps the DPAC model (1) onto a biased generalization of the Glauber-Ising model . In the absence of bias, model (1) reduces to Lushnikov’s model and the the duality transformation (7) maps this model onto Glauber-Ising model .

Let us also consider the dual transformation of $`H^{input}`$ (5):

$`\widehat{H}^{input}`$ $``$ $`{\displaystyle \underset{r=1}{\overset{L}{}}}V^1H_r^{input}V=\gamma {\displaystyle \underset{r=1}{\overset{L}{}}}\left(𝟙\sigma _𝕣^𝕩\sigma _{𝕣+\mathrm{𝟙}}^𝕩\mathrm{}\sigma _{𝕃\mathrm{𝟙}}^𝕩\right),`$ (9)
with periodic boundary conditions, and thus $`\widehat{H}^{DPACI}=\widehat{H}^{DPAC}+\widehat{H}^{input}`$.

## IV The generating function of (the dual of) the DPAC model

In the first part of the paper we study exact properties of DPAC model using generating function of the dual model (8), which are explicitly computed in this section. We recover known results and produce new ones.

We introduce Grassmann numbers $`\eta _m`$ with their usual anticommuting properties: $`\{\eta _m,\eta _n\}=0,m,n`$

Following , we consider the quantity:

$`G^\pm (\{\eta ,j\},t)\widehat{\stackrel{~}{\chi }}|{\displaystyle \underset{j1}{}}\left(1\pm \eta _j\sigma _j^z\right)e^{\widehat{H}^{DPAC}t}|\widehat{P}(0)={\displaystyle \underset{j1}{}}\left(1\pm \eta _j\sigma _j^z\right)(t),`$ (10)
where $`|P(0)`$ symbolically denotes the initial state in the original (1) model with an even number of particles (see footnote 1), and $`|\widehat{P}(0)V^1|P(0)`$. $`G^\pm (\{\eta ,j\},t)`$ is the generating function of the dual model and its derivatives provide correlation function of the dual model (8), e.g $`\sigma _i^z\sigma _j^z(t)=\frac{^2G^\pm (\{\eta ,j\},t)}{\eta _j\eta _i}|_{\{\eta \}=0},(i<j)`$.

It has to be emphasized that at this stage all the correlation functions obtained from the generating function $`G^\pm (\{\eta ,j\},t)`$ are correlators of the dual model (8). For the sequel, we introduce the following notations:

$`b{\displaystyle \frac{h+h^{}+ϵ+ϵ^{}}{2}};c{\displaystyle \frac{hh^{}+ϵϵ^{}}{2}};d{\displaystyle \frac{h^{}h+ϵϵ^{}}{2}};\stackrel{~}{c}`$ $``$ $`c+d=ϵϵ^{};\stackrel{~}{d}cd=hh^{}`$ (11)
It is useful to separate the generating function into two parts. The one that generates the correlators with an even numbers operators $`\sigma ^z`$ denoted $`𝒱^+(\{\eta ,j\},t)`$. The functional, called $`𝒱^{}(\{\eta ,j\},t)`$, generates correlators with an odd number of operators $`\sigma ^z`$, i.e.,

$`𝒱^\pm (\{\eta ,j\},t){\displaystyle \frac{1}{2}}\left(G^+(\{\eta ,j\},t)\pm G^{}(\{\eta ,j\},t)\right)`$ (12)
When the free-fermionic constraint (4) is fulfilled, the equation of motion of the generating function $`𝒱^\pm `$ can be rewritten, using the properties of Grassmann algebra, as a first-order partial differential equation which can be solved by the method of characteristics :

$`𝒱^\pm (\{\eta ,j\},t)={\displaystyle \frac{1}{2}}\left({\displaystyle \underset{j1}{}}\left(1+\eta _j^0\sigma _j^z(t=0)\right)\pm {\displaystyle \underset{j1}{}}\left(1\eta _j^0\sigma _j^z(t=0)\right)\right)\mathrm{exp}\left[{\displaystyle \underset{k_1}{}}{\displaystyle \underset{k_1>k_2}{}}\eta _{k_1}\eta _{k_2}K_{k_1,k_2}^\pm (t)\right],`$ (13)
where $`\eta _k^0e^{bt}_{j=1}^L(F^\pm (t))_{kj}^1\eta _j`$ and (in the thermodynamics limit, $`L\mathrm{}`$),

$`(F^\pm )_{j,k}^1(t)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}𝑑\varphi \mathrm{exp}\left(i(jk)\varphi +\left\{2c\mathrm{cos}\varphi \stackrel{~}{d}e^{i\varphi }\right\}t\right)=\delta ^{\frac{jk}{2}}I_{jk}(\mathrm{}t)`$ (14)
$`K_{j_1,j_2}^\pm (t)=z^{j_2j_1}e^{2bt}{\displaystyle \underset{k_2>k_1}{}}z^{k_2k_1}\left(I_{j_1k_1}(\stackrel{~}{c}t)I_{j_2k_2}(\stackrel{~}{c}t)I_{j_1k_2}(\stackrel{~}{c}t)I_{j_2k_1}(\stackrel{~}{c}t)\right),(j_2>j_1),`$ (15)
where the $`I_\nu (w)`$’s are the usual modified Bessel function of first kind, $`z\frac{b\sqrt{b^2\stackrel{~}{c}^2}}{\stackrel{~}{c}}`$, $`\mathrm{}\sqrt{\stackrel{~}{c}^2\stackrel{~}{d}^2}`$ and $`\delta ^1\left(\frac{\stackrel{~}{c}\stackrel{~}{d}}{\stackrel{~}{c}+\stackrel{~}{d}}\right)`$.

According to the definition of the generating functions (10),(12) the instantaneous correlation functions of the dual model (8) are obtained by taking partial derivatives of $`𝒱^\pm (\{\eta ,j\},t)`$ (13). If one considers the mean-value of an observable $`O`$ of the original DPAC process (1), its dual counterpart is denoted by $`\widehat{O}(t)`$ as: $`\widehat{O}(t)\widehat{\stackrel{~}{\chi }}|\widehat{O}e^{\widehat{H}^{DPAC}t}|\widehat{P}(0)\stackrel{~}{\chi }|V(V^1OV)(V^1e^{H^{DPAC}t}V)V^1|P(0)O(t)`$, and therefore, the dual counterpart of the observable $`O`$ (say $`n_j\frac{1}{2}(𝟙\sigma _𝕛^𝕫)`$) is $`\widehat{O}V^1OV`$ ($`\widehat{n}_j\frac{1}{2}(𝟙\widehat{\sigma }_{𝕛\mathrm{𝟙}}^𝕫\widehat{\sigma }_𝕛^𝕫)`$).

For the unbiased DPAC model (i.e. Lushnikov’s model, where $`h=h^{}`$), we thus recover Glauber’s original result . For the biased DPAC situation, the generating function (13) not only provides the complete solution of the DPAC model but also the complete solution of the biased generalization of Glauber’s model with transition probability $`\omega (\sigma _i\sigma _i)=D\left(1\frac{d}{2D}\sigma _i(\sigma _{i1}+\frac{c}{d}\sigma _{i+1})\right)`$, where the $`\sigma _i=\pm 1`$ are the usual spin variables and we assume the rate $`D=b/2`$. For this model we obtain:$`\sigma _j(t)=\frac{𝒱^{}(\{\eta ,j\})}{\eta _j}|_{\{\eta \}=0}=e^{bt}_k\sigma _k(0)\delta ^{\frac{kj}{2}}I_{kj}(\mathrm{}t)`$. Relabelling the sites which now run from $`\frac{L}{2}+1`$ to $`\frac{L}{2}`$, the long-time behaviour follows as $`\sigma _j(t)\sigma (0)e^{(b\stackrel{~}{c})t}\frac{e^{(b\mathrm{})t}}{\delta ^j\sqrt{\mathrm{}t}}f(\delta )`$, where $`\delta ^1(\frac{\stackrel{~}{c}\stackrel{~}{d}}{\stackrel{~}{c}+\stackrel{~}{d}})`$, $`\mathrm{}\sqrt{\stackrel{~}{c}^2\stackrel{~}{d}^2}`$ and $`f(u)_ke^{iuk}(\sigma _k(0)\sigma (0))`$, which is considered to be an analytic function. This result generalizes recent results  obtained for the biased zero-temperature Glauber’s model.

## V  Density and Correlation functions of the DPAC model

### A  Density and Correlation functions for arbitrary initial states

In the previous section, we have obtained an explicit expression for the generating function of the dual model (see (7)) and we have shown how to compute correlation functions for the dual model (8). In this section we show how to relate correlation functions of the dual model (8) to the correlation function of the original DPAC model (1).

Here we are especially interested in density-density correlation functions: $`n_{j_1}\mathrm{}n_{j_n}(t)=\frac{1}{2^n}(1\widehat{\sigma }_{j_11}^z\widehat{\sigma }_{j_1}^z)(1\widehat{\sigma }_{j_21}^z\widehat{\sigma }_{j_2}^z)\mathrm{}(1\widehat{\sigma }_{j_n1}^z\widehat{\sigma }_{j_n}^z)(t),(j_1<j_2<\mathrm{}<j_n)`$, where the symbol $`\widehat{\sigma }^z`$ means that the mean-value $`\widehat{\sigma }^z(t)`$ is taken with respect to the dual model (8).

In what follows, we need to know how to connect initial correlation functions of the dual model (8) with initial those of the original DPAC model (1): $`\widehat{\sigma }_{j_1}^z\mathrm{}\widehat{\sigma }_{j_{2n}}^z(0)=_{j_1<jj_{2n}}(12n_j(0)),(j_1<j_2<\mathrm{}<j_{2n})`$

In particular, for the density, we have

$`n_r(t)={\displaystyle \frac{1}{2}}\left(1{\displaystyle \frac{^2𝒱^+}{\eta _{r1}_r}}\right)|_{\{\eta \}=0}={\displaystyle \frac{1}{2}}\left[1z+e^{2bt}{\displaystyle \underset{n>0}{}}z^n\left\{I_{n1}(2\stackrel{~}{c}t)I_{n+1}(2\stackrel{~}{c}t)\right\}\right]`$ (16)
$``$ $`{\displaystyle \frac{e^{2bt}}{2}}{\displaystyle \underset{j_2>j_1}{}}{\displaystyle \underset{j_1<jj_2}{}}(12n_j(0))\delta ^{\frac{j_1+j_22r+1}{2}}\left\{I_{j_1+1r}(\mathrm{}t)I_{j_2r}(\mathrm{}t)I_{j_1r}(\mathrm{}t)I_{j_2+1r}(\mathrm{}t)\right\},`$ (17)
where $`\mathrm{}\sqrt{\stackrel{~}{c}^2\stackrel{~}{d}^2}`$, $`\delta ^1\left(\frac{\stackrel{~}{c}\stackrel{~}{d}}{\stackrel{~}{c}+\stackrel{~}{d}}\right)`$ and we have relabelled the sites indices which now run from $`\frac{L}{2}+1`$ to $`\frac{L}{2}`$.

We can also obtain the non-instantaneous two-point correlation functions :

$`n_r(t)n_s(0)={\displaystyle \frac{1}{2}}\left[1z+e^{2bt}{\displaystyle \underset{n>0}{}}z^n\left\{I_{n1}(2\stackrel{~}{c}t)I_{n+1}(2\stackrel{~}{c}t)\right\}\right]`$ (18)
$``$ $`{\displaystyle \frac{e^{2bt}}{2}}{\displaystyle \underset{j_2>j_1}{}}n_s(0){\displaystyle \underset{j_1<jj_2}{}}(12n_j(0))\delta ^{\frac{j_1+j_22r+1}{2}}\left\{I_{j_1+1r}(\mathrm{}t)I_{j_2r}(\mathrm{}t)I_{j_1r}(\mathrm{}t)I_{j_2+1r}(\mathrm{}t)\right\}`$ (19)
To obtain the long-time behaviour of (18), we introduce the Fourier transform of the initial state from which we have substracted the homogeneous part:

$`f(u,v,s)_{j_1,j_2}e^{i(uj_1+vj_2)}\left(n_s(0)_{j_1<jj_2}(12n_j(0))(12\rho (0))^{j_2j_1}\right)\mathrm{\Theta }(j_2j_1)`$, where $`\mathrm{\Theta }(j_2j_1)=1`$ if $`j_2>j_1`$ and vanish otherwise. With the steepest-descendent method, we obtain for $`|r|,|s|L\mathrm{}`$ and $`bt,|\stackrel{~}{c}|t,\mathrm{}t1`$:

$`n_r(t)n_s(0)\rho (t){\displaystyle \frac{e^{2(b\mathrm{})t}}{\delta ^r(\mathrm{}t)^2}}f(\delta ,\delta ,s),`$ (20)
where $`\rho (t)`$ is the translationally invariant and uncorrelated density studied in (21-33). From (20), we see that the decay of the connected correlation functions (18) depends on the bias and is in the form of an exponential, except for the unbiased case (i.e. $`\stackrel{~}{d}=0`$) and when there is no pair-creation (i.e. $`b=\stackrel{~}{c}`$). In the latter case the decay follows the power-law $`(\stackrel{~}{c}t)^2`$. For the case $`\stackrel{~}{d}<0`$, because of the drift, the rightmost sites tend more rapidly to their steady-state. The effect of the initial state appears through the function $`f`$, which is assumed to be analytic.

All the the multipoint correlation functions can be obtained in a similar and systematic way.

### B Density for translationally invariant and uncorrelated initial states

In this subsection we check the results for the density (16) in the translationally invariant DPAC model against previous results obtained for some particular initial states (lattice with initial density $`\rho (0)=0,1/2,1`$), directly by the free-fermionic approach . As original results, we obtain the exact asymptotic behaviour of the density for uncorrelated and translationally invariant states of arbitrary initial density (in the presence as well as in the absence of pair-creation).

As one can check by direct computation, for tranlationally invariant systems, the density in the DPAC model is independent of the bias $`\stackrel{~}{d}`$ and thus the density of Lushnikov’s model coincides with the density of the generalized DPAC model (for the class of translationally invariant systems under consideration here). It has been shown to be a general property of the DPAC model, which holds true for every instantaneous correlation functions . We have seen in (18,20) that this property is lost for non-instantaneous correlation functions.

For translationally invariant systems and with help of the properties of Bessel functions, (16) leads to $`\widehat{\sigma }_r^z\widehat{\sigma }_{r+s}^z(t)=\widehat{\sigma }_0^z\widehat{\sigma }_s^z(t)=z^s+e^{2bt}_{s^{}>0}\left(\widehat{\sigma }_0^z\widehat{\sigma }_s^{}^z(0)z^s^{}\right)\left(I_{ss^{}}(2\stackrel{~}{c}t)I_{s+s^{}}(2\stackrel{~}{c}t)\right)`$, ($`s>0`$). Where, $`z\frac{b\sqrt{b^2\stackrel{~}{c}^2}}{\stackrel{~}{c}}`$. This result coincides with Glauber’s original one .

Let us consider the uncorrelated case where initially $`n_{i_1}\mathrm{}n_{i_n}(0)=(\rho (0))^n`$, ($`i_n>i_{n1}>\mathrm{}>i_1`$). This implies: $`\widehat{\sigma }_s^{}^z\widehat{\sigma }_0^z(0)=(12n_1)\mathrm{}(12n_s^{})(0)=(12\rho (0))^s^{}`$ and therefore,

$`\rho (t)={\displaystyle \frac{1}{2}}\left(1\widehat{\sigma }_0^z\widehat{\sigma }_1^z(t)\right)`$ $`=`$ $`{\displaystyle \frac{(1z)}{2}}{\displaystyle \frac{e^{2bt}}{2}}{\displaystyle \underset{s^{}>0}{}}\left((12\rho (0))^s^{}z^s^{}\right)\left(I_{s^{}1}(2\stackrel{~}{c}t)I_{s^{}+1}(2\stackrel{~}{c}t)\right)`$ (21)
To study expression (21) we notice that it can be rewritten as

$`\rho (t)\rho (\mathrm{})={\displaystyle \frac{e^{2bt}sign\stackrel{~}{c}}{2}}\left\{\left[(|z|y)I_0(2|\stackrel{~}{c}|t)+(z^2y^2)I_1(2|\stackrel{~}{c}|t)\right]+{\displaystyle \underset{n2}{}}\left((1y^2)y^{n1}(1z^2)|z|^{n1}\right)I_n(2|\stackrel{~}{c}|t)\right\}`$ (22)
where $`y(12\rho (0))sign\stackrel{~}{c}`$.

For the case where $`\stackrel{~}{c}=0`$ with $`ϵ=ϵ^{}>0`$ the density (22) reads $`\rho (t)=\frac{1}{2}+(\rho (0)\frac{1}{2})e^{4ϵt}`$. From now we focus on the cases where $`|\stackrel{~}{c}|0`$.

In the absence of pair-creation (i.e. $`ϵ^{}=0`$, $`b=\stackrel{~}{c}=ϵ>0`$), we simply have (imposing $`z=1`$, in (22))

$`\rho (t)=\rho (0)e^{2\stackrel{~}{c}t}\left[I_0(2\stackrel{~}{c}t)+{\displaystyle \underset{n1}{}}(12\rho (0))^{n1}I_n(2\stackrel{~}{c}t)\right],`$ (23)
this result coincides with the result obtained in .

Taking into account the asymptotic behavior of the Bessel functions $`I_n(x)`$ and collecting terms, we arrive at the asymptotic behaviour of the density.

In the massive (non-universal regime) when $`bt>|c|t1,u=L^2/|c|t1`$ and $`|y|,z<10<\rho (0)<1`$, we have:

$`\rho (t)\rho (\mathrm{}){\displaystyle \frac{e^{2(b|\stackrel{~}{c}|)t}sign\stackrel{~}{c}}{32|\stackrel{~}{c}|t\sqrt{\pi |\stackrel{~}{c}|t}}}\left[(|z|y)(13(y+|z|))+{\displaystyle \frac{|z|(|z|+1)(158|z|3z^2)}{(1|z|)^2}}{\displaystyle \frac{y(y+1)(158y3y^2)}{(1y)^2}}\right]`$ (24)
The validity of (24) is restricted to $`|y|<1`$, i.e. to $`0<\rho (0)<1`$, which corresponds to the convergence domain of the geometrical series occuring in computing (24). The cases $`\rho (0)=0`$ and $`\rho (0)=1`$ correspond to $`y=1`$ (for $`\rho (0)=0)`$ and $`y=1`$ (for $`\rho (0)=1)`$. As $`|z|<1`$, the residual summation over the expansion of $`z^{n1}I_n(2|\stackrel{~}{c}|t)`$ can be carried out, leading to

$`\rho (t)\rho (\mathrm{})\{\begin{array}{cccc}\frac{e^{2(b|\stackrel{~}{c}|)t}}{\sqrt{\pi |\stackrel{~}{c}|t}}\hfill & \text{, if }\rho (0)=0\text{ and }\stackrel{~}{c}>0\text{ }\hfill & & \\ \frac{e^{2(b|\stackrel{~}{c}|)t}}{32(\pi |\stackrel{~}{c}|t)^{3/2}}\left(13(|z|1)+\frac{|z|(1+|z|)}{1|z|}(158|z|3z^2)\right)\hfill & \text{, if }\rho (0)=0\text{ and }\stackrel{~}{c}<0\text{ }\hfill & & \\ \frac{e^{2(b|\stackrel{~}{c}|)t}}{32(\pi |\stackrel{~}{c}|t)^{3/2}}\left(13(|z|1)+\frac{|z|(1+|z|)}{1|z|}(158|z|3z^2)\right)\hfill & \text{, if }\rho (0)=1\text{ and }\stackrel{~}{c}>0\text{ }\hfill & & \\ \frac{e^{2(b|\stackrel{~}{c}|)t}}{\sqrt{\pi |\stackrel{~}{c}|t}}\hfill & \text{, if }\rho (0)=1\text{ and }\stackrel{~}{c}<0\text{ }\hfill & & \end{array}`$ (29)
Results (29) coincide with results obtained in .

In the critical case, when there is no pair-creation ($`ϵ^{}=0,b=\stackrel{~}{c}>0,z=1`$), we have for $`bt=\stackrel{~}{c}t1,u=L^2/|\stackrel{~}{c}|t1`$ and $`|y|<1\rho (0)<1`$.

$`\rho (t)`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{\pi \stackrel{~}{c}t}}}\left(1{\displaystyle \frac{1}{16\stackrel{~}{c}t}}\left[(1y)(3y+2)y(1+y){\displaystyle \frac{158y3y^2}{(1y)^2}}\right]\right)+𝒪\left((\stackrel{~}{c}t)^{5/2}\right)`$ (30)
This result is restricted to $`0<\rho (0)<1`$. When initially there are no particles, $`\rho (0)=0,y=1`$, no dynamics take place. The case of an initially full lattice ($`y=1,\rho (0)=1`$) yields: $`\rho (t)_{|_{\rho (0)=1}}=e^{2\stackrel{~}{c}t}I_0(2\stackrel{~}{c}t)=\frac{1}{2\sqrt{\pi \stackrel{~}{c}t}}(1+𝒪((\stackrel{~}{c}t)^1))`$.

A similar asymptotic results would have been obtained for the case where there are only pairs created and no annihilation ($`b=\stackrel{~}{c}=ϵ^{}>0`$ and $`ϵ=0`$). In this case $`z=1`$ and $`\rho (\mathrm{})=1`$. With (22), we obtain in the asymptotic regime ($`ϵ^{}t1`$, for an initially partially filled lattice, $`0\rho (0)<1`$) a critical decay of the density: $`\rho (t)=1\frac{1}{2\sqrt{\pi bt}}\left(1\frac{1}{16bt}\left[(1y)(3y+2)y(1+y)\frac{158y3y^2}{(1y)^2}\right]\right)+𝒪((bt)^{5/2})`$.

So far we have considered asymptotic behaviour for times which were much larger than the typical times of diffusion (we had $`u=2L^2/|\stackrel{~}{c}|t1`$). Now, we study the asymptotic behaviour of the density, both in the massive and critical regimes, for typical times of order of the diffusion time, i.e. $`u=2L^2/|\stackrel{~}{c}|t1`$.

In the massive, non-universal regime, when $`bt>\stackrel{~}{c}t1,u=2L^2/|\stackrel{~}{c}|t1`$, from (22), we obtain ($`y`$, i.e. $`0\rho (0)1`$):

$`\rho (t)\rho (\mathrm{})`$ $`=`$ $`{\displaystyle \frac{e^{2(b|\stackrel{~}{c}|)t}}{4\sqrt{\pi |\stackrel{~}{c}|t}}}\left[|z|(|z|+1)y(y+1)+{\displaystyle \underset{n2}{}}\left(\{(1y^2)y^{n1}(1z^2)|z|^{n1}\}e^{\frac{n^2}{4|\stackrel{~}{c}|t}}\right)+𝒪\left((|\stackrel{~}{c}|t)^1\right)\right]`$ (31)
In the critical regime when $`ϵ^{}=0`$ and $`bt=|\stackrel{~}{c}|t1`$, with $`u=L^2/|\stackrel{~}{c}|t1`$, from (22), we obtain ($`0<\rho (0)1`$):

$`\rho (t)`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{\pi \stackrel{~}{c}t}}}\left[\left(1{\displaystyle \frac{y(y+1)}{2}}\right)+{\displaystyle \frac{1y^2}{2}}{\displaystyle \underset{n2}{}}y^{n1}e^{\frac{n^2}{4\stackrel{~}{c}t}}\right]+𝒪\left({\displaystyle \frac{1}{(\stackrel{~}{c}t)^{3/2}}}\right)`$ (32)
If one considers $`\frac{n^2}{4\pi |\stackrel{~}{c}|t}<\frac{L^2}{4\pi |\stackrel{~}{c}|t}=u/8\pi 1`$, we recover the universal regime (30). A similar study can be performed for the case where there is no pair-annihilation, but only pair-creation ($`ϵ=0`$, $`b=\stackrel{~}{c}=ϵ^{}>0`$, for a partially filled initial lattice: $`\rho (0)<1`$) and, we have $`\rho (t)1\frac{1}{2\sqrt{\pi \stackrel{~}{c}t}}\left[\left(1\frac{y(y+1)}{2}\right)+\frac{1y^2}{2}_{n2}y^{n1}e^{\frac{n^2}{4\stackrel{~}{c}t}}\right]`$.

Let us conclude this subsection with the study of the so-called staggered current for the (bias-)DPAC model. This quantity has been introduced by Grynberg et al.  to measure the flux of particles due to the bias (so it vanishes for Lushnikov’s model): it is defined as $`𝒥(t)h^{}n_m(1n_{m+1})hn_{m+1}(1n_m)(t)`$. In  this quantity was computed for an initially empty lattice. Here we obtain the exact expression of this quantity for translationally invariant (and uncorrelated, but random) initial states with arbitrary initial density $`\rho (0)`$.

For a translational invariant system, the expression of the staggered current is $`𝒥(t)=\stackrel{~}{d}\left(\rho (t)n_{m+1}n_m(t)\right)=\frac{\stackrel{~}{d}}{4}\left(1\widehat{\sigma }_j^z\widehat{\sigma }_{j+2}^z\right)`$. With help of results of this section, we then find:

$`𝒥(t)={\displaystyle \frac{\stackrel{~}{d}}{4}}\left\{\rho (\mathrm{})(\rho (\mathrm{})1)+e^{2bt}{\displaystyle \underset{s^{}>0}{}}\left[(12\rho (0))^s^{}z^s^{}\right]\left(I_{s^{}2}(2\stackrel{~}{c}t)I_{s^{}+2}(2\stackrel{~}{c}t)\right)\right\}`$ (33)
The long-time behavior of this quantity follows similarly as in (24-32).

### C Density for translationally invariant but correlated initial states

In the previous subsection, we have obtained exact results for the ordered DPAC model for a class of translationally invariant and uncorrelated initial state. It is also instructive to consider translationally invariant but correlated initial states of the form $`(\{\sigma _j^x\sigma _{j+1}^x\})|0`$, which according to (7) are transformed into $`V^1(\{\sigma _j^x\sigma _{j+1}^x\})V|0=(\{\sigma _j^x\})|0`$, where $`(O)`$ is a functional of the operator $`O`$. For translationally invariant systems the density is independent of the bias . In the sequel, we focus on two classes of translationally invariant but correlated initial states of the form $`(\{\sigma _j^x\sigma _{j+1}^x\})|0`$, namely:

i)

$`|P(0){\displaystyle \underset{j}{}}\left({\displaystyle \frac{1+\mu }{2}}+{\displaystyle \frac{1\mu }{2}}\sigma _j^x\sigma _{j+1}^x\right)|0`$ (34)
We have the initial correlations: $`\sigma _j^z(0)=\mu `$, $`\sigma _j^z\sigma _{j^{}j}^z(0)=\mu ^2`$ and so the states (34) correspond to uncorrelated states for the dual model. However, via the duality transformation (7), this state is related to a correlated initial state in the original DPAC model (1) . For the density of the DPAC model, with initial states (34), we have (see (16))

$`\rho (t)={\displaystyle \frac{1z}{2}}+{\displaystyle \frac{e^{2bt}}{2}}{\displaystyle \underset{s>0}{}}z^s\left(I_{s1}(2\stackrel{~}{c}t)I_{s+1}(2\stackrel{~}{c}t)\right){\displaystyle \frac{\mu ^2e^{2bt}}{2}}\left(I_0(2\stackrel{~}{c}t)+I_1(2\stackrel{~}{c}t)\right),`$ (35)
When $`\stackrel{~}{c}=0`$, with $`ϵ=ϵ^{}>0`$, this expression reduces to $`\rho (t)=\frac{1\mu ^2e^{4ϵt}}{2}`$. Hereafter, we focus on the more interesting case $`\stackrel{~}{c}0`$.

In the first term of the r.h.s of (35), we recognize the expression of $`\rho (t)\rho (\mathrm{})`$ for an initially uncorrelated state with $`\rho (0)=1/2`$ (21). When there is no creation of pairs of particles (i.e. $`ϵ^{}=0b=\stackrel{~}{c}=ϵ>0,z=1`$),

$`\rho (t)`$ $`=`$ $`{\displaystyle \frac{e^{2\stackrel{~}{c}t}(1\mu ^2)}{2}}\left[I_0(2\stackrel{~}{c}t)+I_1(2\stackrel{~}{c}t)\right]`$ (36)
In fact, such initial (correlated) states have been considered previously by Santos  who computed the density for the diffusion-annihilation version of the DPAC-model (i.e. $`ϵ^{}=0z=1`$) using Jordan-Wigner transformation and the free-fermionic procedures.

Let us study the asymptotic behaviour of the density for this (correlated) initial state (34). To this end, we proceed as in the previous subsection: we begin with the massive, non-universal regime, where $`bt>|c|t1`$, with $`u=L^2/|c|t1`$. In this regime, the density decays as $`\rho (t)\rho (\mathrm{})=\frac{e^{2(b\stackrel{~}{c})t}}{4\sqrt{\pi |\stackrel{~}{c}|t}}\left(\mu ^2(1+sign\stackrel{~}{c})+𝒪((|\stackrel{~}{c}|t)^1)\right)`$. On the other hand, in the critical regime (when $`ϵ^{}=0,b=\stackrel{~}{c}=ϵ>0,\stackrel{~}{c}t1`$, with $`u=2L^2/|\stackrel{~}{c}|t1`$), we have the following power-law decay: $`\rho (t)\frac{1\mu ^2}{2\sqrt{\pi \stackrel{~}{c}t}}\left(1\frac{1}{8\stackrel{~}{c}t}\right)`$

As noticed in , we see that although the initial state is correlated, the long-time behaviour of the density decays as $`t^{\frac{1}{2}}`$, as in the uncorrelated cases. The interesting point for the DPA model (where $`ϵ^{}=0`$) is however that the dynamics, though critical, is no longer universal: the amplitude of the density (i.e. the term $`1\mu ^2`$) depends on the initial state. We infer that in this case initial correlations don’t renormalize the dynamical exponent, yet affect the amplitude in a non-universal manner.

Similarly, a power-law decay is obtained in the case where there is no pair-annihilation (i.e. $`b=\stackrel{~}{c}=ϵ^{}>0`$): $`\rho (t)1\frac{1}{2\sqrt{\pi bt}}`$. Notice however that the asymptotic behavior is independent of $`\mu `$ and so is universal.

ii) Let us now consider another class of translationally-invariant, but correlated, initial states, namely ($`\mu >0`$):

$`|P(0){\displaystyle \frac{1}{\alpha }}{\displaystyle \underset{j}{}}\left({\displaystyle \frac{1+\mu }{2}}+{\displaystyle \frac{1\mu }{2}}\sigma _j^x\sigma _{j+2}^x\right)V|0`$ (37)
for which initial correlation functions are: $`\sigma _{s0}^z\sigma _0^z(0)=e^\beta \mathrm{sinh}\beta \left(1+\frac{\delta _{s,1}}{2}+(1\delta _{s,0}2\delta _{s,1})e^\beta \mathrm{sinh}\beta \right)`$ and $`\sigma _j^z(0)=e^\beta \mathrm{sinh}\beta `$, where $`\beta \mathrm{ln}\frac{1}{\sqrt{\mu }}`$. The initial states are therefore correlated both for the dual model and for the DPAC model.

Notice that here the operators coding correlations act are not nearest-neighbours preventing a direct Jordan-Wigner transformation of the expression (37), i.e. a direct free-fermion approach.

The computation of the density yields:

$`\rho (t)`$ $`=`$ $`{\displaystyle \frac{1z}{2}}+{\displaystyle \frac{e^{2bt}}{2}}{\displaystyle \underset{s>0}{}}z^s\left(I_{s1}(2\stackrel{~}{c}t)I_{s+1}(2\stackrel{~}{c}t)\right){\displaystyle \frac{e^{2bt\beta }}{2}}\mathrm{sinh}\beta (1+e^\beta \mathrm{sinh}\beta )\left(I_0(2\stackrel{~}{c}t)+I_1(2\stackrel{~}{c}t)\right)`$ (38)
$`+`$ $`{\displaystyle \frac{\mathrm{sinh}\beta e^\beta }{2}}\left(2e^\beta \mathrm{sinh}\beta 1\right)e^{2bt}\left(I_0(2\stackrel{~}{c}t)I_2(2\stackrel{~}{c}t)\right)`$ (39)
In the first term of the r.h.s. of (38), we recognize the expression for the uncorrelated density $`\rho (t)\rho (\mathrm{})`$ with an uncorrelated initial state with $`\rho (0)=1/2`$ (21).

When $`\stackrel{~}{c}=0`$ with $`ϵ=ϵ^{}>0`$, (38) reduces to the following expression: $`\rho (t)=\frac{1}{2}e^\beta \mathrm{sinh}\beta (\frac{3}{2}e^\beta \mathrm{sinh}\beta )e^{4ϵt}`$. Below, we focus on the case where $`\stackrel{~}{c}0`$.

When there are no pairs created (i.e. $`ϵ^{}0,b=\stackrel{~}{c}=ϵ>0`$ and thus $`z=1`$), this expression reduces to:

$`\rho (t)=\left[1e^\beta \mathrm{sinh}\beta (1+e^\beta \mathrm{sinh}\beta )\right]{\displaystyle \frac{e^{2\stackrel{~}{c}t}}{2}}\left(I_0(2\stackrel{~}{c}t)+I_1(2\stackrel{~}{c}t)\right)+{\displaystyle \frac{\mathrm{sinh}\beta e^\beta }{2}}\left(2\mathrm{sinh}\beta e^\beta 1\right)e^{2\stackrel{~}{c}t}\left(I_0(2\stackrel{~}{c}t)I_2(2\stackrel{~}{c}t)\right)`$ (40)
Let us now investigate the asymptotic behaviour of the density for this (correlated) initial state (37). This is done as in the previous subsection. We begin with the massive, non-universal regime, where $`bt>|c|t1`$, and ($`b>\stackrel{~}{c}`$ implies $`|z|<1`$) with $`u=L^2/|c|t1`$. In this regime, the density decays as $`\rho (t)\rho (\mathrm{})=\frac{(1+sign\stackrel{~}{c})e^{2(b|\stackrel{~}{c}|)t}}{4\sqrt{\pi |\stackrel{~}{c}|t}}\left[e^\beta \mathrm{sinh}\beta (1+e^\beta \mathrm{sinh}\beta )+𝒪((|\stackrel{~}{c}|t)^1)\right]`$.

On the other hand, in the critical regime (when $`ϵ^{}=0,b=\stackrel{~}{c}=ϵ>0,ct1`$ and $`z=1`$, with $`u=2L^2/|\stackrel{~}{c}|t1`$), we have a power-law decay: $`\rho (t)=\frac{1e^\beta \mathrm{sinh}\beta \left(1+e^\beta \mathrm{sinh}\beta \right)}{2\sqrt{\pi \stackrel{~}{c}t}}+𝒪((ct)^{3/2})`$

As before, in the case where there would be no pair-annihilation (i.e. $`b=\stackrel{~}{c}=ϵ^{}>0`$), we have a power-law decay of the density: $`\rho (t)1\frac{1}{2\sqrt{\pi bt}}`$. However in this case the asymptotic behavior is independent of $`\mu `$ and so is universal.

Notice that in the critical regime, the density decays as a power-law : $`\rho (t)t^{1/2}`$, as for uncorrelated initial cases (30,32). We also remark that the amplitude of the long-time behaviour of the density is non-universal and depends on the initial state through the parameter $`\mu `$ (with $`\beta =\mathrm{ln}\frac{1}{\sqrt{\mu }}`$).

We conclude that for the correlated states under consideration (34,37), initial correlations affect the dynamics through the amplitude of the density, which in absence of pair-creation is no longer universal. We have observed that when there are pairs created, without pair-annihilation, the density decays, for the two classes of correlated states (34,37), as an universal power-law.

### D Comparison with traditional free-fermion methods

In this subsection we discuss and compare the approach devised in this work for the study of the DPAC model with the free-fermion methods .

In most traditional free-fermionic methods, the stochastic Hamiltonian (1) is recasted, via a Jordan-Wigner transformation, in (free-)fermion form. The resulting quadratic (but generally neither hermitian nor normal) Hamiltonian is then reformulated in the Fourier space. As the resulting Hamiltonian is quadratic, the time dependence of the operators is simple. Such an approach is well adapted for translationally invariant systems, and has been extensively used for the study of the DPA model (where $`ϵ^{}=0`$) with homogeneous initial density. Explicit results have in particular been obtained for initially empty and full lattice as well as for initial density $`\rho (0)=1/2`$ (see e.g. ).

Another approach consists in diagonalizing the (free-)fermion version of the DPAC stochastic Hamiltonian through a Bogoliubov-type rotation, dealing, in Fourier space, with a diagonal quadratic form of the so called pseudo-fermions . This method is particularly efficient for the DPAC model in the presence of pair-creation ($`ϵ^{}>0`$), but concrete results only have been obtained for lattice initially empty or completely filled .

Let us also mention that some exact results have also been obtained for particular non-uniform initial distribution for the DPA and DPAC models .

Here, to illustrate the approach devised in this work, in the translationally invariant situation, we have recovered known results for the density for both DPA and DPAC models (23,29,30, 32,36) and have extended these results to the case of arbitrary initial density (22,24,31,33,35). These results should, in principle, be also accessible by the traditional methods described above, but in the case of arbitrary initial density $`0\rho (0)1`$, the computations required are tedious. On the contrary, with the approach devised here, previous results are recovered in a simple and systematic manner. Although the long-time behaviour of the DPA does not depend on $`\rho (0)`$, this not the case for the general DPAC in the presence of pair-creation, where the dynamics is no longer universal but depends on $`\rho (0)`$ (through the amplitude and the power-law, see (22,24,29,31)). Therefore, when $`ϵ>0,ϵ^{}>0`$, the study of the long-time behaviour of the density with respect of the initial density is relevant. An important advantage of the present method is the fact that it is based on an explicit generating function which allow the systematic computation of multi-point correlation functions. In addition, the method employed in this work provides explicit results which are beyond traditional methods: the results (16,18,20) are specific to this approach and are valid for arbitrary initial distribution (uniform, non-uniform, correlated, …). Furthermore, this approach does not make use of the Jordan-Wigner transformation, and so we are not limited to deal with nearest-neighbour correlations (see (37-40)). Another advantage of the present approach is that it can be extended to the case of site-dependent DPAC model and will also be useful in the study of disordered version of this model . In this sense the present approach is complementary to previous ones.

## VI Density and non-instantaneous correlation functions for DPACI model

In this section, we consider the model described by the stochastic Hamiltonian $`H^{DPACI}=H^{DPAC}+H^{input}`$ (according to (4)) where $`H^{DPAC}`$ and $`H^{input}`$ have been defined previously (in (1) and (5), respectively). From a theoretical point of view this model is interesting because it is one case where (non-trivial) nonequilibrium (the detailed balance is violated) steady-states can be computed exactly.

To proceed with this study, following the original approach of Racz , we first establish the equations of motion of the two-point correlation function of the dual counterpart of these models, i.e. $`\widehat{H}^{DPAC}`$ (defined in (8)) and $`\widehat{H}^{input}`$ (defined in (9)), respectively.

We obtain the equations of motion of the correlation functions as $`\frac{d}{dt}\widehat{\sigma }_m^z\widehat{\sigma }_{m+l}^z=[\widehat{\sigma }_m^z\widehat{\sigma }_{m+l}^z,\widehat{H}^{DPAC}+\widehat{H}^{input}]`$ For the translational invariant situation, we adopt the following notation ($`l>0`$): $`r_l(t)\widehat{\sigma }_m^z\widehat{\sigma }_{m+l}^z(t)=\widehat{\sigma }_0^z\widehat{\sigma }_l^z(t)`$. For the non-translational invariant situation, we define $`r_{m+l,m}(t)\widehat{\sigma }_{m+l}^z\widehat{\sigma }_m^z(t)`$.

For the translational invariant system, the equation of motion of the two-point correlation functions is ($`\stackrel{~}{c}0`$):

$`{\displaystyle \frac{d}{dt}}r_l=2(b+\gamma l\mathrm{\Theta }(Ll))r_l(t)+\stackrel{~}{c}\left(r_{l+1}(t)+r_{l1}(t)\right),`$ (41)
with the boundary condition $`r_0=1`$, where $`\mathrm{\Theta }(Ll)=1`$ for $`l<L`$ and $`\mathrm{\Theta }(Ll)=0`$, for $`l=L`$.

For the non-translationally invariant case, we consider the unbiased version of the DPAC model ((1) with $`h=h^{}>0`$, i.e. Lushnikov’s model), and the equation of motion of the correlation functions then reads ($`l>0`$):

$`{\displaystyle \frac{d}{dt}}r_{m+l,m}=2\left[b+\gamma (2m\mathrm{\Theta }(Lm)+l)\right]r_l(t)+c\left(r_{m+l+1,m}(t)+r_{m+l1,m}(t)+r_{m+l,m}+1(t)+r_{m+l,m1}(t)\right),`$ (42)
with the boundary condition $`r_{0,0}=1`$.

i) We start with the translationally invariant problem (41).

The steady state $`\overline{r}_l`$ corresponding to this situation obeys the difference-equation:

$`2(b+\gamma l\mathrm{\Theta }(Ll))\overline{r_l}=\stackrel{~}{c}\left(\overline{r}_{l+1}+\overline{r}_{l1}\right),`$ (43)
To solve (43), we make the Ansatz $`\overline{r}_l=AJ_{l+\alpha }\left(\stackrel{~}{c}/\gamma \right)`$ ( where $`J_\nu (\omega )`$ is the Bessel function of first kind) and take into account the following property of Bessel functions: $`J_{\nu +1}(\omega )+J_{\nu 1}(\omega )=\frac{2\nu }{\omega }J_\nu (\omega )`$.

Therefore: $`\alpha =\frac{b}{\gamma }L\delta _{l,L}`$. Taking into account the boundary condition $`r_0=\overline{r}(0)=1`$, the constant $`A`$ follows as: $`A=\left(J_{\frac{b}{\gamma }}(\stackrel{~}{c}/\gamma )\right)^1`$. The steady-state of the density for translationally invariant systems described by (43) reads: $`\overline{r}_l=\frac{J_{l+\frac{b}{\gamma }L\delta _{l,L}}(\stackrel{~}{c}/\gamma )}{J_{\frac{b}{\gamma }}(\stackrel{~}{c}/\gamma )}`$. To solve the dynamical equation (41), we seek a solution in the form : $`r_l(t)=\overline{r}_l+_l^{}q_{l^{}l}e^{2\lambda _l^{}\stackrel{~}{c}t}`$

With this Ansatz, (41) reduces to the following difference equation: $`q_{l^{}l+1}+q_{l^{}l1}=\frac{2}{\stackrel{~}{c}}\left[b\stackrel{~}{c}\lambda _l^{}\gamma l\mathrm{\Theta }(Ll)\right]q_{l^{}l},`$ which is solved as above, with an Ansatz: $`q_{l^{}l}=B_l^{}J_{l^{}l+\beta }(\stackrel{~}{c}/\gamma )`$. We thus find: $`\beta =\frac{b\stackrel{~}{c}\lambda _l^{}}{\gamma }+2ll^{}L\delta _{l,L}`$, where the spectrum $`\lambda _l^{}`$ is determined by the boundary condition $`r_0=1`$, which implies $`q_l^{}=0`$, and imposes the condition:

$`J_{\left(b\stackrel{~}{c}\lambda _l^{}\right)/\gamma }\left(\stackrel{~}{c}/\gamma \right)=0`$ (44)
The constants $`B_l^{}`$ are fixed by the initial condition: $`_l^{}q_{l^{}l}=r_l(0)\overline{r}_l=_l^{}B_l^{}J_{\left(b\stackrel{~}{c}\lambda _l^{}\right)/\gamma +lL\delta _{l,L}}\left(\stackrel{~}{c}/\gamma \right)`$

Finally, from the the spectrum of $`\{\lambda _l^{}\}`$ defined by (44) and the $`B_l^{}`$ defined above , the solution of the equation (41) follows

$`r_l(t)=\overline{r}_l+{\displaystyle \underset{l^{}}{}}B_l^{}e^{2\lambda _l^{}\stackrel{~}{c}t}J_{\left(b\stackrel{~}{c}\lambda _l^{}\right)/\gamma +lL\delta _{l,L}}\left(\stackrel{~}{c}/\gamma \right),`$ (45)
where the inverse of the relaxation time, say $`\lambda ^{}`$, is determined from the smallest  zero of the Bessel function (44) and the steady-state has been obtained above.

With the definition of the duality transformation (7), as in (16), we obtain for the density of particles $`\rho (t)`$:

$`\rho (t)={\displaystyle \frac{1r_1(t)}{2}}=\rho (\mathrm{}){\displaystyle \underset{l^{}}{}}{\displaystyle \frac{B_l^{}}{2}}e^{2\lambda _l^{}\stackrel{~}{c}t}J_{\left(b\stackrel{~}{c}\lambda _l^{}\right)/\gamma +1}\left(\stackrel{~}{c}/\gamma \right),`$ (46)
where the steady-state is $`\rho (\mathrm{})=\frac{1\overline{r}_1}{2}=\frac{1}{2}\left(1\frac{J_{\left(1+b/\gamma \right)}\left(\stackrel{~}{c}/\gamma \right)}{J_{\left(b/\gamma \right)}\left(\stackrel{~}{c}/\gamma \right)}\right)`$. In the absence of pair-creation ($`ϵ^{}=0b=\stackrel{~}{c}=ϵ>0`$), we recover the results of  (with $`b=\stackrel{~}{c}=ϵ=1`$).

Following reference , we can now study the relaxation-time in the limit of a weak source (i.e. when $`\gamma 1`$) and analyze its interplay with the pair-creation term (described by ($`b|\stackrel{~}{c}|)/2`$). For the DPAC model, where $`ϵ^{}>0`$ and $`\gamma =0`$, we have seen in section V that the density (for translationally invariant systems) decays exponentially fast in time as $`\mathrm{exp}(2(b|\stackrel{~}{c}|)t)`$. In  it has been shown that in absence of pair-creation (i.e. $`b=\stackrel{~}{c}=ϵ>0`$) the density of the DPAI model decays as $`e^{\gamma ^{2/3}t}`$. Here we will study the case where both pair-creation ($`b>\stackrel{~}{c}>0`$) and the source term ($`\gamma 0`$) are present.

As we focus on the smallest value $`\lambda ^{}`$ of the spectrum $`\{|\lambda _l^{}|\}`$, we consider $`\gamma 1`$. In this limit, we obtain an explicit relation which the $`\{\lambda _l^{}\}`$ have to fulfill in terms of the zero’s $`a_l^{}`$ of Airy function ($`Ai(a_l^{})=0`$), namely : $`\frac{\mu (1+\lambda _j)1}{(1\mu \lambda _j)^{1/3}}=\left(\frac{\gamma ^2}{2b^2}\right)^{1/3}|a_j|`$, where $`\mu \frac{\stackrel{~}{c}}{b}`$ (with $`|\mu |1`$) and the $`a_l^{}`$’s are real and negative.

In the presence of pair-creation (i.e. when $`b>|\stackrel{~}{c}|`$), one has to consider the above equation , which should be solved for the (infinite) set of $`a_j`$’s. In so doing, we would obtain the spectrum $`\{\lambda _j\}`$, which in turn provide the inverse of the relaxation-time $`\lambda ^{}min(\{|\lambda _j|\})`$. Here, we prefer to focus on $`\gamma 1`$, and $`\mu 1`$, i.e., we assume that the pair-creation term is also “small” ($`ϵ^{}1`$). Therefore, the spectrum obeys the following equation: $`\lambda _j=\frac{3\mu }{\mu +2}\left[\left(\frac{\gamma ^2}{2b^2}\right)^{1/3}|a_j|+1\mu \right]`$. In the absence of pair-creation (i.e. $`\mu =1`$), we recover the result of : $`\lambda _j=\left(\frac{\gamma ^2}{2b^2}\right)^{1/3}|a_j|`$.

The inverse of relaxation time $`\lambda ^{}`$, follows ($`\gamma 1,ϵ^{}1`$ and $`\mu 1`$)

$`\lambda ^{}={\displaystyle \frac{3\stackrel{~}{c}}{\stackrel{~}{c}+2b}}\left[\left({\displaystyle \frac{\gamma ^2}{2b^2}}\right)^{1/3}|a_1|+{\displaystyle \frac{b\stackrel{~}{c}}{b}}\right]`$ (47)
Thus the long-time behaviour ($`bt>|\stackrel{~}{c}|t1;\gamma 1ϵ^{}1`$ and $`\mu 1`$) of the density reads (see (47)) : $`\rho (t)\rho (\mathrm{})\frac{B^{}}{2}e^{2\lambda ^{}\stackrel{~}{c}t}`$, where $`B^{}`$ corresponds to the term $`B_l^{}`$ for which $`\lambda _l^{}=\lambda ^{}`$.

ii) Let us now pass to the non-translationally invariant situation. In this case, we have to solve the linear differential-difference equation (42). To do this, we follow the same prodedure as above. The steady-state $`\overline{r}_{m+l,m}`$ is inferred from the Ansatz $`\overline{r}_{m+l,m}=A(J_{m+l+\alpha }\left(c/\gamma \right))^2`$ and obtain: $`\alpha =\frac{b}{2\gamma }\frac{l}{2}L\delta _{m,L}`$. Taking into account the boundary condition $`\overline{r}_{m,m}=1`$ we find ($`l>0`$): $`\overline{r}_{m+l,m}=\left(\frac{J_{m+\frac{l}{2}+\frac{b}{2\gamma }L\delta _{m,L}}(c/\gamma )}{J_{m+\frac{b}{2\gamma }L\delta _{m,L}}(c/\gamma )}\right)^2`$. Hereafter for notational simplicity, we will denote $`\lambda _{m^{},l^{}}`$ simply by $`\lambda `$.

To solve the dynamical equation (42), we seek a solution of the form: $`r_{m+l,m}(t)=\overline{r}_{m+l,m}+_{m^{},l^{}}e^{2\lambda ct}q_{m^{}+l^{}ml,m^{}m}`$ and obtain the solution of (42):

$`r_{m+l,m}(t)=\overline{r}_{m+l,l}+{\displaystyle \underset{m^{},l^{}}{}}A_{m^{},l^{}}e^{2\lambda ct}J_{\frac{bc\lambda }{\gamma }+mL\delta _{m,L}}(c/\gamma )J_{\frac{bc\lambda }{\gamma }+m+lL\delta _{m,L}}(c/\gamma ),`$ (48)
where the boundary condition $`q_{m^{}+l^{}m,m^{}m}=0`$ determines the spectrum $`\{\lambda \}`$ through the condition:

$`J_{\frac{bc\lambda }{\gamma }+mL\delta _{m,L}}(c/\gamma )=0`$ (49)
The constants $`A_{m^{},l^{}}`$ follow from the initial conditions.

With (16, 48), we obtain the density of particles ($`m<L`$):

$`n_{m+1}(t)=n_{m+1}(\mathrm{}){\displaystyle \frac{1}{2}}{\displaystyle \underset{m^{},l^{}}{}}A_{m^{},l^{}}e^{2\lambda ct}J_{\frac{bc\lambda }{\gamma }+m}(c/\gamma )J_{\frac{bc\lambda }{\gamma }+m+1}(c/\gamma ),`$ (50)
with

$`n_{m+1}(\mathrm{})={\displaystyle \frac{1}{2}}\left[1\left({\displaystyle \frac{J_{m+\frac{1}{2}+\frac{b}{2\gamma }}(c/\gamma )}{J_{m+\frac{b}{2\gamma }}(c/\gamma )}}\right)^2\right]`$ (51)
iii) From the solution (50) for the density in the non-translationally invariant situation, we can obtain the non-instantaneous two-point correlation functions of the (unbiased) DPACI model.

With $`p1m<L`$, $`n_{m+1}(t)n_p(0)n_{m+1}(\mathrm{})=_{m^{},l^{}}\frac{B_{m^{},l^{}}}{2}e^{2\lambda ct}J_{\frac{bc\lambda }{\gamma }+m}(c/\gamma )J_{\frac{bc\lambda }{\gamma }+m+1}(c/\gamma )`$

The only change from (50) is in the $`B_{m^{},l^{}}`$ which are fixed from the initial conditions.

The spectrum $`\{\lambda \}`$ again follows from equation (49) and provides the inverse of relaxation time $`\lambda ^{}`$, as $`\lambda ^{}=min\left(\{\lambda _{n,j}\}\right)`$.

With $`nmL\delta _{m,L}`$, in the limit $`\gamma |n|1`$, $`ϵ^{}1`$, with $`\mu 1`$, the spectrum satisfy the equation: $`\lambda _{n,j}\frac{b+2\gamma n}{c}\frac{2\gamma }{c(2\mu 1)^{3/2}}\left[2^{1/3}|a_j|+2\mu n\left(\frac{2\gamma }{b}\right)^{1/3}\right]^{3/2}`$, where the $`a_j<0`$ are again the zeros of the Airy function.

It has to be emphasized that steady-states obtained for this DPACI model are nonequilibrium steady-states, in the sense that they violate the detailed-balance condition. In the DPAC model (where $`\gamma =0`$), the steady-states were expressed in terms of the quantity $`z`$ ( with $`\rho (\mathrm{})=\frac{1z}{2}`$) first introduced by Glauber  to describe the steady-state of the time-dependent Ising model (with the prescription of recovering, as steady-state, the thermal equilibrium). Thus, in this sense the steady-state of the DPAC model is reminiscent of the states in equilibrium statistical mechanics. On the contrary, the steady-states of the DPACI model are genuine nonequilibrium steady-states.

## VII Connection with the coagulation model

In this section we translate the results obtained so far for the DPAC-model to another, free-fermionic, model of stochastic hard-core particles on a lattice, namely the diffusion-limited coagulation (DC) model. In this model particles can jump from a site $`r`$ to an adjacent site (provided it was vacant) with rates $`h`$ ($`\mathrm{}+AA+\mathrm{}`$) and $`h^{}`$ ($`A+\mathrm{}\mathrm{}+A`$). Two adjacent particles at site $`r`$ and $`r1`$ can coagulate with rate $`h`$ ($`A+AA+\mathrm{}`$). Two adjacent particles at site $`r`$ and $`r1`$ can also coagulate with rate $`h^{}`$ ($`A+A\mathrm{}+A`$):

It is known that (for initially uncorrelated systems) there is an exact duality between DPA-model (where there no pair-creation: $`ϵ^{}0`$) and the coagulation model (see e.g. and references therein). Here we briefly recall the transformation connecting these two stochastic models and translate the results of the $`DPA`$-model.

Applying the similarity-transformation  $`\stackrel{~}{H}=H^1;BB\mathrm{}B=B^L`$, where $`B\left(\begin{array}{cc}1& 1\\ 0& 2\end{array}\right),`$ the (stochastic-)Hamiltonian $`H`$ of the DPA-model becomes the (stochastic-)Hamiltonian $`\stackrel{~}{H}`$ of the coagulation-model.

For the coagulation model (with initial state $`|P^{}(0)`$), we can express the correlation functions with help of correlation function of the DPA model. The density of the DC then reads:$`n_j(t)_{|P^{}(0)}^{DC}=2n_j(t)_{^1|P^{}(0)}^{DPA}`$. For the two-point correlations functions, we have:$`n_{j+r}(t)n_j(t)_{|P^{}(0)}^{DC}=4n_{j+r}(t)n_j(t)_{^1|P^{}(0)}^{DPA}`$ and $`n_{j+r}(t)n_j(0)_{|P^{}(0)}^{DC}=2n_{j+r}(t)n_j(0)_{^1|P^{}(0)}^{DPA}`$

As an illustration, let us now provide, the density of the $`DC`$ model from the expressions obtained for the DPA-model (where $`ϵ^{}0`$) as (see (16)). The density of the translationally invariant DC model with uncorrelated initial state $`|P^{}(0)=(1\rho (0)\rho (0))^T`$ (of initial density $`\rho `$), is related to the unbiased version of the DPA model with initial density $`\frac{\rho (0)}{2}`$. Thus, in this situation, with (16) , we obtain for the density of DC model:

$`\rho ^{DC}(t)`$ $`=`$ $`e^{2(h+h^{})t}\left\{I_0(2(h+h^{})t)+I_1(2(h+h^{})t){\displaystyle \underset{s>0}{}}(1\rho (0))^s\left(I_{s1}(2(h+h^{})t)I_{s+1}(2(h+h^{})t)\right)\right\}`$ (52)
As for the DPA-model, for translationally invariant systems, the instantaneous correlation functions do not depend on the bias $`hh^{}`$. We could also extend in a similar manner all the results obtained in the previous sections V to the DC model. The discussion of the long-time behaviour of $`\rho ^{DC}(t)`$ follows directly from that of the DPA model given in subsection V.B.

## VIII Summary and Conclusion

This work provides a complete and unified solution of some related one-dimensional nonequilibrium models.

In a first part we have solved the diffusion-limited pair-annihilation-creation (DPAC) model (under the free-fermionic constraint) by a novel approach, alternative and complementary to the traditional free-fermion methods. Combining domain-wall duality  with a recently developped technique , we were able to compute the generating function of the dual of the DPAC model (section IV), i.e. to obtain correlation functions (section V). In this paper, we take advantage of this approach to recover results previously obtained (see (23,29,30,36)) by traditional methods, extend these (see(22,24,31,32,33,35)) and derive original ones (see(16,18,20,38,40)), in a simple manner. The method is very promising for site-dependent and disordered systems . We focused on the density, the staggered-current, and the two-point correlation functions (instantaneous and non-instantaneous) (see (16,18,20)). In particular we studied explicitly the density (in absence as well as in presence of pair-creation) for translationally invariant systems with arbitrary initial density and have discussed the asymptotic behaviour of these quantities in different regimes. We obtained explicit expression of the staggered-current in arbitrary translationally-invariant (uncorrelated but random) initial conditions (see (33)). We investigated the densities for two classes of correlated (but tranlationally invariant) systems. We have computed explicitly the non-instantaneous two-point correlation functions for initially translationally invariant and uncorrelated (but random) states: these quantities depend on an eventual bias (contrary to instantaneous correlation functions). As a by-product we provided the solution of a biased generalization of Glauber’s model, which is the dual counterpart of the DPAC model.

In a second part we studied the density and non-instantaneous two-point correlation functions of the diffusion-limited pair-annihilation-creation model in the presence of source of particles (DPACI model). Via the domain-wall duality transformation, we solved exactly the equations of motion of the correlation functions of the dual model. We obtained the exact density and steady-states for translationally invariant as well as non translationally-invariant systems. We also computed the non-instantaneous two-point correlation functions. In particular for the translationally invariant situation in the absence of the pair-creation term, we recovered, the density first computed by Racz . We extended these results to the case where both pair-creation term ($`ϵ^{}>0`$) and source term were present ($`\gamma >0`$) and obtained the relaxation-time for the case when both $`\gamma 1`$ and $`ϵ^{}1`$ (see (47)).

In the last part of this work, we have shown how to translate the results obtained for the DPA model to the diffusion-limited coagulation (DC) model.

## IX Acknowledgments

The support of the Swiss National Fonds is gratefully acknowledged.