# Richness Estimation with Species Identity Error

Jai-Hua Yen
Department of Agronomy
National Taiwan University
Taipei 10617, Taiwan
<EMAIL_ADDRESS>
&Chun-Huo Chiu
Department of Agronomy
National Taiwan University
Taipei 10617, Taiwan
<EMAIL_ADDRESS>

###### Abstract

Richness estimation of an interesting area is always a challenge statistical
work due to small sample size or species identity error. In the literatures,
most richness estimators were only proposed to tackle the underestimation of
the size-limited sample. However, species identity error almost occurs in each
species survey and seriously reduces the accuracy of observed, singleton, and
doubleton richness in turns to influence the behavior of richness estimator.
Therefore, to estimate the true richness, the biased collected data due to
species identity error should be modified before processing the richness
estimation work. In the manuscript, we propose a new approach to correct the
bias of richness estimation due to species identity error. First, a species
list inventory from a subplot obtained by the investigator was used to
estimate the species identity error rate. Then, we can correct the biased
observed, singleton, and doubleton richness of the raw sampling data from the
interesting area. Finally, the richness estimators proposed in the literatures
could be supplied to get the more correct estimates based on adjusted observed
data. To investigate the behavior of the proposed method, we performed
simulations by generating data sets from various species models with different
species identity error rates. For the purpose of illustration, the real data
was supplied to demonstrate our proposed approach. A presence/absence of weeds
species was surveyed in the organic farmland located at Soft Bridge County in
the North of Taiwan.

_K_ eywords Biodiversity $\cdot$ Singleton $\cdot$ Doubleton $\cdot$ Sampling
error

##  1 Introduction

Long-term biodiversity monitoring is the basis for ecological research and
promotion of organic agriculture. In recent years, more and more non-
professional citizen scientists have participated in the projects of
monitoring diversity, so the possibility of species identity errors may
increase dramatically in the collected data. Therefore, correcting the impact
of species identity error becomes an important statistical issue.

Species richness is the most intuitive and widely used as biodiversity index
due to its ecological intuitive concept and simplest form. However, due to the
sampling limitation of time or other resources, completely species inventories
in the wild field are almost unattainable goals. Therefore, the observed
richness in the sample always underestimates the true species richness in the
assemblage. In the literatures, among the discussed estimation approaches of
species richness, the non-parametric methods are widely used in practical
application, which include first order Jackknife approach, second order
Jackknife approach by Burnham and Overton (1978) and Chao1 (or Chao2) lower
bound estimator by Chao (1984) and Chao (1987). They all use the observed rare
species in the sample (i.e. singletons and doubletons) to estimate the unseen
richness in the sample. However, species identity error almost occurred in
each survey especially in vegetation sampling, and it was ignored before and
recently discussed in the literatures by Vittoz and Guisan (2007), Burg et al.
(2015), and Morrison (2015). This identity error may seriously make observed
richness biased and in turn the estimation of true richness will be seriously
biased. Therefore, without error adjustment, the species richness estimation
will be inaccurate based on original sampling data. In this manuscript, we
have proposed a modify approach to revise the biased sampling data caused by
species identity error. From the results of simulation study in session 3 show
that our adjusting approach can revise the biased observed richness, singleton
and doubleton richness. Also, the richness estimators based on the revised
data effectively correct the bias caused by the species identity error.

##  2 Methodology

In this article, we choose Chao2 lower bound estimator for incidence data as
our species richness estimator. Since we assume that species identity error
exists in the process of sampling, adjustment of richness estimator should be
considered.

First, we need to estimate the mean species identity error rate of observer or
investigator. Plant inventories from subplot of the area which the survey is
conducted. We assume that the number of species ($S_{sub}$) and the categories
of species in the subplot are known by the experiment designer but unknown by
the observer who goes conducting inventories. After conducting inventories, we
have the information that the number of observed species belongs to the
subplot ($S_{sub,e}$) and the number of observed species does not exist in the
subplot ($f_{sub,0}$).$\ X_{i}$ represents the record status of the survey of
species $i$. When $X_{i}=1$, species $i$ has been recorded. When $X_{i}=0$,
species $i$ has not been recorded. We assume the species identity error ($e$)
is a random variable follows the distribution of $F(e)$ with mean $\bar{e}$.
$r$ denotes the mean probability that a species is misidentified into another
species which belongs to the sampling plot. $f_{sub,0}$ equals to the number
of species which is misidentified and recorded as species do not exist in the
subplot. Also, if plant inventories of the subplot are correct, then
$S_{sub,e}$ should be equal to $S_{sub}$ species. However, when species
identity error occurs, $S_{sub,e}$ may not be equal to $S_{sub}$ species. When
the $i$-th species is misidentified and other species are not misidentified to
the $i$-th species, $i$-th species is not recorded. After that, we have the
equations:

$\displaystyle E\left(f_{sub,0}\right)$ $\displaystyle=\int{S_{sub}\times
e\times\left(1-r\right)}dF\left(e\right)$ (1) $\displaystyle\approx
S_{sub}\times\bar{e}\times\left(1-r\right),$

and

$\displaystyle E\left(S_{sub,e}\right)$
$\displaystyle=S_{sub}-\sum_{i=1}^{S_{sub}}E\left[I\left(X_{i}=0\right)\right]$
(2) $\displaystyle\approx
S_{sub}-S_{sub}\int{e\times\left(1-\frac{e}{\frac{S_{sub}}{r}-1}\right)^{S_{sub}-1}dF\left(e\right)}$
$\displaystyle\approx
S_{sub}-S_{sub}\times\bar{e}\times\left(1-\frac{\bar{e}\times
r}{S_{sub}-r}\right)^{S_{sub}-1}.$

By solving those two equations, we have the estimate of $\bar{e}$ and $r$
which are denoted by $\widehat{\bar{e}}$ and $\hat{r}$.

Second, the sampled observed, singleton, and doubleton richness should be
adjusted after sampling in the plot. The true observed, singleton, and
doubleton richness are denoted by $S_{obs}$, $Q_{1}$, and $Q_{2}$,
respectively. The sampled observed, singleton, and doubleton richness without
adjustment are denoted by $S_{obs,e}$, $Q_{1e}$, and $Q_{2e}$, respectively.
When species identity error occurs, the sampled observed richness is formed by
the observed species which do not misidentified and observed species which
misidentified as species do not exist in the plot. Thus, we have the expected
sampled observed richness:

$E\left(S_{obs,e}\right)\approx
E\left\\{S_{obs}\left[\left(1-e\right)+e\times\left(1-r\right)\right]\right\\}.$

Next, we have the expected observed richness adjustment:

$S_{obs,a}=\frac{S_{obs,e}}{1-\widehat{\bar{e}}\times\hat{r}}.$ (3)

When species identity error occurs, the possibilities of sampled singleton
species are as follows: (1) singleton species which do not misidentified, and
other species would not be misidentified as the singleton species at the same
time, and (2) singleton species which misidentified as species do not exist in
the plot, and other species would not be misidentified as the singleton
species at the same time. Thus, we have the expected sampled singleton
richness:

$\displaystyle E\left(Q_{1e}\right)$ $\displaystyle\approx
E\left\\{Q_{1}\left[\left(1-e\right)+e\times\left(1-r\right)\right]\times\left(1-\frac{e}{\frac{S_{obs}}{r}-1}\right)^{S_{obs}-1}\right\\}$
$\displaystyle\approx
E\left\\{Q_{1}\left[\left(1-e\right)+e\times\left(1-r\right)\right]\times
exp{\left(-e\times r\right)}\right\\}.$

Similarly, when species identity error occurs, the possibilities of sampled
doubleton species are as follows: (1) doubleton species which do not
misidentified, and other species would not be misidentified as the singleton
species at the same time, (2) doubleton species which misidentified as species
do not exist in the plot, and other species would not be misidentified as the
singleton species at the same time, and (3) when a singleton species
misidentified to a singleton species, the doubleton richness increases by one
unit, and other species would not be misidentified as the doubleton species
which is formed by singleton species at the same time. Accordingly, we have
the expected sampled doubleton richness:

$\displaystyle E\left(Q_{2e}\right)$ $\displaystyle\approx
E\left\\{Q_{2}\left[\left(1-e\right)+e\times\left(1-r\right)\right]\times
exp{\left(-e\times r\right)}\right\\}$ $\displaystyle+E\left\\{Q_{1}\times
e\times r\times\left(1-\frac{1}{T}\right)\times\frac{Q_{1}}{S_{obs,a}}\times
exp{\left(-e\times r\right)}\right\\}$

where $T$ denotes the number of sampling unit. By solving the two equations
above, we have the singleton and doubleton richness adjustment:

$Q_{1a}=\frac{Q_{1e}}{(1-\widehat{\bar{e}}\times\hat{r})exp{\left(-\widehat{\bar{e}}\times\hat{r}\right)}},$
(4)

and

$Q_{2a}=\frac{Q_{2e}-Q_{1a}\times\widehat{\bar{e}}\times\hat{r}\times\left(1-\frac{1}{T}\right)\times\frac{Q_{1a}}{S_{obs,a}}\times
exp{\left(-\widehat{\bar{e}}\times\hat{r}\right)}}{(1-\widehat{\bar{e}}\times\hat{r})\times
exp{\left(-\widehat{\bar{e}}\times\hat{r}\right)}}.$ (5)

However, the estimation of traditional Chao2 estimator will be inaccurate even
though $Q_{1a}$ and $Q_{2a}$ are asymptoticly unbiased. It causes the value of
$\frac{Q_{1a}^{2}}{2Q_{2a}}$ overestimated. Hence, we choose first-order
Jackknife and Chao2 richness estimator as the theoretical foundation of
deriving the adjusted richness estimator. We propose an adjusted richness
estimator by Taylor series expansion of
$E\left(\frac{Q_{1}^{2}}{2Q_{2}}\right)$ by the mean $Q_{1}$ and $Q_{2}$. Then
we get the difference between $\frac{[E(Q1)]^{2}}{2E(2Q2)}$ and
$E\left(\frac{Q_{1}^{2}}{2Q_{2}}\right)$ to have the adjust term:

$E\left(\frac{Q_{1}^{2}}{2Q_{2}}\right)\approx\
\frac{\left[E\left(Q_{1}\right)\right]^{2}}{E\left(2Q_{2}\right)}+\frac{V\hat{a}r\left(Q_{1}\right)}{2E\left(Q_{2}\right)}-\frac{E\left(Q_{1}\right)C\hat{o}v\left(Q_{1},Q_{2}\right)}{\left[E\left(Q_{2}\right)\right]^{2}}+\frac{\left[E\left(Q_{1}\right)\right]^{2}V\hat{a}r\left(Q_{2}\right)}{2\left[E\left(Q_{2}\right)\right]^{3}}$

where $C\hat{o}v\left(Q_{1},Q_{2}\right)=-\frac{Q_{1}Q_{2}}{\hat{S}}$,
$V\hat{a}r\left(Q_{i}\right)=Q_{i}\left(1-\frac{Q_{i}}{\hat{S}}\right)$.
Therefore, we have the adjusted richness estimator:

${\hat{S}}_{adj}=S_{obs,a}+\frac{T-1}{T}max\left\\{\left(\frac{Q_{1a}^{2}}{2Q_{2a}}-\frac{Q_{1a}}{2Q_{2a}}-\frac{Q_{1a}^{2}}{2Q_{2a}^{2}}\right),\
0\right\\}$ (6)

When $0\leq Q_{2a}\leq 1$, by simulation studies, the adjusted richness
estimator will be:

${\hat{S}}_{adj}=S_{obs,a}+\frac{T-1}{T}Q_{1a}$ (7)

##  3 Result

###  3.1 Simulation Results

To test the performance of the adjusted richness estimator, we presented the
simulation results by several species detection models and different settings
of number of sampling units. We fixed $S_{sub}=40$ and $S=100$. 500 simulation
data sets were generated and 200 bootstrapping trials were conducted by each
simulation data. The bootstrapping method is regenerating $S_{obs,a}$,
$Q_{1a}$, and $Q_{2a}$ by binomial distribution independently in order to
increase the estimated standard error while the traditional bootstrapping
method usually underestimates the standard error in this case. In true method,
the estimation of species richness used the traditional Chao2 estimator by the
data without species identity error. In observed method, the estimation of
species richness used the traditional Chao2 estimator by the data with species
identity error. In adjusted method, the estimation of species richness used
the adjusted richness estimator by the data with species identity error.

When species identity error occurs, the estimate of species richness by
observed method will be underestimated, which causes larger bias. The large
bias still exists even though the increase of the number of sampling units.
Since adjusted method slightly overestimated species richness when the species
identity error rate is large, it reduces a great quantity of bias. The
variation of observed method is lower, and it remains the same by different
species identity error rate. The adjusted method has a higher variation. When
species identity error rate is larger, the variation of adjusted method is
larger. By evaluating both bias and variation, the observed method has a
larger RMSE (Root Mean Square Error) due to its larger bias. The adjusted
method has about half RMSE of the observed method when the number of sampling
unit is large.

Table 1: Comparison of species richness estimator for incidence data based on
500 simulation data sets and 200 bootstrapping trials under random uniform (0,
1) model, with $\bar{p}=0.51$, $CV=0.53$, $S=100$, $S_{sub}=40$, $T=5$, and
$r=0.91$.

$\bar{e}$ | $\widehat{\bar{e}}$ | Method | $S_{obs}$ | $Q_{1}$ | $Q_{2}$ | $\hat{S}$ | Bias | s.e. | $\hat{s.e.}$ | RMSE |
---|---|---|---|---|---|---|---|---|---|---|---
0 | 0 | True | 85.2 | 15.3 | 17.3 | 91.37 | -8.63 | 4.82 | 4.19 | 9.89 |
0.053 | 0.058 | Observed | 81.5 | 13.9 | 15.8 | 87.22 | -12.78 | 5.46 | 4.06 | 13.9 |
|  | Adjusted | 86.3 | 15.6 | 17.5 | 92.05 | $-7.95^{\star}$ | 7.17 | 8.33 | $10.71^{\dagger}$ |
0.097 | 0.098 | Observed | 78.3 | 13.2 | 14.8 | 83.72 | -16.28 | 5.29 | 3.95 | 17.12 |
|  | Adjusted | 86.3 | 15.9 | 17.5 | 92.2 | $-7.8^{\star}$ | 7.92 | 9.4 | $11.12^{\dagger}$ |
0.15 | 0.157 | Observed | 74 | 11.7 | 13.4 | 78.86 | -21.14 | 5.24 | 3.75 | 21.78 |
|  | Adjusted | 86.8 | 16 | 17.6 | 92.89 | $-7.11^{\star}$ | 10.33 | 10.2 | $12.54^{\dagger}$ |
0.199 | 0.209 | Observed | 70.7 | 10.3 | 12.7 | 74.71 | -25.29 | 5.01 | 3.34 | 25.78 |
|  | Adjusted | 88.3 | 15.8 | 18.5 | 94.34 | $-5.66^{\star}$ | 14.05 | 11.12 | $15.15^{\dagger}$ |

⋆ Denotes the smaller bias. † Denotes the smaller RMSE.

Table 2: Comparison of species richness estimator for incidence data based on
500 simulation data sets and 200 bootstrapping trials under random uniform (0,
1) model, with $\bar{p}=0.51$, $CV=0.53$, $S=100$, $S_{sub}=40$, $T=20$, and
$r=0.91$.

$\bar{e}$ | $\widehat{\bar{e}}$ | Method | $S_{obs}$ | $Q_{1}$ | $Q_{2}$ | $\hat{S}$ | Bias | s.e. | $\hat{s.e.}$ | RMSE |
---|---|---|---|---|---|---|---|---|---|---|---
|  | True | 95.3 | 4.1 | 3.9 | 98.8 | -1.2 | 4.9 | 4.25 | 5.06 |
0.053 | 0.055 | Observed | 91.2 | 3.9 | 3.6 | 94.8 | -5.2 | 5.46 | 4.45 | 7.53 |
|  | Adjusted | 96.1 | 4.3 | 4 | 97.85 | $-2.15^{\star}$ | 5.26 | 5.39 | $5.68^{\dagger}$ |
0.097 | 0.095 | Observed | 87.3 | 3.3 | 3.5 | 90.1 | -9.9 | 5.15 | 3.76 | 11.15 |
|  | Adjusted | 95.8 | 4 | 4.1 | 97.1 | $-2.9^{\star}$ | 6.52 | 5.72 | $7.14^{\dagger}$ |
0.15 | 0.151 | Observed | 82.9 | 3.1 | 2.9 | 85.61 | -14.39 | 5.21 | 3.79 | 15.31 |
|  | Adjusted | 96.7 | 4.1 | 3.9 | 97.94 | $-2.06^{\star}$ | 8.94 | 6.23 | $9.17^{\dagger}$ |
0.199 | 0.21 | Observed | 79.2 | 2.9 | 2.7 | 81.79 | -18.21 | 5.25 | 3.66 | 18.95 |
|  | Adjusted | 98.8 | 4.4 | 4 | 100.5 | $0.46^{\star}$ | 11.52 | 7.04 | $11.53^{\dagger}$ |

⋆ Denotes the smaller bias. † Denotes the smaller RMSE.

Table 3: Comparison of species richness estimator for incidence data based on
500 simulation data sets and 200 bootstrapping trials under $(0.8\times
Uniform\left(0.1,\ 0.3\right)+0.2\times Uniform\left(0.4,\ 1\right))$, with
$\bar{p}=0.29$, $CV=0.7$, $S=100$, $S_{sub}=40$, $T=5$, and $r=0.67$.

$\bar{e}$ | $\widehat{\bar{e}}$ | Method | $S_{obs}$ | $Q_{1}$ | $Q_{2}$ | $\hat{S}$ | Bias | s.e. | $\hat{s.e.}$ | RMSE |
---|---|---|---|---|---|---|---|---|---|---|---
|  | True | 72 | 32.4 | 19.8 | 94.98 | -5.02 | 11.38 | 10.76 | 12.44 |
0.053 | 0.056 | Observed | 69.7 | 30.4 | 19 | 90.91 | -9.09 | 11.13 | 10.22 | 14.37 |
|  | Adjusted | 72.5 | 32.9 | 19.9 | 94.7 | $-5.3^{\star}$ | 12.56 | 12.98 | $13.63^{\dagger}$ |
0.097 | 0.1 | Observed | 67.3 | 28.8 | 18.3 | 87.32 | -12.68 | 11.12 | 9.91 | $16.87^{\dagger}$ |
|  | Adjusted | 72.3 | 33.1 | 19.8 | 95.78 | $-4.22^{\star}$ | 21.14 | 15.12 | 21.56 |
0.15 | 0.155 | Observed | 64.7 | 26.4 | 17.7 | 82.27 | -17.73 | 11.77 | 9.06 | $21.28^{\dagger}$ |
|  | Adjusted | 72.7 | 33.1 | 20.1 | 96.26 | $-3.74^{\star}$ | 21.81 | 17.28 | 22.13 |
0.199 | 0.203 | Observed | 63.1 | 24.9 | 17.2 | 78.81 | -21.19 | 9.08 | 8.36 | 23.06 |
|  | Adjusted | 73.9 | 33.9 | 20.3 | 98.02 | $-1.98^{\star}$ | 22.62 | 19.58 | $22.71^{\dagger}$ |

⋆ Denotes the smaller bias. † Denotes the smaller RMSE.

Table 4: Comparison of species richness estimator for incidence data based on
500 simulation data sets and 200 bootstrapping trials under $(0.8\times
Uniform\left(0.1,\ 0.3\right)+0.2\times Uniform\left(0.4,\ 1\right))$, with
$\bar{p}=0.29$, $CV=0.7$, $S=100$, $S_{sub}=40$, $T=20$, and $r=0.67$.

$\bar{e}$ | $\widehat{\bar{e}}$ | Method | $S_{obs}$ | $Q_{1}$ | $Q_{2}$ | $\hat{S}$ | Bias | s.e. | $\hat{s.e.}$ | RMSE |
---|---|---|---|---|---|---|---|---|---|---|---
|  | True | 97.8 | 7 | 11.9 | 100.25 | 0.25 | 2.56 | 2.43 | 2.57 |
0.053 | 0.056 | Observed | 94.7 | 6.6 | 11.1 | 97.08 | -2.92 | 2.98 | 2.37 | $4.17^{\dagger}$ |
|  | Adjusted | 98.5 | 7.1 | 12 | 100.62 | $0.62^{\star}$ | 4.55 | 5.8 | 4.59 |
0.097 | 0.102 | Observed | 91.5 | 6.2 | 10.4 | 93.78 | -6.22 | 3.72 | 2.34 | 7.25 |
|  | Adjusted | 98.6 | 7.2 | 12 | 100.76 | $0.76^{\star}$ | 6.24 | 6.97 | $6.29^{\dagger}$ |
0.15 | 0.151 | Observed | 88.2 | 5.8 | 9.8 | 90.42 | -9.58 | 3.69 | 2.31 | 10.27 |
|  | Adjusted | 98.5 | 7.2 | 12.1 | 100.62 | $0.62^{\star}$ | 7.5 | 7.5 | $7.53^{\dagger}$ |
0.199 | 0.204 | Observed | 85.4 | 5.4 | 9.1 | 87.45 | -12.55 | 4.2 | 2.3 | 13.24 |
|  | Adjusted | 99.9 | 7.3 | 12.2 | 102.08 | $2.08^{\star}$ | 9.64 | 7.98 | $9.86^{\dagger}$ |

⋆ Denotes the smaller bias. † Denotes the smaller RMSE.

###  3.2 Real Data Analysis

The data set was collected of weed species from organic farmland located at
Soft Bridge county in the North of Taiwan. There are 12 transect lines with
length 20m each were conducted. Only the incidence (detection or non-
detection) of species in each transect line was recorded. Before richness
estimation, a subplot occupied by 40 known weed species was treated as the
testing of the degree of investigator’s skill. Compare these 40 weed species
list with the inventories of the investigator, we have $S_{sub}=40$,
$S_{sub,e}=35$, and $f_{sub,0}=1$. Therefore, we have the estimate of
$\widehat{\bar{e}}=0.14$ and $\hat{r}=0.82$ based on equations (1) and (2).
Many of the misidentified species were misidentified as species which did not
exist in the plot. The summary of the frequency counts of weed species is in
Table 5. The result using our adjusted estimator is in Table 6. By simulation
studies, the error rate is high in this case. Hence, the estimate of species
richness using row data directly underestimates and the adjusted estimator
should be applied to get the accurate estimate of species richness.

Table 5: Summary of the data set of weed species frequency counts at Soft Bridge county in the North of Taiwan, with $T=12$. Frequency | $Q_{1}$ | $Q_{2}$ | $Q_{3}$ | $Q_{4}$ | $Q_{5}$ | $Q_{6}$ | $Q_{7}$ | $Q_{8}$ | $Q_{9}$ | $Q_{10}$ | $Q_{11}$ | $Q_{12}$
---|---|---|---|---|---|---|---|---|---|---|---|---
Counts | 18 | 9 | 12 | 8 | 6 | 4 | 1 | 4 | 3 | 3 | 2 | 3
Table 6: Species richness adjustment for data set of weed species from Soft Bridge county in the North of Taiwan in farmland, with $T=12$, $\hat{r}=0.82$, and $\widehat{\bar{e}}=0.14$. Method | $S_{obs}$ | $Q_{1}$ | $Q_{2}$ | $\hat{S}$ | $\hat{s.e.}$
---|---|---|---|---|---
Observed | 74.0 | 19.0 | 9.0 | 92.4 | 11.27
Adjusted | 83.6 | 24.1 | 10.6 | 105.4 | 18.68

##  4 Discussion and Conclusion

Species richness is the simplest and most popular measure of biodiversity. The
approach of estimating species richness is widely discussed due to its
application in many ecological or agricultural issues mentioned by Carvalheiro
et al. (2011) and Garibaldi et al. (2013). In the manuscript, we demonstrated
the effect of species identity error while sampling in estimating species
richness. When the mean probability that a species is misidentified into
another species which belongs to the sampling plot is high, the observed
richness and singleton richness will be seriously negative biased which
implying most richness estimators’ serious underestimation even though
increasing sampling units. Our simulations show that the adjusted richness
estimator removes a large proportion of the negative bias under different
settings of sampling units, species identity error, and species detection
model. We suggest that the adjusted richness estimator for incidence data
should be applied to estimate species richness of the target region since
species identity error occurs almost in every investigation of species.

##  5 Acknowledgements

The research was supported by the Taiwan National Science Council under
Project 107-2118-M-002-001-MY2 and Council of Agriculture under Project
107AS-1.2.7-ST-a6.

## References

  * [1] Burg, S., Rixen, C., Stöckli, V. & Wipf, S. (2015).  Observation bias and its causes in botanical surveys on high-alpine summits.  Journal of Vegetation Science 26, 191–200.
  * [2] Burnham, K. P., & Overton, W. S. (1978).  Estimation of the size of a closed population when capture probabilities vary among animals.  Biometrika 65(3), 625-633.
  * [3] Carvalheiro, L. G., Veldtman, R., Shenkute, A. G., Tesfay, G. B., Pirk, C. W. W., Donaldson, J. S., & Nicolson, S. W. (2011).  Natural and within-farmland biodiversity enhances crop productivity.  Ecology letters 14(3), 625-633.
  * [4] Chao, A. (1984).  Non-parametric estimation of the number of classes in a population.  Scandinavian Journal of Statistics 11, 265-270.
  * [5] Chao, A. (1987).  Estimating the population size for capture-recapture data with unequal catchability.  Biometrics 43, 783-791.
  * [6] Garibaldi, L. A., Steffan-Dewenter, I., Winfree, R., Aizen, M. A., Bommarco, R., Cunningham, S. A., … & Bartomeus, I. (2013).  Wild pollinators enhance fruit set of crops regardless of honey bee abundance.  Science 339(6127), 1608-1611.
  * [7] Morrison, LW. (2015).  Observer error in vegetation surveys: a review.  Journal of Plant Ecology 9, 367-379.
  * [8] Vittoz, P. & Guisan, A. (2007).  How reliable is the monitoring of permanent vegetation plots? A test with multiple observers.  Journal of Vegetation Science 18, 416-422.


functions up to some order in $\epsilon$, then the Feynman integral
$A_{g}(x_{1},x_{2},x_{3})$ is approximated by the sum of the simpler Feynman
integrals. From (5) we obtain the condition that in each graph of the
approximation the edge $01$ acquires a weight such that the scaling $N_{g}$ of
$g$ is preserved, see (6). Figures 12 and 13 are important examples for such
approximations to various order in $\epsilon$ in four dimensions. We may
replace $g$ in $G$ by the approximation if the order in $\epsilon$ suffices
and if no singular (by power-counting) subgraph of $G$ uses edges in $g$.

The latter condition can be relaxed if the power-counting of the singularity
is preserved by each graph in the approximation. As an example of this (rather
rare) situation, consider the first approximation of the bubble three-star in
Figure 12. Assume the left hand side is substituted into an ambient graph with
an edge of weight 1 from the left to the top vertex. In this case one obtains
a logarithmic singularity in the shape of Dunce’s cap (spanned by the left,
top, and central vertex). The cap contains the bubble, the central vertex, and
the top edge of the graph $g$ on the left hand side in Figure 12. The weight
of this chain is $3-(\lambda+1)/\lambda=1-\epsilon/2\lambda$, see (30). The
left-top edge on the sole graph in the first approximation on the right hand
side reproduces this weight. So, the first approximation in Figure 12 can be
used in the presence of a left-top edge of weight 1 in the ambient graph. The
other approximations in Figure 12 cannot be used in this case because there
exist graphs in the approximation with left-top edges of weight 1 ($\neq
1-\epsilon/2\lambda$).

Some of the identities in Figures 12 and 13 generalize to weights
$(1+a_{e}\epsilon)/\lambda$, see (27), which makes them significantly more
powerful. In HyperlogProcedures [34] many more three-vertex approximations are
implemented. A systematic and comprehensive search for three-vertex
approximations is still missing.

Approximations of three-vertex cuts are an important tool for the calculation
of renormalization functions in four-dimensional $\phi^{4}$ theory. There
exist analogous three-vertex cut approximations in six-dimensional $\phi^{3}$
theory. As a rule of thumb, however, IBP (see Section 8.2) is more important
in $\phi^{3}$ theory whereas three-vertex cut approximations are more
important in $\phi^{4}$ theory.

##  10\. Subtraction of subdivergences by rerouting

Figure 14. Linear combinations of subgraphs with lifted singularities. The
shaded areas connect to an ambient graph to represent a graphical function.
All but the leftmost graphs have two-chains of edges which simplify by (30).
The central vertices of the two-chains must be internal in the full graphical
function.

When using approximations (Section 9) for the calculation of a graphical
function, one can only handle a certain number of orders in the regulator
$\epsilon$. In this context it can be very helpful to reduce the pole order of
a graphical function by subtracting subdivergences. Subtraction of
subdivergences has already been used to handle divergences for appending
edges, see Section 5. Here, we use a purely graphical method which was
developed by F. Brown and D. Kreimer in [14]. We change the routing of some
edges in a graph $G$ in such a way that the singularity structure of the
original graph is maintained. Further, we want to generate a two-chain of
edges in the subtraction graph(s), so that we can use (30) to eliminate the
internal vertex of the two-chain. We subtract the rerouted graph(s) from $G$
to lift the singular substructure in $G$. It can be necessary to subtract (or
add) several graphs to $G$ to fully cancel a subdivergence, see Figure 14. A
given graphical function $f_{G}(z)$ can thus be expressed as a sum of graphs
with a reduced subdivergence plus the subtraction graphs with a reduced
vertex. Both terms are simplified, which may allow one to compute $f_{G}(z)$,
see Section 12.

A particularly accessible case occurs when one needs to know $f_{G}(z)$ only
to order $\epsilon^{-1}$ (note that $G$ may have a high pole order due to many
nested subdivergences). In this case a complete regularization by rerouting
has low degree $\epsilon^{0}$ and can therefore be ignored. One may thus
replace $G$ by the sum of the subtraction graphs.

If a period $P_{G}$ (or a graphical function $f_{G}(z)$) needs to be
calculated to order $\epsilon^{0}$, a complete regularization by rerouting
allows one to work in integer dimensions. Then, the period $P_{G}$ can be
completed without generating many edges that vanish in the limit $\epsilon\to
0$, see Sections 7 and 9.

Using the renormalization Hopf algebra by A. Connes and D. Kreimer [19, 28] it
is a purely combinatorial task to check if a linear combination of graphs is
finite. The reduced coproduct is given by

(44) $\Delta^{\prime}(G)=\sum_{g\subset G}g\otimes G/g,$

which extends linearly to sums of graphs. In (44) the sum is over all disjoint
unions $g$ of 1PI subgraphs which are logarithmically divergent by power-
counting. The graph $G/g$ is obtained from $G$ by contracting all edges of
$g$. A linear combination of graphs $\Gamma$ is finite if and only if
$\Delta^{\prime}(\Gamma)=0$. To show that $\Delta^{\prime}(\Gamma)=0$ one can
use the fact that tadpoles vanish in dimensionally regularized massless
theories (one needs this e.g. to prove that the second line is finite).

As an example, consider the first line $\Gamma$ in Figure 14. For both graphs
in $\Gamma$ the only 1PI subdivergent graph is the bubble. Both graphs have
the same contraction, so $\Delta^{\prime}(\Gamma)=0$.

After rerouting, approximations can only be used simultaneously for all graphs
in the regularized sum. The approximation is limited to the part of the
ambient graph which is not affected by the rerouting.

In HyperlogProcedures [34], subtraction by rerouting has only been implemented
for simple cases. For the eight loop result of the $\phi^{4}$ gamma function
it was necessary to reroute some more complex graphs by hand. A general
algorithm for rerouting would be desirable in the future. Note that by three-
valence, six-dimensional $\phi^{3}$ theory has no overlapping divergences
which leads to a significantly simpler rerouting algorithm.

##  11\. HyperInt

If all attempts fail to calculate a graphical function, then a last resort is
parametric integration.

Graphical functions have parametric representations where the integration is
over edge-variables instead of vertices [22]. A method of brute force
parametric integration has been devised by F. Brown [8, 9]. The method was
developed and implemented in Maple by E. Panzer (HyperInt [29]).

The method is restricted to graphical functions whose parametric
representation is ‘linearly reducible’. In practice, further restrictions come
from severe memory and time consumption in many calculations. The complexity
of a parametric integration heavily depends on the number of edges, the order
in $\epsilon$, the topology of the graph and the weights of the edges. It
proved very efficient to parametrically calculate families of graphical
functions with weights $(n_{e}+a_{e}\epsilon)/\lambda$, see (27), where
$n_{e}\in{\mathbb{Z}}$ is fixed and $a_{e}\in{\mathbb{R}}$ is a free
parameter. The parameters $a_{e}$ can be handled by HyperInt with moderate
loss in speed and extra memory demand. Knowledge of graphical functions with
parameters can be useful for the calculation of large sets of relevant
graphical functions and Feynman periods. Note that from rather small graphical
function a whole series of larger graphical functions can be constructed with
the techniques of the previous sections. Small graphical functions can thus be
considered as germs which considerably enrich the set of available results.

Parametric integration in four dimensions is limited to small graphs and low
orders in $\epsilon$ (typically three internal vertices with three orders in
$\epsilon$, which already needs some tweaks by the author). In six dimensions
the restrictions are even more severe, so that mostly graphs with two internal
vertices were computed to low orders in $\epsilon$. This difficulty results
from a higher power in the denominator of the parametric integrand and it is a
main source for the slow progress in calculating the $\phi^{3}$ beta function
to loop order six. There exist ideas to bypass the difficulty (most
prominently by analytically pre-calculating the first integrations using
general Dodgson identities [9]). More intense work on the level of parametric
integration however has low priority to the author of this article because it
stands somewhat outside the theory of graphical functions.

##  12\. A worked example

Figure 15. The reduction of the graphical function $f_{G}(z)$ uses the graphs
$G_{1},\ldots,G_{14}$ in intermediate steps. Edges with no labels have weight
$1$.  Figure 16. The graphical function $f_{G}(z)$ in Figure 15 can be
expressed in terms of the graphical functions $F^{3}_{a,b,c}(z)$ (order
$\epsilon^{2}$), $F^{4}_{a,b,c,d}(z)$ (order $\epsilon^{3}$),
$F^{5}_{a,b,c,d,e}(z)$ (order $\epsilon^{1}$), and $F^{6}_{a,b,c,d,e,f}(z)$
(order $\epsilon^{1}$), see (12). The graphical function
$F^{5}_{a,b,c,d,e}(z)$ can be calculated from $\widetilde{F}^{51}_{a,b,c}(z)$
and $\widetilde{F}^{52}_{a,b,c}(z)$. Edges with no labels have weight $1$.
(The graphical functions are represented by their graphs.)

In this section we show how to combine the tools of the previous sections to
calculate the graphical function of the graph $G$ on the top left of Figure 15
to order $\epsilon^{0}$. The graph $G$ has a triple subdivergence (Dunce’s cap
and a bubble) giving rise to a low degree $\epsilon^{-3}$ in $f_{G}(z)$. The
graph $G$ is not constructible, so we want to approximate Dunce’s cap by the
identities in Figure 13. We can use the second identity as it covers four
orders in $\epsilon$ (from $\epsilon^{-2}$ to $\epsilon^{1}$). However, the
second identity in Figure 13 has a rather complicated first term where one
loses one internal vertex but still has to calculate four orders in $\epsilon$
(after the substitution in $G$ one needs the coefficients for $\epsilon^{-2}$
to $\epsilon^{1}$). In the first equation this term is absent and the
remaining terms either feature less orders in $\epsilon$ or they have lost
both internal vertices. To be able to use the first identity we regularize $G$
by subtracting the second graph $G_{1}$ in Figure 15. In the subtraction, the
rightmost bubble is regularized (see Section 10 and the first line in Figure
14). For the difference $f_{G}(z)-f_{G_{1}}(z)$ we only need to calculate the
three coefficients of $\epsilon^{-2}$ to $\epsilon^{1}$. Notice that the
approximation does not use the rerouted edge. We obtain

$\displaystyle
f_{G}(z)-f_{G_{1}}(z)=-2(f_{G_{3}}(z)+f_{G_{4}}(z)-f_{G_{5}}(z)-f_{G_{6}}(z)-f_{G_{7}}(z)+f_{G_{8}}(z))$
$\displaystyle\qquad+\frac{2}{\epsilon^{2}(1-\epsilon)(1-\frac{3\epsilon}{2})\Gamma(\lambda)^{2}}(f_{G_{9}}(z)+f_{G_{10}}(z)-f_{G_{11}}(z)-f_{G_{12}}(z)-f_{G_{13}}(z)+f_{G_{14}}(z))+\mathcal{O}(\epsilon).$

The only singularity in the graphical functions of $G_{3}$, $G_{4}$, and
$G_{5}$ is the bubble. They thus have pole order $-1$ and we can approximate
the bubble using the first identity in Figure 12. Using the graphical
functions in Figure 16 we obtain

$\displaystyle f_{G_{3}}(z)$
$\displaystyle=\frac{2F^{6}_{-\frac{1}{2},\frac{1}{2},-\frac{1}{2},0,0,0}(z)}{\epsilon(1-\epsilon)\Gamma(\lambda)}+\mathcal{O}(\epsilon),\qquad
f_{G_{4}}(z)=\frac{2F^{6}_{-\frac{1}{2},-\frac{1}{2},\frac{1}{2},0,0,0}(z)}{\epsilon(1-\epsilon)\Gamma(\lambda)}+\mathcal{O}(\epsilon),$
$\displaystyle f_{G_{5}}(z)$
$\displaystyle=\frac{2F^{5}_{-\frac{1}{2},\frac{1}{2},-\frac{1}{2},0,0}(z)}{\epsilon(1-\epsilon)\Gamma(\lambda)}+\mathcal{O}(\epsilon).$

In the graphs $G_{6}$, $G_{7}$, $G_{8}$, $G_{12}$, $G_{13}$, and $G_{14}$ we
reduce the two-chain of edges using (30). This gives

$\displaystyle f_{G_{6}}(z)$
$\displaystyle=f_{G_{7}}(z)=\frac{2F^{6}_{-\frac{1}{2},0,0,0,0,0}(z)}{\epsilon(1-\epsilon)\Gamma(\lambda)},\qquad
f_{G_{8}}(z)=\frac{2F^{5}_{-\frac{1}{2},0,0,0,0}(z)}{\epsilon(1-\epsilon)\Gamma(\lambda)},$
$\displaystyle f_{G_{12}}(z)$
$\displaystyle=f_{G_{13}}(z)=\frac{2F^{4}_{-\frac{1}{2},0,0,-1}(z)}{\epsilon(1-\epsilon)\Gamma(\lambda)},\qquad
f_{G_{14}}(z)=\frac{2F^{4}_{-\frac{3}{2},0,0,0}(z)}{\epsilon(1-\epsilon)\Gamma(\lambda)}.$

For the graphs $G_{9}$ and $G_{11}$ we use the third approximation in Figure
12. In some of the resulting graphs we swap the external vertices 0 and 1
which transforms $z$ to $1-z$, see (29). This gives

$\displaystyle f_{G_{9}}(z)=$ $\displaystyle-
f_{G_{10}}(z)-F^{3}_{0,-1,0}(1-z)-\frac{\epsilon}{2(1-\epsilon)}(F^{5}_{0,-1,0,0,0}+F^{6}_{0,-1,0,0,0,0}+F^{6}_{0,0,-1,0,0,0})(1-z)$
$\displaystyle+\,\Big{(}\frac{2}{\epsilon(1-\epsilon)\Gamma(\lambda)}+\zeta(3)\epsilon^{2}\Big{)}(F^{4}_{-\frac{1}{2},0,0,-1}+F^{4}_{0,-\frac{1}{2},0,-1}+F^{4}_{0,0,-\frac{1}{2},-1})(z)+\mathcal{O}(\epsilon^{3}),$
$\displaystyle f_{G_{11}}(z)=$ $\displaystyle-
f_{G_{11}}(z)-F^{3}_{-1,0,0}(1-z)-\frac{\epsilon}{2(1-\epsilon)}(F^{5}_{-1,0,0,0,0}+2F^{6}_{-1,0,0,0,0,0})(1-z)$
$\displaystyle+\,\Big{(}\frac{2}{\epsilon(1-\epsilon)\Gamma(\lambda)}+\zeta(3)\epsilon^{2}\Big{)}(F^{4}_{-\frac{3}{2},0,0,0}+2F^{4}_{-1,-\frac{1}{2},0,0})(z)+\mathcal{O}(\epsilon^{3}).$

In the rerouted graph $G_{1}$ we eliminate the two-chain using (30),

$f_{G_{1}}(z)=\frac{2f_{G_{2}}(z)}{\epsilon(1-\epsilon)\Gamma(\lambda)}.$

In the reduced graph $G_{2}$ we use the second identity in Figure 13, yielding

$\displaystyle f_{G_{2}}(z)=$
$\displaystyle\,\frac{3}{\epsilon(1-\epsilon)\Gamma(\lambda)}F^{3}_{-\frac{1}{2},0,-\frac{1}{2}}(z)$
$\displaystyle+\frac{1}{4(1-\epsilon)(1-\frac{3\epsilon}{2})\Gamma(\lambda)}(2F^{6}_{-\frac{1}{2},0,0,-3,-\frac{3}{2},4}-F^{5}_{1,0,0,-1,-1})(z)$
$\displaystyle-\frac{1}{\epsilon^{2}(1-\epsilon)(1-\frac{3\epsilon}{2})\Gamma(\lambda)^{2}}(2F^{4}_{1,-\frac{3}{2},-1,0}-F^{4}_{\frac{3}{2},-\frac{3}{2},-\frac{3}{2},0})(1-z)+\mathcal{O}(\epsilon^{2}).$

Altogether, we are able to express the graphical function $f_{G}(z)$ in terms
of the graphical functions $F^{3}_{a,b,c}(z)$, $F^{4}_{a,b,c,d}(z)$,
$F^{5}_{a,b,c,d,e}(z)$, and $F^{6}_{a,b,c,d,e,f}(z)$ in Figure 16,

$\displaystyle
f_{G}(z)=\frac{(8F^{6}_{-\frac{1}{2},0,0,0,0,0}-4F^{6}_{-\frac{1}{2},\frac{1}{2},-\frac{1}{2},0,0,0}-4F^{6}_{-\frac{1}{2},-\frac{1}{2},\frac{1}{2},0,0,0}-4F^{5}_{-\frac{1}{2},0,0,0,0}+4F^{5}_{-\frac{1}{2},\frac{1}{2},-\frac{1}{2},0,0})(z)}{\epsilon(1-\epsilon)\Gamma(\lambda)}$
$\displaystyle\quad+\,\frac{(F^{6}_{-1,0,0,0,0,0}-F^{6}_{0,-1,0,0,0,0}-F^{6}_{0,0,-1,0,0,0}+\frac{1}{2}F^{5}_{-1,0,0,0,0}-F^{5}_{0,-1,0,0,0})(1-z)}{\epsilon(1-\epsilon)^{2}(1-\frac{3\epsilon}{2})\Gamma(\lambda)^{2}}$
$\displaystyle\quad+\,\frac{(F^{6}_{-\frac{1}{2},0,0,-3,-\frac{3}{2},4}-\frac{1}{2}F^{5}_{1,0,0,-1,-1})(z)}{\epsilon(1-\epsilon)^{2}(1-\frac{3\epsilon}{2})\Gamma(\lambda)^{2}}+\frac{(F^{3}_{-1,0,0}-2F^{3}_{0,-1,0})(z)}{\epsilon^{2}(1-\epsilon)(1-\frac{3\epsilon}{2})\Gamma(\lambda)^{2}}+\frac{6F^{3}_{-\frac{1}{2},0,-\frac{1}{2}}(1-z)}{\epsilon^{2}(1-\epsilon)^{2}\Gamma(\lambda)^{2}}$
$\displaystyle\quad+\,\frac{4}{\epsilon^{3}(1-\epsilon)^{2}(1-\frac{3\epsilon}{2})\Gamma(\lambda)^{3}}\Big{(}\Big{(}1+\frac{\zeta(3)\epsilon^{3}}{2}\Big{)}(F^{4}_{0,-\frac{1}{2},0,-1}+F^{4}_{0,0,-\frac{1}{2},-1}-F^{4}_{-1,-\frac{1}{2},0,0})(z)$
(45)
$\displaystyle\qquad-\,\Big{(}1-\frac{\zeta(3)\epsilon^{3}}{2}\Big{)}(F^{4}_{-\frac{1}{2},0,0,-1}-\frac{1}{2}F^{4}_{-\frac{3}{2},0,0,0})(z)-(F^{4}_{1,-\frac{3}{2},-1,0}-\frac{1}{2}F^{4}_{\frac{3}{2},-\frac{3}{2},-\frac{3}{2},0})(1-z)\Big{)}.$

We need the expansions of the graphical functions $F^{3}_{a,b,c}(z)$,
$F^{4}_{a,b,c,d}(z)$, and $F^{6}_{a,b,c,d,e,f}(z)$ to orders $\epsilon^{2}$,
$\epsilon^{3}$, and $\epsilon^{1}$, respectively. This can be achieved by
parametric integration using (an optimized version of) HyperInt, see [29] and
Section 11.

The graphical function $F^{5}_{a,b,c,d,e}(z)$ needs to be expanded to order
$\epsilon^{1}$. With its ten edges it is too complex to use HyperInt directly.
We use completion (see Section 7) and swap the vertices
$0\leftrightarrow\infty$ and $1\leftrightarrow z$ to derive
$F^{5}_{a,b,c,d,e}(z)$ from the graphical function
$\widetilde{F}^{5}_{a,b,c,d,e}(z)$ in Figure 16 (where we have deleted edges
between external vertices). In $\widetilde{F}^{5}_{a,b,c,d,e}(z)$ we have
three edges of weights that vanish in the limit $\epsilon\to 0$. Because we
need $\widetilde{F}^{5}_{a,b,c,d,e}(z)$ only to order $\epsilon^{1}$ we can
expand the integrand in $\epsilon$, see Section 9.1. This allows us to express
$\widetilde{F}^{5}_{a,b,c,d,e}(z)$ in terms of $\widetilde{F}^{51}_{a,b,c}(z)$
and $\widetilde{F}^{52}_{a,b,c}(z)$. Both $\widetilde{F}^{51}_{a,b,c}(z)$ and
$\widetilde{F}^{52}_{a,b,c}(z)$ only have eight edges. The first two orders of
their graphical functions can easily be calculated with HyperInt.

##  13\. Limitations and future developments

At the time of writing, in the calculation of the eight loop $\phi^{4}$ beta
function results for 19000 of 233490 Feynman periods are missing. The
calculation of the six loop beta function in $\phi^{3}$ theory lacks 1100 of
40788 Feynman periods. Computations are hampered by limited access to
computing facilities. Nevertheless, the author expects that a six loop result
for the $\phi^{3}$ beta function will be available in the near future. Files
of intermediate results are in [34].

Difficulties in the calculation of renormalization functions are threefold:

  1. (1)

computing time,

  2. (2)

memory consumption,

  3. (3)

inaccessible Feynman periods.

Computing time is the least severe limitation. Still, a large number of
Feynman periods needs to be calculated and each calculation can be quite time
consuming. In most cases the majority of computing time is used for finding a
way how to do the calculation by using all implemented formulae for non-
constructible Feynman periods. One has to try all internal vertices for the
external vertex $z$ and for each choice of $z$ one needs to use a large set of
transformations in the attempt to calculate the graphical function. Once a
chain of transformations is found, the actual calculation typically is much
faster. An exception is the case when vast amounts of memory are used.

Memory consumption is typically rather mild (a few GB). An exception occurs if
the alphabet of GSVHs uses letters that lead to numbers which are not Euler
sums. This happens quite frequently at six loops in $\phi^{3}$ theory. A
prominent type of such numbers are iterated integrals [18] in the alphabet

$X=\Big{\\{}-1,0,\frac{1}{2},1,2\Big{\\}}.$

The alphabet $X$ is a transformed (real) subset of the alphabet of sixth roots
of unity $\exp(2k\pi\mathrm{i}/6)$, $k=0,1,\ldots,5$. Although
HyperlogProcedures can handle these numbers (using the ‘Schnetz basis’), the
${\mathbb{Q}}$-basis of these number is quite large at high weights and
transformations into the $f$-alphabet (that is used by HyperlogProcedures for
transformation into a ${\mathbb{Q}}$-basis [10, 11]) can be slow. Moreover,
taking limits of such GSVHs at $z=1$ and at $z=\infty$ can be time and memory
consuming.

Hitherto inaccessible Feynman periods can be calculated by using some hand-
crafted reduction chains, see Section 12, or by generating results for
graphical functions using parametric integration, see Section 11. Both
approaches are cumbersome and tedious.

Systematic approaches to increase the power of the graphical functions method
are the following subsections.

Moreover, the application of the graphical functions method to Gauge Theories
[24] seems possible and desirable.

###  13.1. Extra identities

The list of identities presented in this article is quite incomplete. Many
identities were found in the search for a solution of a particular Feynman
period. It seems possible that many important techniques for calculating
graphical functions are yet to be found.

###  13.2. Recursive use of known identities

In HyperlogProcedures [34] identities are used to simplify graphical
functions. In some cases, transformations into graphs of similar complexity
were accepted because the new graphs may be amenable to reductions. In the
Maple implementation HyperlogProcedures this immediately leads to very slow
searches for chains of transformations that solve a given graphical function.
The search for a reduction chain is independent of the actual calculation. It
can be programmed recursively in C$++$ with much better efficiency. In the
case of convergent graphical functions in six dimensions this was successfully
done by M. Borinsky in [6]. It became evident that reduction chains often use
detours over seemingly more complex graphical functions and are hence missed
by HyperlogProcedures. In particular, the IBP method (see Section 8.2) which
is very powerful in $\geq 6$ dimensions can (by its complexity) only be
handled in this approach.

The recursive calculation of graphical functions and Feynman periods in
integer dimensions [6] can be modified to compute the $\epsilon$-coefficients
of regularized integrals. It seems efficient to introduce two types of edge
weights: general weights $(n_{e}+a_{e}\epsilon)/\lambda$ with
$n_{e}\in{\mathbb{Z}}$, $a_{e}\in{\mathbb{R}}$, as in (27) and special weights
$(n_{e}-\frac{\epsilon}{2})/\lambda$ which, after external differentiation or
integration, may become weight 1 edges that can be reduced by the appending an
edge algorithm of Section 4.

Very desirable would also be a systematic implementation of rerouting, see
Section 10. This method for the subtraction of subdivergences (and the
reduction of the relevant orders in $\epsilon$) has proved very efficient.
Only the simplest cases are implemented in HyperlogProcedures.

###  13.3. Deep learning

Even in a C$++$ implementation the recursive use of known identities quickly
becomes too big for exhaustive results. One needs to introduce cutoffs in the
size of admissible graphs or needs to limit the set of identities which is
used. The setup seems quite accessible to methods of deep learning. The power
of deep learning in this context is not yet known.

## References

  * [1] P.A. Baikov, K.G. Chetyrkin and J.H. Kühn, Five-loop running of the QCD coupling constant, Phys. Rev. Lett. 118 8, 082002 (2017).
  * [2] P.A. Baikov, K.G. Chetyrkin, J.H. Kühn, J. Rittinger, Vector correlator in massless QCD at order $\mathcal{O}(\alpha_{s}^{4})$ and the QED $\beta$-function at five Loop, JHEP No. 7 017 (2012).
  * [3] M. Borinsky, Tropical monte carlo quadrature for Feynman integrals, to appear in Annales de l’Institut Henri Poincaré D, arXiv:2008.12310 [math-ph] (2020).
  * [4] M. Borinsky, J.A. Gracey, M.V. Kompaniets, O. Schnetz, Five loop renormalization of $\phi^{3}$ theory with applications to the Lee-Yang edge singularity and percolation theory, Phys. Rev. D 103, 116024 (2021).
  * [5] M. Borinsky, O. Schnetz, Graphical functions in even dimensions, Comm. in Number Theory and Physics 16, No. 3, 515-614 (2022).
  * [6] M. Borinsky, O. Schnetz, Recursive computation of Feynman periods, JHEP No. 08, 291 (2022).
  * [7] F.C.S. Brown, Single-valued multiple polylogarithms in one variable, C.R. Acad. Sci. Paris, Ser. I 338, 527-532 (2004).
  * [8] F.C.S. Brown, The massless higher-loop two-point function, Comm. Math. Phys. 287 925-958 (2009).
  * [9] F.C.S. Brown, On the periods of some Feynman integrals, arXiv:0910.0114 [math.AG] (2009).
  * [10] F.C.S. Brown, Mixed Tate Motives over ${\mathbb{Z}}$, Annals of Mathematics 175, No. 2, 949-976 (2012).
  * [11] F.C.S. Brown, On the decomposition of motivic multiple zeta values, in Galois-Teichmüller theory and arithmetic geometry, vol. 68 of Adv. Studies in Pure Math. (Tokyo), Math. Soc. Japan, 31-58 (2012).
  * [12] F.C.S. Brown, Feynman amplitudes, coaction principle, and cosmic Galois group, Comm. in Number Theory and Physics 11, No. 3, 453-555 (2017).
  * [13] F.C.S. Brown, Notes on motivic periods, Comm. in Number Theory and Physics 11, No. 3, 557-655 (2017).
  * [14] F.C.S. Brown, D. Kreimer Angles, Scales and Parametric Renormalization, Lett. Math. Phys. 103, 933–1007 (2013).
  * [15] F.C.S. Brown, O. Schnetz, A K3 in $\phi^{4}$, Duke Mathematical Journal, Vol. 161, No. 10, 1817-1862 (2012).
  * [16] S. Caron-Huot, L.J. Dixon, F. Dulat, M. von Hippel, A.J. McLeod, G. Papathanasiou, The Cosmic Galois Group and Extended Steinmann Relations for Planar $N=4$ SYM Amplitudes, J. High Energ. Phys. 61 (2019).
  * [17] F. Chavez, C. Duhr, Three-mass triangle integrals and single-valued polylogarithms, JHEP 1211, 114-144 (2012).
  * [18] K. Chen, Algebras of Iterated Path Integrals and Fundamental Groups, Transactions of the American Mathematical Society, Vol. 156, 359-379 (1971).
  * [19] A. Connes, D. Kreimer Hopf Algebras, Renormalization and Noncommutative Geometry, Comm. Math. Phys. 199, 203–242 (1998).
  * [20] J. Drummond, C. Duhr, P. Heslop, J. Pennington, V.A. Smirnov, Leading singularities and off-shell conformal integrals, JHEP No. 8, 133-190 (2013).
  * [21] M. Golz, Graphical functions in parametric space, Master thesis, Humboldt University Berlin, available at http://www2.mathematik.hu-berlin.de/$\sim$kreimer/wp-content/uploads/GolzMasterThesis.pdf (2015).
  * [22] M. Golz, E. Panzer, O. Schnetz, Graphical functions in parametric space, Lett. Math. Phys. 107, No. 6, 1177-1182 (2017).
  * [23] F. Herzog, B. Ruijl, T. Ueda, J.A.M. Vermaseren, A. Vogt, The five-loop beta function of Yang-Mills theory with fermions, JHEP No. 2, 90 (2017).
  * [24] J.C. Itzykson, J.B. Zuber, Quantum Field Theory, Mc-Graw-Hill, (1980).
  * [25] S. Jeffries, K.A. Yeats, A degree preserving delta wye transformation with applications to 6-regular graphs and Feynman periods, arXiv:2110.07764 [math.CO] (2021).
  * [26] M.V. Kompaniets, E. Panzer, Minimally subtracted six loop renormalization of $O(n)$-symmetric $\phi^{4}$ theory and critical exponents, Phys. Rev. D96, No. 3, 036016 (2017).
  * [27] M.V. Kompaniets, A. Pikelner, Critical exponents from five-loop scalar theory renormalization near six-dimensions, Phys. Lett. B 817 (2021).
  * [28] D. Kreimer, On Overlapping Divergences, Commun. Math. Phys. 204, 669-689 (1999).
  * [29] E. Panzer, Algorithms for the symbolic integration of hyperlogarithms with applications to Feynman integrals, Computer Physics Comm. 188, 148-166 (2015).
  * [30] E. Panzer, O. Schnetz, The Galois coaction on $\phi^{4}$ periods, Comm. in Number Theory and Physics 11, No. 3, 657-705 (2017).
  * [31] O. Schnetz, Graphical functions and single-valued multiple polylogarithms, Comm. in Number Theory and Physics 8, No. 4, 589–675 (2014).
  * [32] O. Schnetz, Quantum field theory over ${\mathbb{F}}_{q}$, Electron. J. Comb. 18N1:P102 (2011).
  * [33] O. Schnetz, Numbers and Functions in Quantum Field Theory, Phys. Rev. D 97, 085018 (2018).
  * [34] O. Schnetz, HyperlogProcedures, Version 0.6, Maple package available on the homepage of the author, https://www.math.fau.de/person/oliver-schnetz/ (2022).
  * [35] O. Schnetz, The Galois coaction on the electron anomalous magnetic moment, Comm. in Number Theory and Physics 12, No. 2, 335–354 (2018).
  * [36] O. Schnetz, Geometries in perturbative quantum field theory, Comm. in Number Theory and Physics 15, No. 4, 743 – 791 (2021).
  * [37] O. Schnetz, Generalized single-valued hyperlogarithms, arXiv:2111.11246 [hep-th], submitted to Comm. in Number Theory and Physics (2021).
  * [38] D. Zagier, The dilogarithm function in geometry and number theory, in: Number theory and related topics. Papers presented at the Ramanujan Colloquium, Bombay 1988, Studies in Mathematics No. 12, TIFR and Oxford University Press, pp. 231-249 (1989) and J. Math. Phys. Sci. No. 22, 131-145 (1988).
  * [39] J. Zinn-Justin, Quantum Field Theory and Critical Phenomena, International Series of Monographs on Physics 113, Clarendon Press, Oxford (2002).


# Two-species reaction-diffusion system in the presence of random velocity
fluctuations

Michal Hnatič1,2,3, Matej Kecer1 Tomáš Lučivjanský1

###### Abstract

We study random velocity effects on a two-species reaction-diffusion system
consisting of three reaction processes
$\textit{A}+\textit{A}\rightarrow(\emptyset,A),$
$\textit{A}+\textit{B}\rightarrow\textit{A}$. Using the field-theoretic
perturbative renormalization group we analyze this system in the vicinity of
its upper critical dimension $d_{c}=2$. Velocity ensemble is generated by
means of stochastic Navier-Stokes equations. In particular, we investigate the
effect of thermal fluctuations on reaction kinetics. The overall analysis is
performed to the one-loop approximation and possible macroscopic regimes are
identified.

1 Institute of Physics, Faculty of Science, P. J. Šafárik University, Park
Angelinum 9, 040 01 Košice, Slovakia
2 Institute of Experimental Physics, Slovak Academy of Sciences, Watsonova 47,
040 01 Košice, Slovakia
3 Joint Institute for Nuclear Research, 141980 Dubna, Russia

##  1 Introduction

Diffusion-limited reactions constitute prominent models in non-linear
statistical physics [1]. Theoretical study of such systems attracted a lot of
attention in the past [2, 3]. A straightforward approach to theoretical
analysis of such systems is based on kinetic rate equations, which might be
regarded as a simple mean-field-like approximation [3, 4]. However, reaction
systems are known to exhibit non-trivial behavior, especially in low space
dimensions [5], where density fluctuations become especially pronounced. There
the kinetic rate equations approach is not adequate and more sophisticated
approaches are called for. In this paper, we study a multi-species reaction-
diffusion system [6, 7, 8, 9], which consists of the following three reaction
processes

$\displaystyle A+A$ $\displaystyle\rightarrow\begin{cases}A\hskip
14.22636pt\text{coalescence },\\\ \emptyset\hskip
17.07182pt\text{annihilation},\end{cases}$ (1) $\displaystyle A+B$
$\displaystyle\rightarrow A\hskip 25.6073pt\text{trapping},$

where coalescence process occurs with probability $p\,(0\leq p\leq 1)$, and
annihilation process with a complementary probability $1-p$. The model becomes
even more intricate when additional effects are taken into account. Their
investigation is especially important, as they naturally arise in many
practical circumstances. For instance, the majority of chemical reactions in
typical experimental settings occur in some fluid environment. Various aspects
of such a problem have already been studied recently [7, 8, 9]. Here, our aim
is to investigate the influence of thermal fluctuations of a surrounding
environment on the kinetics of reaction-scheme (1). We model the environment
as a fluid at a constant temperature using a well-known approach based on the
stochastic Navier-Stokes equation [10, 11].

A powerful tool for analyzing the asymptotic behavior of stochastic systems is
provided by the renormalization group (RG) method [12, 13]. It allows us to
determine the long-time and large-scale – or infrared (IR) – asymptotic
regimes of the system and also is a very efficient tool for the calculation of
various universal physical quantities, e.g. critical exponents. The aim of
this paper is to address the possible IR behavior of the reaction-diffusion
process (1) under the influence of advecting velocity fluctuations and to
determine their IR regimes.

The paper is organized as follows. In Sec. 2 we give a field-theoretic
formulation of the model and specify the main ingredients of the perturbation
theory. Sec. 3 is devoted to the analysis of ultraviolet divergences and
renormalization of the model in one-loop order of perturbation scheme. The
analysis of fixed points (FP) and their regions of stability are discussed in
Sec. 4. Conclusions are drawn in Sec. 5.

##  2 Field-theoretic formulation of the model

The field theory for the reaction-diffusion system described by the scheme (1)
can be constructed from the master equation by means of Doi-Peliti formalism
[4, 14, 15, 16]. For brevity, we omit the derivation as it can be easily found
elsewhere (see e.g., [3]). We start our analysis with the field-theoretic
action for the reaction scheme (1) augmented with diffusion processes

$\displaystyle{\mathcal{S}}_{r}[\Psi]$
$\displaystyle=\psi^{\dagger}_{A}(-\partial_{t}+\nu_{0}u_{A0}\partial^{2})\psi_{A}^{~{}}+\psi^{\dagger}_{B}(-\partial_{t}+\nu_{0}u_{B0}\partial^{2})\psi_{B}^{~{}}-\nu_{0}u_{A0}\lambda_{0}\psi_{A}^{\dagger}\psi_{A}^{2}$
$\displaystyle-\nu_{0}u_{A0}\lambda_{0}\psi_{A}^{\dagger
2}\psi_{A}^{2}-\lambda_{0}^{\prime}Q\nu_{0}u_{A0}\psi_{B}^{\dagger}\psi_{A}^{~{}}\psi_{B}^{~{}}-\nu_{0}u_{A0}\lambda_{0}^{\prime}\psi_{A}^{\dagger}\psi_{B}^{\dagger}\psi_{A}^{~{}}\psi_{B}^{~{}},$
(2)

where
$\Psi\equiv\\{\psi_{A},\psi_{A}^{\dagger},\psi_{B},\psi_{B}^{\dagger}\\}$ are
bosonic-like coherent fields arising in taking a continuum limit in the Doi-
Peliti approach [3], $\partial^{2}=\partial_{i}\partial_{i}$ denotes Laplace
operator in $d$-dimensions and diffusion parameters are expressed through
respective Prandtl numbers $u_{A0}$, $u_{B0}$ and viscosity $\nu_{0}$ (see
below Eq. (7)). The parameters $\lambda_{0},\lambda^{\prime}_{0}$ denote
reaction constants, and parameter $Q=1/(2-p)$ is related to the probability of
whether annihilation or coagulation process takes place. In this work, we
employ the RG method, which introduces two different kinds of variables - bare
(unrenormalized) quantities and their renormalized counterparts. Therefore we
denote the former ones with the subscript “0”, whereas the latter will be
written without the subscript “0”.

Reaction process (1) is an example of a genuine non-equilibrium system and,
therefore, we have to specify its initial conditions. We choose them in the
following form

${\mathcal{S}}_{init}[\Psi]=(a_{0}\psi_{A}^{\dagger}+b_{0}\psi_{B}^{\dagger})\delta(t),$
(3)

where $a_{0},b_{0}$ are appropriately rescaled initial average densities [4,
7].

In writing actions (2) and (3) we have employed a condensed notation, in which
integrations over space and time variables in the expressions for action
functionals are implied. For instance, the first term in the action (2)
corresponds to

$-\psi^{\dagger}_{A}\partial_{t}\psi_{A}=-\int\mathrm{d}x\,\psi^{\dagger}_{A}(x)\partial_{t}\psi_{A}(x),$
(4)

where we have written coordinates compactly as $x=(t,\bm{x})$ and integration
measure as $\mathrm{d}x=\mathrm{d}t\mathrm{d}^{d}x$.

The aim of this paper is to study the case where chemical particles are
advected within the fluid environment with random fluctuations. We introduce
advection processes into the formalism by the inclusion of convective
derivative [17]. This corresponds to the replacement of the time derivative as
follows

$\partial_{t}\rightarrow\partial_{t}+\bm{v}\cdot{\bm{\nabla}}=\partial_{t}+v_{j}\partial_{j},$
(5)

where summation over repeated indices is implied in the last term. Let us
stress that the advection for both particle types is considered to be passive,
i.e., the velocity field itself is not affected by the particles or reactions
processes, respectively. Corresponding advective terms to the action (2) take
the form

${\mathcal{S}}_{adv}[\Psi,\bm{v}]=-\psi_{A}^{\dagger}v_{j}\partial_{j}\psi_{A}^{~{}}-\psi_{B}^{\dagger}v_{j}\partial_{j}\psi_{B}^{~{}}.$
(6)

To finalize the model construction we need to specify velocity field $\bm{v}$.
Here, we assume that velocity field $\bm{v}(t,\bm{x})$ is a random variable
with zero mean, whose dynamics is governed by stochastic Navier-Stokes
equation [10, 18].

$\partial_{t}v_{i}+(v_{j}\partial_{j})v_{i}=\nu_{0}\partial^{2}v_{i}-\partial_{i}P+f_{i},$
(7)

where $P=P(x)$ is the pressure field, and $f_{i}=f_{i}(x)$ denotes $i$-th
component of an external random force $\bm{f}$. Following earlier works [10,
18, 19] we assume the force $\bm{f}$ is a random Gaussian variable with zero
mean and correlation function of the prescribed form

$\langle
f_{i}(t,\bm{x})f_{j}(0,\bm{0})\rangle=\int\frac{\mathrm{d}^{d}k}{(2\pi)^{d}}D_{ij}(t,\bm{k})\mathrm{e}^{i\bm{k}\cdot\bm{x}}.$
(8)

We consider the case of an incompressible fluid, which implies transversality
of the field $\bm{v}$ ($\partial_{i}v_{i}=0$). Using this condition it is
possible to express pressure in terms of velocity field [18]. This is
equivalent to work in transversal space by taking the following replacement
for velocity field $\bm{v}$ in the momentum representation

$v_{i}(\bm{k})\rightarrow P_{ij}(\bm{k})v_{j}(\bm{k}),$ (9)

where $P_{ij}(\bm{k})=\delta_{ij}-k_{i}k_{j}/k^{2}$ with $(k=|\bm{k}|)$ is
transverse projection operator.

The incompressibility condition implies that the kernel $D_{ij}$ in the
momentum representation is proportional to transverse projector
$P_{ij}(\bm{k})$. In fact, it can be readily shown that for incompressible
medium $D_{ij}\sim\delta_{ij}$ is sufficient. However, we follow the
traditional notation in previous works and keep $P_{ij}$ in the expression for
kernel $D_{ij}$. Using a specific choice for the momentum dependence of
$D_{ij}$ term it is possible to generate fluctuations of the velocity field
near thermal equilibrium.

These considerations finally lead to

$D_{ij}(t,\bm{k})=\delta(t)D_{0}k^{2}P_{ij}(\bm{k}),$ (10)

where $\delta=\delta(t)$ is Dirac delta function. It can be shown that delta
correlations in time of the kernel $D_{ij}$ ensures that the present model
possesses the Galilean symmetry [11, 20].

In hindsight, this particular form (10) is convenient for the application of
RG method, because both velocity fluctuations and reaction processes of the
original reaction-diffusion system become simultaneously marginal in the
critical space dimension $d=d_{c}=2$. The stochastic problem (7)-(10) can be
recast into a field theory with the doubled set of fields
$\Phi=\\{\bm{v},\bm{v}^{\prime}\\}$ described by the De Dominicis-Janssen
action functional [12, 11],

${\mathcal{S}}_{v}[\Phi]=\frac{1}{2}v_{i}^{\prime}D_{ij}v_{j}^{\prime}+v_{i}^{\prime}\left(-\partial_{t}v_{i}-v_{j}\partial_{j}v_{i}+\nu_{0}\partial^{2}v_{i}\right),$
(11)

where response field $v_{i}^{\prime}$ is incompressible, and again condensed
notation in the sense of Eq. (4) is assumed. Let us note that quadratic term
in the response field $\bm{v}^{\prime}$ in the action (11) actually stands for

$v^{\prime}_{i}D_{ij}{v^{\prime}}_{j}=\int\mathrm{d}x\int\mathrm{d}x^{\prime}v^{\prime}_{i}(x)D_{ij}(x-x^{\prime})v^{\prime}_{k}(x^{\prime}),$
(12)

where $D_{ij}$ corresponds to the inverse Fourier transform of the kernel
(10).

The sum of action functionals (2), (3), (6), and (11), respectively, then
gives us a final field-theoretic action

${\mathcal{S}}={\mathcal{S}}_{r}+{\mathcal{S}}_{v}+{\mathcal{S}}_{adv}+{\mathcal{S}}_{init}.$
(13)

Expectation values of some physical observable $A=A(t,\bm{x})$ can be, in
principle, calculated as a functional integral [4, 12]

$\langle
A(t,\bm{x})\rangle=\mathcal{N}^{-1}\int\mathcal{D}\Psi\mathcal{D}\Phi\,A(t,\bm{x})\mathrm{e}^{S},$
(14)

where $\mathcal{N}$ is a normalization constant.

In what follows we analyze field-theoretic action (13) using the field-
theoretic renormalization group. This technique was employed in the past on
similar problems as well [3, 21, 22, 23, 24, 25, 26]. We apply it here in a
perturbation setting, which is based on expressing Green functions as a series
in coupling constants of a theory. The perturbation theory of the model is
then constructed using well-known Feynman diagrammatic rules [3, 12, 13]. The
part of the action (13) quadratic in fields determines the bare propagators,
which in frequency-momentum representation take form

$\displaystyle\langle\psi_{A}^{~{}}\psi_{A}^{\dagger}\rangle_{0}$
$\displaystyle=\frac{1}{-i\omega+\nu_{0}u_{A0}k^{2}},$
$\displaystyle\langle\psi_{B}^{~{}}\psi_{B}^{\dagger}\rangle_{0}$
$\displaystyle=\frac{1}{-i\omega+\nu_{0}u_{B0}k^{2}},$ (15)
$\displaystyle\langle v_{i}v_{j}\rangle_{0}$
$\displaystyle=\frac{D_{0}k^{2}P_{ij}(\bm{k})}{\omega^{2}+\nu_{0}^{2}k^{4}},$
$\displaystyle\langle v_{i}v^{\prime}_{j}\rangle_{0}$
$\displaystyle=\frac{P_{ij}(\bm{k})}{-i\omega+\nu_{0}k^{2}}.$ (16)

The nonlinear terms determine interaction vertices with associated vertex
factors [12]. They can be calculated with the help of the formula

$V_{N}(x_{1},\ldots,x_{N};\varphi)=\frac{\delta^{N}{\mathcal{S}}_{\text{int}}}{\delta\varphi(x_{1})\ldots\delta\varphi(x_{N})},\quad\varphi\in\\{\psi_{A},\psi_{A}^{\dagger},\psi_{B},\psi_{B}^{\dagger},\bm{v},\bm{v}^{\prime}\\},$
(17)

where ${\mathcal{S}}_{\text{int}}$ corresponds to the non-linear terms of the
action (13). In a straightforward manner, we get the following bare vertices
without an inclusion of velocity field

$\displaystyle V_{\psi_{A}^{\dagger}\psi_{A}^{~{}}\psi_{A}^{~{}}}$
$\displaystyle=-2\lambda_{0}\nu_{0}u_{A0},$ $\displaystyle
V_{\psi_{B}^{\dagger}\psi_{B}^{~{}}\psi_{A}^{~{}}}$
$\displaystyle=-\lambda^{\prime}_{0}\nu_{0}u_{A0}Q,$ $\displaystyle
V_{\psi_{A}^{\dagger}\psi_{A}^{\dagger}\psi_{A}^{~{}}\psi_{A}^{~{}}}$
$\displaystyle=-4\lambda_{0}\nu_{0}u_{A0},$ $\displaystyle
V_{\psi_{A}^{\dagger}\psi_{B}^{\dagger}\psi_{A}^{~{}}\psi_{B}^{~{}}}$
$\displaystyle=-\lambda^{\prime}_{0}\nu_{0}u_{A0}.$ (18)

On the other hand, there are three additional vertices that include the
velocity field

$V_{\psi_{A}^{\dagger}(\bm{k})\psi_{A}^{~{}}v_{j}}=ik_{j},\quad
V_{\psi_{B}^{\dagger}(\bm{k})\psi_{B}^{~{}}v_{j}}=ik_{j},\quad
V_{v^{\prime}_{i}(\bm{k})v_{l}v_{j}}=i(k_{l}\delta_{ij}+k_{j}\delta_{il}).$
(19)

First, two describe advection processes and the latter vertex is responsible
for interactions between velocity fluctuations. Also, we have explicitly
written, the momentum of which field enters a given interaction vertex. For
instance, in expression for the vertex factor
$V_{v^{\prime}_{i}(\bm{k})v_{l}v_{j}}$ the momentum $k_{j}$ is carried by the
response field $v_{i}^{\prime}$ [11, 12].

##  3 Renormalization of the model

The analysis of UV divergences starts with a determination of the canonical
dimensions for model parameters. In dynamical models, there are two
independent scales that need to be considered [3, 12]. These are frequency and
momentum scales (time and length). Then any quantity $F$ is characterized with
both frequency dimension $d_{F}^{\omega}$ and a momentum dimension
$d_{F}^{k}$, respectively. Canonical dimensions are determined from
normalization conditions

$d_{k}^{k}=-d_{x}^{k}=1,\ d_{\omega}^{\omega}=-d_{t}^{\omega}=1,\
d_{k}^{\omega}=d_{\omega}^{k}=0,$ (20)

and the fact that the action functional has to be a dimensionless quantity
[12]. The total canonical dimension of any $F$ is then given as
$d_{F}=d_{F}^{k}+2d_{F}^{\omega}$ (because of
$\partial_{t}\propto\mathbf{\partial}^{2}$ proportionality in quadratic part
of the action functional). Canonical dimensions of all the fields and
parameters of model (13) are listed in Tab. 1.

$F$ |  $\psi_{A}$, $\psi_{B}$ |  $\psi^{\dagger}_{A}$, $\psi^{\dagger}_{B}$ | $\bm{v}$ | $\bm{v}^{\prime}$ |  $\lambda_{0}$, $\lambda^{\prime}_{0}$, $g_{0}$ | $Q$ |  $a_{0}$, $b_{0}$ | $\nu_{0}$ | $D_{0}$ |  $u_{A0}$, $u_{B0}$
---|---|---|---|---|---|---|---|---|---|---
$d_{F}^{k}$ | $d$ | $0$ | $-1$ | $d+1$ | $2-d$ | $0$ | $d$ | $-2$ | $-d-4$ | $0$
$d_{F}^{\omega}$ | $0$ | $0$ | $1$ | $-1$ | $0$ | $0$ | $0$ | $1$ | $3$ | $0$
$d_{F}$ | $d$ | $0$ | $1$ | $d-1$ | $2-d$ | $0$ | $d$ | $0$ | $2-d$ | $0$
Table 1: Canonical dimensions of fields and parameters.

There are altogether five charges (coupling constants) of the theory

$g_{0}=\frac{D_{0}}{\nu_{0}^{3}},\ u_{A0},\ u_{B0},\ \lambda_{0},\
\lambda^{\prime}_{0}.$ (21)

In the space dimension $d=2$, all of these charges become simultaneously
dimensionless and the model becomes logarithmic. Therefore this dimension is
identified as an upper critical dimension $d_{c}$ of the model. In dimensional
regularisation, the UV divergences manifest themselves as poles in expansion
parameter $\varepsilon=2-d$, whereas the IR divergences are regulated by the
sharp cutoff at $k=m$, which is an analog of the inverse of integral
turbulence scale $L=1/m$. Let us note that the latter divergences do not
affect renormalization constants [12].

Probably the most economical way to renormalize the translationally invariant
model is through the renormalization of its one-particle irreducible (1PI)
Green functions. This is a restricted class of Feynman diagrams that consists
of such diagrams that remain connected even after one internal line is cut off
[3, 12]. An arbitrary one-particle irreducible (1PI) Green’s function will be
denoted as $\Gamma_{\\{\varphi\\}}=\langle\varphi\ldots\varphi\rangle_{1PI}$,
where $\varphi\in\Psi\cup\Phi$ denotes an arbitrary field from the full set of
fields of the model (13). Its total canonical dimension is given by a general
formula [12, 13]

$d_{\Gamma}=d+2-\sum_{\varphi}N_{\varphi}d_{\varphi},$ (22)

where the sum runs through all the types of fields $\varphi$, $N_{\varphi}$
denotes the number of times the given field appears in the particular 1PI
function and $d_{\varphi}$ is its canonical dimension. Following the standard
approach [12] the task is to identify superficial divergences in 1PI functions
and construct renormalized action, in which introduced additional counter-
terms ensure the removal of these divergences in the given order of
perturbation theory.

The UV divergences, which require further treatment, are identified with those
1PI Green functions, which possess a non-negative formal index of divergence
$\delta_{\Gamma}=d_{\Gamma}|_{\varepsilon=0}$. However, for the present case,
this statement is to be adjusted based on the following considerations. First,
the 1PI functions not involving any of the response functions
$\psi_{B}^{\dagger},\psi_{A}^{\dagger},\bm{v}^{\prime}$ as external fields
vanish as they necessarily contain closed cycles of causal propagators [12].
Since vertex factor $V_{\bm{v}^{\prime}\bm{v}\bm{v}}$ is proportional to the
momentum carried by field $\bm{v}^{\prime}$ (see the corresponding expression
in (19)), every instance of $\bm{v}^{\prime}$ appearing as external field
lowers the overall index of divergence. Thus the real index of divergence is
defined as

$\tilde{\delta}_{\Gamma}=\delta_{\Gamma}-N_{\bm{v}^{\prime}}.$ (23)

Second, the number of counter-terms is further reduced because of the
invariance property of generating functional of model (13) with respect to
Galilean transformations. This symmetry implies that the function
$\langle\bm{v}^{\prime}\bm{v}\bm{v}\rangle_{1PI}$ does not diverge (for
further discussions on the subject see e.g. [11, 12, 27]). Taking these into
account along with available diagrammatic elements and transversality of the
velocity field we can identify the following irreducible functions with
superficial UV divergences

$\displaystyle\langle\bm{v}^{\prime}\bm{v}^{\prime}\rangle_{1PI},$
$\displaystyle\langle\psi_{A}^{\dagger}\psi_{A}^{~{}}\psi_{A}^{~{}}\rangle_{1PI},$
$\displaystyle\langle\bm{v}^{\prime}\bm{v}\rangle_{1PI},$
$\displaystyle\langle\psi_{B}^{\dagger}\psi_{B}^{~{}}\psi_{A}^{~{}}\rangle_{1PI},$
$\displaystyle\langle\psi_{A}^{\dagger}\psi_{A}^{~{}}\rangle_{1PI},$
$\displaystyle\langle\psi_{A}^{\dagger}\psi_{A}^{\dagger}\psi_{A}^{~{}}\psi_{A}^{~{}}\rangle_{1PI},$
$\displaystyle\langle\psi_{B}^{\dagger}\psi_{B}^{~{}}\rangle_{1PI},$
$\displaystyle\langle\psi_{B}^{\dagger}\psi_{A}^{\dagger}\psi_{B}^{~{}}\psi_{A}^{~{}}\rangle_{1PI}$
. (24)

All of these have the form that is already present in the bare action
functional. This implies that the model is multiplicatively renormalizable.

The total renormalized action takes the form

$\displaystyle S_{R}$
$\displaystyle=\psi^{\dagger}_{A}\left(-\partial_{t}+Z_{1}u_{A}\nu\partial^{2}\right)\psi_{A}^{~{}}+\psi^{\dagger}_{B}\left(-\partial_{t}+Z_{2}u_{B}\nu\partial^{2}\right)\psi_{B}^{~{}}+\frac{v_{i}^{\prime}\mu^{\epsilon}Z_{3}D_{ij}v_{j}^{\prime}}{2}$
$\displaystyle+v_{i}^{\prime}(-\partial_{t}+Z_{4}\nu\partial^{2})v_{i}-u_{A}\nu\lambda\mu^{\epsilon}Z_{5}\left[\psi_{A}^{\dagger}+\psi_{A}^{\dagger
2}\right]\psi_{A}^{2}-\lambda^{\prime}u_{A}\nu\mu^{\epsilon}Z_{6}\left[Q\psi_{A}^{~{}}+\psi_{A}^{\dagger}\psi_{A}^{~{}}\right]\psi_{B}^{\dagger}\psi_{B}^{~{}}$
$\displaystyle-
v_{i}^{\prime}(\bm{v}\cdot{\bm{\nabla}})v_{i}-\left[\psi_{A}^{\dagger}(\bm{v}\cdot{\bm{\nabla}})\psi_{A}^{~{}}+\psi_{B}^{\dagger}(\bm{v}\cdot{\bm{\nabla}})\psi_{B}\right]+\delta(t)\left(\psi_{A}^{\dagger}\
a_{0}+\psi_{B}^{\dagger}\ b_{0}\right),$ (25)

and was obtained from the bare action (13) by introducing the following
renormalization of fields and parameters of the model

$\displaystyle\varphi$ $\displaystyle\rightarrow Z_{\varphi}\varphi,$
$\displaystyle u_{A0}$ $\displaystyle\rightarrow Z_{u_{A}}u_{A},$
$\displaystyle u_{B0}$ $\displaystyle\rightarrow Z_{u_{B}}u_{B},$
$\displaystyle\nu_{0}$ $\displaystyle\rightarrow Z_{\nu}\nu,$ $\displaystyle
g_{0}$ $\displaystyle\rightarrow\mu^{\varepsilon}Z_{g}g,$
$\displaystyle\lambda_{0}$
$\displaystyle\rightarrow\mu^{\varepsilon}Z_{\lambda}\lambda,$
$\displaystyle\lambda^{\prime}_{0}$
$\displaystyle\rightarrow\mu^{\varepsilon}Z_{\lambda^{\prime}}\lambda^{\prime}.$
(26)

Here,
$\varphi\in\\{\psi_{A},\psi_{A}^{\dagger},\psi_{B},\psi_{B}^{\dagger},\bm{v},\bm{v}^{\prime}\\}$,
$\mu$ is an arbitrary momentum scale and $Z_{F}$ denotes the corresponding
renormalization constant.

By direct inspection, we get relations between renormalization constants in
the renormalized action (25) and RG constants (26)

$\displaystyle Z_{u_{A}}$ $\displaystyle=Z_{1}Z_{4}^{-1},$ $\displaystyle
Z_{u_{B}}$ $\displaystyle=Z_{2}Z_{4}^{-1},$ $\displaystyle Z_{\lambda}$
$\displaystyle=Z_{5}Z_{1}^{-1},$ $\displaystyle Z_{\lambda}^{\prime}$
$\displaystyle=Z_{6}Z_{1}^{-1},$ $\displaystyle Z_{g}$
$\displaystyle=Z_{3}Z_{4}^{-3},$ $\displaystyle Z_{\nu}$
$\displaystyle=Z_{4},$ $\displaystyle Z_{\Psi}$ $\displaystyle=Z_{Q}=1.$ (27)

The explicit form of RG constants $Z_{1}-Z_{6}$ is calculated from the one-
loop 1PI Feynman diagrams using dimensional regularisation and a minimal
subtraction scheme. The final expressions read

$\displaystyle Z_{1}=1-\frac{\hat{g}}{4u_{A}(u_{A}+1)\varepsilon},$
$\displaystyle Z_{2}$
$\displaystyle=1-\frac{\hat{g}}{4u_{B}(u_{B}+1)\varepsilon},$ $\displaystyle
Z_{3}$ $\displaystyle=Z_{4}=1-\frac{\hat{g}}{16\varepsilon},$ $\displaystyle
Z_{5}=1+\frac{\hat{\lambda}}{\varepsilon},$ $\displaystyle Z_{6}$
$\displaystyle=1+\frac{\hat{\lambda}^{\prime}u_{A}}{(u_{A}+u_{B})\varepsilon},$
(28)

where $\hat{F}\equiv FS_{d}/(2\pi)^{d}$, $S_{d}=2\pi^{d/2}/\Gamma(d/2)$ is the
area of unit $d$-dimensional sphere, and $\Gamma(x)$ is Euler’s gamma
function.

##  4 RG functions and scaling regimes

Once the calculation of RG constants is successfully accomplished, it is
possible to analyze the asymptotic behavior of the system. A fundamental
equation that governs the behavior of renormalized Green functions is
expressed with the help of the RG operator, which can be expressed in the
present case as

$D_{RG}=\mu\partial_{\mu}+\sum_{e}\beta_{e}\partial_{e}-\gamma_{\nu}\nu\partial_{\nu},$
(29)

where the given sum runs through all charges of the theory
$e=\\{g,u_{A},u_{B},\lambda^{\prime},\lambda\\}$. The coefficient functions
are defined as

$\beta_{e}=\mu\frac{\partial
e}{\partial\mu}\biggr{|}_{0},\quad\gamma_{F}=\frac{\partial\ln
Z_{F}}{\partial\ln\mu}\biggr{|}_{0},$ (30)

for any parameter, $F$ and $|_{0}$ means that bare parameters are held
constant during evaluation. For the model (13), we have altogether five beta
functions

$\displaystyle\beta_{g}=$ $\displaystyle-g(\varepsilon+\gamma_{g}),$
$\displaystyle\beta_{u_{A}}$ $\displaystyle=-u_{A}\gamma_{u_{A}},$
$\displaystyle\beta_{u_{B}}$ $\displaystyle=-u_{B}\gamma_{u_{B}},$
$\displaystyle\beta_{\lambda}=$
$\displaystyle-\lambda(\varepsilon+\gamma_{\lambda}),$
$\displaystyle\beta_{\lambda^{\prime}}$
$\displaystyle=-\lambda^{\prime}(\varepsilon+\gamma_{\lambda^{\prime}}),$ (31)

with corresponding anomalous dimensions (by definition Eq. (30))

$\displaystyle\gamma_{g}=$ $\displaystyle-\frac{\hat{g}}{8},$
$\displaystyle\gamma_{u_{i}}$
$\displaystyle=\hat{g}\biggl{(}\frac{1}{4u_{i}(1+u_{i})}-\frac{1}{16}\biggr{)};i\in\\{A,B\\},$
$\displaystyle\gamma_{\lambda}=$
$\displaystyle-\hat{\lambda}-\frac{\hat{g}}{4u_{A}(1+u_{A})},$
$\displaystyle\gamma_{\lambda^{\prime}}$
$\displaystyle=-\hat{\lambda^{\prime}}\frac{u_{A}}{u_{A}+u_{B}}-\frac{\hat{g}}{4u_{A}(1+u_{A})},$
(32)

where the higher order corrections $\hat{g}^{2},\hat{\lambda}^{2}$ are
neglected in the one-loop approximation.

The long-time asymptotic behavior of the model is governed by the IR stable
fixed points (FP) [12, 13] of beta functions. These are such points
$e^{*}=(g^{*},u_{A}^{*},u_{B}^{*},\lambda^{*},\lambda^{\prime*})$ in coupling
constant space that satisfy

$\beta_{g}(e^{*})=\beta_{u_{A}}(e^{*})=\beta_{u_{B}}(e^{*})=\beta_{\lambda}(e^{*})=\beta_{\lambda^{\prime}}(e^{*})=0.$
(33)

IR stability is determined by the eigenvalues of the matrix of the first
derivatives

$\Omega_{ij}=\frac{\partial\beta_{i}}{\partial e_{j}}\biggl{|}_{e^{*}},$ (34)

where index $i$ and charge $e_{j}$ belong to the set
$\\{g,u_{A},u_{B},\lambda^{\prime},\lambda\\}$. For IR stable regimes the
eigenvalues of the matrix (34) have to have positive real parts. We have found
eight FPs, however, only two of them are IR stable (see Tab. 2). These are

  1. 1.

Gaussian fixed point (FP1): $g^{*}=0$, $u_{A}^{*}=$ arbitrary, $u_{B}^{*}=$
arbitrary, $\lambda^{*}=0$, $\lambda^{\prime*}=0$. IR stable for
$\varepsilon<0$.

  2. 2.

Thermal fixed point (FP8): $g^{*}=8\varepsilon$,
$u_{A}^{*}=u_{B}^{*}=(\sqrt{17}-1)/2$, $\lambda^{*}=\varepsilon/2$,
$\lambda^{\prime*}=\varepsilon$. IR stable for $\varepsilon>0$.

Let us note that in non-trivial (thermal) FP both velocity fluctuations and
reaction interactions are simultaneously IR-relevant. RG predicts also an FP
for which only reaction processes are relevant (FP4). However, even though it
would have been stable without the velocity field [16], it can never be truly
IR stable in the presence of thermal fluctuations, which are inevitable in
practice. A similar conclusion was obtained in the past for different
reaction-diffusion model [24]. On the borderline of the two regimes, i.e. for
the case $\varepsilon=0$, couplings of the theory become marginally
irrelevant, and logarithmic corrections are expected to appear in expressions
for Green’s functions. Based on the standard analysis [12] we predict that
these corrections will be different from the ones realized if the velocity
field was not present (and FP4 would be stable). Therefore the behavior in
two-dimensional systems is also expected to be affected by the presence of
velocity field fluctuations. The proof of this statement is deferred to future
work.

Fixed point | $\hat{g}^{*}$ | $u_{A}^{*},u_{B}^{*}$ | $\hat{\lambda}^{*}$ | $\hat{\lambda}^{\prime*}$ | Stability region
---|---|---|---|---|---
FP1 | $0$ | arbitrary | $0$ | $0$ | $\varepsilon<0$
FP2 | $0$ | arbitrary | $0$ | $\varepsilon(u_{A}+u_{B})/u_{A}$ | unstable
FP3 | $0$ | arbitrary | $\varepsilon$ | $0$ | unstable
FP4 | $0$ | arbitrary | $\varepsilon$ | $\varepsilon(u_{A}+u_{B})/u_{A}$ | unstable
FP5 | $8\varepsilon$ | $(\sqrt{17}-1)/2$ | $0$ | $0$ | unstable
FP6 | $8\varepsilon$ | $(\sqrt{17}-1)/2$ | $0$ | $\varepsilon$ | unstable
FP7 | $8\varepsilon$ | $(\sqrt{17}-1)/2$ | $\varepsilon/2$ | $0$ | unstable
FP8 | $8\varepsilon$ | $(\sqrt{17}-1)/2$ | $\varepsilon/2$ | $\varepsilon$ | $\varepsilon>0$
Table 2: RG fixed points and their region of stability.

##  5 Conclusion

We have investigated the influence of thermal fluctuations on a reaction-
diffusion system with reactions
$\textit{A}+\textit{A}\rightarrow(\emptyset,A),$
$\textit{A}+\textit{B}\rightarrow\textit{A}$. Using the field-theoretic
formulation of the model we have analyzed possible macroscopic behavior
utilizing the renormalization group approach. In particular, we have
renormalized the model to the one-loop order of the perturbation scheme. The
RG analysis revealed the existence of two IR-stable FPs which govern the long-
time behavior of the system.

## Conflict of Interest

The authors declare that they have no conflicts of interest.

## Acknowledgment

The work was supported by VEGA grant No. 1/0535/21 of the Ministry of
Education, Science, Research and Sport of the Slovak Republic.

Conflict of Interest: The authors declare that they have no conflicts of
interest.

## References

  * [1] P. L. Krapivsky, S. Redner, E. Ben-Naim, A kinetic view of statistical physics, Cambridge University Press (2010).
  * [2] G. Ódor, Rev. Mod. Phys. 76, 663 (2004).
  * [3] U. C. Täuber, Critical Dynamics: A Field Theory Approach to Equilibrium and Non-Equilibrium Scaling Behavior, Cambridge University Press, Cambridge (2014).
  * [4] U. C. Täuber, M. Howard, B. P. Vollmayr-Lee, J. Phys. A: Math. Gen. 38, R79 (2005).
  * [5] A. A. Ovchinnikov, S. F. Timashev, A. A. Belyy, Kinetics of Diffusion controlled chemical processes, Nova Science Pub Incorporated (1989).
  * [6] R. Rajesh and O. Zaboronski, Phys. Rev. E 70, 036111 (2004).
  * [7] B. Vollmayr-Lee, J. Hanson, R. S. McIsaac, J. D. Hellerick, J. Phys. A: Math. Theor. 51, 034002 (2018).
  * [8] J. D. Hellerick, R. C. Rhoades, B. P. Vollmayr-Lee, Phys. Rev. E 101, 042112 (2020).
  * [9] D. Shapoval, V. Blavatska, M. Dudka, J. Phys. A: Math. Theor. 55, 455002 (2022).
  * [10] D. Forster, D. R. Nelson, M. J. Stephen, Phys. Rev. A 16, 732 (1977).
  * [11] L. Ts. Adzhemyan, N. V. Antonov, A. N. Vasiliev, Field theoretic renormalization group in fully developed turbulence, CRC Press (1999).
  * [12] A. N. Vasil’ev, The Field Theoretic Renormalization Group in Critical Behavior Theory and Stochastic Dynamics, [in Russian], PIYaF, St. Petersburg (1998); English trans., Chapman and Hall/CRC, Boca Raton, Fla (2004).
  * [13] J. Zinn-Justin, Quantum Field Theory and Critical Phenomena, Claredon Press 1989, 4th edn. Oxfod University Press (2002).
  * [14] M. Doi, J. Phys. A 9, 1465 (1976).
  * [15] L. Peliti, J. Phys. 46, 1469 (1985).
  * [16] B. P. Lee, J. Phys. A: Math. Gen. 27, 2633 (1994).
  * [17] L. D. Landau and E. M. Lifshitz, Fluid Mechanics, Pergamon Press, (1959).
  * [18] U. Frisch, Turbulence: The Legacy of A.N. Kolmogorov, Cambridge Univ. Pr., Cambridge, (1995).
  * [19] D. Forster, D. R. Nelson, M. J. Stephen, Phys. Rev. Lett. 36, 867 (1976).
  * [20] M. Hnatič, J. Honkonen, T. Lučivjanský, Symmetry 11 (10), 1193 (2019).
  * [21] J.-M. Park, M. W. Deem, Phys. Rev. E 57, 3618 (1998).
  * [22] M. J. E. Richardson and J. Cardy, J. Phys. A 32, 4035 (1999).
  * [23] M. W. Deem, J.-M. Park, Phys. Rev. E 58, 3223 (1998).
  * [24] M. Hnatich, J. Honkonen, Phys. Rev. E 61, 3904 (2000).
  * [25] J. Honkonen, Acta Phys. Slovaca 52, 533 (2002).
  * [26] M. Hnatich, J. Honkonen, T. Lučivjanský, Theor. Math. Phys. 169, 1481 (2011).
  * [27] M. Hnatič, J. Honkonen, T. Lučivjanský, Symmetry 11, 1193 (2019).