Yinsheng Li

Yang Chen

Yining Hu

Ahmed Oukili

Limin Luo

Wufan Chen

Christine Toumoulin
email: christine.toumoulin@univ-rennes1.fr

Strategy of CT sinogram inpainting based on sinusoid-like curve decomposition and eigenvector-guided interpolation

Keywords: computed tomography (CT), sinogram completion, detector truncation, sparse sampling

published or not. The documents may come

Strategy of computed tomography sinogram inpainting

based on sinusoid-like curve decomposition and eigenvector-guided interpolation.

Introduction

In To X-ray computed tomography (CT), sparse sampling and detector gaps are two typical cases of sinogram incompleteness [START_REF] Kak | Principles of Computerized Tomographic imaging[END_REF][START_REF] Prince | Constrained sinogram restoration for limited-angle tomography[END_REF][START_REF] Joseph | View sampling requirements in fan beam computed tomography[END_REF][START_REF] Ohnesorge | Efficient correction for CT image artifacts caused by objects extending outside the scan field of view[END_REF][START_REF] Maltz | CT truncation artifact removal using water-equivalent thicknesses derived from truncated projection data[END_REF][START_REF] Xu | Statistical projection completion in X-ray CT using consistency conditions[END_REF][START_REF] Chen | Mathematical models for local non-texture inpaintings[END_REF]. Sparse sampling with limited projections of views has wide applications in fast scanning (~4-10s) such as cardiac CT and neonate CT examination [START_REF] Prince | Constrained sinogram restoration for limited-angle tomography[END_REF][START_REF] Joseph | View sampling requirements in fan beam computed tomography[END_REF]. Detector gaps in sinogram are generally caused by physical damage of CT detectors along some specific projection views [START_REF] Ohnesorge | Efficient correction for CT image artifacts caused by objects extending outside the scan field of view[END_REF][START_REF] Maltz | CT truncation artifact removal using water-equivalent thicknesses derived from truncated projection data[END_REF][START_REF] Xu | Statistical projection completion in X-ray CT using consistency conditions[END_REF]. However, sinogram data missing from sparse sampling and detector gaps often lead to degraded CT images with severe streak or ring artifacts. suppress these artifacts, many approaches have been proposed in the last decade. They are generally classified as inpainting-like and statistical methods. Typical inpainting-like techniques include TV (total variation) inpainting [START_REF] Chen | Mathematical models for local non-texture inpaintings[END_REF][START_REF] Xue | Metal artifact reduction in dual energy CT by sinogram segmentation based on active contour model and TV inpainting[END_REF], partial differential equation (PDE) based inpainting [START_REF] Gu | X-ray CT metal artifacts reduction through curvature based sinogram inpainting[END_REF][START_REF] Kostler | Adaptive variational sinogram interpolation of sparsely sampled CT data[END_REF] and directional interpolation [START_REF] Bertram | Directional view interpolation for compensation of sparse angular sampling in cone-beam CT[END_REF]. In some of these methods such as in [START_REF] Constantino | Sinogram recovery for sparse angle tomography using sinusoidal Hough transform[END_REF][START_REF] Zamyatin | Extension of the reconstruction field of view and truncation correction using sinogram decomposition[END_REF][START_REF] Chityala | Artifact reduction in truncated CT using sinogram completion[END_REF], the inherent property of sinusoid-like curve (S-curve) composition is exploited to restore the incomplete sinograms. The second category of methods, which has been higly developped in the literature, makes use of statistical methods with prior constraints to suppress these artifacts. Some approaches such as maximum likelihood expectation maximization (MLEM) [START_REF] Oehler | Statistical image reconstruction for inconsistent CT projection data[END_REF][START_REF] Zbijewski | Statistical reconstruction for x-ray CT systems with non-continuous detectors[END_REF], Bayesian based EM optimization [START_REF] Lemmens | Suppression of metal artifacts in CT using a reconstruction procedure that combines MAP and projection completion[END_REF] and TV-based compressed sensing (CS) reconstructions [START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained total-variation minimization[END_REF][START_REF] Duan | Few-view projection reconstruction with an iterative reconstruction re-projection algorithm a TV constraint[END_REF][START_REF] Ritschl | A new approach to limited angle tomography using the compressed sensing framework[END_REF][START_REF] Tang | Performance Comparison Between Total Variation (TV)-Based Compressed Sensing and Statistical Iterative Reconstruction Algorithms[END_REF], were found to be efficient in improving the reconstructed images in the presence of sparse sampling and detector gaps. It was noted in [START_REF] Zbijewski | Statistical reconstruction for x-ray CT systems with non-continuous detectors[END_REF] that statistical reconstruction allows suppression of the artifacts brought on by the asymmetric detector gaps, but remains inefficient to remove the artifacts brought on by the symmetric detector gaps. In [START_REF] Zamyatin | Extension of the reconstruction field of view and truncation correction using sinogram decomposition[END_REF][START_REF] Chityala | Artifact reduction in truncated CT using sinogram completion[END_REF], the intensities of the unmeasured points were set to the minimum intensities of the known intensities along the decomposed S-curves. This method was shown to be effective in the case of external truncation because the high attenuation tissues are always located within the internal field of view [START_REF] Chityala | Artifact reduction in truncated CT using sinogram completion[END_REF]. Nevertheless, this method might not work in the cases of detector gaps and sparse sampling, in which the missing sinogram data correspond to tissues or structures for which there is a high attenuation. In

Method: am decomposition

the 2D fan-beam CT geometry as illustrated in Fig. 1: , O Q isocenter of FOV and X-ray source this paper, with an aim to provide effective sinogram completion for sparse sampling and detector gaps, we proposed a new strategy of inpainting, which includes three steps:S-curve decomposition, sinusoid fitting and eigenvector-guided interpolation, respectively. The first step is performed to determine the family of S-curves that goes through each unmeasured sinogram point [START_REF] Constantino | Sinogram recovery for sparse angle tomography using sinusoidal Hough transform[END_REF][START_REF] Zamyatin | Extension of the reconstruction field of view and truncation correction using sinogram decomposition[END_REF][START_REF] Chityala | Artifact reduction in truncated CT using sinogram completion[END_REF]. In the second step, a sinusoid fitting method is applied to retrieve the S-curves that go through each sinogram point with unmeasured intensity. Finally, to estimate the unmeasured sinogram intensities, an eigenvector-guided interpolation is realized to determine the S-curve to be used for interpolation. This process is described in section 2. In section 3, sparse sampling and symmetric detector gaps under 2D fan-beam CT geometry are then simulated to validate the effectiveness of the proposed method and Compute Unified Device Architecture (CUDA) based parallelization is performed to accelerate the sinusoid fitting and interpolation operations [START_REF] Kharlamov | Image Denoising Tech. Rep[END_REF][START_REF][END_REF][START_REF]Accelerating MATLAB with CUDA Using MEX Files[END_REF][START_REF]GPU Acceleration in MATLAB[END_REF]. Section 4 provides a qualitative and quantitative evaluation of the method as well as a comparison with two other inpainting methods and a CS iterative reconstruction algorithm. Results show that better performances can be reached with our approach in the cases of sparse sampling and detector gaps. Section 5, finally,concludes on the relevance of the method.

A. Geometry and sinogr

Let consider the notations listed in Table 1 and If R de te source-to-isocenter distance, the circul trajectory of Q can be given by: no s the ar ur ( ) ( cos , sin

) T R R β β β = y (1)
The projection of each X-ray can be jointly defined by the project ion view β and the fan angle γ . T he fan-b eam transformation is given by: 0 ( , ) ( ( ) )

β γ μ β θ +∞ = + ∫ ur r g t

dt y

(2) where, θ r denotes the unit ray directional vector and is calculated as [cos(

),sin( )]

β γ β γ + + . [0, ] V β ∈ and [ , 1] D D γ ∈ -
-, with V and D denoting the total projection nu ec sion, Eq. ( 2), t deno s the integral path along each ray, which leads to the sinogram intensity at point ( , ) mber and the half det tor dimen respectively. In te β γ . From [START_REF] Zamyatin | Extension of the reconstruction field of view and truncation correction using sinogram decomposition[END_REF], we know that each point in the reconstructed image uniquely corresponds to one specifi rve in sinogram space, and the sinogram data c S-cu g to the superposition of all the S-curves.

2 2 2 2 2 2 | ' | | | cos( ) cos( ) | ' | cos( ) | ' | | ' | | ' | sin( ) | | ' ' | | ' ' | | ' ' OX OX r QX R r QX R OX XX r QX QX XX QX QX XX QX QX XX ϕ ϕ ϕ β ϕ β ϕ β ϕ β ϕ ⎫ ⎫ = -= - = - - ⎪ ⎪ ⎪ ⎪ = - ⇒ = - ⇒ ⎬ ⎬ ⎪ ⎪ = + = + ⎪ ⎪ ⎭ ⎭ = + =
uuuuur uuur uuuur uuuur uuuur uuuur uuur uuur uuur

2 2 2 2 ( sin( )) ( cos( )) 2 cos( ) r R r r R R r ϕ ϕ ϕ ϕ β ϕ β ϕ - + - - = + - - β ϕ (6)
rom Eq. ( 6), we can obtain the relation between

and β , r ϕ , R γ F : 2 2 2 2 sin( | '| r XX ϕ β ) sin( ) sin arcsin | | 2 cos( ) 2 cos( ) r QX r R Rr r R Rr ϕ ϕ ϕ ϕ ϕ ϕ β ϕ γ γ β ϕ β ϕ - = = ⇒ = + - - + - - uuur (7)
ere, each sinogram point

- uuuur ( , ) β γ H
is from the path of the X-ray beam MN , and the point X moves from M to N as ϕ varies from 2 ϕ to 1 ϕ . Then from the set of all the po in ints MN , we can termine a set S-curves ( , ) S de of β γ going thro h po t ( , ) ug in β γ in the sinogram space (as illust d in Fig. 3):

rate 2 2 sin( ) ( , ) ( , ) | arcsin 2 cos( ) r S r R Rr ϕ ϕ ϕ β ϕ β γ β γ γ β ϕ ⎧ ⎫ - ⎪ ⎪ ⊆ = ⎨ ⎬ + - - ⎪ ⎪ ⎩ ⎭ ) ) ) ) ) (8)
where, [0, 2 ]

β π ∈ ) and [ ] 1 2 , ϕ ϕ ϕ ∈
.

B. The method:

S-point ( , )

β γ lies within a group of S-curves ( , )

S β γ

Each unmeasured , so the interpolation of each S-curve going through the poin ) t ( , β γ will lead to one estimation. Theref all the S-curves going through the point ( , ) ore β γ should be jointl idered to estimate the intensity at ( , ) y cons β γ . We thus developed thus an eigenvecto ed interpolation to estimate the intensities of the S-points erform this interpolation, we r-guid . To p set a patch (as the red rectangle in Fig. 3) P centered at the current unmeasured point. Then, taking into account the local sinogram structure, the tensor product matrix T for patch P is calculated in Eq.( 9) as the convolution of a kernel K(•) and the tensor production of the partial differentia on matrix of P : ti [START_REF] Gu | X-ray CT metal artifacts reduction through curvature based sinogram inpainting[END_REF] where, ( ) ) uuuuuu r e will also not damage the region consistency because the eigenvalues of T tend to be equal in this ca

2 2 11 1 T T T 2 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) T x x y x y y x y ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ∂ ⎛ ⎞ ⎛ ⎞ ∂ ∂ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ = ∇ ∇ = = = ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ∂ ∂ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ∂ ⎝ ⎠ ⎝ ⎠ ⎝
( ) ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ( ) ( ) x y ∂ ∂ ∂ ∂ ∂ ∂ P x P x x y ∂ ∂ ∂ ∂ ( ) / x ∂ ∂ P x directi
We applied the Jacobi rotation method [27 o calculate the eigenvalues of th s ] t e symmetric ock matrix [START_REF] Bertram | Directional view interpolation for compensation of sparse angular sampling in cone-beam CT[END_REF] bl T to and identify the eigenvector associated with the smallest eigenvalue. The least-square method is used find the interpolated S-curve with the closest direction to the vector S ( , β γ uuuuuuuuu r e . In this least-square method, we define each feature line in the sinogram coordinate as: In summary, we outline the steps in the proposed approach as follows:

Step1: Compute a family of S-curv

Experiment Settings

periments were conducted on fan-beam clinical data to test the proposed approach. The 2D dataset was acquired on a Siemens Somatom scanner. Shown in Fig. 6 (a), the original complete sinogram data is from one pelvis scanning of a 56-year-old male patient. The scanning protocol is 100 mAs, 120kVp, 6 mm slice-thickness, 256 radial samples and 360 angular samples. The source circle radius and the detector-center distance are 816 mm and 408mm. The target image size is 512 512 × . The acquired sinogram was normalized, corrected for attenuation and scatter and then rebinned to parallel format. We applied a Filtered Back Projection (FBP) reconstruction algorithm with a Ramp kernel to build the image (Fig. 6 (d)) from the original complete sinogram (Fig. 6 (a)). As illustrated in Fig. 6 (b) and (c), sparse sampling and detector gaps are respectively simulated. In the case of sparse view sampling, projection data contains 90 views on the 360. Fig. 6 (c) illustrates the simulated detector gaps, in which two 10 detector gaps are symmetrically located on each side of the center-of-detector. Fig. 6 (e) and (f) show the reconstructed images from these two sinograms. We can see in Fig. 6 that the data missing from sparse sampling and detector gaps significantly degrades the reconstruction by introducing severe cyclic and streak artifacts in the reconstructed images. We can note that the data missing from sparse sampling and detector gaps leads to degraded reconstructions with severe cyclic and streak artifacts.

We then compare then our approach with two other inpainting methods which are linear interpolation and TV inpainting. For all the inpainted sinograms, the FBP algorithm is applied with a Ramp filter to obtain the reconstructed images.

Linear interpolation estimates the unmeasured sinogram data using a linear weighting that is inversely proportional to the spatial distance between the known data and unmeasured data. Based on [START_REF] Kothe | Edge and junction detection with an improved structure tensor[END_REF], 1-D linear interpolation along both the projection view and the detector ordinate are applied to interpolate the missing data for sparse sampling and detector gaps, respectively.

TV inpainting estimates the unknown data by computing the numerical solution to the partial differential equation (PDE) [START_REF] Chen | Mathematical models for local non-texture inpaintings[END_REF]. The PDE equation includes a positive shifting parameter to avoid zero denominator, and a time parameter to determine the maximum iteration to inpaint all the missing points.

We also compare the reconstruction results from the proposed inpainting with the TV-based CS statistical algorithm in [START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained total-variation minimization[END_REF][START_REF] Duan | Few-view projection reconstruction with an iterative reconstruction re-projection algorithm a TV constraint[END_REF][START_REF] Ritschl | A new approach to limited angle tomography using the compressed sensing framework[END_REF][START_REF] Tang | Performance Comparison Between Total Variation (TV)-Based Compressed Sensing and Statistical Iterative Reconstruction Algorithms[END_REF]. In this method, the artifact suppressed CT image μ is obtained by minimizing the energy function that combines the data consistency term in sinogram space with the TV term in image space:

2 l 1 l ( ) ( ) 1 ˆarg max( ( )) 2 T TV g A D g A U β μ μ μ = - - + μ (11) ( ) ( ) 2 2 2 ( ) TV h v j j j U ε μ μ μ = Δ + Δ + ∑ (12)
where, g is the measured projection sinogram data, is the system matrix and A D is the diagonal covariance matrix. Matrixes and A D can be determined well from the fan-beam CT geometry. h j μ Δ and v j μ Δ denote the linear operators corresponding to the horizontal and vertical first order differences at pixel j, respectively. ε is a small positive number introduced to avoid a zero denominator. Parameter β controls the trade-off between the data consistency and the TV terms.

In the experiments, the parameters involved in the different methods are set to minimize the mean square error (MSE) on the sinogram. The parameter β in TV-based CS reconstruction is set using the rule of MSE minimization on the reconstructed images. The MSE value for the inpainted sinogram ( g MSE ) and the reconstructed images ( MSE ) are respectively given by Eqs. ( 13) and ( 14):

μ 2 1 ( o g j j g ) j MSE g g N = - ∑ (13) 2 1 ( o j j j MSE N μ μ ) μ μ = - ∑ (14)
where, g N and N μ denote the total element numbers in the sinograms and the reconstructed images, respectively. g and o g denote the inpainted sinograms and the original complete sonogram respectively, In our approach, to limit the computational costs, we choose equidistantly discretized S-curves in the estimation of each of the unmeasured points. Patch size P has to be suitably set to characterize local sinogram textures. It has to be larger than the original gaps to guarantee the original known intensity information is used for inpainting. Nevertheless a too large Patch may include many non-local distributed structures, and might cause a failure in characterizing the local intensities. Additionally, a large P increases the computationnal costs of the eigenvalues and eigenvectors. Therefore, the size definition of P should consequently take into account the practical distribution of the missing data in sinogram. We also found the proposed approach was not sensitive to the size of P . No visual difference could be detected if we changed the patch size from 7x7 to 15x15 for the case of sparse sampling.

A GPU implementation using the Compute Unified Device Architecture (CUDA) framework was applied to accelerate the process of fitting and interpolation of all the S-curves on each unmeasured points. Routinely, the CUDA computation model is broken down into a host (CPU) and the device (GPU), both of which are connected by PCI-E bus. Threads, blocks and grids make up of the physical structure of a CUDA multithreading frame [START_REF] Kharlamov | Image Denoising Tech. Rep[END_REF][START_REF][END_REF][START_REF]Accelerating MATLAB with CUDA Using MEX Files[END_REF][START_REF]GPU Acceleration in MATLAB[END_REF]. Once the data enter device memory and are ready to be processed, kernel functions are called to fit and interpolate the S-curves simultaneously. The algorithms have been written in C language using the MATLAB release R2006b and NVIDIA CUDA libraries and were run on a PC with an Inter Xeon processor, GHZ, 8G RAM and GPU (NVIDIA GTX465).

2 ×

Results

Experiments on a sparse sampling sinogram

Parameter setting for the different methods are listed in Table . 2. Fig. 7 (a)-(c) show the inpainted sinograms from linear interpolation, TV inpainting and the proposed inpainting method. The delineated ROIs in Fig. 7 (a)-(c) are zoomed in Fig. 7 (a1)-(c1). Fig. 7 (d1) provides the corresponding zoomed ROI of the original reference sinogram of Fig. 6 (a). These enlarged regions show that our proposed method can lead to a better restoration of the sinogram textures. Fig. 8 (a)-(c) display the reconstructed images from the inpainted sinograms of Fig. 7 (a)-(c). We can still note that the proposed inpainting provide a better reconstruction with more effective artifact-suppression and structure-preservation than with the other methods. The iterative TV-based CS reconstruction in Fig. 8 (d) also demonstrates a good artifact-suppression property. Comparing with the TV-based CS reconstruction, the proposed approach provides a reconstruction with comparable artifact-suppression, and a better preservation of fine structures. Fig. 9 plots the profiles along the 245 th rows in the reconstructed images of Fig. 8. This graph confirms the performance of our approach by showing better agreement with the reference image (Fig. 6 (d)) than is found with the other methods. The MSE computed for both the inpainted sinograms and the reconstructed images are given in Table . 3 for each method. The figures confirms the supremacy of our approach with a MSE that is the lowest in the set of all methods.

Method Parameter settings

Linear interpolation [START_REF] Kothe | Edge and junction detection with an improved structure tensor[END_REF] 5 points in each symmetric side are used in the interpolation TV inpainting [START_REF] Chen | Mathematical models for local non-texture inpaintings[END_REF] Inpainting time: 300; Positive shifting parameter : Table 3. MSE of the inpainted sinograms and reconstructed images of Fig. 7 and Fig. 8 respectively

Method

Parameter settings

Linear interpolation [START_REF] Kothe | Edge and junction detection with an improved structure tensor[END_REF] 5 points in each symmetric side are used in interpolation TV inpainting [START_REF] Chen | Mathematical models for local non-texture inpaintings[END_REF] Inpainting times: 300; Positive shifting parameter :

Computation costs:

The computation cost (in CPU seconds) for different methods are listed in Table 6. As to the methods of linear interpolation, TV inpainting, and the proposed method, the total computation time includes both the inpainting and FBP reconstructions. For TV-based CS reconstruction, the computation time relates just to the reconstruction process. We can see that more computation costs are required for the TV-based CS method and our proposed method. Moreover, when considering the sparse sampling and detector gap, this computational costs increased with the number of unmeasured data to be estimated.

Conclusion

In Ho this paper, a new sinogram inpainting strategy is proposed to suppress the artifacts brought by sparse sampling and detector gaps in CT. This strategy includes S-curve based sinogram discretization, sinusoid fitting and eigenvector-guided interpolation. A CUDA-based parallelization is applied to accelerate the processing. Experiments on real fan-beam scanning data show that the proposed completion strategy can lead to improved CT imaging in the two cases of sparse view sampling and symmetric detector gaps.

wever, from the above inpainting and reconstruction results from Fig. 7-Fig. 12, we can still observe some error estimation in the inpainted sinograms from the proposed approach. One reason for this might be that the interpolation must rely on the remaining known data along each decomposed S-curve to complete the S-curves with missing data. Therefore in the cases with too much missing data (very low sampling ratios or very large detector gaps), performance of the method would be lowered. Another limit is the computation time that still remains relative high even when using CUDA acceleration. Further work will be devoted to the estimation of more accurate intensity information, the extension of the method to 3D cone-beam CT scan, and reducing the computation time.

  source rotates around the object along the circular trajectory with the angle view e 1. Parameter def nition for the 2D fan-beam CT geometr OV of region within the Q β . ( ) X μ denotes the intens y at position it X in the target 2D attenuation distribution.

Fig. 1 .Fig. 2 .

 12 Fig.1. Fan beam geometry and sinogram decomposition

Fig. 3 .

 3 Fig.3. A family of S-curves going though the point

e

  to identity matrix because the sinograms are assumed to be free of noise. We compute the eigenvalu and eigenvector of T . Then for each S-point ( , ) β γ , according to the tensor theory in[START_REF] Kothe | Edge and junction detection with an improved structure tensor[END_REF] and as illustrated in Fig.4, the eigen ctor L ve of T represents the direction of the largest intensity va , while the eigenvector S um eigenvalue of T denotes the direction of the smallest intensity variation, this latter being able to used to represent the local edge direction. To preserve the original local sinogram textures, we choose the interpolated S-curve whose direction is the closest to the vector S among all the interpolated S-curves going through ( , )β γ , and we set the intensity at the unmeasured ( , ) point β γ to the interpolated intensity at point ( , ) β γ in hosen S-curve. If the point ( , ) the c β γ belongs to a local ogeneous region other than an edge region, this choice of the S-curve with the b tch to the direction S

Fig. 4 .Fig. 5 .

 45 Fig.4. Eigenvalue analysis of one patch P in the sinogram. β and γ characterize the pr tion the detector ojec view and coordinate, respectively. ) L ( , β γ uuuuuur e and ) S ( , β γ ur uuuuu e correspond to the dire tion of the largest and smallest intensity variation, respectively.c

)

  es ( , S β γ for a unmeasured sinogram po , ) int (β γ according to (8). Eqs. (3)-Step2: For each S-curve in the S-curve family S( , ) β γ , estimate the positions and intensities of the red points by sinusoid fitting and cubic spline interpolation, using all the measured points in each unmeasu S-curve. o eigenval c estimatio Step3: C mpute the tensor product matrix T for a patch P using Eq. (9), and then calculate the ues and eigenvectors of T by Jacobi rotation. Step4: A cording to the direction of eigenvector with the minimum eigenvalue, determine the ultimate n of this current S-point ( , ) timation of unmeasured points is performed sequentially from the upper left to the lower right epeat Step1-Step4 recursively until all the unmeasured sinogram points are estimated.

Fig. 5

 5 Fig. 5 illustrates the process of estimating one specific unmeasured point.

Fig. 6 .

 6 Fig.6. Reconstruction of incomplete sinograms. (a), original complete sinogram; (b), simulated sinogram with a sparse sampling: 90 projection views on 360 are available; (c), simulated sinogram containing 2 symmetric detector gaps: two 10 row detector gaps are located at each side of the center-of-detector; (d), image reconstructed from the complete sinogram in (a); (e), image reconstructed from the sparse sampling sinogram in (c); (f), image reconstructed from the sinogram including detector gaps in (d).We can note that the data missing from sparse sampling and detector gaps leads to degraded reconstructions with severe cyclic and streak artifacts.

  μ and o μ denote the reconstructed images from the inpainted sinogram and the reference reconstructed image (Fig. 6(d)) from the original complete sinogram .g o g

Table 2 .Fig. 7 .

 27 Fig.7. Inpainting results in the experiment of sparse sampling after application of (a), the linear interpolation; (b), the TV inpainting method in [7]; (c), our proposed method; (a1)-(c1): Zoomed ROI from (a)-(c); (d1) is the zoomed region from the original complete reference sinogram in Fig.6 (a).

Fig. 8 .

 8 Fig.8. Reconstruction results in the experiment of sparse sampling. (a)-(c), from the corresponding inpainted sinograms in Fig.7 (a)-(c) and (d) from the TV-based CS method.We can note that the proposed inpainting can lead to reconstruction with better artifact-suppression and structure-preservation than the other inpainting methods, and reconstruction comparable to the TV-based CS reconstruction can be obtained from the inpainted sinogram using the proposed method.

Fig. 10 .

 10 Fig.10. Inpainting results in the experiment of symmetric detector gaps. Sinogram restoration using (a), the linear interpolation; (b), the TV inpainting method in [7]; (c), the proposed method; (a1)-(c1): zoomed ROIs from (a)-(c); (d1) is the zoomed region from the original complete reference sinogram in Fig.6 (a).

Table 5 .

 5 Fig.11. Reconstruction results in the experiment of symmetric detector gaps. (a)-(c), from the inpainted sinograms in Fig.10 (a)-(c); (d), from the TV-based CS reconstruction.We can note that the proposed inpainting can lead to reconstruction (Fig.11 (c)) with better artifact-suppression and structure-preservation than other inpainting methods and the iterative CS reconstruction.

  Fig.12. Illustration of the profiles along the 100 th row in the reconstructed images in Fig.11

  ⎠

					P x		P x P x	P x
	T K P P	K	P x	P x	P x	K	P x P x	P x	1 T

interpolation +FBP TV Inpainting +FBP Proposed inpainting +FBP Unparallelized proposed method +FBP TV-based CS reconstruction

  Table 6 also shows that CUDA parallelization led to nearly 40-fold acceleration.

		Linear
	Sparse sampling	1.02	22.16	452.68	18021.42	14298.47
	Detector gaps	0.98	13.96	750.71	28800.28	14298.47

Table 6 .

 6 CPU computation cost (in CPU seconds) for different methods.

Acknowledgment

This research was supported by National Basic Research Program of China under grant (2010CB732503), National Natural Science Foundation under grant (81000636), and the Project supported by Natural Science Foundation of Jiangsu Province (BK2009012).

Experiment on detector gaps

 [START_REF] Zbijewski | Statistical reconstruction for x-ray CT systems with non-continuous detectors[END_REF], Zbijewski et al. tested the application of the iterative reconstruction algorithm in the case of an incomplete CT sinogram with non-continuous detectors. The authors pointed out that, although the method was able to reconstruct artifact-free images from sinograms containing asymmetric detector gaps, the iterative algorithm was inefficient in suppressing artifacts when sinograms contained symmetric detector gaps. We applied our proposed inpainting to suppress the artifacts in Fig. 6 (f), which were caused by symmetric detector gaps in the sinogram and compared the results with those obtained from the 2 inpainting methods and the TV-based CS reconstruction algorithm. The parameter settings for the different methods are provided in Table. [START_REF] Ohnesorge | Efficient correction for CT image artifacts caused by objects extending outside the scan field of view[END_REF].

Fig. 10 (a)-(c) show the completed sinograms from linear interpolation, TV inpainting and our proposed inpainting method and the zoomed ROIs in Fig. 10 (a)-(c) are displayed in Fig. 10 (a1)-(c1). Fig. 10 (d1) displays the corresponding zoomed region in the original complete sinogram. Keeping the original complete sinogram (Fig. 6 (a)) as the reference, we can observe in Fig. 10 that linear interpolation and TV inpainting fail to restore the original sinogram textures, while our proposed method allows a better restoration of sinogram textures. Fig. 11 (a)-(c) depict the reconstructed images from the corresponding inpainted sonograms in Fig. 10 (a)-(c). We can observe that, with respect to the reference reconstruction in Fig. 6 (d), the proposed inpainting method still give a reconstruction with less artifacts than the other inpainting methods. Fig. 11 (d) shows the images from the TV-based CS reconstruction algorithm. We can note that the TV constraint leads to smooth structures without significantly improved artifact suppression. Compared with the TV-based CS reconstruction, the proposed approach is more efficient (Fig. 11 (c)) in suppressing ring artifacts, and preserving original structures.

.12 plots the profiles along the 100 th row in the reconstructed image of Fig. 11. As previously observed with sparse sampling sinograms, the proposed approach behaves in the same way and provides a reconstruction result that is in better conformity with the reference image than are the results obtained with the other methods. MSE computed on the inpainted sinogram and the reconstructed images (and given in Table 5) also confirms these results.


Ekaterina Vl Bulinskaya

ASYMPTOTIC BEHAVIOR OF LOCAL PARTICLES NUMBERS IN BRANCHING RANDOM WALK

Keywords: walk, Bellman-Harris process with particles of six types, Yaglom type conditional limit theorems, Kolmogorov's equations, 2010 AMS classification: 60F05

Critical catalytic branching random walk on an integer lattice Z d is investigated for all d ∈ N. The branching may occur at the origin only and the start point is arbitrary. The asymptotic behavior, as time grows to infinity, is determined for the mean local particles numbers. The same problem is solved for the probability of particles presence at a fixed lattice point. Moreover, the Yaglom type limit theorem is established for the local number of particles. Our analysis involves construction of an auxiliary Bellman-Harris branching process with six types of particles. The proofs employ the asymptotic properties of the (improper) c.d.f. of hitting times with taboo. The latter notion was recently introduced by the author for a non-branching random walk on Z d .

Keywords and phrases: critical branching random

Introduction

Catalytic branching random walk (CBRW) on d-dimensional integer lattice is a model of particles population evolution. We recall its main features. Each particle independently of others may perform random walk on Z d and produce offsprings at the source of branching located w.l.g. at the origin. Symmetric branching random walk (SBRW) on Z d studied earlier, e.g., in [START_REF] Albeverio | Asymptotics of branching symmetric random walk on the lattice with a single source[END_REF], [START_REF] Bogachev | The moment analysis of a branching random walk on a lattice with a single source[END_REF] and [START_REF] Yarovaya | A limit theorem for critical branching random walk on Z d with a single source[END_REF] is a particular case of CBRW on Z d (see [START_REF] Yarovaya | Criteria of the exponential growth of particles numbers in models of branching random walks[END_REF]).

The model under consideration was proposed in [START_REF] Vatutin | Catalytic branching random walk and queueing systems with random number of independent servers[END_REF] for d = 1 and studied for other d ∈ N in [START_REF] Bulinskaya | Catalytic Branching Random Walk on a Two-Dimensional Lattice[END_REF], [START_REF] Bulinskaya | Catalytic Branching Random Walk on Three-Dimensional Lattice[END_REF], [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF] and [START_REF] Yarovaya | Criteria of the exponential growth of particles numbers in models of branching random walks[END_REF]. The analysis of CBRW in [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF] and [START_REF] Yarovaya | Criteria of the exponential growth of particles numbers in models of branching random walks[END_REF] has shown that similarly to many kinds of branching processes (see [START_REF] Sewastianow | [END_REF]) CBRW on Z d is classified as supercritical, critical or subcritical. According to [START_REF] Yarovaya | Criteria of the exponential growth of particles numbers in models of branching random walks[END_REF], the exponential growth (as time tends to infinity) of total number of particles in population and local numbers of particles as well is characteristic for the supercritical CBRW on Z d . The term local refers to the (number of) particles located at a lattice point.

Quite different situation occurs for critical CBRW which is the main object of study in this paper. For example, for d = 1 or d = 2 the particles population degenerates with probability 1 but survives with strictly positive probability for d ≥ 3 (see [START_REF] Bulinskaya | Catalytic Branching Random Walk on Three-Dimensional Lattice[END_REF], [START_REF] Bulinskaya | Limit Distributions for the Number of Particles in Branching Random Walks[END_REF], [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF] and [START_REF] Vatutin | Catalytic branching random walk and queueing systems with random number of independent servers[END_REF]). Moreover, the total number of particles conditioned on non-degeneracy has non-trivial discrete limit distribution, different for d < 3 and d ≥ 3 (see [START_REF] Bulinskaya | Catalytic Branching Random Walk on a Two-Dimensional Lattice[END_REF], [START_REF] Bulinskaya | Catalytic Branching Random Walk on Three-Dimensional Lattice[END_REF] and [START_REF] Vatutin | Catalytic branching random walk and queueing systems with random number of independent servers[END_REF]). Thus, in the model of critical

Main results

Now we dwell on the definition of a critical CBRW on Z d . At the initial time t = 0 there is a single particle on the lattice located at a point x ∈ Z d . If x = 0, the particle performs a continuous time random walk until the time of the first hitting the origin. The random walk outside the origin is symmetric, homogeneous, irreducible (i.e. a particle passes from an arbitrary u ∈ Z d to any υ ∈ Z d with positive probability within a finite time) and has a finite variance of jumps. Accordingly, we assume this random walk be specified by an infinitesimal matrix A = (a(u, υ)) u, υ∈Z d such that a(u, υ) = a(υ, u), a(u, υ) = a(0, υ -u) := a(υ -u), u, υ ∈ Z d , υ∈Z d a(υ) = 0 where a(0) < 0 and a(υ) ≥ 0 if υ = 0, υ∈Z d υ 2 a(υ) < ∞. If x = 0 or the particle has just hit the origin it spends there an exponentially distributed time (with parameter 1). Afterwards, it either dies with probability α ∈ (0, 1) producing before the death a random number of offsprings ξ or leaves the source of branching with probability 1 -α. In the latter case the intensity of transition from the origin to a point υ = 0 is given by a(0, υ) = -(1 -α) a(υ) a(0) .

At the origin the branching is determined by a probability generating function

f (s) := Es ξ = ∞ k=0 f k s k , s ∈ [0, 1].
In [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF] CBRW on Z d is called critical if the following relations hold

αf ′ (1) + (1 -α)(1 -h d ) = 1 and σ 2 := α f ′′ (1) < ∞. (1)
Here h d is the probability of the event that a particle leaving the origin will never return there. By the recurrence of a random walk on Z and Z 2 one has h 1 = h 2 = 0. It is well known that h d ∈ (0, 1) for d ≥ 3. Newborn particles are located at the origin at the birth moment. They evolve according to the scheme described above independently of each other as well as of the parent particles. The number of particles located at a point y ∈ Z d at time t ≥ 0 is denoted by µ(t; y).

The goal of the paper is three-fold. Firstly, we find the asymptotic behavior (as t → ∞) of the mean number of particles m(t; x, y) := E x µ(t; y) located at a point y ∈ Z d , y = 0, at time t ≥ 0 (everywhere the index x means that our CBRW starts at x ∈ Z d ). Secondly, we retrieve the asymptotic behavior of the probability q(t; x, y) := P x (µ(t; y) > 0) of the presence of particles at the point y at time t. Thirdly, we establish a limit theorem for properly normalized local numbers µ(t; y) conditioned on µ(t; y) > 0 as t → ∞.

To formulate the main results of the paper we introduce some more notation. Let p(t; x, y) be the transition probability from x to y within time t ≥ 0 for a random walk on Z d generated by matrix A. Set

G λ (x, y) := ∞ 0 e -λt p(t; x, y) dt, λ > 0, x, y ∈ Z d .
Note that the Green's function G 0 (x, y) := lim λ→0+ G λ (x, y) is well-defined and takes finite values for d ≥ 3 by virtue of the transience of our random walk on Z d , d ≥ 3. One can check (see [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF]) that h d = (aG 0 (0, 0)) -1 , d ∈ N, where a := -a(0).

As shown in [START_REF] Yarovaya | Branching Random Walks in Inhomogeneous Medium[END_REF], Theorem 2.1.1 (see also [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF]), for any fixed x, y ∈ Z d , one has

p(t; x, y) ∼ γ d t d/2 , p(t; 0, 0) -p(t; x, y) ∼ γd (y -x) t 1+d/2 , t → ∞, (2)
where

γ d := (2π) d |det φ ′′ θθ (0)| -1/2 , φ(θ) := z∈Z d a(z, 0) cos(z, θ), θ ∈ [-π, π] d , φ ′′ θθ (0) = ∂ 2 φ(θ) ∂θ i ∂θ j θ=0 i,j∈{1,...,d} , γd (z) := 1 2(2π) d R d (υ, z) 2 e (φ ′′ θθ (0)υ,υ)/2 dυ, z ∈ Z d ,
and (•, •) stands for the scalar product in R d . In particular, it follows that the value

m d := 1 -(1 -α)a -1 + 2(1 -α)a -1 G -2 0 (0, 0) ∞ 0 tp(t; 0, 0) dt is finite for d ≥ 5. Set also q(s, t; x, y) := 1 -E x s µ(t;y) , s ∈ [0, 1], t ≥ 0, x, y ∈ Z d . For d = 2 we use the function J(s; y) := α ∞ 0 (f (1 -q(s, u; 0, y)) -1 + q(s, u; 0, y)) du, s ∈ [0, 1], y ∈ Z d .
The main results are contained in the following three theorems. For the sake of completeness their statements include the case y = 0 studied earlier in [START_REF] Bulinskaya | Catalytic Branching Random Walk on a Two-Dimensional Lattice[END_REF]- [START_REF] Bulinskaya | Limit Distributions Arising in Branching Random Walks on Integer Lattices[END_REF], [START_REF] Hu | Branching random walk in Z 4 with branching at the origin only[END_REF], [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF], [START_REF] Topchii | Individuals at the origin in the critical catalytic branching random walk[END_REF] and [START_REF] Vatutin | Catalytic branching random walk and queueing systems with random number of independent servers[END_REF].

Theorem 1 Let x, y ∈ Z d . The following relations are valid being different for y = 0 and y = 0, namely, as t → ∞,

m(t; x, y) ∼ γ 1 √ t , m(t; x, 0) ∼ γ 1 a (1 -α) √ t , d = 1, m(t; x, y) ∼ γ 2 t , m(t; x, 0) ∼ γ 2 a (1 -α)t , d = 2, m(t; x, y) ∼ G 0 (x, 0)G 0 (0, y) 2πγ 3 √ t , m(t; x, 0) ∼ aG 0 (x, 0)G 0 (0, 0) 2πγ 3 (1 -α) √ t , d = 3, m(t; x, y) ∼ G 0 (x, 0)G 0 (0, y) γ 4 ln t , m(t; x, 0) ∼ aG 0 (x, 0)G 0 (0, 0) γ 4 (1 -α) ln t , d = 4, m(t; x, y) → (1 -α)G 0 (x, 0)G 0 (0, y) a G 2 0 (0, 0) m d , m(t; x, 0) → G 0 (x, 0) G 0 (0, 0)m d , d ≥ 5.
Theorem 2 For x, y ∈ Z d and t → ∞ the following formulae hold true

q(t; x, y) ∼ 2(1 -α) σ 2 γ 1 a √ t ln t , d = 1, q(t; x, y) ∼ γ 2 t 1 - a 1 -α J(0; y) , y = 0, d = 2, q(t; x, 0) ∼ γ 2 a (1 -α)t (1 -J(0; 0)), d = 2, q(t; x, y) ∼ 4πγ 3 (1 -α)G 0 (x, 0) σ 2 a G 3 0 (0, 0) √ t ln t , d = 3, q(t; x, y) ∼ 3γ 4 (1 -α)G 0 (x, 0) ln t σ 2 a G 3 0 (0, 0) t , d = 4, q(t; x, y) ∼ 2 m d G 0 (x, 0) σ 2 G 0 (0, 0)t , d ≥ 5,
where for d = 2 and s ∈ [0, 1] the strict inequalities J(s; y) < (1 -s)(1 -α)/a, y = 0, and J(s; 0) < 1 -s are valid.

Theorem 3 Given x, y ∈ Z d , λ ∈ [0, ∞) and s ∈ [0, 1], one has, as t → ∞, lim t→∞ E x exp - λ µ(t; y) E x (µ(t; y)|µ(t; y) > 0) µ(t; y) > 0 = 1 λ + 1 , d = 1, d = 3 or d ≥ 5, lim t→∞ E x s µ(t;y) µ(t; y) > 0 = (1 -α)s -a(J(0; y) -J(s; y)) 1 -α -aJ(0; y) , y = 0, d = 2, lim t→∞ E x s µ(t;0) µ(t; 0) > 0 = s -(J(0; 0) -J(s; 0)) 1 -J(0; 0) , d = 2, lim t→∞ E x exp - λ µ(t; y) E x (µ(t; y)|µ(t; y) > 0) µ(t; y) > 0 = 1 3 + 2 3 • 2 2 + 3λ , d = 4.
Observe that the normalizing factor E x (µ(t; y)|µ(t; y) > 0) arising in Theorem 3 is exactly m(t; x, y)/q(t; x, y) and the asymptotic behavior of the functions m(t; x, y) and q(t; x, y) is given by Theorems 1 and 2, respectively.

To establish Theorem 1 it is useful to invoke the forward and backward Kolmogorov's differential equations (considered in appropriate Banach spaces) for mean numbers of particles at different points of the lattice and also the resulting integral equations (see [START_REF] Yarovaya | Branching Random Walks in Inhomogeneous Medium[END_REF]). As for Theorems 2 and 3, note that for proving results in [START_REF] Bulinskaya | Catalytic Branching Random Walk on a Two-Dimensional Lattice[END_REF], [START_REF] Hu | Branching random walk in Z 4 with branching at the origin only[END_REF], [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF], [START_REF] Topchii | Individuals at the origin in the critical catalytic branching random walk[END_REF] and [START_REF] Vatutin | Catalytic branching random walk and queueing systems with random number of independent servers[END_REF] concerning the number of particles at the origin the method of introduction of an auxiliary Bellman-Harris branching process with particles of two types was efficient. However, for proving Theorems 2 and 3 we have to involve a Bellman-Harris branching process with particles of six types. To apply the latter method we attend to a new notion of the hitting time with taboo in the framework of a (non-branching) random walk on Z d . More precisely, we use our recent results (see [START_REF] Bulinskaya | The Hitting Times with Taboo for a Random Walk on an Integer Lattice[END_REF]) on the asymptotic behavior of the tail of the (improper) cumulative distribution function of this time. Due to that one can employ the theorems by V.A.Vatutin for Bellman-Harris branching processes with particles of several types (see, e.g., [START_REF] Vatutin | Discrete Limit Distributions of the Number of Particles in a Multitype Age-Dependent Branching Processes[END_REF]- [START_REF] Vatutin | On a Class of Critical Bellman-Harris Branching Processes with Several Types of Particles[END_REF]). Afterwards we have to deal with sophisticated analytic estimates of the solutions of the parametric integral equations (see, e.g., [START_REF] Hu | Branching random walk in Z 4 with branching at the origin only[END_REF], [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF], [START_REF] Topchii | Individuals at the origin in the critical catalytic branching random walk[END_REF] and [START_REF] Vatutin | Catalytic branching random walk and queueing systems with random number of independent servers[END_REF]).

Proof of Theorem 1

Let us recall some useful results employed within this section. According to [START_REF] Gikhman | The Theory of Stochastic Processes II[END_REF], Ch.3, Sec.2, the transition probabilities p(t; x, y), t ≥ 0, x, y ∈ Z d , of the random walk generated by matrix A satisfy the backward Kolmogorov's equations

d p(t; x, y) d t = (Ap(t; •, y)) (x), p(0; x, y) = δ y (x). (3)
Here (Ap(t; •, y))(x) = z∈Z d a(x, z)p(t; z, y) and δ y (•) is a column vector in the space l 2 (Z d ) with zero components except for the component 1 indexed by y. In a similar way, the backward Kolmogorov's equations for m(t; x, y), t ≥ 0, x, y ∈ Z d , (see, e.g., Theorem 2.1 in [START_REF] Yarovaya | Criteria of the exponential growth of particles numbers in models of branching random walks[END_REF]) take the form d m(t; x, y) d t = Am(t; •, y) (x) + β c (∆ 0 m(t; •, y)) (x), m(0; x, y) = δ y (x),

where A = (a(u, υ)) u,υ∈Z d := A+(a -1 (1 -α) -1) ∆ 0 A, ∆ 0 := δ 0 δ T 0 (T stands for transposition) and β c := (1 -α)a -1 G -1 0 (0, 0). Here we follow the notation of [START_REF] Yarovaya | Criteria of the exponential growth of particles numbers in models of branching random walks[END_REF].

Lemma 1 For each y ∈ Z d , the function m(t; y, y) is non-increasing in t.

Proof. The monotonicity of m(•; y, y) for SBRW on Z d was established in Lemma 3.3.5 of [START_REF] Yarovaya | Branching Random Walks in Inhomogeneous Medium[END_REF]. The key step of its proof was to use self-adjointness of the operator H := A + β c ∆ 0 where β c := G -1 0 (0, 0). For CBRW the analog of H is the non self-adjoint operator H := A + β c ∆ 0 . However, Lemma 3.1 in [START_REF] Yarovaya | Criteria of the exponential growth of particles numbers in models of branching random walks[END_REF] permits to pass to (self-adjoint) symmetrization of H and then apply Lemma 3.3.5 in [START_REF] Yarovaya | Branching Random Walks in Inhomogeneous Medium[END_REF]. Further argument is similar to the proof of Theorem in [START_REF] Yarovaya | The monotonicity of the probability of return into the source in models of branching random walks[END_REF]. Equation (4) was obtained by differentiating at s = 1 the following backward Kolmogorov's equation for the generating function F (s, t; x, y) := E x s µ(t;y) , s ∈ [0, 1], t ≥ 0, x, y ∈ Z d , (see [START_REF] Yarovaya | Criteria of the exponential growth of particles numbers in models of branching random walks[END_REF])

∂F (s, t; x, y) ∂t = AF (s, t; •, y) (x) + ∆ 0 f (F (s, t; •, y)) (x), F (s, 0; x, y) = s δy(x) . ( 5
)
Here

f (s) := α(f (s) -s), s ∈ [0, 1]
, is an infinitesimal generating function of the number of offsprings of a parent particle. We will employ (5) in Section 5.

In Lemma 2 we derive a counterpart of the forward Kolmogorov's equation for the fuction F (s, t; x, y), s ∈ [0, 1], t ≥ 0, x, y ∈ Z d , and, as a consequence, the forward Kolmogorov's equation for m(t; x, y). Recall that A * denotes an adjoint operator for A and

A * m(t; x, •) (y) = z∈Z d m(t; x, z)a(z, y). Lemma 2 For s ∈ [0, 1], t ≥ 0, x, y ∈ Z d , the following relation holds true ∂F (s, t; x, y) ∂t = (s -1)
z∈Z d , z =y a(z, y)E x s µ(t;y) µ(t; z) + (s -1)a(y, y)E x s µ(t;y)-1 µ(t; y)

+ δ 0 (y)f(s)E x s µ(t;y)-1 µ(t; y), F (s, 0; x, y) = s δx(y) . (6)
Moreover, one has

d m(t; x, y) d t = A * m(t; x, •) (y) + β c (∆ 0 m(t; x, •)) (y), m(0; x, y) = δ x (y). ( 7
)
Proof. As usual in derivation of forward Kolmogorov's equations, we consider all possible evolutions of the particles population within the time interval [t, t + h) and let h → 0+. To justify arising passages to the limit we involve the Lebesgue theorem on dominated convergence and useful estimates for transition probabilities (see proof of Lemma 3 in [START_REF] Gikhman | The Theory of Stochastic Processes II[END_REF], Ch.3, Sec.2). We also benefit from finiteness of the mean total number of particles M(t; x) := E x z∈Z d µ(t; z) for each x ∈ Z d and t ≥ 0. The latter observation is true since the last function belonging to l ∞ (Z d ) is a solution of the linear differential equation in ( 4) with the initial condition M(0; x) = 1 for all x (instead of δ y (x) in ( 4)), see [START_REF] Yarovaya | Criteria of the exponential growth of particles numbers in models of branching random walks[END_REF]. Equation ( 7) is an immediate consequence of ( 6) due to formula m(t; x, y) = ∂ s F (s, t; x, y)| s=1 . We also take into account that f ′ (1) = β c in view of [START_REF] Albeverio | Asymptotics of branching symmetric random walk on the lattice with a single source[END_REF].

Consider equations ( 4) and ( 7) as inhomogeneous ones for differential equation (3) in Banach space l ∞ (Z d ). Applying the variation of constant formula (see [START_REF] Daletsky Yu | Stability of Solution of Differential Equations in Banach Space[END_REF], Ch.2, Sec.1) we infer that

m(t; x, y) = p(t; x, y) + 1 - a 1 -α t 0 p(t -u; x, 0)m ′ (u; 0, y) du + a β c 1 -α t 0 p(t -u; x, 0)m(u; 0, y) du, (8)
m(t; x, y) = p(t; x, y) + 1 -α a -1 t 0 m(t -u; x, 0)p ′ (u; 0, y) du + β c t 0 m(t -u; x, 0)p(u; 0, y) du. (9)
An analogous result for SBRW on Z d can be found in [START_REF] Yarovaya | Branching Random Walks in Inhomogeneous Medium[END_REF], Theorem 1.4.1. Now we can give Proof of Theorem 1. To find the asymptotic behavior of m(t; x, y), t → ∞, x, y ∈ Z d , y = 0, we estimate each of the summands in the right-hand sides of ( 8) and ( 9) when x = 0 and x = 0, respectively, as t → ∞. Namely, we will show that, for d = 1 and d = 2, the main contribution to the asymptotic behavior of the right-hand side of ( 8), as well as of ( 9), is due to the first summand. However, for d ≥ 3, the asymptotic behavior of the right-hand sides of ( 8) and ( 9) is determined only by the third summands. It is worth mentioning that, for d = 1 and d = 2, the third summands in ( 8) and ( 9) vanish in view of equality β c = 0. Let x = 0. The asymptotic behavior of the first summand in the right-hand side of ( 9) is given by [START_REF] Bogachev | The moment analysis of a branching random walk on a lattice with a single source[END_REF]. The estimate of the second summand could be obtained on account of Lemma 6 in [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF] and, in particular, relation [START_REF] Vatutin | Catalytic branching random walk and queueing systems with random number of independent servers[END_REF]. However, to avoid verifying the bounded variation of the functions p(t; x, y) and m(t; x, y) in variable t we choose another approach consisting in direct estimation of the second summand. Recall that representation (2.1.15) in [START_REF] Yarovaya | Branching Random Walks in Inhomogeneous Medium[END_REF] entails the inequalities p ′ (t; 0, 0) ≤ 0, p ′ (t; 0, 0) ≤ p ′ (t; 0, y), p ′′ (t; 0, 0) ≥ 0 and p ′′ (t; 0, 0) -p ′′ (t; 0, y) ≥ 0, t ≥ 0. Then by virtue of (2) as well as the classical results on differentiating the asymptotic formulae (see, e.g., [START_REF] Bruijn | Asymptotic methods in analysis[END_REF], Ch.7, Sec.3), for d ∈ N, one has

p ′ (t; 0, 0) ∼ - d γ d 2 t d/2+1 , p ′ (t; 0, 0) -p ′ (t; 0, y) ∼ - (d + 2)γ d (y) 2t d/2+2 , t → ∞.
Whence taking into account Lemma 5.1.2 in [START_REF] Yarovaya | Branching Random Walks in Inhomogeneous Medium[END_REF] ("lemma on convolutions") and the already proved assertion of Theorem 1 for x = 0 and y = 0 we deduce that, as t → ∞,

t 0 m(t -u; 0, 0)p ′ (u; 0, y) du = t 0 m(t -u; 0, 0) (p ′ (u; 0, y) -p ′ (u; 0, 0)) du + t 0 m(t -u; 0, 0)p ′ (u; 0, 0) du = m(t; 0, 0) -m(t; 0, 0) + o(m(t; 0, 0)) = o(m(t; 0, 0)). ( 10
)
Combining relations (2), ( 9) and ( 10) we establish Theorem 1 for d = 1 and d = 2 when x = 0. The statement of Theorem 1 for d ≥ 3 and x = 0 follows from formulae (2), ( 9) and ( 10) by Lemma 5.1.2 in [START_REF] Yarovaya | Branching Random Walks in Inhomogeneous Medium[END_REF] and in view of Theorem 1 for the known case x = y = 0. Let x = 0. Similarly to the case x = 0 we see that

t 0 p(t -u; x, 0)m ′ (u; 0, y) du = t 0 m(t -u; 0, y)p ′ (u; x, 0) du = o(m(t; 0, y)), t → ∞. (11)
Thus, the combination of ( 2), ( 8) and [START_REF] Hu | Branching random walk in Z 4 with branching at the origin only[END_REF] proves Theorem 1 for d = 1 or d = 2 and x = 0. For d ≥ 3 and x = 0 we estimate the third summand in [START_REF] Bulinskaya | The Hitting Times with Taboo for a Random Walk on an Integer Lattice[END_REF] with the help of Lemma 5.1.2 in [START_REF] Yarovaya | Branching Random Walks in Inhomogeneous Medium[END_REF], relation [START_REF] Bogachev | The moment analysis of a branching random walk on a lattice with a single source[END_REF] and the assertion of Theorem 1 for d ≥ 3 and x = 0 established above.

Auxiliary Bellman-Harris branching process

Let us briefly describe a Bellman-Harris branching process with particles of six types. It is initiated by a single particle of type i = 1, . . . , 6. The parent particle has a random life-length with a cumulative distribution function (c.d.f.) G i (t), t ≥ 0. When dying the particle produces offsprings according to a generating function f i ( s ), s = (s 1 , . . . , s 6 ) ∈ [0, 1] 6 . The new particles of type j = 1, . . . , 6 evolve independently with the life-length distribution G j (t) and an offspring generating function f j ( s ). Let M := ∂ s j f i | s=(1,...,1) i,j=1,...,6 be the mean matrix of the process. The Bellman-Harris branching process is called critical indecomposable if the Perron root of M (i.e. eigenvalue having the maximal modulus) equals 1 and for some integer n all elements of M n are positive (see, e.g., [START_REF] Sewastianow | [END_REF], Ch.4, Sec.6 and 7). Denote the number of particles of type j existing at time t by Z j (t), t ≥ 0, j = 1, . . . , 6. Set

F i (t; s ) = E i 6 j=1 s Z j (t) j , i = 1, . . . , 6, t ≥ 0, s ∈ [0, 1] 6
, where the index i means that the parent particle is of type i. In other words, F i (t; s ) is a generating function of the numbers of particles of all types existing at time t given that the process is initiated by a single particle of type i. Before demonstrating how an auxiliary Bellman-Harris process can be constructed in the framework of CBRW on Z d we have to introduce some notation. Recall that in [START_REF] Bulinskaya | The Hitting Times with Taboo for a Random Walk on an Integer Lattice[END_REF] a new notion of a hitting time with taboo was proposed for a (non-branching) random walk on Z d generated by matrix A. Namely, let τ - y,z , y, z ∈ Z d , y = z, be the time spent by the particle (performing the random walk) after leaving the starting point until the first hitting y if particle's trajectory does not pass z. Otherwise (if particle's trajectory passes point z before the first hitting y), τ - y,z = ∞. Denote by H - x,y,z (t), t ≥ 0, the improper c.d.f. of τ - y,z given that the starting point of the random walk is x ∈ Z d .

Return to CBRW on Z d . In this section we assume that CBRW may start at the origin or at a fixed point y = 0. We divide the particles population existing at time t ≥ 0 into seven groups. The particles located at time t at the origin (respectively, at y) form the first (respectively, second) group having cardinality µ(t; 0) (respectively, µ(t; y)). Next consider at time t a family of particles labeled by a collection (u, v, w) of lattice points, its cardinality being µ u,v,w (t). It consists of the particles which have left u at least once within time interval [0, t], upon the last leaving u have yet reached neither v nor w but eventually will hit v before possible hitting w. Our third group corresponds to (u, v, w) = (0, y, 0), the fourth to (y, 0, y), the fifth to (0, 0, y) and the sixth to (y, y, 0). The seventh group comprises the rest of particles not included into the above six groups. Note that the last group consists of the particles having infinite life-length since after time t they will not hit the origin any more. So, after time t these particles will not produce any offsprings and have no influence on the numbers of particles in other six groups.

Now we can introduce an auxiliary Bellman-Harris process and use it for the study of CBRW on Z d . Consider a six-dimensional Bellman-Harris process having the following c.d.f. G i and generating function f i , i = 1, . . . , 6,

G 1 (t) = 1 -e -t , f 1 ( s ) = αf (s 1 ) + (1 -α)H - 0,y,0 (0)s 2 + (1 -α)(H - 0,y,0 (∞) -H - 0,y,0 (0))s 3 +(1 -α)H - 0,0,y (∞)s 5 + (1 -α)(1 -H - 0,y,0 (∞) -H - 0,0,y (∞)), G 2 (t) = 1 -e -at , f 2 ( s ) = H - y,0,y (0)s 1 + (H - y,0,y (∞) -H - y,0,y (0))s 4 +H - y,y,0 (∞)s 6 + (1 -H - y,0,y (∞) -H - y,y,0 (∞)), G 3 (t) = H - 0,y,0 (t) -H - 0,y,0 (0) H - 0,y,0 (∞) -H - 0,y,0 (0) , f 3 ( s ) = s 2 , G 4 (t) = H - y,0,y (t) -H - y,0,y (0) H - y,0,y (∞) -H - y,0,y (0) , f 4 ( s ) = s 1 , G 5 (t) = H - 0,0,y (t) H - 0,0,y (∞) , f 5 ( s ) = s 1 , G 6 (t) = H - y,y,0 (t) H - y,y,0 (∞) , f 6 ( s ) = s 2 ,
where H - x,y,z (∞) := lim t→∞ H - x,y,z (t). The symmetry and homogeneity of the random walk generated by matrix A imply identities H - 0,y,0 ≡ H - y,0,y and H - 0,0,y ≡ H - y,y,0 , whence G 3 ≡ G 4 and G 5 ≡ G 6 . It is not difficult to see that for the branching process constructed in this way one has (µ(t; 0), µ(t; y), µ 0,y,0 (t), µ y,0,y (t), µ 0,0,y (t), µ y,y,0 (t)) Law = (Z 1 (t), . . . , Z 6 (t)), t ≥ 0. Observe that the introduced Bellman-Harris branching process with particles of six types is critical indecomposable. Indeed, it is an easy computation task to check that all entries of M 6 are positive. Furthermore, if H - 0,y,0 (0) = 0 (that is a(0, y) > 0) then already all entries of M 4 are positive. Hence, the constructed process is indecomposable. To verify its criticality note that in view of Theorem 3 in [START_REF] Bulinskaya | The Hitting Times with Taboo for a Random Walk on an Integer Lattice[END_REF] one can rewrite the first relation in (1) as follows

αf ′ (1) = 1 -(1 -α) H - 0,0,y (∞) + (H - 0,y,0 (∞)) 2 1 -H - 0,0,y (∞) . ( 12
)
Then by inspecting the explicit expression for the characteristic polynomial of the mean matrix M we deduce that it has the form

det (M -κI) = κ 2 (κ -1)R(κ)
where I is a unit matrix, κ ∈ C and

R(κ) := κ 3 + κ 2 (1 -αf ′ (1)) + κ 1 -αf ′ (1) -(2 -α)H - 0,0,y (∞) -(1 -α)(H - 0,y,0 (0)) 2 + (1 -α)(H - 0,y,0 (∞) -H - 0,y,0 (0)) 2 -(1 -α)(H - 0,0,y (∞)) 2 .
The polynomial R(κ) has no real roots greater than 1 because R(1) > 0 and R ′ (κ) > 0 for κ ≥ 1. In fact, due to identity [START_REF] Sewastianow | [END_REF] we obtain the representation with strictly positive summands

R(1) = (1 -H - 0,0,y (∞))((1 -α)H - 0,0,y (∞) + 1) + (1 -α)H - 0,y,0 (∞)(H - 0,y,0 (∞) -H - 0,y,0 (0)) + (1 -α)H - 0,y,0 (∞)(H - 0,y,0 (∞) -H - 0,y,0 (0)(1 -H - 0,0,y (∞))) 1 -H - 0,0,y (∞) + (1 -α)(H - 0,y,0 (∞)) 2 1 -H - 0,0,y (∞)
.

Moreover, if κ ≥ 1 then R ′ (κ) = 3κ 2 + 2κ(1 -αf ′ (1)) + 1 -αf ′ (1) -(2 -α)H - 0,0,y (∞) -(1 -α)(H - 0,y,0 (0)) 2 > 3 -2H - 0,0,y (∞) -H - 0,y,0 (0) > 0.
Thus, the greatest positive real root of the characteristic polynomial of M is 1. Hence, by the Frobenius theorem (see, e.g., Theorem 2 in [START_REF] Sewastianow | [END_REF], Ch.4, Sec.5) 1 is the Perron root of M. So, the auxiliary Bellman-Harris process is critical.

Denote by v = (v 1 , . . . , v 6 ) and u = (u 1 , . . . , u 6 ) the left and right positive eigenvectors corresponding to the Perron root of M such that ( u, 1 ) = 1 and ( v, u ) = 1 where 1 = (1, . . . , 1) ∈ R 6 . Taking into account [START_REF] Sewastianow | [END_REF] we rewrite the components of u and υ in the convenient form

u 1 = u 4 = u 5 = 1 -H - 0,0,y (∞) U , u 2 = u 3 = u 6 = H - 0,y,0 (∞) U , (13)
v 1 = U V , v 2 = U(1 -α)H - 0,y,0 (∞) V (1 -H - 0,0,y (∞)) , v 3 = U(1 -α)(H - 0,y,0 (∞) -H - 0,y,0 (0)) V , (14)
v 4 = v 2 (H - 0,y,0 (∞) -H - 0,y,0 (0)), v 5 = U(1 -α)H - 0,0,y (∞) V , v 6 = v 2 H - 0,0,y (∞) (15)
where the auxiliary variables U and V are defined by way of

U := 3(1 -H - 0,0,y (∞) + H - 0,y,0 (∞)), V := 3-2αf ′ (1)-(2-α)H - 0,0,y (∞)+(1-α)((H - 0,y,0 (∞)-H - 0,y,0 (0)) 2 -(H - 0,y,0 (0)) 2 -(H - 0,0,y (∞)) 2 ).
Using decomposition f (1-x) = 1-f ′ (1)x+f ′′ (1)x 2 /2+o(x 2 ), x → 0+, along with formulae ( 12)-( 15) and the definition of f ( s ) = (f 1 ( s ), . . . , f 6 ( s )), it is not difficult to verify by standard calculations that

x -v, 1 -f 1 -ux ∼ Bx 2 , x → 0+, where B := σ 2 1 -H - 0,0,y (∞) 2 2UV . (16)
In the next two lemmas we apply theorems proved in papers [START_REF] Vatutin | Discrete Limit Distributions of the Number of Particles in a Multitype Age-Dependent Branching Processes[END_REF]- [START_REF] Vatutin | On a Class of Critical Bellman-Harris Branching Processes with Several Types of Particles[END_REF] to the constructed six-dimensional Bellman-Harris branching process and then reformulate the obtained results for CBRW on Z d when d ≥ 5. Common to these theorems are the conditions of criticality and indecomposability of the Bellman-Harris process which were established above. Another common condition on the behavior of the function x -( v, 1 -f ( 1 -ux)) is fulfilled due to [START_REF] Vatutin | Discrete Limit Distributions of the Number of Particles in a Multitype Age-Dependent Branching Processes[END_REF]. However, various Vatutin's theorems involve different assumptions on the order of asymptotic decrease of the tails of G k (•), k = 1, . . . , 6. It is worth to mention that such asymptotic behavior was established in [START_REF] Bulinskaya | The Hitting Times with Taboo for a Random Walk on an Integer Lattice[END_REF], Theorem 3. Namely, our result for d ≤ 5 corresponds to condition of Theorem 1 in [START_REF] Vatutin | Discrete Limit Distributions of the Number of Particles in a Multitype Age-Dependent Branching Processes[END_REF] whereas the cases d = 6 and d ≥ 7 meet the respective conditions of Theorem 3 in [START_REF] Vatutin | On a Class of Critical Bellman-Harris Branching Processes with Several Types of Particles[END_REF] and Theorem 2 in [START_REF] Vatutin | Limit Theorem for a Critical Multitype Bellman-Harris Branching Process with Infinite Second Moments[END_REF]. Lemma 3 Given y ∈ Z 5 , y = 0, for CBRW on Z 5 one has q(t; 0, y) = o t -3/4 , q(t; y, y) = o t -3/4 , t → ∞.

Proof. To apply Theorem 1 in [START_REF] Vatutin | Discrete Limit Distributions of the Number of Particles in a Multitype Age-Dependent Branching Processes[END_REF] to the six-dimensional Bellman-Harris process constructed above for CBRW on Z 5 we verify the conditions of that theorem. According to the definition of G(•) = (G 1 (•), . . . , G 6 (•)) and by Theorem 3 in [START_REF] Bulinskaya | The Hitting Times with Taboo for a Random Walk on an Integer Lattice[END_REF] for d = 5, the variable β in condition 2) of Theorem 1 in [START_REF] Vatutin | Discrete Limit Distributions of the Number of Particles in a Multitype Age-Dependent Branching Processes[END_REF] is equal to 3/2 whereas the function L 1 (t) in the same condition tends to a constant, as t → ∞. The validity of condition 3) of Theorem 1 in [START_REF] Vatutin | Discrete Limit Distributions of the Number of Particles in a Multitype Age-Dependent Branching Processes[END_REF] is implied by Theorem 1 in [START_REF] Vatutin | Limit Theorems for Critical Markov Branching Processes with Several Types of Particles and Infinite Second Moments[END_REF] (for our process the function L 1 (t) in this theorem tends to 1/B, as t → ∞, in view of ( 16)) combined with the definition of G(•) and Theorem 3 in [START_REF] Bulinskaya | The Hitting Times with Taboo for a Random Walk on an Integer Lattice[END_REF] for d = 5. Thus, we may employ Theorem 1 in [START_REF] Vatutin | Discrete Limit Distributions of the Number of Particles in a Multitype Age-Dependent Branching Processes[END_REF]. Taking into account Theorem 3 in [START_REF] Bulinskaya | The Hitting Times with Taboo for a Random Walk on an Integer Lattice[END_REF] for d = 5 and formulae ( 13)- [START_REF] Vatutin | Critical Bellman-Harris Branching Processes Starting with a Large Number of Particles[END_REF] we deduce from Theorem 1 in [START_REF] Vatutin | Discrete Limit Distributions of the Number of Particles in a Multitype Age-Dependent Branching Processes[END_REF] that lim t→∞ E i (s

Z 2 (t) 2
| Z(t) = 0 ) = 1 for each s 2 ∈ [0, 1] and i = 1, 2 (as usual, Z(t) = (Z 1 (t), . . . , Z 6 (t)) and 0 = (0, . . . , 0) ∈ R 6 ). Setting s 2 = 0 in the last relation one has P i (Z 2 (t) > 0) = o(P i ( Z(t) = 0 )), as t → ∞. Moreover, examining the proof of Theorem 1 in [START_REF] Vatutin | Discrete Limit Distributions of the Number of Particles in a Multitype Age-Dependent Branching Processes[END_REF] we can show that for our Bellman-Harris process the slowly varying function L * (x) in the assertion of that theorem turns equivalent to 1/ √ B, as x → 0+. Consequently, the indicated in [START_REF] Vatutin | Discrete Limit Distributions of the Number of Particles in a Multitype Age-Dependent Branching Processes[END_REF] formula (0.4) can be sharpened in our case, namely, the function P i ( Z(t) = 0 ) has an order of decreasing t -3/4 , as t → ∞. Whence by the connection between CBRW on Z 5 and the auxiliary Bellman-Harris process we complete the proof.

Lemma 4

In the framework of CBRW on Z d with d ≥ 6 the following relations hold true for y ∈ Z d , y = 0,

q(t; 0, y) ∼ 2 m d σ 2 t , q(t; y, y) ∼ 2 m d G 0 (0, y) σ 2 G 0 (0, 0) t , t → ∞.
Proof. Let us apply Theorem 3 in [START_REF] Vatutin | On a Class of Critical Bellman-Harris Branching Processes with Several Types of Particles[END_REF] to our Bellman-Harris branching process when d = 6. To this end we verify whether all the conditions of Theorem 3 in [START_REF] Vatutin | On a Class of Critical Bellman-Harris Branching Processes with Several Types of Particles[END_REF] are satisfied. In view of ( 16) relation ( 6) in [START_REF] Vatutin | On a Class of Critical Bellman-Harris Branching Processes with Several Types of Particles[END_REF] is valid for our process and the function L 1 (n) in ( 6) tends to 1/B, as n → ∞. Equality [START_REF] Bulinskaya | Limit Distributions for the Number of Particles in Branching Random Walks[END_REF] in [START_REF] Vatutin | On a Class of Critical Bellman-Harris Branching Processes with Several Types of Particles[END_REF] is also satisfied due to [START_REF] Bulinskaya | Limit Distributions Arising in Branching Random Walks on Integer Lattices[END_REF] in [START_REF] Vatutin | On a Class of Critical Bellman-Harris Branching Processes with Several Types of Particles[END_REF] in view of the definition of G(•) and Theorem 3 in [START_REF] Bulinskaya | The Hitting Times with Taboo for a Random Walk on an Integer Lattice[END_REF] for d = 6. Now we may apply Theorem 3 in [START_REF] Vatutin | On a Class of Critical Bellman-Harris Branching Processes with Several Types of Particles[END_REF].

In particular, it follows that for each i = 1, 2 the expressions lim t→∞ P i (Z 1 (t) = 0| Z(t) = 0 ) and lim t→∞ P i (Z 2 (t) = 0| Z(t) = 0 ) coincide, are positive and strictly less than 1. Consequently,

P i (Z 1 (t) > 0) ∼ P i (Z 2 (t) > 0), as t → ∞.
The asymptotic behavior of q(t; 0, 0) = P 1 (Z 1 (t) > 0) and q(t; y, 0) = P 2 (Z 1 (t) > 0) can be found in [START_REF] Bulinskaya | Limit Distributions Arising in Branching Random Walks on Integer Lattices[END_REF], Lemmas 2 and 4. Thus, Lemma 4 is proved for d = 6.

For d ≥ 7 we will employ Theorem 2 in [START_REF] Vatutin | Limit Theorem for a Critical Multitype Bellman-Harris Branching Process with Infinite Second Moments[END_REF]. Condition [START_REF] Bulinskaya | Limit Distributions Arising in Branching Random Walks on Integer Lattices[END_REF] of that theorem is valid due to Theorem 1 in [START_REF] Vatutin | Limit Theorems for Critical Markov Branching Processes with Several Types of Particles and Infinite Second Moments[END_REF] (for our process, the function L 1 (t) in this theorem tends to 1/B, as t → ∞) by virtue of the definition of G(•) and Theorem 3 in [START_REF] Bulinskaya | The Hitting Times with Taboo for a Random Walk on an Integer Lattice[END_REF] for d ≥ 7. The definition of G k (•) and Theorem 3 in [START_REF] Bulinskaya | The Hitting Times with Taboo for a Random Walk on an Integer Lattice[END_REF] for d ≥ 7 also imply that ∞ 0 t dG k < ∞ for each k = 1, . . . , 6. So, all the conditions of Theorem 2 in [START_REF] Vatutin | Limit Theorem for a Critical Multitype Bellman-Harris Branching Process with Infinite Second Moments[END_REF] are satisfied and it follows that lim t→∞ P i (Z k (t) = 0| Z(t) = 0 ) = 0 for each k = 1, . . . , 6 and i = 1, 2. Hence,

P i (Z 1 (t) > 0) ∼ P i (Z 2 (t) > 0), t → ∞.
Notably, the asymptotic behavior of q(t; 0, 0) = P 1 (Z 1 (t) > 0) and q(t; y, 0) = P 2 (Z 1 (t) > 0) can be found in [START_REF] Bulinskaya | Limit Distributions Arising in Branching Random Walks on Integer Lattices[END_REF], Lemmas 2 and 4. Lemma 4 is proved for d ≥ 7.

Concluding this section we derive an integral equation in function q(•; 0, y), y = 0, which is a counterpart of equation (2.6) in [START_REF] Hu | Branching random walk in Z 4 with branching at the origin only[END_REF] for q(•; 0, 0). Our integral equation will be essentially used for proving Theorem 2 when d = 4. Before formulating the corresponding statement we have to introduce some more notation. Let τ z be the time spent by a particle performing a random walk generated by matrix A until the first hitting a point z ∈ Z d . In a similar way, τ - z is the time spent by the particle after leaving the starting point of the random walk until the first hitting the point z. If the starting point of the random walk is z then the first hitting z means the first return to z. Denote by H x,z (t) and H - x,z (t), t ≥ 0, the (improper) c.d.f. of τ z and τ - z , respectively, given that the starting point of the ran-

dom walk is x ∈ Z d . Obviously, H x,z (t) = G 2 * H - x,z (t) for t ≥ 0 and x, z ∈ Z d . Set also K(t) := αf ′ (1)G 1 (t) + (1 -α)G 1 * H - 0,0 (t) and h(s) := α(f (1 -s) -1 + f ′ (1)s), s ∈ [0, 1]
. Note that the function K(t) and the function K d (t), d ∈ N, arising in [START_REF] Hu | Branching random walk in Z 4 with branching at the origin only[END_REF] and [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF], coincide for each t ≥ 0. Thus, Lemma 2.3 in [START_REF] Hu | Branching random walk in Z 4 with branching at the origin only[END_REF] and Lemma 11 in [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF] in which the asymptotic properties of c.d.f. K d (t) and its density k d (t) are established, as t → ∞, may be applied to our function K.

Lemma 5 For y ∈ Z d , y = 0, one has

q(t; 0, y) = (1 -α)G 1 * H - 0,y (t) -H 0,y (t) + q(•; 0, y) * K(t) -h(q(•; 0, y)) * G 1 (t). (17)
Proof. Recall integral equations (see, e.g., [START_REF] Sewastianow | [END_REF], Ch.8, Sec.1) for probability generating functions F (t; s ) := (F 1 (t; s ), . . . , F 6 (t; s )) of a six-dimensional Bellman-Harris process

F i (t; s ) = s i (1 -G i (t)) + t 0 f i F (t -u; s ) dG i (u), t ≥ 0, s i ∈ [0, 1], i = 1, . . . , 6.
By setting here s = (1, 0, 1, 1, 1, 1) and substituting the explicit formulae for G j , j = 3, 4, 5, 6, and f i , i = 1, . . . , 6, we get six integral equations in functions F i (t) := F i (t; (1, 0, 1, 1, 1, 1)), t ≥ 0, i = 1, . . . , 6. Substituting the fourth and the sixth ones into the second equation and solving the obtained renewal equation in F 2 (•) we find

F 2 (t) = G 2 * 1 -H - 0,y,0 (t) -H - 0,0,y * ∞ k=0 H * k 0,0,y (t) + F 1 * G 2 * H - 0,y,0 (0) + H - 0,y,0 (•) -H - 0,y,0 (0) * ∞ k=0 H * k 0,0,y (t)
where H x,z,r (t) := G 2 * H - x,z,r (t), t ≥ 0, x, z, r ∈ Z d , z = r. Now we substitute the last equation as well as the third and the fifth equations in functions F i into the first one. After some algebraic transformations we obtain the following non-linear integral equation in function

F 1 F 1 (t) = 1 -αG 1 * (1 -f (F 1 (t))) -(1 -α)G 1 * (1 -F 1 (•)) * H - 0,0 (t) -(1 -α)G 1 * (H - 0,y (t) -H 0,y (t)) (18)
provided that the following two equalities are valid

H - 0,0 (t) = H - 0,0,y (t) + ∞ k=0 H - 0,y,0 * H * k y,y,0 * H y,0,y (t), H - 0,y (t) = H - 0,y,0 * ∞ k=0 H * k 0,0,y (t),
for each t ≥ 0. The first of them is true since any trajectory from 0 to 0 of a particle performing a random walk on Z d either passes y exactly k times, k = 1, 2, . . . , or does not hit y until the first returning to 0. Similar argument justifies the second equality as well. Recall that due to the connection between the CBRW on Z d and the constructed Bellman-Harris process one has q(t; 0, y)

= P 1 (Z 2 (t) > 0) = 1 -F 1 (t)
. Hence, rewriting [START_REF] Vatutin | Limit Theorem for a Critical Multitype Bellman-Harris Branching Process with Infinite Second Moments[END_REF] as an equation in q(t; 0, y) we come to [START_REF] Vatutin | Limit Theorems for Critical Markov Branching Processes with Several Types of Particles and Infinite Second Moments[END_REF].

Proofs of Theorems 2 and 3

First of all, we derive some integral equations to be treated in this section. Consider equation ( 5) as inhomogeneous one for differential equation ( 4) in Banach space l ∞ (Z d ). By the variation of constant formula we infer (for a similar deduction see [START_REF] Bulinskaya | Catalytic Branching Random Walk on Three-Dimensional Lattice[END_REF]) that

q(s, t; x, y) = (1 -s)m(t; x, y) - t 0 m(t -u; x, 0)h(q(s, u; 0, y)) du (19)
where q(s, t; x, y) = 1 -F (s, t; x, y), s ∈ [0, 1], t ≥ 0, x, y ∈ Z d . Substituting x = 0 in the last equation we come to an integral equation in function q(s, t; 0, y) q(s, t; 0, y) = (1 -s)m(t; 0, y) -t 0 m(t -u; 0, 0)h(q(s, u; 0, y)) du.

Note that q(0, t; x, y) is equal to q(t; x, y). Thus, on account of ( 19) one has q(t; x, y) = m(t; x, y) -t 0 m(t -u; x, 0)h(q(u; 0, y)) du.

Substituting x = 0 in [START_REF] Yarovaya | A limit theorem for critical branching random walk on Z d with a single source[END_REF] we derive an integral equation in function q(t; 0, y) q(t; 0, y) = m(t; 0, y) -t 0 m(t -u; 0, 0)h(q(u; 0, y)) du. [START_REF] Vatutin | Catalytic branching random walk and queueing systems with random number of independent servers[END_REF], Theorem 2 in [START_REF] Bulinskaya | Catalytic Branching Random Walk on a Two-Dimensional Lattice[END_REF], Theorem 4 in [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF] (item 3) and Theorem 4 in [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF] (item 4). Moreover, by virtue of Lemma 5 the proof of Theorem 2 for x = 0 and d = 4 is similar to that of Theorem 1.1 in [START_REF] Hu | Branching random walk in Z 4 with branching at the origin only[END_REF]. So, we give only a few comments on the proof of Theorem 2 for x = 0 and d ≤ 5.

If d = 1 then the equality ∞ 0 h(q(u; 0, y)) du = (1 -α)a -1 is valid. Furthermore, in view of (2), ( 8), [START_REF] Daletsky Yu | Stability of Solution of Differential Equations in Banach Space[END_REF] and Theorem 5 in [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF] one gets the useful estimate

m(t; 0, y) -(1 -α)a -1 m(t; 0, 0) = O t -3/2 , t → ∞.
When d = 2 one can check the strict inequality J(0; y) = ∞ 0 h(q(u; 0, y)) du < (1 -α)a -1 . However, if d = 3 then ∞ 0 h(q(u; 0, y)) du = (1 -α)a -1 G 0 (0, y)G -1 0 (0, 0) and

m(t; 0, y) -(1 -α)a -1 G 0 (0, y)G -1 0 (0, 0) m(t; 0, 0) = O(t -1 ), t → ∞,
by virtue of ( 2), ( 9) and Theorem 5 in [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF]. For d = 4 the first summand in [START_REF] Vatutin | Limit Theorems for Critical Markov Branching Processes with Several Types of Particles and Infinite Second Moments[END_REF] is o(t -1 ), t → ∞, by Lemma 3 in [START_REF] Bulinskaya | Limit Distributions Arising in Branching Random Walks on Integer Lattices[END_REF] and it does not contribute to the (main term of) asymptotic behavior of q(t; 0, y). As for d = 5, one has ∞ 0 h(q(u; 0, y)) du = (1 -α)a -1 G 0 (0, y)G -1 0 (0, 0) and

m(t; 0, y) -(1 -α)a -1 G 0 (0, y)G -1 0 (0, 0) m(t; 0, 0) = O t -3/2 , t → ∞,
in view of ( 2), ( 9), Theorem 5 and Corollary 1 in [START_REF] Topchii | Catalytic Branching Random Walk in Z d with branching at the origin only[END_REF]. Thus, Theorem 2 is proved for x = 0. Turn to Theorem 3 when x = 0. The proof of Theorem 3 for x = 0 is similar to those of Theorem 4 in [START_REF] Topchii | Individuals at the origin in the critical catalytic branching random walk[END_REF], Theorem 2 in [START_REF] Bulinskaya | Catalytic Branching Random Walk on a Two-Dimensional Lattice[END_REF] and Theorem 4 in [START_REF] Bulinskaya | Limit Distributions Arising in Branching Random Walks on Integer Lattices[END_REF] for d = 1, 3, d = 2 and d ≥ 5, respectively. Note only that the constant c * arising in the proof of Theorem 3 for x = 0 in contrast to its counterpart in Theorem 4 in [START_REF] Topchii | Individuals at the origin in the critical catalytic branching random walk[END_REF] is equal to σ 2 γ 2 1 a/(2(1-α)) and σ 2 aG 3 0 (0, 0)G 0 (0, y)/(8π 2 γ 2 3 (1-α)) when d = 1 and d = 3, respectively. At last, the constant c * d appearing in Theorem 4 in [START_REF] Bulinskaya | Limit Distributions Arising in Branching Random Walks on Integer Lattices[END_REF] equals (1 -α)G 0 (0, y)σ 2 / (2a G 0 (0, 0)m 2 d ) in the case of Theorem 3 for x = 0 and d ≥ 5. Since the limit theorem for µ(t; 0) when d = 4 was established by another approach, namely the moment method, we give the detailed proof of the limit theorem for µ(t; y) when d = 4. So, to complete the proof of Theorem 3 for x = 0 we dwell on the case d = 4 in detail.

Set s(t) := s(t; λ) = exp{-λ ln 2 t/(c * t)} where c * := σ 2 aG 3 0 (0, 0)G 0 (0, y)/(3γ 2 4 (1-α)), t > 0 and λ ≥ 0. By Theorems 1 and 2 for x = 0 and d = 4 we see that E 0 ( µ(t; y)| µ(t; y) > 0) = m(t; 0, y) q(t; 0, y)

∼ c * t ln 2 t , t → ∞. (23)
The inequality 1 -e -z ≤ z for z ≥ 0 yields q(s(t), u;

0, y) = E 0 1 -exp - λ ln 2 t µ(u; y) c * t ≤ λ ln 2 t c * t E 0 µ(u; y) = λ ln 2 t c * t m(u; 0, y)
where u ≥ 0 and t > 0. By virtue of this estimate combined with Theorem 1 and the inequality h(z) ≤ σ 2 z 2 (being true for z ≥ 0 small enough) one has for t large enough t/ ln 3 t 0 m(t -u; 0, 0)h(q(s(t), u; 0, y))du

≤ σ 2 λ 2 ln 4 t c * 2 t 2 t/ ln 3 t 0 m 2 (u; 0, y)m(t -u; 0, 0)du = ρ 1 (t; λ) ln t t . (24)
Here ρ 1 ∈ U and U is the class of all bounded functions ρ(t; λ) vanishing as t → ∞ uniformly in λ ∈ [0, b], whatever positive b is taken. In a similar way, we obtain where ρ 3 ∈ U. After changing the variable u = tυ and using Theorems 1 and 2 for x = 0 and d = 4 we get

I(t; λ) = 3 2q(t; 0, y) 1-1/ ln 2 t 1/ ln 3 t q 2 (s(t), tυ; 0, y) dυ (1 + ρ 4 (t; λ)), ρ 4 ∈ U. (27)
In the last integral the function q(s(t; λ), tυ; 0, y) can be replaced by q(s(tυ; λυ), tυ; 0, y). Indeed, as 1 -e -z ≤ z for z ≥ 0, we have |q(s(t; λ), tυ; 0, y) -q(s(tυ; λυ), tυ; 0, y)|

= E 0 exp - λυ ln 2 (tυ) c * tυ µ(tυ; y) -exp - λ ln 2 t c * t µ(tυ; y) ≤ E 0 1 -exp - λ(-2 ln t ln υ -ln 2 υ) c * t µ(tυ; y) ≤ λ(-2 ln t ln υ -ln 2 υ) c * t m(tυ; 0, y).
Since functions z ln z and z ln 2 z are bounded for z ∈ (0, 1), by virtue of Theorems 1 and 2 for x = 0 along with relation (26) we see that uniformly in υ ∈ [1/ ln 3 t, 1 -1/ ln 2 t] and 0 ≤ λ ≤ Λ with an arbitrary positive Λ q(s(t; λ), tυ; 0, y) q(tυ; 0, y) -q(s(tυ; λυ), tυ; 0, y) q(tυ; 0, y) → 0, t → ∞.

Set ϕ(t; λ) := q(s(t; λ), t; 0, y)/(λq(t; 0, y)), t > 0, λ ≥ 0. Then dividing both sides of (20) by λq(t; 0, y) and using ( 24)-(28) along with Theorem 2 for x = 0 and relation 1 -e -z ∼ z, z → 0, we obtain ϕ(t; λ) = 1 + ρ 5 (t; λ) -3λ 2 The argument similar to the proof of Theorem 4 in [START_REF] Topchii | Individuals at the origin in the critical catalytic branching random walk[END_REF] establishes that

lim t→∞ ϕ(t; λ) = ϕ(λ) = 2 3λ + 2 , 0 < λ ≤ Λ 0 , (29)
where Λ 0 is some positive number and ϕ(λ) is the unique solution of the equation

ϕ(λ) = 1 - 3 2 λ 0 ϕ 2 (w) dw, λ ≥ 0.
Invoking the definition of ϕ(t; λ) we rewrite relation (29) by way of

lim t→∞ E 0 exp - λ ln 2 t µ(t; y) c * t µ(t; y) > 0 = 1 -λ lim t→∞ ϕ(t; λ) = 1 3 + 2 3 • 2 3λ + 2 (30)
for 0 < λ ≤ Λ 0 . Since both the Laplace transform of a non-negative random variable and the function 1/3 + 2/3 • 2/(3λ + 2) are analytic and bounded in the domain {λ : Re λ > 0} ⊂ C, by the uniqueness theorem for analytic functions relation (30) is valid for each λ with Re λ > 0 (for an analogous deduction see, e.g., [START_REF] Vatutin | Critical Bellman-Harris Branching Processes Starting with a Large Number of Particles[END_REF]). Combining ( 23) and (30) we complete the proof of Theorem 3 for x = 0 and d = 4. Thus, Theorem 3 is proved for x = 0. Next we prove Theorems 2 and 3 when x = 0. As a preliminary we derive some more integral equations. In the framework of CBRW on Z d , the parent particle can either hit the point 0 or not within time interval [0, t]. In the latter case at time t there is a single particle on Z d located at the point y or outside it. Consequently, E x s µ(t;y) = E x s µ(t;y) I(τ 0 ≤ t) + E x s µ(t;y) I(τ 0 > t, µ(t; y) = 1) + E x s µ(t;y) I(τ 0 > t, µ(t; y) = 0) = E x s µ(t;y) I(τ 0 ≤ t) + sP x (τ 0 > t, µ(t; y) = 1) + P x (τ 0 > t, µ(t; y) = 0) (31)

where I(•) stands for the indicator of a set. Evidently, the first summand in (31) can be rewritten in the form q(s, t -u; 0, y) dH y,0 (u).

In particular, for s = 0 one has q(t; x, y) = H x,y,0 * ∞ k=0 H * k y,y,0 * (1 -G 2 (t)) + t 0 q(t -u; 0, y) dH x,0 (u) when x = y, (37) q(t; y, y) = ∞ k=0 H * k y,y,0 * (1 -G 2 (t)) + t 0 q(t -u; 0, y) dH y,0 (u).

(38)

Now we have the tools for proving Theorems 2 and 3 for x = 0. To establish Theorem 2 for x = 0 and d = 2 we employ equations (37) and (38). It is not difficult to see that the first summands in the right side of (37) and (38) are equal to p(t; x, y) -t 0 p(t -u; 0, y) dH x,0 (u) for x = y and x = y, respectively. The latter expression can be rewritten as follows p(t; x, y) -t 0 p(t -u; 0, y)dH x,0 (u) = p(t; x, y) -p(t; x, 0) + t 0 (p(t -u; 0, 0) -p(t -u; 0, y))dH x,0 (u) (39) due to the obvious relation p(t; x, 0) = t 0 p(t -u; 0, 0) dH x,0 (u). The asymptotic behavior of the first summand at the right-hand side of (39) is given by formula (2) whereas the asymptotic behavior of the second summand in (39) can be found with the help of relation (2), Lemma 3 in [START_REF] Bulinskaya | Limit Distributions Arising in Branching Random Walks on Integer Lattices[END_REF] and Lemma 5.1.2 in [START_REF] Yarovaya | Branching Random Walks in Inhomogeneous Medium[END_REF]. Finally, the first summands in (37) and (38) are O(t -3/2 ) when d = 1 and O(t -d/2 ) when d ≥ 3. Hence, the first summands in (37) and (38) are o(q(t; 0, y)), as t → ∞, by Theorem 2 for x = 0. Moreover, on account of Lemma 3 in [START_REF] Bulinskaya | Limit Distributions Arising in Branching Random Walks on Integer Lattices[END_REF] and Lemma 5.1.2 in [START_REF] Yarovaya | Branching Random Walks in Inhomogeneous Medium[END_REF] we reveal that the last summands in (37) and (38) are equivalent to q(t; 0, y) and q(t; 0, y) G 0 (x, 0)G -1 0 (0, 0) for d = 1 and d ≥ 3, respectively, as t → ∞. Hence Theorem 2 is proved for x = 0 and d = 2. As for Theorem 2 when x = 0 and d = 2 as well as Theorem 3 for x = 0, we only note that their proofs bear on analysis of equations ( 19) and [START_REF] Yarovaya | A limit theorem for critical branching random walk on Z d with a single source[END_REF]. Since the proofs are similar to that of Theorem 5 in [START_REF] Bulinskaya | Limit Distributions Arising in Branching Random Walks on Integer Lattices[END_REF], they are omitted. So, Theorems 2 and 3 are proved completely.

( 22 )

 22 Now let us prove Theorems 2 and 3 for x = 0. Since their proofs depend on d ∈ N essentially, we have to consider the cases d = 1, d = 2, d = 3, d = 4 and d ≥ 5 separately. Evidently, Theorem 2 for x = 0 and d ≥ 6 is implied by Lemma 4. Due to Lemmas 1-3 and equation (22) the proof of Theorem 2 for x = 0 in the respective cases d = 1, d = 2, d = 3 and d = 5 mainly follows the scheme proving, respectively, Theorem 2 in

E 0 E 0 E

 00 x s µ(t;y) I(τ 0 ≤ t) = {τ 0 ≤t} s µ(t;y) d P x = {τ 0 ≤t} E x s µ(t;y) τ 0 d Px = t x s µ(t;y) τ 0 = u dH x,0 (u) = t 0 s µ(t-u;y) dH x,0 (u). (32)It is easily seen that the probability at the second summand in (31) can be represented as followsP x (τ 0 > t, µ(t; y) = 1) = H x,y,0 * ∞ k=0 H * k y,y,0 * (1 -G 2 (t)) when x = y,(33)P y (τ 0 > t, µ(t; y) = 1) = ∞ k=0 H * k y,y,0 * (1 -G 2 (t)). (34)It also turns convenient to write the third summand in (31) in the formP x (τ 0 > t, µ(t; y) = 0) = 1 -H x,0 (t) -H x,y,0 * ∞ k=0 H * k y,y,0 * (1 -G 2 (t)) if x = y, (35) P y (τ 0 > t, µ(t; y) = 0) = 1 -H y,0 (t) -∞ k=0 H * k y,y,0 * (1 -G 2 (t)).(36)Combining relations (31)-(36) we come to the desired integral equations q(s, t; x, y)= (1 -s)H x,y,0 * ∞ k=0 H * k y,y,0 * (1 -G 2 (t)) + t 0q(s, t -u; 0, y) dH x,0 (u) if x = y, q(s, t; y, y) = (1 -s) ∞ k=0 H * k y,y,0 * (1 -G 2 (t)) + t 0

  It is not difficult to show that uniformly in u ∈ [t/ ln 3 t, t -t/ ln 2 t]

		t
		t-t/ ln 2 t	m(t -u; 0, 0)h(q(s(t), u; 0, y))du
	≤	σ 2 λ 2 ln 4 t c * 2 t 2	t t-t/ ln 2 t	m 2 (u; 0, y)m(t -u; 0, 0)du =	ρ 2 (t; λ) ln t t	(25)
	for ρ 2 ∈ U. ln u ∼ ln t, ln(t -u) ∼ ln t, t → ∞.	(26)
	These facts, Theorem 1 and the relation h(z) ∼ σ 2 z 2 /2, z → 0, allow us to claim that
					t-t/ ln 2 t
	I(t; λ) : =		t/ ln 3 t	h(q(s(t), u; 0, y))m(t -u; 0, 0)du
		=	σ 2 m(t; 0, 0) 2	t-t/ ln 2 t t/ ln 3 t

q 2 (s(t), u; 0, y)du (1 + ρ 3 (t; λ))

The author is grateful to Associate Professor E.B.Yarovaya for permanent attention and to Professor V.A.Vatutin for valuable remarks. Special thanks are to Professors I.Kourkova and G.Pagès for invitation to LPMA UPMC (Paris-6) where this work was started.

 [START_REF] Bogachev | The moment analysis of a branching random walk on a lattice with a single source[END_REF]The work is partially supported by the RFBR grant 10-01-00266.

1

Faculty of Mathematics and Mechanics, Lomonosov Moscow State University, Moscow 119991, Russia

E-mail address: bulinskaya@yandex.ru