11institutetext: École nationale Supérieure en Informatique, Laboratoire LMCS,
Oued-Smar, Algiers, Algeria, 11email: {e_guerrout, r_mahiou, fr_boukabene,
<EMAIL_ADDRESS>
# Image Segmentation Methods for Non-destructive testing Applications

EL-Hachemi Guerrout  Ramdane Mahiou  Randa Boukabene  Assia Ouali

###### Abstract

In this paper, we present new image segmentation methods based on hidden
Markov random fields (HMRFs) and cuckoo search (CS) variants. HMRFs model the
segmentation problem as a minimization of an energy function. CS algorithm is
one of the recent powerful optimization techniques. Therefore, five variants
of the CS algorithm are used to compute a solution. Through tests, we conduct
a study to choose the CS variant with parameters that give good results
(execution time and quality of segmentation). CS variants are evaluated and
compared with non-destructive testing (NDT) images using a misclassification
error (ME) criterion.

###### Keywords:

I

mage segmentation; Non-Destructive Testing (NDT) images; Hidden Markov Random
Fields; Cuckoo Search Variants; Misclassification Error.

##  1 INTRODUCTION

Non Destructive Testing (NDT) is a set of valuable techniques. They have many
uses. For example, they are used in manufacturing to inspect materials and
extract defective regions without causing damages [1, 7, 9, 5]. NDT methods
rely on image analysis where segmentation image is the most important.

Hidden Markov random fields (HMRFs) [3, 2, 17, 4] proved their efficiency and
robustness in image segmentation. These latter take into consideration the
global and local properties of the image. The global properties include the
mean and standard deviation values of different pixel classes. The local
properties are characterized by the interactions of neighboring pixels.
Looking for the segmented image is seen as an optimization problem. Therefore,
optimization techniques are used to compute a solution.

The cuckoo search algorithm is one of the latest powerful nature-inspired
meta-heuristic algorithms [15]. It is inspired by the behavior of some cuckoo
types that put their eggs in the nests of other birds (known as brood
parasitism).

In this paper, we shed light on the combinations of HMRFs with five variants
of the cuckoo search algorithm that are: standard CS (SCS) [15, 16, 14],
improved CS (ICS) [11], auto adaptive modified CS (AACS) [8], modified CS
(MCS) [12] and novel modified (NMCS) [13]. The quality of the solution is very
sensitive to the choice of parameters. For that reason, we conduct an
evaluative study in order to choose parameters that give a good segmentation.
Misclassification error (ME) criterion [10] measures the difference between
the segmentation result and the ground truth (GT). Tests are carried out on
NDT (Non-Destructive Testing) image dataset [10].

The rest of the paper is organized as follows. In section 2, we present
segmentation methods that are based on HMRF and CS variants. Section 3 is
devoted to some experimental results. Finally, we conclude in section 4.

##  2 HMRFs and cuckoo search variants

In image segmentation using HMRFs both image to segment and segmented image
are seen as realizations of Markov random variables.

Let $y=({y}_{1},\dots,{y}_{s},{\dots},{y}_{M})$ be the image to segment,
$y_{s}$ is the pixel value of the site $s$ that takes its values in the gray
level space $E_{y}=\\{0,\ldots,255\\}$ and configurations set
$\Omega_{y}=E_{y}^{M}$.

Let $x=({x}_{1},\dots,{x}_{s},{\dots},{x}_{M})$ be the segmented image into
$K$ classes, $x_{s}$ is the class of the site $s$ that takes its values in the
discrete space $E_{x}=\left\\{1,{\dots},K\right\\}$ and configurations set
$\Omega_{x}=E_{x}^{M}$.

Let $\mu=(\mu_{1},\dots,\mu_{j},\dots,\mu_{K})$ be the means and
$\sigma=(\sigma_{1},\dots,\sigma_{j},\dots,\sigma_{K})$ be the standard
deviations of $K$ classes in the segmented image
$x=(x_{1},\dots,x_{s},\dots,x_{M})$ i.e.,

${}\begin{array}[]{l}\begin{cases}\mu_{j}={\frac{1}{|S_{j}|}\sum_{s\in
S_{j}}y_{s}}\\\ \\\ \sigma_{j}=\sqrt{\frac{1}{|S_{j}|}\sum_{s\in
S_{j}}(y_{s}-\mu_{j})^{2}}\\\ \\\ S_{j}=\\{s\ |\
x_{s}=j\\}\end{cases}\end{array}$ (1)

Image segmentation using HMRFs is modeled in [4, 6] as a minimization of an
energy function, which is presented in (2). We can always compute $x$ through
$\mu$ by classifying $y_{s}$ into the nearest mean $\mu_{j}$ i.e., $x_{s}=j$
if the nearest mean to $y_{s}$ is $\mu_{j}$. Thus, we look for $\mu^{*}$
instead of $x^{*}$. The configuration set of $\mu$ is $\Omega_{\mu}=[0\dots
255]^{K}$.

$\begin{array}[]{l}\begin{cases}{\mu}^{*}=\operatorname*{arg}_{\mu\in\Omega_{\mu}}{\mathit{\min}}\left\\{\Psi(\mu)\right\\}\\\
\\\ \Psi(\mu)=\sum_{j=1}^{K}\sum\limits_{s\in
S_{j}}{[\ln(\sigma_{j})+\frac{(y_{s}-\mu_{j})^{2}}{2\sigma_{j}^{2}}]}+\frac{B}{T}\sum_{c_{2}=\\{s,t\\}}{(1-2\delta(x_{s},x_{t}))}\end{cases}\end{array}$
(2)

where $B$ is a constant, $T$ is a control parameter called temperature,
$\delta$ is Kronecker’s delta and $S_{j}$, $\mu_{j}$ and $\sigma_{j}$ are
defined in (1). When $B>0$, the most likely segmentation corresponds to the
constitution of large homogeneous regions. The size of these regions is
controlled by the parameter $B$.

To apply unconstrained optimization techniques, we redefine the function
$\Psi(\mu)$ for $\mu\in\mathbb{R}^{K}$ instead of $\mu\in[0\dots 255]^{K}$.
Therefore, the new function $\Psi(\mu)$ becomes as follows:

$\begin{array}[]{l}\Psi(\mu)=\begin{cases}\sum_{j=1}^{K}\sum\limits_{s\in
S_{j}}{[\ln(\sigma_{j})+\frac{(y_{s}-\mu_{j})^{2}}{2\sigma_{j}^{2}}]}&\hbox{\multirowsetup
$\mbox{if}\quad\mu\in[0\dots 255]^{K}$}\\\
+\frac{B}{T}\sum_{c_{2}=\\{s,t\\}}{(1-2\delta(x_{s},x_{t}))}&\\\ \\\
+\infty&\mbox{otherwise}\end{cases}\end{array}$ (3)

To solve the minimization problem expressed in 3, we have used Cuckoo Search
(CS) variants. The principal steps are set out below.

Each egg in a host nest represents a solution
$\mu^{i,t}=(\mu^{i,t}_{1},\dots,\mu^{i,t}_{j},\dots,\mu^{i,t}_{K})$ at time
$t$. Let $c^{i,t}=(c^{i,t}_{1},,\dots,c^{i,t}_{j},\dots,c^{i,t}_{K})$ be
cuckoo egg generated at time $t$. Let $n$ be the number of available host
nests (or different solutions). The initial population
$\\{\mu^{i,0}\\}_{i=1,\dots,n}$ is generated by random initialization. Let
$\mbox{best}^{t}=(\mbox{best}^{t}_{1},\dots,\mbox{best}^{t}_{j},\dots,\mbox{best}^{t}_{K})$
be the best solution at time $t$.

$\mbox{best}^{t}:=\arg\min\limits_{\mu^{i,t}_{i=0,\dots,n}}{\Psi(\mu^{i,t})}$
(4)

Cuckoo Search (CS) algorithm is based on three rules defined in [15, 16] and
proceeds as follows:

  1. 1.

Generating new cuckoos $\\{c^{i,t}\\}_{i=1,\dots,n}$ can be performed as
follows (Lévy flight):

$c^{i,t}:=\mu^{i,t}+\alpha\times\mbox{\mbox{step}}\otimes(\mu^{i,t}-\mbox{best}^{t})\otimes\mbox{randn}(1,K)$
(5)

where $\otimes$ means the entry-wise product of two vectors. The
$\mbox{randn}(1,K)$ returns $K$ random numbers from a normal distribution with
mean 0 and variance 1. $\alpha>0$ denotes a scaling factor, which is related
to the problem under study. In our implementation, the step length (step) is
calculated by using Mantegna’s algorithm (as defined in [14]).

The nests are updated as follows:

$\mu^{i,t}:=\begin{cases}c^{i,t}&\mbox{if }\Psi(c^{i,t})\leq\Psi(\mu^{i,t})\\\
\\\ \mu^{i,t}&\mbox{otherwise}\end{cases}$ (6)

  2. 2.

A fraction $p_{a}$ of deficient nests are abandoned. New nests
$\\{v^{i,t}\\}_{i=1,\dots,n}$ are generated by biased random walks (crossover
operator) as follows:

$v^{i,t}:=\mu^{i,t}+\mbox{rand()}\otimes(\mu^{r1,t}-\mu^{r2,t})\otimes[H(p_{a}-\mbox{rand()})]$
(7)

where $\mu^{r1,t}$ and $\mu^{r2,t}$ are two different solutions selected
randomly by random permutation. $H$ is Heaviside function. rand() returns a
single uniformly distributed random number in the interval (0,1). $[\ ]$ is a
$K$ vector.

  3. 3.

The best nests will be carried over to the next generations ($t+1$) as
follows:

$\mu^{i,t+1}:=\begin{cases}v^{i,t}&\mbox{if
}\Psi(v^{i,t})\leq\Psi(\mu^{i,t})\\\ \\\
\mu^{i,t}&\mbox{otherwise}\end{cases}$ (8)

Each variant of CS (SCS, ICS, NMCS, MCS and AACS) manages in a specific way
$p_{a}$, $\alpha$ and step. The standard CS (SCS) uses constant values for the
parameters $p_{a}$ and $\alpha$. As for the improved CS (ICS), the parameters
$p_{a}$ and $\alpha$ are calculated at the time $t$ as expressed in (9) and
(10):

$p_{a}(t)={p_{a}}_{\mbox{max}}-\frac{t}{\mbox{NI}}({p_{a}}_{\mbox{max}}-{p_{a}}_{\mbox{min}})$
(9)

$\alpha(t)=\alpha_{\mbox{max}}\exp\left(\frac{t\times\ln\left(\frac{\alpha_{\mbox{min}}}{\alpha_{\mbox{max}}}\right)}{\mbox{NI}}\right)$
(10)

where NI is the maximum number of iterations.

Each time, the novel modified CS (NMCS) computes the parameters $p_{a}$ and
$\alpha$ as expressed in (11) and (12).

$\begin{array}[]{l}p_{a}(t)=\begin{cases}0.1&\mbox{if}\qquad
p_{a}(0)/\xi(t)+p_{a}(0)\times\theta(t)<0.1\\\ \\\ 0.85&\mbox{if}\qquad
p_{a}(0)/\xi(t)+p_{a}(0)\times\theta(t)>0.85\\\ \\\
p_{a}(0)/\xi(t)+p_{a}(0)\times\theta(t)&\mbox{otherwise}\end{cases}\end{array}$
(11)

$\alpha(t)=\alpha(0)/\xi(t)+\alpha(0)\times\theta(t)$ (12)

where $\xi(t)$ and $\theta(t)$ are the scale conversion factors computed using
the speed factor ($\omega_{s}(t)$) and the aggregation factor
($\omega_{a}(t)$) as expressed in (13).

${}\begin{array}[]{l}\begin{cases}\xi(t)=10^{10\times\tan(\arctan(0.1)\times(2\times\omega_{s}(t)-1))}\\\
\\\ \theta(t)=10^{10\times\tan(\arctan(0.1)\times(2\times\omega_{a}(t)-1))}\\\
\\\ \omega_{s}(t)=\frac{\Psi(\mbox{best}^{t})}{\Psi(\mbox{best}^{t-1})}\\\ \\\
\omega_{a}(t)=\frac{\Psi(\mbox{best}^{t})}{\frac{1}{n}\sum_{i=1}^{n}(\Psi(\mu^{i,t}))}\end{cases}\end{array}$
(13)

The Auto adaptive CS (AACS) variant uses a procedure parameter setting that
computes for every nest $i$ at the time $t$ a self-adaptive parameter
$p_{a}^{i}(t)$ (see [8] for more detail). AACS builds a new nests
$\\{v^{i,t}\\}_{i=1,\dots,n}$ using the mutation, crossover and selection
operators as follows:

$\begin{array}[]{l}\mbox{If}\quad\mbox{rand()}>=(1-t/\mbox{NI})\qquad\mbox{then}\\\
\quad
v^{i,t}:=\begin{cases}\mbox{best}^{t}+\varphi\otimes(\mu^{r1,t}-\mu^{r2,t}+\mu^{r3,t}-\mu^{r4,t})&\mbox{if}\quad\mbox{rand()}<p^{i}_{a}(t)\\\
\\\ \mu^{i,t}&\mbox{otherwise}\end{cases}\\\ \mbox{Else}\\\ \quad
v^{i,t}:=\begin{cases}\mu^{r1,t}+\varphi\otimes(\mu^{r2,t}-\mu^{r3,t})&\mbox{if}\quad\mbox{rand()}<p^{i}_{a}(t)\\\
\\\ \mu^{i,t}&\mbox{otherwise}\end{cases}\\\ \mbox{End If}\end{array}$ (14)

where $\mu^{r1,t}$, $\mu^{r2,t}$, $\mu^{r3,t}$ and $\mu^{r4,t}$ are four
different solutions selected randomly by random permutation. rand() returns a
single uniformly distributed random number in the interval (0,1). $\varphi$ is
the scale factor that represents a Gaussian distribution with mean $0.5$ and
standard deviation $0.1$.

In Algorithm 1, we summarize four combinations: HMRF-SCS, HMRF-ICS, HMRF-NMCS,
and HMRF-AACS.

The objective function is
$\Psi(\mu),\mu=(\mu_{1},\dots,\mu_{K})\in\mathbb{R}^{K}$

Generate the initial population of $n$ host nests
$\\{\mu^{i,0}\\}_{i=1,\dots,n}$

while  _( $t<\mbox{NI}$) _ do

Compute the best solution $\mbox{best}^{t}$ using (4)

if  _ICS_ then

compute $\alpha(t)$ using (10)

end if

if  _NMCS_ then

compute $\alpha(t)$ using (12)

end if

Generating new cuckoos $\\{c^{i,t}\\}_{i=1,\dots,n}$ using (5)

Update nests $\\{\mu^{i,t}\\}_{i=1,\dots,n}$ using (6)

if  _ICS_ then

compute $p_{a}(t)$ using (9)

end if

if  _NMCS_ then

compute $p_{a}(t)$ using (11)

end if

if  _AACS_ then

compute $p^{i}_{a}(t)$ using the procedure parameter setting

end if

A fraction $p_{a}$ of deficient nests are abandoned

if  _AACS_ then

New nests $\\{v^{i,t}\\}_{i=1,\dots,n}$ are built using (14)

else

New nests $\\{v^{i,t}\\}_{i=1,\dots,n}$ are built using (7)

end if

Keep best nests to the next generations ($t+1$) using (8)

$t:=t+1$

end while

Algorithm 1 HMRF-SCS, HMRF-ICS, HMRF-NMCS, and HMRF-AACS algorithms.

The modified CS (MCS) divides nests into two groups. The first group is the
top nests with high fitness ($25\%$ of nests). The second group is the
abandoned nest with low fitness ($75\%$ of nests). The HMRF-MCS combination is
summarized in Algorithm 2.

The objective function is
$\Psi(\mu),\mu=(\mu_{1},\dots,\mu_{K})\in\mathbb{R}^{K}$

Generate the initial population of $n$ host nests
$\\{\mu^{i,0}\\}_{i=1,\dots,n}$

while  _( $t<\mbox{NI}$) _ do

Split eggs into $N1=$top nests (25%) and $N2=$abandoned nests (75%)

foreach  _i $\in$ N2_ do

Get new cuckoos using 5 with $\alpha=\frac{1}{\sqrt{t}}$

Replace the nest $i$ by the new one.

end foreach

foreach  _i $\in$ N1_ do

pick a nest $j\in N1$ randomly

if  _$\mu^{i,t}=\mu^{j,t}$_ then

Get new cuckoos using 5 with $\alpha=\frac{1}{t^{2}}$

Choose a random nest $l$ from all nests

Update the nest $l$ by the new one using 6

else

dx $=|\mu^{i,t}-\mu^{j,t}|/\phi$ with $\phi=(1+\sqrt{5})/2$

Move distance dx from the worst nest to the best nest to generate a new one

Choose a random nest $l$ from all nests

Update the nest $l$ by the new one using 6

end if

end foreach

Update $t:=t+1$

end while

Algorithm 2 HMRF-MCS algorithm.

##  3 Experimental Results

To have a meaningful study and comparison, we have performed our tests in two
steps. In the first step, we have conducted an evaluative study for each CS
variant to select parameters that give a good segmentation. Each algorithm is
tested upon 25 images of the NDT dataset [10]. The values of parameter $n$
(number of nests) tested are $5,10,15,20,25,\mbox{ and }30$. The values of
parameter $NI$ (number of iterations) tested are $50\mbox{ and }100$. The
values of parameter $T$ (temperature) tested are $2\mbox{ and }3$. The other
parameters ($\alpha,p_{a},\phi,\varphi,\mbox{ and }\mbox{ step}$) are used as
recommended in the basic articles [15, 11, 8, 12, 13]. In the second step, we
carried out a comparative study of different combinations with the best
parameters selected in the first step. The comparison focuses on the execution
time and misclassification error results.

###  3.1 Misclassification error criterion (ME)

In the case of binary images, ME measures the difference between segmented
image and ground truth (GT). The binary images are constituted by a foreground
and a background. ME provides the percentage of misclassification pixels as
follows[10]:

$ME=1-\frac{|B_{O}\bigcap B_{T}|+|F_{O}\bigcap F_{T}|}{|B_{O}|+|F_{O}|}$ (15)

where $F_{O}$ and $B_{O}$ indicate foreground and background of ground truth
(GT), $F_{T}$ and $B_{T}$ indicate foreground and background of the segmented
image. ME equals zero that means the perfect match between segmented image and
its ground-truth.

###  3.2 Tests context

Algorithms are implemented in MATLAB 2017a on a computer with Intel Core i7
1.8 GHz CPU, 8G RAM, and Microsoft Windows 10.

###  3.3 Results

After the evaluative study, we have selected for each algorithm the parameters
that give the best results, which are shown in Table 1. Mean ME values and
mean execution time are calculated using all NDT images. In Table 1, the best
results are given in bold type.

Methods | Parameters | ME | Time(s)
---|---|---|---
SCS | n=20, NI=100, T=2 | 0.067 | 13.634
ICS | n=20, NI=100, T=2 | 0.046 | 13.516
AACS | n=25, NI=100, T=2 | 0.055 | 26.313
MCS | n=5, NI=100, T=3 | 0.058 | 1.329
NMCS | n=5, NI=50, T=2 | 0.066 | 1.633
Table 1: Parameters, mean ME values, and mean execution time.

Figures 1 and 2 provide a comparison between algorithms on ME and execution
time respectively compared to NDT images. The algorithms are executed with the
parameters presented in Table 1. Figure 3 gives the visual results, which
provide an integrated view with ME.

Figure 1: Misclassification error of each algorithm compared to images. Figure 2: The execution time of each algorithm with the best parameters compared to images. No | Image | GT | SCS | ICS | AACS | MCS | NMCS
---|---|---|---|---|---|---|---
1 |  |  |  |  |  |  |
2 |  |  |  |  |  |  |
3 |  |  |  |  |  |  |
4 |  |  |  |  |  |  |
5 |  |  |  |  |  |  |
6 |  |  |  |  |  |  |
7 |  |  |  |  |  |  |
8 |  |  |  |  |  |  |
9 |  |  |  |  |  |  |
10 |  |  |  |  |  |  |
11 |  |  |  |  |  |  |
12 |  |  |  |  |  |  |
13 |  |  |  |  |  |  |
14 |  |  |  |  |  |  |
15 |  |  |  |  |  |  |
16 |  |  |  |  |  |  |
17 |  |  |  |  |  |  |
18 |  |  |  |  |  |  |
19 |  |  |  |  |  |  |
20 |  |  |  |  |  |  |
21 |  |  |  |  |  |  |
22 |  |  |  |  |  |  |
23 |  |  |  |  |  |  |
24 |  |  |  |  |  |  |
25 |  |  |  |  |  |  |
Figure 3: Visual results of each algorithm.

##  4 Conclusion and discussion

In this paper, we have presented new approaches that combine hidden Markov
random fields and cuckoo search variants to perform segmentation. Then, we
have conducted an evaluative study for each CS variant in order to select good
parameters. After that, we have carried out a comparative study between
different combinations. Each algorithm is tested upon 25 images where ground
truth (GT) segmentation is known. The comparison focuses on the execution time
and misclassification error (ME). HMRC-ICS combination shows the best results
$ME=0.046$ and 21 segmented images visually close to the ground truth. On the
other hand, HMRF-MCS presents an interesting execution time and ME not far
from the best results. To make a fair choice, prospectively we will invest in
a statistical study of the parameters: step, $\alpha$, $\phi$, $\varphi$, and
$p_{a}$.

## References

  * [1] Cartz, L.: Nondestructive testing (1995)
  * [2] Deng, H., Clausi, D.A.: Unsupervised image segmentation using a simple MRF model with a new implementation scheme. Pattern recognition 37(12), 2323–2335 (2004)
  * [3] Geman, S., Geman, D.: Stochastic relaxation, Gibbs distributions, and the bayesian restoration of images. Pattern Analysis and Machine Intelligence, IEEE Transactions on (6), 721–741 (1984)
  * [4] Guerrout, E.H., Ait-Aoudia, S., Michelucci, D., Mahiou, R.: Hidden Markov random field model and Broyden-Fletcher-Goldfarb-Shanno algorithm for brain image. Journal of Experimental & Theoretical Artificial Intelligence 30(3), 415–427 (2018)
  * [5] Guerrout, E.H., Mahiou, R., Ait-Aoudia, S.: Hidden Markov random fields and swarm particles: a winning combination in image segmentation. IERI Procedia 10, 19–24 (2014)
  * [6] Guerrout, E., Mahiou, R., Ait-Aoudia, S.: Hidden markov random fields and particle swarm combination for brain image segmentation. Int. Arab J. Inf. Technol. 15(3), 462–468 (2018)
  * [7] Hellier, C.: Handbook of nondestructive evaluation (2001)
  * [8] Li, X., Yin, M.: Modified cuckoo search algorithm with self adaptive parameter method. Information Sciences 298, 80–97 (2015)
  * [9] Mix, P.E.: Introduction to nondestructive testing: a training guide. John Wiley & Sons (2005)
  * [10] Sezgin, M., Sankur, B.: Survey over image thresholding techniques and quantitative performance evaluation. Journal of Electronic imaging 13(1), 146–166 (2004)
  * [11] Valian, E., Mohanna, S., Tavakoli, S.: Improved cuckoo search algorithm for global optimization. International Journal of Communications and Information Technology 1(1), 31–44 (2011)
  * [12] Walton, S., Hassan, O., Morgan, K., Brown, M.: Modified cuckoo search: a new gradient free optimisation algorithm. Chaos, Solitons & Fractals 44(9), 710–718 (2011)
  * [13] Yang, B., Miao, J., Fan, Z., Long, J., Liu, X.: Modified cuckoo search algorithm for the optimal placement of actuators problem. Applied Soft Computing 67, 48–60 (2018)
  * [14] Yang, X.S.: Nature-inspired metaheuristic algorithms. Luniver press (2010)
  * [15] Yang, X.S., Deb, S.: Cuckoo search via Lévy flights. In: 2009 World congress on nature & biologically inspired computing (NaBIC). pp. 210–214. IEEE (2009)
  * [16] Yang, X.S., Deb, S.: Engineering optimisation by cuckoo search. arXiv preprint arXiv:1005.2908 (2010)
  * [17] Zhang, Y., Brady, M., Smith, S.: Segmentation of brain MR images through a hidden Markov random field model and the expectation-maximization algorithm. IEEE transactions on medical imaging 20(1), 45–57 (2001)


1000

# Lattice pure gauge compact QED in the Landau gauge: the photon propagator,
the phase structure and the presence of Dirac strings

Lee C. Loveridge1,2<EMAIL_ADDRESS>Orlando Oliveira1<EMAIL_ADDRESS>Paulo J. Silva1<EMAIL_ADDRESS>
1 CFisUC, Department of Physics, University of Coimbra, 3004-516 Coimbra,
Portugal

2 Los Angeles Pierce College, 6201 Winnetka Ave., Woodland Hills CA 91371, USA

###### Abstract

In this work we investigate the lattice Landau gauge photon propagator
together with the average number of Dirac strings in the compact formulation
of QED for the pure gauge version of the theory as a function of the coupling
constant. Their $\beta$ dependence show that these two quantities can be used
to identify the confinement-deconfinement transition and that the nature of
this transition is first order. Our results show that in the confined phase
the propagator is always finite, the theory has a mass gap and the number of
Dirac strings present in the configuration is two orders of magnitude larger
than in the deconfined phase. Furthermore, in the deconfined phase where
$\beta\geq 1.0125$ the theory becomes massless, there are essentially no Dirac
strings and the photon propagator diverges when the limit $p\rightarrow 0^{+}$
is taken. Our results illustrate the importance of the topological structures
in the dynamics of the two phases.

##  I Introduction and Motivation

The regularization of QED using a hypercubic lattice Wilson:1974sk  takes as
fundamental fields the link variables

$U_{\mu}(x)\equiv\exp\\{i\,e\,a\,A_{\mu}(x+a\,\hat{e}_{\mu}/2)\\}\ ,$ (1)

where $A_{\mu}$ is the bare continuum photon field, $a$ is the lattice
spacing, $e$ is the bare coupling constant and $\hat{e}_{\mu}$ is the unit
vector along direction $\mu$. The link variables $U_{\mu}(x)$ are defined on a
compact manifold, they cover the unit circle centered around the origin of the
complex plane. On the other hand, the continuum photon field spans the real
numbers. The phase diagram of the regularized compact formulation of QED has
two phases that are distinguished by the value of the bare coupling constant
or, equivalently, by $\beta=1/e^{2}$.

For low values of $\beta$, i.e in the strong coupling limit, the static
potential between fermions fields grows linearly with the distance between the
fermion sources and the theory is confining. For $\beta\gtrsim 1$, i.e. in the
weak coupling limit, the static potential becomes essentially constant at
large distance separations, the theory is no longer confining, and the results
of the lattice formulation of QED approach those of the perturbative solution
of the continuum theory after taking the thermodynamic limit Banks:1977cc ;
Glimm:1977gz ; Fradkin:1978th ; Creutz:1979zg ; Guth:1979gz ; Lautrup:1980xr ;
Frohlich:1982gf ; Berg:1983is ; Kogut:1987cd ; Coddington:1987yz ;
Nakamura:1991ww ; Jersak:1996mj ; Gubarev:2000eu ; Arnold:2002jk ;
Panero:2005iu ; Loveridge:2021qzs .

The confinement mechanism seems to be related to the topological structure of
the Abelian gauge group $U(1)$ that, for the 3D case, is associated with the
presence of Dirac monopoles. These classical configurations are absent in the
deconfined phase but are observed in the confined phase Polyakov:1975rs ;
Polyakov:1976fu ; Chernodub:2002gp ; Bakker:2000kg ; Greensite:2011zz .
Furthermore, in four dimensions, the confined phase has a mass gap
Loveridge:2021qzs . In the strong coupling limit compact QED has a mechanism
that generates a photon mass gap111The photon propagator being finite in the
confined phase, means that its functional form should be of the type
$Z(p^{2})/(p^{2}+M^{2}(p^{2}))$ and that the running photon mass
$M^{2}(p^{2})$ does not vanish in the infrared region. In the main text when
we do refer to a photon mass we have in mind the function $M(p^{2})$. Recall
that in the deconfined phase the propagator is compatible with a divergent
behaviour that implies a vanishing $M^{2}(p^{2})$ at zero momentum. It is in
this sense that the deconfined phase is understood as being a massless theory.
and the theory is confining.

The presence of a mass gap is also observed in the simulation of non-Abelian
gauge theories. Indeed, the generation of mass scales was observed in non-
Abelian gauge theories as in QCD and in the pure Yang-Mills SU(2) gauge theory
Alkofer:2003jj ; Cucchieri:2007md ; Fischer:2008uz ; Aguilar:2008xm ;
Boucaud:2008ky ; Dudal:2008sp ; Bogolubsky:2009dc ; Dudal:2010tf ;
Oliveira:2012eh ; Duarte:2016iko ; Cucchieri:2016jwg ; Cyrol:2016tym ;
Dudal:2018cli ; Huber:2018ned ; Oliveira:2018lln ; Gracey:2019xom ;
Dudal:2019gvn ; Huber:2020keu ; Fischer:2020xnb ; Falcao:2020vyr ;
Aguilar:2021lke ; Gao:2021wun ; Lo:2021qkw ; Hayashi:2021nnj .

The understanding of the confinement mechanism for non-Abelian gauge theories
and its possible connection with the generation of mass gaps is still a
fundamental open problem for Particle Physics. Hopefully, the comparison of
the compact QED formulation with the non-Abelian case will bring further
insight into the confinement problem.

The dynamical properties of compact QED depend on the value of $\beta$. This
dependence translates into deep changes in the properties of the QED Green
functions, i.e. in the propagators and vertices of the theory. Indeed, the
computation of the photon propagator in the Landau gauge for pure gauge
lattice compact QED at low $\beta$ and at high $\beta$ illustrates those
differences Berg:1983is ; Nakamura:1990dq ; Nakamura:1991ww ; Durr:2002jc ;
Loveridge:2021qzs . For the confined low $\beta$ phase, the photon propagator
seems to be described by a Yukawa-type of propagator. The theory has a mass
gap and the topological structures, measured by the average number of Dirac
strings that are observed in the gauge configurations, are orders of magnitude
higher than in the large $\beta$ phase where confinement is not present.
Furthermore, in the deconfined high $\beta$ phase the photon propagator seems
to approach the continuum like $1/p^{2}$ behaviour of a free field theory as
we go towards higher volumes.

In the current work we revisit the Landau gauge photon propagator computed
using the lattice formulation of pure gauge compact QED in four dimensions and
study its $\beta$ dependence. The investigations reported here complement the
recently published paper Loveridge:2021qzs  that illustrates how different the
photon propagator is in the confined and deconfined phases.

Our study shows that in the confined phase the photon propagator is finite
over the full range of momenta, i.e. it has a mass gap, while in the
deconfined phase it seems to diverge as $p\rightarrow 0$. For the deconfined
phase the zero momentum propagator increases with $\beta$. The average number
of Dirac strings present in the gauge configurations is also phase dependent.
The number of Dirac strings is much larger in the confined phase, in
comparison with the deconfined phase where the Dirac strings are almost non-
existent. Moreover, the average number of Dirac strings is strongly correlated
with the properties of the photon propagator and this correlation suggests
that the origin of the photon mass observed in the confined phase is due to
these topological structures. Although we do not compute the photon mass as a
function of the coupling constant, as this would require modelling the
propagator, in the confined phase the propagator is finite for all momenta a
result that can be viewed as due to the generation of a finite non-vanishing
running mass for the photon field. In the same sense, the divergence of the
photon propagator at zero momentum in the deconfined phase can be read as an
indication that the theory becomes massless.

Our results show that the transition between the two phases can be clearly
identified by looking either at the photon propagator data or at the average
number of Dirac strings. From either of these two quantities it is possible to
estimate the corresponding critical value of $\beta_{c}$ where the phase
transition occurs. In this sense, the analysis of the propagator and/or the
number of Dirac strings can be used as order parameters for the confinement-
deconfinement transition for pure gauge compact QED. Further, the photon
propagator and the Dirac strings data show that the transition to the
deconfined phase is of first order in good agreement with the literature
Fradkin:1978th ; Guth:1979gz ; Lautrup:1980xr ; Frohlich:1982gf ; Kogut:1987cd
; Arnold:2000hf .

The understanding of the phase diagram of compact QED is important per se.
Historically QED is the first quantum field theory to be studied and has
furnished a number of results and methods that were later generalized to other
types of theories. Also QED is a part of the Standard Model of Particle
Physics and the interplay of electrodynamics with strong interactions is
currently a subject of research. Furthermore, it is essential to have a good
theoretical understanding of QED to perceive the Higgs sector of the Standard
Model and of Abelian Higgs models, that are also relevant for Condensed Matter
systems. QED has a fundamental role in Condensed Matter Physics and Atomic
Physics and $U(1)$ gauge theories are also being used to explore the
possibility of performing realistic simulations with quantum computers.

The paper is organized as follows. In Sec. II we define the theoretical setup
of our study, giving details on the gauge fixing and on the computation of the
photon propagator. In Sec III we report on the photon propagator results for
various $\beta$ and how the propagator evolves with $\beta$. In Sec. IV the
Dirac string content and distributions in the different phases of the theory
is analyzed. The correlation of these topological structures with the photon
propagator data is also discussed there. Issues related to the performance of
the sampling with the hybrid Monte Carlo method are also touched. Finally, in
Sec. V we summarize and conclude.

##  II Pure Gauge Compact QED

The simulations considered herein refer to the compact version of QED defined
over an hypercubic lattice and described by the Wilson action. In Euclidean
space the Wilson action reads

$S_{W}(U)=\beta\sum_{x}\sum_{1\leqslant\mu,\nu\leqslant
4}\bigg{\\{}1-\Re\,\left[U_{\mu\nu}(x)\right]\bigg{\\}}\ ,$ (2)

where the plaquette operator is given by

$U_{\mu\nu}(x)=U_{\mu}(x)\,U_{\nu}(x+a\,\hat{e}_{\mu})\,U^{\dagger}_{\mu}(x+a\,\hat{e}_{\nu})\,U^{\dagger}_{\nu}(x)$
(3)

that, in the continuum limit, is

$U_{\mu\nu}(x)=\exp\left\\{i\,e\,\oint_{C}A_{\mu}(z)\,dz_{\mu}\right\\}$ (4)

with $C$ being any closed curve that contains in its interior the point $x$
and whose points are infinitesimally close to $x$. On an hypercubic lattice
whose lattice spacing is $a$, the link variables are related to the photon
field $A_{\mu}$ by Eq. (1). It follows from the definition given in Eq. (3)
that the exponential term is the change of the photon field around a plaquette
centered at $x+a\,(\hat{e}_{\mu}+\hat{e}_{\nu})/2$. Indeed, writing for the
plaquette

$U_{\mu\nu}(x)=\exp\left\\{i\,e\,a\,\Big{(}\,\Delta
A_{\mu\nu}(x)\,\Big{)}\right\\}\ ,$ (5)

given that $-\pi\leqslant e\,a\,A_{\mu}\leqslant\pi$ and that $-\pi\leqslant
e\,a\leavevmode\nobreak\ \Delta A_{\mu\nu}(x)\leqslant\pi$, it follows that

$\displaystyle\Delta A_{\mu\nu}(x)$ $\displaystyle=$ $\displaystyle
A_{\mu}\left(x+\frac{a}{2}\,\hat{e}_{\mu}\right)+A_{\nu}\left(x+a\,\hat{e}_{\mu}+\frac{a}{2}\,\hat{e}_{\nu}\right)$
(6) $\displaystyle\quad-
A_{\mu}\left(x+a\,\hat{e}_{\nu}+\frac{a}{2}\,\hat{e}_{\mu}\right)-A_{\mu}\left(x+\frac{a}{2}\,\hat{e}_{\nu}\right)$
$\displaystyle\qquad\quad+\,\frac{2\,\pi\,m_{\mu\nu}(x)}{e\,a}$

where the integer field $m_{\mu\nu}(x)$ measures the number of Dirac strings
that cross the plaquette associated with $U_{\mu\nu}(x)$. In a simulation, for
a given gauge configuration, the field $m_{\mu\nu}(x)$ can be computed by
combining links and plaquettes as described in DeGrand:1980eq ;
Loveridge:2021qzs . As was shown in Loveridge:2021qzs , the confined and
deconfined phases can be distinguished by looking at $m_{\mu\nu}(x)$ or to its
average over the lattice.

In the simulations reported below, the links were sampled with the Wilson
action (2), relying on the hybrid Monte Carlo method (HMC) Duane:1987de  that
was implemented with the Chroma library that requires the QDP++ library
Edwards:2004sx . It is well known that the HMC has problems simulating compact
QED Grosch:1985cz ; Kerler:1995va , as the various topological sectors are not
properly sampled. The topological freezing of the HMC becomes more important
close to the transition between the confined and deconfined region and also in
the deconfined phase.

In an effort to estimate the bias introduced by the sampling method, close and
above the transition between the two phases, i.e. for $\beta\geq\beta_{c}$,
several simulations were performed using the same $\beta$ and taking different
starting points to initiate the Monte Carlo. We also considered long runs for
some of the $\beta$’s. Our results show that for $\beta\leq\beta_{c}$ the
different Markov chains result in the same photon propagator. However, for
$\beta\geq\beta_{c}$ the photon propagators differ between the different
simulations. A correlation between the photon propagator and the average value
of Dirac strings is observed; see the results of Sec IV.1 for details. For
$\beta\geq\beta_{c}$ the propagators for the different simulations with the
same $\beta$ value overlap in the UV region but show different functional
forms in the infrared region. We also also find that for $\beta\geq\beta_{c}$
and for the different simulations with the same $\beta$, the photon propagator
resulting from the configurations with the smaller number of Dirac strings is
closer to the perturbative propagator. If for the deconfined high $\beta$
phase the photon propagator obtained with the hybrid Monte Carlo method shows
a dependence on the initial configuration in the infrared region, the overall
analysis of the $\beta$ dependence is robust and does not change qualitatively
by taking either of the photon propagators computed for $\beta>\beta_{c}$. On
the other hand, for the confined low $\beta$ phase our results suggest that
the HMC is robust. For the deconfined high $\beta$ phase our simulations
confirm the sampling problems of the hybrid Monte Carlo method previously
identified in the literature.

After importance sampling, the links were rotated towards the Landau gauge as
in Loveridge:2021qzs . On a first stage, we rely on the linear definition for
the photon field given by

$e\,a\,A_{\mu}\left(x+\frac{a}{2}\hat{e}_{\mu}\right)=\frac{U_{\mu}(x)-U^{\dagger}_{\mu}(x)}{2\,i}$
(7)

and maximize the functional

$F[U;g]=\frac{1}{V\,D}\sum_{x,\mu}\,\Re\left[\,g(x)\,U_{\mu}(x)\,g^{\dagger}(x+a\,\hat{e}_{\mu})\,\right]$
(8)

over the gauge orbit. In Eq. (8) $g(x)\in U(1)$, $V$ is the total number of
lattice points and $D$ the Euclidean spacetime dimension. In a second stage,
the photon field is computed with the logarithmic definition

$e\,a\,A_{\mu}\left(x+\frac{a}{2}\hat{e}_{\mu}\right)=-i\,\ln\Big{(}U_{\mu}(x)\Big{)}\
,$ (9)

that provides an exact definition, up to machine precision, and does not rely
on the use of a small lattice spacing. Then, the Landau gauge condition is
achieved by maximizing, over the gauge orbits, the functional
Ilgenfritz:2010gu

$\widetilde{F}[U;g]=\frac{1}{V\,D}\sum_{x,\mu}\,\bigg{\\{}1-a^{2}e^{2}\left[A^{(g)}_{\mu}\left(x+\frac{a}{2}\hat{e}_{\mu}\right)\right]^{2}\bigg{\\}}$
(10)

where the field $e\,a\,A^{(g)}$ is the photon field as given by Eq. (9) using
the links $U_{\mu}(x)$ gauge transformed by $g(x)$. The gauge fixing towards
the Landau gauge is, in all cases, monitored computing the lattice version of
$\partial\cdot A$ and, for each stage, the gauge fixing was stopped for an
averaged value over the lattice of $|\partial\cdot A|^{2}<10^{-15}$. Further
details on the gauge fixing can be found in Loveridge:2021qzs .

The optimizing functions that define the Landau gauge on the lattice have
multiple maxima that lead to different Landau gauge configurations, the Gribov
copies. In all the simulations we ignored the various maxima and the Landau
gauge fixing was performed starting the iterative process with $g(x)=1$ and
performing a single maximization for each gauge configuration obtained with
the Monte Carlo sampling of the Wilson action.

From the link variables, the momentum space photon field is computed with the
defnition

$A_{\mu}(p)=~\sum_{x}\,e^{-ip\cdot(x+\frac{a}{2}\hat{e}_{\mu})}A_{\mu}\left(x+\frac{a}{2}\hat{e}_{\mu}\right)\
,$ (11)

using the spacetime photon field given in (9). The Landau gauge propagator is
defined as

$\langle A_{\mu}(p_{1})\leavevmode\nobreak\
A_{\mu}(p_{2})\rangle=V\,\delta(p_{1}+p_{2})\,D_{\mu\nu}(p_{1})$ (12)

where $\langle\cdots\rangle$ stands for the average over the gauge fields. The
analysis of the propagator is performed assuming that the propagator has the
same tensor structure as in the continuum

$D_{\mu\nu}(p)=\left(\delta_{\mu\nu}-\frac{p_{\mu}p_{\nu}}{p^{2}}\right)\,D(\hat{p})\
.$ (13)

The form factor $D(\hat{p})$ is a function of the tree level improved momenta

$\displaystyle\hat{p}=\frac{2}{a}\sin\left(\frac{\pi}{L}~\,n_{\mu}\right)\ ,$
$\displaystyle\quad\quad
n_{\mu}=-\frac{L}{2},\,-\frac{L}{2}+1,\,\dots,\,0,\,1,\,\dots,\,\frac{L}{2}-1$
(14)

where $L$ is the number of lattice points in each side of the hypercubic
lattice. The rationale for using $\hat{p}$ instead of the naive lattice
momenta $p=2\pi n_{\mu}/aL$ is to reduce the finite space effects in the
propagator Leinweber:1998uu ; Catumba:2021hcx , a procedure developed for
asymptotically free field theories. The lattice data for the propagator shown
here satisfies the conical and cylindrical momentum cuts Leinweber:1998uu  for
$a\,\hat{p}>\Lambda_{\text{IR}}$. These momentum selections were set to
further suppress the finite space effects and to produce a well defined curve
from the lattice data. For $a\,\hat{p}\leqslant\Lambda_{\text{IR}}$ we follow
Dudal:2018cli  and consider all the momenta accessed in the simulations. In
all cases we use $\Lambda_{\text{IR}}=0.4$.

Assuming that the propagator tensor structure is as given in Eq. (13), then

$D(\hat{p})=\left\\{\begin{array}[]{l @{\hspace{1cm}}
l}\frac{1}{3}\sum_{\mu=1}^{4}D_{\mu\mu}(p),\hfil\hskip 28.45274pt&\hat{p}\neq
0\ ,\\\ \hfil\hskip 28.45274pt&\\\
\frac{1}{4}\sum_{\mu=1}^{4}D_{\mu\mu}(p),\hfil\hskip 28.45274pt&\hat{p}=0\
.\end{array}\right.$ (15)

The propagators for the different $\beta$ values were computed on a $48^{4}$
hypercubic lattice using the last (in the Markov chain) 200 gauge
configurations. The statistical errors are evaluated with the bootstrap method
for a confidence level of 67.5%, except for the fits where we use Gaussian
error propagation.

Figure 1: Dimensionless photon propagator as a function of $\beta$ for the
momenta selected by the cuts. Note the differences in the vertical scale. In
the bottom plots the vertical scale is logarithmic. The left graphs report the
bare photon propagator, while the right graphs show the renormalized photon
propagator as described in the main text. For the renormalized propagator the
inserted plots include only the high momenta, with the exception of the bottom
one that also shows the low momenta propagator for the larger $\beta$ values
considered.

##  III The photon propagator for various $\beta$

Figure 2: Dimensionless zero momentum bare (upper plot) and renormalized
(lower plot) photon propagator for various $\beta$. Note the logarithmic scale
for the vertical axis. The points associated with the same $\beta$ but
different color are the outcome of different sampling histories (see text for
further details).

In this work we measure the photon propagator for $\beta$ ranging from 0.2 up
to 2.0, covering in detail the region $\beta\approx 1$ where the transition
between the two phases is expected to occur.

The bare photon propagator form factor $D(\hat{p})$ is reported in Fig. 1 as a
function of the dimensionless improved momentum. The various plots have
different vertical scales and, in particular, for the higher $\beta$ the
vertical scale is logarithmic. In Fig. 2 (upper plot) the bare $D(0)$ is shown
as a function of $\beta$. For the deconfined phase the plots also include the
results of different simulations that where performed starting the Markov
chain differently. The data just referenced is bare lattice data and,
therefore, its comparison should be done with care. The bottom plot in Fig. 2
refers to renormalized propagators; more on renormalization later. This bottom
plot reports only the data associated with the orange points in the upper
(unrenormalized) plot.

As Fig. 2 shows, for large $\beta$ there is some dependence of $D(0)$ on the
starting configuration used to start the Markov chain; see the points marked
with different colors for the the same $\beta$. However, the overall result,
i.e. the sharp transition from a low $D(0)$ to a large $D(0)$, is not spoiled
by the problems that can be ascribed to the sampling method.

The data in Fig. 1 includes only the results of the runs that are reported in
Fig. 2 in orange. The bare data in Fig. 1 suggests that the photon propagator
is enhanced at low momenta, when compared to the higher momenta, for all
$\beta$ values, with the enhancement of $D(0)$ increasing with $\beta$.
Moreover, comparing the various plots there seems to be a change in the
functional form of $D(\hat{p})$ at $\beta\approx 1.0125$, with $D(\hat{p})$
becoming much steeper in the low momentum region for larger $\beta$. For
$\beta\gtrsim 1.0125$ the propagator becomes steeper and steeper as $\beta$
increases. This behaviour hints that above this $\beta$ value, the photon
propagator is divergent in the thermodynamical limit and, hopefully, recovers
the $1/p^{2}$ perturbative behaviour of a free field theory. The analysis of
the volume dependence of the $\beta=1.2$ Landau gauge photon propagator data
of Loveridge:2021qzs  supports this claim.

The bare lattice data for $D(0)$ in Fig. 2 also shows a sharp variation when
$\beta$ goes from $\sim 0.8$ to $\sim 1.1$ that is associated with the
transition to the deconfined continuum-like phase. It is possible to identify
two baselines for $D(0)$ in this Fig. for $\beta\lesssim 0.8$ and for
$\beta\gtrsim 1.1$. Note however, that for large $\beta$ the curve associated
with $D(0)$ does not seem to be so well defined due to problems with the
sampling method. Despite these problems the presence of a much higher typical
value for $D(0)$ for $\beta>\beta_{c}$ survives to the tests that were
performed.

In order to compare the propagators for the different $\beta$ values, one
should renormalize the lattice data. The problem for pure gauge QED being that
there is no clear way to set the scale for its lattice formulation, which
makes it impossible to choose a given momentum scale to renormalize the photon
field and, therefore, one has to look for alternatives. As a tentative method
to overcome the problem of the scale setting, we choose to fix the high
momentum behaviour of the photon propagator to renormalize the theory assuming
that $D(\hat{p})$ is independent of $\beta$ in the ultraviolet region. In
order to use such a procedure one has to assume that in all the simulations
the asymptotic ultraviolet region is accessed for all $\beta$ values. As the
results show, see the insertions in Fig. 1, this assumption seems to be
validated by the lattice data.

As a first step towards renormalizing the lattice data, for each $\beta$ value
the data for momenta in the range $a\,\hat{p}\in[3.5\,,\,4]$ is fitted to the
following functional form

$D_{\text{Fit}}(x)=a+\frac{b}{x^{2}}\ .$ (16)

Typical values for the $\chi^{2}/d.o.f.$ for the fits are below 1.6, with 0.62
and 2.2 being the smallest and highest observed figures, respectively. Then,
for each $\beta$, the corresponding fit is used to set $D(a\mu)=1$ at
$a\mu=3.8$. Recall that $a\hat{p}=4$ is the highest dimensionless momentum
accessed in the simulation and our choice of $a\mu$ is the middle point of the
fitting range.

Similarly as observed for the bare lattice data, the renormalized propagators
show a drastic change in the low momentum region for $\beta\gtrsim 1.0125$.
The data in Fig. 2 shows a transition from low $\beta$ and low $D(0)$, that
starts starts at $\beta\approx 0.8$ and ends at $\beta\approx 1.1$, with
$D(0)$ taking relatively large values for $\beta\gtrsim 1.1$ (note the log
scale on the vertical axis). Note also that again for $\beta\approx 1.1$ and
above there seems to be some dispersion on the values of $D(0)$ that are, in
principle, due to the use of the hybrid Monte Carlo method as sampling method.

The photon propagator data shows that this two point correlation function,
through the form factor $D(\hat{p})$, can be used to distinguish the phases of
pure gauge compact QED. It follows that the dynamics of the gauge fields in
the confined and deconfined phases are rather different.

Our conclusion for the photon propagator is inline with similar studies for
QCD where the gluon propagator was studied as a function of the temperature
Aouane:2011fv ; Silva:2013maa . Indeed, these studies shows that the gluon
propagator can be used to distinguish the confined and deconfined phases and
some authors suggested that $D(0)$ can be used as an order parameter for the
deconfinement transition. Furthermore, for the deconfined region the gluon
propagator can also be used to distinguish the topological sectors that are
associated with the center symmetry in the pure gauge Yang-Mills theory
Silva:2016onh .

##  IV Dirac strings and the Confinement and Deconfinement Phases

Figure 3: Dirac string distribution on a $32^{4}$ lattice for the confined
(upper plot) and deconfined (bottom plot) phases. The value $m_{\mu\nu}(x)=0$
was removed from the histogram as its contribution is much larger than those
reported in the plots. See text for details. Figure 4: Ensemble average of the
mean density of Dirac strings as a function of $\beta$.

The integer field $m_{\mu\nu}(x)$ is related to the topology of the gauge
fields as it measures the presence of Dirac strings in the gauge
configurations. For the two classes of gauge configurations considered in
Loveridge:2021qzs , those in the confined phase ($\beta=0.8$) and those in the
deconfined phase ($\beta=1.2$), $m_{\mu\nu}(x)$ was measured and in Fig. 9 of
Loveridge:2021qzs  the average value over the lattice of $|m_{\mu\nu}|$ for
the Landau gauge configurations for $\beta=0.8$ and $\beta=1.2$ is
illustrated. For these $\beta$’s the average value over the lattice of
$|m_{\mu\nu}|$ differs by two orders of magnitude, with the larger value being
associated with the configurations in the confined phase ($\beta=0.8$). This
was taken as an indication that Dirac strings play a major role in the
dynamics of the confined phase and are at the origin of the observed mass gap.

In Fig. 3 we report the typical distribution of the integer field $m_{\mu\nu}$
for $32^{4}$ gauge configurations in either of the phases. The upper plot
refers to a configuration in the confined phase, while the lower plot is for a
configuration in the deconfined phase. In both cases the averages over the
lattice of $m_{\mu\nu}$ are compatible with zero. In Fig. 3 we do not show the
bar associated with $m_{\mu\nu}=0$ that clearly dominates the distribution.
Indeed, for the configuration in the confined ($\beta=0.8$) phase the number
of vanishing $m_{\mu\nu}$ is about 16 times the number of $m_{\mu\nu}=\pm 1$,
while for deconfined phase this ratio increases to about 695. In the
deconfined phase, the number of vanishing $m$ is larger by a factor of $\sim
43$ relative to the confined case and the distribution of $m\neq 0$ becomes
asymmetric. These observations explains the pattern of the non-gauge invariant
quantity

$m=\frac{1}{6V}\sum_{x,\,\mu<\nu}\,|m_{\mu\nu}(x)|$ (17)

reported in Loveridge:2021qzs , see their Fig. 9, and why $m$ decreases
drastically in the deconfined phase in comparison with the confined phase.

Our new simulations offer an opportunity to monitor how the number of Dirac
strings evolves with $\beta$ and its connection with the generation of the
photon mass gap. We remind the reader that the mass gap is identified with a
finite $D(0)$ but no attempt is made to measure its value.

Fig 4 reports on the ensemble averages of $m$, named $\langle m\rangle$, as a
function of $\beta$. This quantity takes larger values in the confined phase
and smaller $\beta$ values, it has a sudden drop around $\beta\sim 1$ and
fluctuates at rather small values for $\beta\gtrsim 1$. To illustrate its
quantitative behaviour, the quantity $\langle m\rangle$ goes from
$0.1764749(47)$ at $\beta=0.2$ to $1.206(18)\times 10^{-3}$ for $\beta=1.2$
that represents a reduction by a factor of $\sim 146$.

The comparison of Figs. 4 and 2 establishes a correlation between $D(0)$ and
$\langle m\rangle$. Indeed, $D(0)$ follows the opposite behaviour of $\langle
m\rangle$ with a fast change and an almost infinite slope, taking place at the
same $\beta$ values. The sudden increase of $D(0)$, that can be translated
into a sudden decrease of the photon mass gap, occurs exactly when $\langle
m\rangle$ drops and becomes close to zero. This behaviour suggests, once more,
that the Dirac strings have a primordial role in the confined phase. The
correlation between the evolution with $\beta$ of $D(0)$ and $\langle
m\rangle$ suggests that the Dirac strings are also at the origin of the photon
mass that is non-vanishing only in the confined phase.

###  IV.1 The hybrid Monte Carlo sampling and the photon propagator

As mentioned previously, the hybrid Monte Carlo method seems to have problems
in sampling the Wilson action for compact QED. To check for possible bias, for
a number of $\beta$s several simulations were run using different gauge
configurations to start the Markov chain, and then compare the corresponding
propagators and the associated topological structures as measured by $\langle
m\rangle$.

Figure 5: Bare dimensionless data for $D(p^{2})$ and for several $\beta$ from
simulations using a different starting point in the Markov chain. Note the
logarithmic scale for the plots associated with $\beta>1$.

Figure 6: Bare dimensionless data for $p^{2}D(p^{2})$ and for several $\beta$
from simulations using a different starting point in the Markov chain.

In Fig. 5 we report on the photon propagator computed at several $\beta$ using
the results of the Monte Carlo simulations initiated with different starting
points. For the smallest $\beta$ value, that is below the confinement-
deconfinement transition, the photon propagators coming from the two
simulations are equal within one standard deviation, both evaluations of the
propagator are finite at zero momentum and the corresponding $\langle
m\rangle$’s are compatible within errors. This suggests that below the
confinement-deconfinement transition the hybrid Monte Carlo method is robust
and does a proper sampling of the compact QED action.

On the other hand, the simulations performed in the deconfined phase, i.e. for
$\beta>1$, show that the evaluation of the photon propagator with the hybrid
Monte Carlo method is (i) sensitive to the initialization of the Markov chain,
although in all cases the propagator increases and seems to diverge as one
approaches the $p\rightarrow 0^{+}$ limit, and (ii) for the ensembles with the
same $\beta$ those with smaller average number of Dirac strings seem to have a
propagator that is closer to the free field propagator $1/p^{2}$ behaviour in
the infrared region. This is better seen in Fig. 6 where the photon dressing
function $p^{2}D(p^{2})$ is plotted. For a free field type of propagator the
dressing function should be constant.

As Figs. 5 and 6 show, in the deconfined phase the configurations with the
smallest $\langle m\rangle$ have a dressing function that is essentially flat
and it is in this sense that the propagator is closer to a free field theory
propagator. The deviations from a $1/p^{2}$ propagator are larger in the
infrared region. Furthermore, these Figs. also show that it is closer to the
confinement-deconfinement transition that the hybrid Monte Carlo performs
worst, as this is where the propagators differ the most. See, in particular,
the propagators computed from the various simulations at $\beta=1.0125$.

##  V Summary and Conclusions

In this work the lattice Landau gauge photon propagator together with the
average number of Dirac strings is studied in the compact formulation of QED
for the pure gauge version of the theory. Following the procedure implemented
in Loveridge:2021qzs , the confined phase and deconfined phase propagators are
computed. Our results show a correlation between the two quantities that
clearly identify the confinement-deconfinement transition. The analysis of the
propagator data suggests that it is for $\beta\geq 1.0125$ that photon
propagator becomes closer to a free propagator.

The simulations also show that the functional form of the photon propagator is
different in each phase. In the confined phase the propagator is finite over
the full range of momenta, an indication that the theory has a mass gap for
low $\beta$. For the deconfined phase the photon propagator becomes compatible
with a divergent propagator in the infrared region, suggesting that it
reproduces a free field like propagator in the thermodynamic limit. Note that
in the present work we do not investigate the thermodynamic limit of the
theory but this conclusion comes from combining the new results with those of
Loveridge:2021qzs .

The differences observed in the photon propagator as a function of $\beta$
correlated well with the average number of Dirac strings present in the gauge
configurations, i.e. with the topological structure of the gauge group.
Indeed, the simulations show that the presence of a large number of Dirac
strings results in a theory with a mass gap that vanishes in the deconfined
phase. It is in the deconfined phase where there are essentially no Dirac
strings found in the gauge configurations.

The results discussed here also show that both the photon propagator and the
number of Dirac strings can be used to distinguish the two phases of the
lattice compact QED formulation. Furthermore, from the analysis of the data
for these two quantities, it seems that the transition between the confined
and deconfined phase is first order.

## Acknowledgments

This work was partly supported by the FCT – Fundação para a Ciência e a
Tecnologia, I.P., under Projects Nos. UIDB/04564/2020 and UIDP/04564/2020. P.
J. S. acknowledges financial support from FCT (Portugal) under Contract No.
CEECIND/00488/2017. The authors acknowledge the Laboratory for Advanced
Computing at the University of Coimbra (http://www.uc.pt/lca) for providing
access to the HPC resource Navigator.

## References

  * (1) K. G. Wilson, doi:10.1103/PhysRevD.10.2445
  * (2) T. Banks, R. Myerson and J. B. Kogut, Nucl. Phys. B 129, 493-510 (1977) doi:10.1016/0550-3213(77)90129-8
  * (3) J. Glimm and A. M. Jaffe, Commun. Math. Phys. 56, 195 (1977) doi:10.1007/BF01614208
  * (4) E. H. Fradkin and L. Susskind, Phys. Rev. D 17, 2637 (1978) doi:10.1103/PhysRevD.17.2637
  * (5) M. Creutz, L. Jacobs and C. Rebbi, doi:10.1103/PhysRevD.20.1915
  * (6) A. H. Guth, Phys. Rev. D 21, 2291 (1980) doi:10.1103/PhysRevD.21.2291
  * (7) B. E. Lautrup and M. Nauenberg, doi:10.1016/0370-2693(80)90400-1
  * (8) J. Fröhlich and T. Spencer, Commun. Math. Phys. 83, 411-454 (1982) doi:10.1007/BF01213610
  * (9) B. Berg and C. Panagiotakopoulos, Phys. Rev. Lett. 52, 94 (1984) doi:10.1103/PhysRevLett.52.94
  * (10) J. B. Kogut, E. Dagotto and A. Kocic, Phys. Rev. Lett. 60, 772 (1988) doi:10.1103/PhysRevLett.60.772
  * (11) G. Arnold, T. Lippert, K. Schilling and T. Neuhaus, Nucl. Phys. B Proc. Suppl. 94, 651-656 (2001) doi:10.1016/S0920-5632(01)01001-5 [arXiv:hep-lat/0011058 [hep-lat]].
  * (12) J. Jersak, C. B. Lang and T. Neuhaus, Phys. Rev. D 54, 6909-6922 (1996) doi:10.1103/PhysRevD.54.6909 [arXiv:hep-lat/9606013 [hep-lat]].
  * (13) M. Panero, JHEP 05, 066 (2005) doi:10.1088/1126-6708/2005/05/066 [arXiv:hep-lat/0503024 [hep-lat]].
  * (14) L. C. Loveridge, O. Oliveira and P. J. Silva, Phys. Rev. D 103, no.9, 094519 (2021) doi:10.1103/PhysRevD.103.094519 [arXiv:2102.01387 [hep-lat]].
  * (15) F. V. Gubarev, L. Stodolsky and V. I. Zakharov, Phys. Rev. Lett. 86, 2220-2222 (2001) doi:10.1103/PhysRevLett.86.2220 [arXiv:hep-ph/0010057 [hep-ph]].
  * (16) G. Arnold, B. Bunk, T. Lippert and K. Schilling, Nucl. Phys. B Proc. Suppl. 119, 864-866 (2003) doi:10.1016/S0920-5632(03)01704-3 [arXiv:hep-lat/0210010 [hep-lat]].
  * (17) P. Coddington, A. Hey, J. Mandula and M. Ogilvie, Phys. Lett. B 197, 191-194 (1987) doi:10.1016/0370-2693(87)90366-2
  * (18) A. Nakamura and M. Plewnia, Phys. Lett. B 255, 274-278 (1991) doi:10.1016/0370-2693(91)90247-N
  * (19) S. Durr and P. de Forcrand, Phys. Rev. D 66, 094504 (2002) doi:10.1103/PhysRevD.66.094504 [arXiv:hep-lat/0206022 [hep-lat]].
  * (20) A. M. Polyakov, Phys. Lett. B 59, 82-84 (1975) doi:10.1016/0370-2693(75)90162-8
  * (21) A. M. Polyakov, Nucl. Phys. B 120, 429-458 (1977) doi:10.1016/0550-3213(77)90086-4
  * (22) M. N. Chernodub, E. M. Ilgenfritz and A. Schiller, Phys. Rev. D 67, 034502 (2003) doi:10.1103/PhysRevD.67.034502 [arXiv:hep-lat/0208013 [hep-lat]].
  * (23) B. L. G. Bakker, M. N. Chernodub and A. I. Veselov, Phys. Lett. B 502, 338-344 (2001) doi:10.1016/S0370-2693(01)00180-0 [arXiv:hep-lat/0011062 [hep-lat]].
  * (24) J. Greensite, Lect. Notes Phys. 821, 1-211 (2011) doi:10.1007/978-3-642-14382-3
  * (25) D. Dudal, O. Oliveira and P. J. Silva, Annals Phys. 397, 351-364 (2018) doi:10.1016/j.aop.2018.08.019 [arXiv:1803.02281 [hep-lat]].
  * (26) A. G. Duarte, O. Oliveira and P. J. Silva, Phys. Rev. D 94, no.1, 014502 (2016) doi:10.1103/PhysRevD.94.014502 [arXiv:1605.00594 [hep-lat]].
  * (27) A. Cucchieri, D. Dudal, T. Mendes and N. Vandersickel, Phys. Rev. D 93, no.9, 094513 (2016) doi:10.1103/PhysRevD.93.094513 [arXiv:1602.01646 [hep-lat]].
  * (28) O. Oliveira and P. J. Silva, Phys. Rev. D 86, 114513 (2012) doi:10.1103/PhysRevD.86.114513 [arXiv:1207.3029 [hep-lat]].
  * (29) D. Dudal, O. Oliveira and N. Vandersickel, Phys. Rev. D 81, 074505 (2010) doi:10.1103/PhysRevD.81.074505 [arXiv:1002.2374 [hep-lat]].
  * (30) I. L. Bogolubsky, E. M. Ilgenfritz, M. Muller-Preussker and A. Sternbeck, Phys. Lett. B 676, 69-73 (2009) doi:10.1016/j.physletb.2009.04.076 [arXiv:0901.0736 [hep-lat]].
  * (31) A. Cucchieri and T. Mendes, PoS LATTICE2007, 297 (2007) doi:10.22323/1.042.0297 [arXiv:0710.0412 [hep-lat]].
  * (32) M. Q. Huber, Phys. Rev. D 101, 114009 (2020) doi:10.1103/PhysRevD.101.114009 [arXiv:2003.13703 [hep-ph]].
  * (33) M. Q. Huber, Phys. Rept. 879, 1-92 (2020) doi:10.1016/j.physrep.2020.04.004 [arXiv:1808.05227 [hep-ph]].
  * (34) A. K. Cyrol, L. Fister, M. Mitter, J. M. Pawlowski and N. Strodthoff, Phys. Rev. D 94, no.5, 054005 (2016) doi:10.1103/PhysRevD.94.054005 [arXiv:1605.01856 [hep-ph]].
  * (35) C. S. Fischer, A. Maas and J. M. Pawlowski, Annals Phys. 324, 2408-2437 (2009) doi:10.1016/j.aop.2009.07.009 [arXiv:0810.1987 [hep-ph]].
  * (36) A. C. Aguilar, D. Binosi and J. Papavassiliou, Phys. Rev. D 78, 025010 (2008) doi:10.1103/PhysRevD.78.025010 [arXiv:0802.1870 [hep-ph]].
  * (37) P. Boucaud, J. P. Leroy, A. Le Yaouanc, J. Micheli, O. Pene and J. Rodriguez-Quintero, JHEP 06, 099 (2008) doi:10.1088/1126-6708/2008/06/099 [arXiv:0803.2161 [hep-ph]].
  * (38) D. Dudal, J. A. Gracey, S. P. Sorella, N. Vandersickel and H. Verschelde, Phys. Rev. D 78, 065047 (2008) doi:10.1103/PhysRevD.78.065047 [arXiv:0806.4348 [hep-th]].
  * (39) R. Alkofer, W. Detmold, C. S. Fischer and P. Maris, Phys. Rev. D 70, 014014 (2004) doi:10.1103/PhysRevD.70.014014 [arXiv:hep-ph/0309077 [hep-ph]].
  * (40) J. A. Gracey, M. Peláez, U. Reinosa and M. Tissier, Phys. Rev. D 100, no.3, 034023 (2019) doi:10.1103/PhysRevD.100.034023 [arXiv:1905.07262 [hep-th]].
  * (41) O. Oliveira, P. J. Silva, J. I. Skullerud and A. Sternbeck, Phys. Rev. D 99, no.9, 094506 (2019) doi:10.1103/PhysRevD.99.094506 [arXiv:1809.02541 [hep-lat]].
  * (42) A. F. Falcão, O. Oliveira and P. J. Silva, Phys. Rev. D 102, no.11, 114518 (2020) doi:10.1103/PhysRevD.102.114518 [arXiv:2008.02614 [hep-lat]].
  * (43) A. C. Aguilar, F. De Soto, M. N. Ferreira, J. Papavassiliou and J. Rodríguez-Quintero, Phys. Lett. B 818, 136352 (2021) doi:10.1016/j.physletb.2021.136352 [arXiv:2102.04959 [hep-ph]].
  * (44) F. Gao, J. Papavassiliou and J. M. Pawlowski, Phys. Rev. D 103, no.9, 094013 (2021) doi:10.1103/PhysRevD.103.094013 [arXiv:2102.13053 [hep-ph]].
  * (45) P. M. Lo, K. Redlich and C. Sasaki, Phys. Rev. D 103, no.7, 074026 (2021) doi:10.1103/PhysRevD.103.074026 [arXiv:2101.12663 [hep-ph]].
  * (46) Y. Hayashi and K. I. Kondo, Phys. Rev. D 103, no.11, L111504 (2021) doi:10.1103/PhysRevD.103.L111504 [arXiv:2103.14322 [hep-th]].
  * (47) C. S. Fischer and M. Q. Huber, Phys. Rev. D 102, no.9, 094005 (2020) doi:10.1103/PhysRevD.102.094005 [arXiv:2007.11505 [hep-ph]].
  * (48) D. Dudal, O. Oliveira, M. Roelfs and P. Silva, Nucl. Phys. B 952, 114912 (2020) doi:10.1016/j.nuclphysb.2019.114912 [arXiv:1901.05348 [hep-lat]].
  * (49) A. Nakamura and R. Sinclair, Phys. Lett. B 243, 396-402 (1990) doi:10.1016/0370-2693(90)91403-X
  * (50) T. A. DeGrand and D. Toussaint, doi:10.1103/PhysRevD.22.2478
  * (51) S. Duane, A. D. Kennedy, B. J. Pendleton and D. Roweth, Phys. Lett. B 195, 216-222 (1987) doi:10.1016/0370-2693(87)91197-X
  * (52) R. G. Edwards et al. [SciDAC, LHPC and UKQCD], Nucl. Phys. B Proc. Suppl. 140, 832 (2005) doi:10.1016/j.nuclphysbps.2004.11.254 [arXiv:hep-lat/0409003 [hep-lat]].
  * (53) V. Grosch, K. Jansen, J. Jersak, C. B. Lang, T. Neuhaus and C. Rebbi, Phys. Lett. B 162, 171-175 (1985) doi:10.1016/0370-2693(85)91081-0
  * (54) W. Kerler, C. Rebbi and A. Weber, Phys. Lett. B 348, 565-570 (1995) doi:10.1016/0370-2693(95)00188-Q [arXiv:hep-lat/9501023 [hep-lat]].
  * (55) E. M. Ilgenfritz, C. Menz, M. Muller-Preussker, A. Schiller and A. Sternbeck, Phys. Rev. D 83, 054506 (2011) doi:10.1103/PhysRevD.83.054506 [arXiv:1010.5120 [hep-lat]].
  * (56) D. B. Leinweber et al. [UKQCD], Phys. Rev. D 60, 094507 (1999) [erratum: Phys. Rev. D 61, 079901 (2000)] doi:10.1103/PhysRevD.60.094507 [arXiv:hep-lat/9811027 [hep-lat]].
  * (57) G. T. R. Catumba, O. Oliveira and P. J. Silva, Phys. Rev. D 103, no.7, 074501 (2021) doi:10.1103/PhysRevD.103.074501 [arXiv:2101.04978 [hep-lat]].
  * (58) M. Pippig, SIAM J. Sci. Comput. 35, C213 (2013)
  * (59) http://www.uc.pt/lca
  * (60) J. Jersak, [arXiv:hep-lat/0010014 [hep-lat]].
  * (61) V. G. Bornyakov, V. K. Mitrjushkin, M. Muller-Preussker and F. Pahl, Phys. Lett. B 317, 596-603 (1993) doi:10.1016/0370-2693(93)91378-Z [arXiv:hep-lat/9307010 [hep-lat]].
  * (62) I. L. Bogolubsky, V. K. Mitrjushkin, M. Muller-Preussker and P. Peter, Phys. Lett. B 458, 102-108 (1999) doi:10.1016/S0370-2693(99)00575-4 [arXiv:hep-lat/9904001 [hep-lat]].
  * (63) I. L. Bogolubsky, L. Del Debbio and V. K. Mitrjushkin, Phys. Lett. B 463, 109-116 (1999) doi:10.1016/S0370-2693(99)00968-5 [arXiv:hep-lat/9903015 [hep-lat]].
  * (64) L. Giusti, M. L. Paciello, S. Petrarca, B. Taglienti and M. Testa, Phys. Lett. B 432, 196-202 (1998) doi:10.1016/S0370-2693(98)00631-5 [arXiv:hep-lat/9803021 [hep-lat]].
  * (65) A. Cucchieri and T. Mendes, [arXiv:hep-lat/9902024 [hep-lat]].
  * (66) A. Cucchieri and F. Karsch, Nucl. Phys. B Proc. Suppl. 83, 357-359 (2000) doi:10.1016/S0920-5632(00)91672-4 [arXiv:hep-lat/9909011 [hep-lat]].
  * (67) H. Nakajima, S. Furui and A. Yamaguchi, [arXiv:hep-lat/0007001 [hep-lat]].
  * (68) I. L. Bogolubsky and V. K. Mitrjushkin, [arXiv:hep-lat/0204006 [hep-lat]].
  * (69) T. A. DeGrand and D. Toussaint, Phys. Rev. D 24, 466 (1981) doi:10.1103/PhysRevD.24.466
  * (70) G. Cella, U. M. Heller, V. K. Mitrjushkin and A. Vicere, Phys. Rev. D 56, 3896-3902 (1997) doi:10.1103/PhysRevD.56.3896 [arXiv:hep-lat/9704012 [hep-lat]].
  * (71) R. Aouane, V. G. Bornyakov, E. M. Ilgenfritz, V. K. Mitrjushkin, M. Muller-Preussker and A. Sternbeck, Phys. Rev. D 85, 034501 (2012) doi:10.1103/PhysRevD.85.034501 [arXiv:1108.1735 [hep-lat]].
  * (72) P. J. Silva, O. Oliveira, P. Bicudo and N. Cardoso, Phys. Rev. D 89, no.7, 074503 (2014) doi:10.1103/PhysRevD.89.074503 [arXiv:1310.5629 [hep-lat]].
  * (73) P. J. Silva and O. Oliveira, Phys. Rev. D 93, no.11, 114509 (2016) doi:10.1103/PhysRevD.93.114509 [arXiv:1601.01594 [hep-lat]].
  * (74) A. Cucchieri, T. Mendes and A. R. Taurines, Phys. Rev. D 71, 051902 (2005) doi:10.1103/PhysRevD.71.051902 [arXiv:hep-lat/0406020 [hep-lat]].
  * (75) P. J. Silva and O. Oliveira, AIP Conf. Proc. 892, no.1, 529 (2007) doi:10.1063/1.2714468
  * (76) D. Dudal, O. Oliveira and P. J. Silva, Phys. Rev. D 89, no.1, 014010 (2014) doi:10.1103/PhysRevD.89.014010 [arXiv:1310.4069 [hep-lat]].
  * (77) O. Oliveira, A. G. Duarte, D. Dudal and P. J. Silva, Few Body Syst. 58, no.2, 99 (2017) doi:10.1007/s00601-017-1269-3 [arXiv:1611.03642 [hep-lat]].
  * (78) J. Horak, J. Papavassiliou, J. M. Pawlowski and N. Wink, [arXiv:2103.16175 [hep-th]].
  * (79) S. W. Li, P. Lowdon, O. Oliveira and P. J. Silva, Phys. Lett. B 803, 135329 (2020) doi:10.1016/j.physletb.2020.135329 [arXiv:1907.10073 [hep-th]].
  * (80) D. Binosi and R. A. Tripolt, Phys. Lett. B 801, 135171 (2020) doi:10.1016/j.physletb.2019.135171 [arXiv:1904.08172 [hep-ph]].
  * (81) R. A. Tripolt, J. Weyrich, L. von Smekal and J. Wambach, Phys. Rev. D 98, no.9, 094002 (2018) doi:10.1103/PhysRevD.98.094002 [arXiv:1807.11708 [hep-ph]].
  * (82) A. K. Cyrol, J. M. Pawlowski, A. Rothkopf and N. Wink, SciPost Phys. 5, no.6, 065 (2018) doi:10.21468/SciPostPhys.5.6.065 [arXiv:1804.00945 [hep-ph]].


# On the use of PN Ranging with High-rate Spectrally-efficient Modulations

B. Ripani  Politecnico di Torino, Corso Duca degli Abruzzi 24, Torino, Italy
{barbara.ripani, roberto.garello, gabriel.maiolinicapez,
<EMAIL_ADDRESS>Corresponding Author  A. Modenini  European Space
Agency, Keplerlaan 1, Noordwijk, The Netherlands
<EMAIL_ADDRESS>R. Garello  Politecnico di Torino, Corso Duca degli
Abruzzi 24, Torino, Italy
{barbara.ripani, roberto.garello, gabriel.maiolinicapez,
<EMAIL_ADDRESS>G. Maiolini Capez  Politecnico di Torino, Corso
Duca degli Abruzzi 24, Torino, Italy
{barbara.ripani, roberto.garello, gabriel.maiolinicapez,
<EMAIL_ADDRESS>G. Montorsi  Politecnico di Torino, Corso Duca
degli Abruzzi 24, Torino, Italy
{barbara.ripani, roberto.garello, gabriel.maiolinicapez,
<EMAIL_ADDRESS>
###### Abstract

In this paper, we study the feasibility of coupling the PN ranging with
filtered high-order modulations, and investigate the simultaneous demodulation
of a high-rate telemetry stream while tracking the PN ranging sequence.
Accordingly, we design a receiver scheme that is able to perform a parallel
cancellation, in closed-loop, of the ranging and the telemetry signal
reciprocally. From our analysis, we find that the non-constant envelope
property of the modulation causes an additional jitter on the PN ranging
timing estimation that, on the other hand, can be limited by properly sizing
the receiver loop bandwidth.

Our study proves that the use of filtered high-order modulations combined with
PN ranging outperforms the state-of-the-art in terms of spectral efficiency
and achievable data rate, while having comparable ranging performance.

Keywords: PN Ranging, GMSK, high-order modulations, SRS, Chip tracking loop.

##  1 Introduction

Space agencies are currently replacing transparent ranging systems with
regenerative schemes to meet future Space Research (SR) mission needs. For
this purpose, the Consultative Committee for Space Data Systems (CCSDS)
defined the Pseudo-Noise (PN, [1, 2]) ranging standard, a state-of-the-art
positioning technique. This standard, already implemented in the ESA missions
BepiColombo and Solar Orbiter, demonstrated an accuracy as high as a few
centimetres during its validation [3, 4]. Additionally, the possibility of
coupling the PN ranging with a Gaussian Minimum Shift Keying (GMSK, [5])
modulated signal allows simultaneous reception of high-rate telemetry streams
and orbit determination.

However, Near Earth SR missions transmitting telemetry in X-Band, in agreement
with frequency assignment regulations [6], are constrained to 10 MHz of
bandwidth. Thus, when using GMSK, these missions are limited to a data rate of
8.7 Mbps, penalizing payload data generation (see the case of GAIA mission
[7]). To exceed this bound, it is clear that high-order modulation schemes
must be adopted. However, their compatibility with PN ranging (to the best of
the authors’ knowledge) was never investigated.

In this paper, we study the coupling of PN ranging with filtered (amplitude)
Phase-Shift Keying (PSK/APSK) signals, by taking the CCSDS standard 131.2 [8,
9] as a reference. Namely, we define a receiver scheme that is able to
simultaneously demodulate the PSK/APSK stream, with modulation orders as high
as 64, and track the PN ranging sequence. This is done by performing a
parallel closed-loop cancellation of the PN ranging sequence over the
telemetry symbols and vice-versa. To evaluate the performance of the system we
designed, we assess both the Bit Error Rate (BER) and the ranging timing
jitter, revealing the presence of a noise floor due to the non-constant
envelope of PSK/APSK modulations. By modeling this phenomenon we derive a
closed-form expression that constitutes an upper-bound for the ranging jitter.

In this work, we find the performance of the proposed scheme to be comparable
to the classical GMSK approach, while enabling much higher data rates, thus
paving the way towards a new generation of SR missions with more ambitious
scientific objectives.

##  2 System Model

We consider the transmission system described by the block diagram reported in
Figure 1. The complex base-band modulated telemetry signal, with unitary
power, is defined as

$x_{\textrm{TM}}(t)=\sum_{k}{a_{k}p(t-kT)}\,,$ (1)

where $p(t)$ is the Square Root Raised Cosine (SRRC) shaping pulse, and
$\\{a_{k}\\}$ the sequence of telemetry symbols, each transmitted at symbol
time $kT$, and belonging to a complex PSK/APSK constellation. Instead, the
phase-modulated ranging signal $r_{\textrm{RG}}(t)$ is defined as
$r_{\textrm{RG}}(t)=e^{j\Phi_{\textrm{RG}}(t)}$ and its phase can be expressed
as

$\Phi_{\textrm{RG}}(t)=m_{\textrm{RG}}\sum_{k}{c_{k}h(t-kT_{c})}\,,$ (2)

where $m_{\textrm{RG}}$ is the ranging modulation index, $h(t)$ the sinusoidal
chip shaping pulse defined in [1], and $\\{c_{k}\\}$ the sequence of chips
transmitted at chip time $kT_{c}$. Note that, in line with CCSDS
recommendation 2.4.22A [10], we imposed $T_{c}\neq nT$, with
$n\in\mathbb{Z}^{+}$.

Figure 1: Transmitter block diagram.

Therefore, by considering a transmission over the Additive White Gaussian
Noise (AWGN) channel, the received signal reads

$y(t)=x_{\textrm{TM}}(t)\cdot r_{\textrm{RG}}(t-\tau_{\textrm{RG}})+w(t)\,,$
(3)

where $w(t)$ is white Gaussian noise with spectral density equal to $N_{0}$
and $\tau_{\textrm{RG}}\in[-\frac{T_{c}}{2},\frac{T_{c}}{2})$ is an arbitrary
ranging timing delay.

##  3 Receiver for PSK/APSK-modulated signals and PN Ranging

We implement a modified version of the receiver defined in [11] that performs
parallel closed-loop cancellation of telemetry and ranging signals from the
received stream. The receiver’s block diagram is shown in Figure 2.

The ranging component is removed from the received signal by using a locally
generated replica of the ranging signal. Mathematically,

Figure 2: Receiver block diagram.

$y(t)r^{*}_{\textrm{RG}}(t-\hat{\tau}_{\textrm{RG}})=x_{\textrm{TM}}(t)+w(t)r^{*}_{\textrm{RG}}(t-\hat{\tau}_{\textrm{RG}})\,,$
(4)

being $r^{*}_{\textrm{RG}}(t)$ the ranging signal’s complex conjugate, locally
generated using the estimated time delay $\hat{\tau}_{\textrm{RG}}$. It is
easy to see that, being $r_{\textrm{RG}}(t)$ a phasor, the resultant process
$w(t)r^{*}_{\textrm{RG}}(t-\hat{\tau}_{\textrm{RG}})$ is still white Gaussian,
and $x_{\textrm{TM}}(t)$ is perfectly recovered as long as
$\hat{\tau}_{\textrm{RG}}=\tau_{\textrm{RG}}$, since
$|r_{\textrm{RG}}(t-\tau_{\textrm{RG}})|^{2}=1$. The telemetry stream thus
obtained is demodulated to extract data, and then re-modulated to obtain the
sequence of symbols. With it, we generate the complex conjugate
$x^{*}_{\textrm{TM}}(t)$, necessary to sub-optimally cancel the telemetry
component from the received signal. We can write

$y(t)x^{*}_{\textrm{TM}}(t)=r_{\textrm{RG}}(t-\tau_{\textrm{RG}})|x_{\textrm{TM}}(t)|^{2}+\mathcal{W}(t)\,,$
(5)

where $\mathcal{W}(t)=x^{*}_{\textrm{TM}}(t)w(t)$ can be proved111It is
pointed out that the signal $x_{{\textrm{TM}}}(t)$ is actually a
cyclostationary process, thus making $\mathcal{W}(t)$ cyclostationary as well.
However, without loss of generality, in this paper we consider their
stationary statistics by simply referring to the cyclic auto-correlation at
null cyclic frequency [12]. to be a white process with power spectral density
equal to $N_{0}$. Since $x_{\textrm{TM}}(t)$ is SRRC filtered,
$|x_{\textrm{TM}}(t)|^{2}\neq 1$, and thus perfect cancellation cannot be
achieved. However, taking into account that
$\mathbb{E}[|x_{\textrm{TM}}(t)|^{2}]=1$, sub-optimal cancellation can be
performed by averaging the samples of the telemetry signal’s time-varying
envelope, as shown in Section 4.

Once recovered, the ranging sequence is input to the Chip Tracking Loop (CTL),
whose block diagram is depicted in Figure 3. Since the ranging chip sequence
resembles a clock signal, we implement the CTL as a modified version of the
Data Transition Tracking Loop (DTTL) [2, 13]. In this case, the mid-phase
integrator operates according to an Integrate and Dump approach, integrating
the PN sequence over a time interval $T_{c}$ between two consecutive chips.
Then, the mid-phase integrator output is multiplied by the PN transition
sequence to adjust the sign of the error at the input of the loop filter.

Finally, the estimated CTL timing $\hat{\tau}_{\textrm{RG}}$ is fed back to
the local PN generator in Figure 2, thus closing the receiver loop.

Figure 3: Chip tracking loop block diagram.

##  4 Ranging Timing Jitter

In this section we analyze the timing jitter of the receiver scheme presented
in Section 3.

We define $\mathcal{P}(t)=|x_{\textrm{TM}}(t)|^{2}-1$ as a zero-mean signal
with variance $\sigma^{2}_{P}=\mathbb{E}[|x_{\textrm{TM}}(t)|^{2}]$, which
depends on the adopted telemetry modulation and is equal to the values
reported in Table 1.

Table 1: Telemetry power variance $\sigma^{2}_{P}$ for different PSK/APSK modulations and roll-off factors. Roll-off | QPSK | 8-PSK | 16-APSK | 32-APSK | 64-APSK
---|---|---|---|---|---
0.2 | 0.240 | 0.243 | 0.435 | 0.560 | 0.483
0.25 | 0.224 | 0.226 | 0.426 | 0.554 | 0.474
0.3 | 0.209 | 0.211 | 0.418 | 0.551 | 0.467
0.35 | 0.197 | 0.199 | 0.412 | 0.549 | 0.462

Thus, Equation (5) can be re-written as

$y(t)x^{*}_{\textrm{TM}}(t)=r_{\textrm{RG}}(t-\tau_{\textrm{RG}})+\mathcal{W}(t)+r_{\textrm{RG}}(t-\tau_{\textrm{RG}})\mathcal{P}(t)\,,$
(6)

where the last term can be understood as additional noise arising from the
non-constant envelope of the considered modulation formats. To assess its
impact, we analyze the output of the CTL. As Figure 3 shows, the CTL extracts
the imaginary component of the ranging signal, which reads

$\alpha\Phi_{\textrm{RG}}(t-\tau_{\textrm{RG}})+\mathcal{N}(t)+\mathcal{D}(t)\,,$
(7)

where $\mathcal{N}(t)$ is the imaginary component of the Gaussian noise,
$\mathcal{D}(t)=\alpha\Phi_{\textrm{RG}}(t-\tau_{\textrm{RG}})\mathcal{P}(t)$,
and222Considering sine-shaped chip pulses, applying the Jacobi-Anger expansion
we have $r_{\textrm{RG}}(t)\approx
J_{0}(m_{\textrm{RG}})+j2J_{1}(m_{\textrm{RG}})\Phi_{\textrm{RG}}(t)$ [11].
$\alpha\approx 2J_{1}(m_{\textrm{RG}})$.

Therefore, we define
$n\triangleq\int_{0}^{T_{c}}{\mathcal{N}(t)\;\textrm{d}t}$ and
$\zeta\triangleq\int_{0}^{T_{c}}{\mathcal{D}(t)}\;\textrm{d}t$, i.e, the
output of the mid-phase integrator when $\mathcal{N}(t)$ and $\mathcal{D}(t)$
are input, respectively. While it is easy to see that
$\mathbb{E}[|n|^{2}]=T_{c}N_{0}$, a closed-form solution for
$\mathbb{E}[|\zeta|^{2}]$ is not easy to obtain but we derive its upper-bound.
In fact, by considering that the auto-correlation of $\mathcal{D}(t)$ is
always lower than $\mathbb{E}[|\mathcal{D}(t)|^{2}]$

$\mathbb{E}[|\zeta|^{2}]\leq\int_{0}^{T_{c}}{\int_{0}^{T_{c}}{\mathbb{E}[|\mathcal{D}(t)|^{2}]\;\textrm{d}t_{1}\textrm{d}t_{2}}}=P\sigma^{2}_{P}T_{c}^{2}\,,\
$ (8)

with $P=\mathbb{E}[|\alpha\Phi_{\textrm{RG}}(t)|^{2}]$, i.e., the useful power
of the received ranging signal at the input of the CTL.

We assume $\zeta$ to be a white process. Thus, we derive a linearized model of
the CTL, as reported in Figure 4, in which both the discretized noise terms at
the $k$-th time ($n_{k}$ and $\zeta_{k}$) are added to the phase detector
output $\varepsilon_{k}$.

Figure 4: Linearized CTL model.

The loop error is defined as

$\eta_{k}=\varepsilon_{k}+\frac{n_{k}}{K_{\varepsilon}}+\frac{\zeta_{k}}{K_{\varepsilon}}\,,$
(9)

where $K_{\varepsilon}=2\sqrt{2P}$ is the mid-phase integrator gain [2].
Considering that $n_{k}$ and $\zeta_{k}$ are uncorrelated, the loop error
variance is
$\sigma^{2}_{\eta}=(\mathbb{E}[|n_{k}|^{2}]+\mathbb{E}[|\zeta_{k}|^{2}])/K^{2}_{\varepsilon}$.
Thus, by using Equation (8), and filtering the loop error $\eta_{k}$ through
the Loop Filter with transfer function $H(f)$, the ranging timing variance is
upper-bounded by

$\sigma^{2}_{\tau}\leq
T_{c}\int_{0}^{1/T_{c}}{\left|\frac{K_{\varepsilon}H(f)}{1+K_{\varepsilon}H(f)}\right|^{2}\cdot
T_{c}\left(\frac{N_{0}}{8P}+\frac{\sigma^{2}_{P}T_{c}}{8}\right)df}\,.$ (10)

Recalling the definition of the noise bandwidth $B_{L}$ [14] from Equation
(10) we obtain

${\sigma}^{2}_{\tau}\leq
B_{L}T_{c}^{2}\left(\frac{N_{0}}{8P}+\frac{\sigma^{2}_{P}T_{c}}{8}\right)\,.$
(11)

It is interesting to observe that the first noise RHS term in Equation (11)
coincides with the theoretical PN ranging jitter [2], and is reduced as
$P/N_{0}B_{L}$ grows. However, due to the modulation’s non-constant envelope,
the second RHS term appears. As it is directly proportional to the product
$\sigma^{2}_{P}B_{L}$, the only way to reduce this noise term is by narrowing
the noise loop bandwidth $B_{L}$, i.e., by averaging the samples of the
varying envelope.

##  5 Numerical Results

We simulate the transmission of a SRRC filtered telemetry signal, with channel
symbol rate $R=1/T=4.2$ Msymbol/s, together with a PN ranging signal of rate
$R_{c}=1/T_{c}=3$ Mchip/s.

First, we analyze the occupied bandwidth in an AWGN channel (as 99-percent of
the signal power): in absence of ranging, the telemetry signal has
approximately 4.8 MHz of bandwidth. Instead, in the presence of the PN signal,
the spectral occupancy changes as a function of the ranging modulation index
$m_{\textrm{RG}}$ as shown in Table 2. For comparison, we include the
bandwidth of a GMSK modulated signal (with the same bitrate of a QPSK
modulated signal). It can be seen that, even in presence of ranging, PSK/APSK
modulations are more spectrally-efficient than GMSK. Thus, when subject to
spectral regrowth due to power amplification (not considered in this paper),
the two schemes will have similar bandwidths.

Table 2: Occupied bandwidth for PSK/APSK (SRRC filtered with roll-off 0.35) and GMSK modulated telemetry, over the AWGN channel, as a function of the ranging modulation index $m_{\textrm{RG}}$. $m_{\textrm{RG}}$ | PSK/APSK | GMSK
---|---|---
0.111 | 4.87 MHz | 7.22 MHz
0.222 | 5.04 MHz | 7.31 MHz
0.444 | 6.38 MHz | 7.72 MHz
0.666 | 7.31 MHz | 8.23 MHz

Next, we analyze the ranging jitter in the ideal condition of perfectly
demodulated telemetry symbols. To do so, we implement and simulate a genie-
aided receiver that is provided with the correct sequence of symbols
$\\{a_{k}\\}$. We evaluate the timing jitter $\sigma^{2}_{\tau}$, normalized
to $T_{c}^{2}$, as a function of the signal-to-noise ratio $P/N_{0}B_{L}$ and
for different values of $B_{L}$. Figure 5 shows the simulation results for
SRRC filtered (roll-off = 0.2) QPSK (Figure 5(a)) and 64-APSK (Figure 5(b))
modulations. For comparison, the upper-bound of Equation (11) is shown
together with the theoretical PN ranging jitter derived in [2]. As predicted,
the jitter has a floor that decreases with $B_{L}$, and thus the performance
gets closer to the theoretical jitter. In particular, a value $B_{L}=150$ Hz
provides almost ideal performance, without limiting the timing synchronization
dynamic. Additionally, Figure 5 shows that as the modulation order decreases,
the upper-bound is tighter to the actual jitter.

(a) QPSK

(b) 64-APSK

Figure 5: Impact of SRRC-filtered modulations on timing jitter in a genie-
aided receiver (roll-off = 0.2).

Finally, we evaluate the end-to-end performance of the receiver scheme of
Figure 2 from the telemetry and ranging points of view, simulating the
complete communication chain. Figure 6 and 7 present the ranging timing jitter
and the telemetry’s BER as a function of $P/N_{0}B_{L}$ and $E_{b}/N_{0}$
(being $E_{b}$ the energy per bit), respectively, for all SRRC-filtered
PSK/APSK modulations (roll-off = 0.2). For the two modulation sets the same
trend is observed: despite an initial deviation (more pronounced as the
modulation order grows), the simulated ranging jitter curves (Figure 6(a) and
7(a)) coincide with the genie-aided curve, which is free of telemetry losses.
Although the telemetry and ranging performances are inherently intertwined,
the worsened jitter performance for low values of $P/N_{0}B_{L}$, with respect
to ideal cancellation, is not reflected in the BER. In fact, the end-to-end
BER is comparable to the modulations’ theoretical performance in an AWGN
channel in the absence of ranging. (Figure 6(b) and 7(b)).

(a) Ranging timing jitter

(b) Telemetry BER

Figure 6: End-to-end ranging and telemetry performance for QPSK and 8-PSK
modulations ($B_{L}=1.5$ kHz).

(a) Ranging Timing Jitter

(b) Telemetry BER

Figure 7: End-to-end ranging and telemetry performance for 16-APSK, 32-APSK,
and 64-APSK modulations ($B_{L}=1.5$ kHz).

##  6 Conclusions

In this paper, we investigated the possibility of combining the PN ranging
with high-order, spectrally-efficient PSK/APKS modulations. At the receiver
side, we designed, implemented and simulated a receiver scheme for
simultaneous demodulation of high-order telemetry symbols and tracking of the
received ranging sequence for orbit determination. The obtained results show
that it is possible to achieve ranging performance comparable to current
state-of-the-art receiver architecture by simply lowering the noise loop
bandwidth $B_{L}$. Furthermore, the introduction of such modulations for
combined telemetry and ranging more than doubles/triples the current data rate
bound of 8.7 Mbps.

## Disclaimer

The view expressed herein can in no way be taken to reflect the official
opinion of the European Space Agency.

## References

  * [1] CCSDS 414.1-B-2, _Pseudo-Noise PN Ranging Systems_ , 2014, Available at http://public.ccsds.org/.
  * [2] CCSDS 414.0-G-2, _Pseudo-Noise PN Ranging Systems_ , 2014, Available at http://public.ccsds.org/.
  * [3] P. Cappuccio, V. Notaro, L. Iess, S. Asmar, J. Border, S. Ciarcia Sr, A. Di Ruscio, E. Montagnon, J. De Vicente, M. Mercolino _et al._ , “First results from cruise tests of the mercury orbiter radio science experiment (MORE) of ESA’s BepiColombo mission,” _AGUFM_ , vol. 2019, pp. P34A–04, 2019.
  * [4] M. Mascarello, G. Sessler, E. Vassallo, G. R. Iglesias, K. Kewin, M. Montagna, and L. Manso, “The Solar Orbiter X-Band TT&C New Features: GMSK with PN Regenerative Ranging and DDOR Semaphores Implementation,” in _2019 8th International Workshop on Tracking, Telemetry and Command Systems for Space Applications (TTC)_ , 2019, pp. 1–6.
  * [5] K. Murota and K. Hirade, “GMSK modulation for digital mobile radio telephony,” _IEEE Transactions on communications_ , vol. 29, no. 7, pp. 1044–1050, 1981.
  * [6] Space Frequency Coordination Group, “Use of the 8450-8500 MHz band for space research, category A,” 1998, REC SFCG 5-1R5.
  * [7] European Space Agency, “GAIA flight model payload data handling unit delivered and integrated,” Available at https://sci.esa.int/web/gaia/-/51777-09-gaia-flight-model-payload-data-handling-unit-delivered-and-integrated .
  * [8] CCSDS 130.11-G-1, _SCCC—Summary of definition and performance_ , 2019, Available at http://public.ccsds.org/.
  * [9] CCSDS 131.2-B-1, _Flexible advanced coding and modulation scheme for high rate telemetry applications_ , 2012, Available at http://public.ccsds.org/.
  * [10] CCSDS 401.0-B-31, _Radio Frequency and modulation systems - Part 1 Earth Stations and Spacecrafts_ , 2020, Available at http://public.ccsds.org/.
  * [11] CCSDS 413.1-G-1, _Simultaneous transmission of GMSK telemetry and PN Ranging_ , 2017, Available at http://public.ccsds.org/.
  * [12] A. Papoulis, _Probability, Random Variables and Stochastic Processes_ , New York, NY, 1991.
  * [13] Jet Propulsion Laboratory, California Institute of Technology, “Autonomous software-defined radio receivers for deep space applications,” JPL, Tech. Rep., 2006.
  * [14] F. M. Gardner, _Phaselock techniques_. John Wiley & Sons, 2005.