text:
\section{Introduction}
\label{intro}
\sloppy
Laser-plasma accelerators can produce ultra short electron bunches with duration in the range from femtosecond to picosecond. Both quasi-monoenergetic high-energy electrons (hundreds of MeV)~\cite{mirzaie_effect_2018,faure_experiments_2008,thaury_shock_2015} and high-charge, low-energy (few MeV) electrons with broadband spectrum~\cite{guillaume_physics_2015} have been obtained with 100 TW class lasers. The peculiarity of these sources is the ultra-high peak dose rate in the pulse (up to \SI{e10}{Gy/s}) that is orders of magnitude higher than conventional Linac (\({\sim\SI{e2}{Gy/s}}\))~\cite{ashraf_dosimetry_2020}. The impact of electrons delivered in such short pulses of ultra-high dose rates on living cells or tissues still need to be extensively explored. Recent results obtained with FLASH-RT~\cite{favaudon_ultrahigh_2014,vozenin_advantage_2019,vozenin_biological_2019} and laser-accelerated protons~\cite{raschke_ultra-short_2016,bayart_fast_2019} indicate that the temporal structure of high dose rate, pulsed irradiation may have an impact on the radiobiological response and ask for a deeper understanding of the phenomena triggered by high dose rate irradiation.\\ To date, Laser-Accelerated Electrons (LAEs) produced by commercially available J-class lasers running at 10 Hz have been used for radiation biology studies~\cite{andreassi_radiobiological_2016,laschinsky_radiobiological_2012,oppelt_comparison_2015}. The main limit of these beams is the shot-to-shot pointing fluctuation in the order of the beam divergence, that hinders the reach of the stability standards required in clinics. In this context, low-energy electrons (few MeV) driven by kHz lasers may represent an interesting alternative for radiobiology applications. Such laser-plasma accelerators have been developed in recent years~\cite{he_high_2013,salehi_mev_2017,guenot_relativistic_2017,faure_review_2019} and, to our knowledge, no radiobiology study in these conditions is reported in literature at the time of writing. The key asset of these beams is the high repetition rate, which boosts the mean beam current and allows averaging of the shot-to-shot fluctuations by integration of a large number of shots. From a dosimetric point of view, this translates into higher mean dose rates (Gy/s), compared to those obtained with {J-class} lasers (Gy/min), and into a higher stability of the dose distribution at the sample. \\In this article we present the dosimetric characterisation of a low-energy kHz laser-driven electron beam for radiation biology applications and report preliminary results on \textit{in-vitro} irradiation of HCT116 colorectal cancer cells. The dose distribution at the irradiated sample was measured at each irradiation with absolutely calibrated radiochromic films, which have been shown to be dose-rate independent over a wide range of dose rate up to \SI{e12}{Gy/s}~\cite{karsch_dose_2012,jaccard_high_2017,bazalova-carter_comparison_2015}. The use of an IBA Razor Nano Chamber provided a real time dose monitoring during irradiations and evaluation of the dose uncertainty in the post-processing analysis.

\section{Electron source and set-up}
\label{source and set up}
\sloppy
The experiment was performed with the \textit{Salle Noire} laser system of the \textit{Laboratoire d'Optique Appliquée}, delivering 3.5 fs pulses at 1 kHz repetition rate with \SI{3}{mJ} on target. The laser can generate low energy electrons in a reproducible way in a range from sub-MeV up to few MeV with charges of \(\sim\)\SI{10}{pC} and few pC respectively, depending on the  gas density and profile~\cite{gustas_high-charge_2018}. In the experiment, the laser pulse was focused with a parabola (\SI{50}{mm} focal length) onto an N\(_2\) gas-jet. The gas-jet density was optimised so as to maximise the accelerated charge per bunch, which corresponds to the production of an electron beam featuring a low-energy quasi-thermal spectrum up to \(\sim\)\SI{2}{MeV} and a divergence of \SI{70}{mrad}.
\begin{figure}[]

	\includegraphics[width=1\columnwidth]{figures/set-up_salleNoire.pdf}

	\caption{Schematic draft (horizontal cut, top view) of the experimental set-up.}
	\label{set-up_sallenoire}
\end{figure}
Fig.~\ref{spectra} shows the electron spectrum taken in four different days. As shown, the spectral shape can be reproduced from one day to another with good precision. The measured charge (total charge at the source) was \(9.8\)~pC/shot with an instability of \(\pm0.7\)~pC/shot (standard deviation).\\A draft of the experimental setup is shown in Fig.~\ref{set-up_sallenoire}. A tube was inserted in the vacuum chamber to place the irradiation site (in air) closer to the electron source, compatibly with the beam diagnostics. A \SI{100}{\micro\meter} thick vacuum-air Mylar window was placed \SI{38}{\centi\meter} far from the source. At this location, the transverse dimension of the beam was \SI{3}{\centi\meter} (FWHM), which allowed to irradiate a circular spot of  \SI{1}{\centi\meter} diameter with a sufficiently homogeneous dose. The cell sample was positioned vertically in air, \SI{2}{\centi\meter} after the Mylar window, with a \SI{800}{\micro\meter} thick plastic holder aligned to the beam-axis. \begin{figure}[]
	\centering

	\includegraphics[width=0.9\columnwidth]{figures/spectra_4days.pdf}

	\caption{Spectrum of the electron beam taken in four different days, obtained by integrating 15 shots. The lower limit of the spectrometer is \SI{400}{keV}.}
	\label{spectra}
\end{figure}
A customised metallic support allowed to hold in place the cell holder and to place it at the same depth in the tube at each irradiation. The support was also designed to house an IBA Razor Nano Chamber (RNC), which was inserted behind the cell sample. The chamber (\SI{2}{\milli\meter} outer electrode diameter, \SI{0.003}{\centi\meter\squared} active volume) was used as on-line dose monitor during irradiations and as day-to-day dose reference during the optimisation of the interaction conditions between the laser and the gas-jet prior to an experimental run. A radiochromic EBT3 film was placed on the front side (in beam's view) of the cell holder at each irradiation for dosimetry. A motorised set-up, located in the vacuum chamber between the source and the Mylar window, was used for beam imaging and spectroscopy and was removed during irradiations. It contains a removable dipole (\SI{0.058}{T}) and a phosphor screen coupled with a CCD camera. Imaging of the beam transverse section, performed by removing the dipole magnet from the beam axis, allowed to align the electron beam to the cell sample and to measure the total accelerated charge per bunch.  The details of these diagnostics can be found in D. Gustas (2019)~\cite{gustas_high-charge_2018}.
\section{Dosimetry}
\sloppy
The dose was measured with an EBT3 radiochromic film placed on the cell plastic support at each irradiation. The films were scanned 24~hours after irradiation, as recommended by Devic et al. (2016)~\cite{devic_reference_2016}, with an EPSON V800 flatbed scanner and a resolution of 300 DPI. The dose was measured within a circular ROI of \SI{10}{mm} diameter aligned with the cell sample. The EBT3 film batch was calibrated at the Elekta Synergy linear accelerator of the \textit{University Hospital of Lausanne} (CHUV, Switzerland) with \SI{4}{MeV} electrons at a dose rate of \SI{6}{Gy/min}. The calibration curve is shown in Fig.~\ref{calib ebt3 sallenoire}, together with the five degree polynomial expression used to fit the experimental points.
\begin{figure}[]
	\centering

	\includegraphics[width=1\columnwidth]{figures/calib_ebt3_chuv.pdf}

	\caption{Calibration of the EBT3 film batch (red channel) with the \SI{4}{MeV} electron beam generated by the Synergy Elekta LINAC available at the University Hospital of Lausanne (CHUV, Switzerland). The experimental values (blue points) are fitted with a five degree polynomial (black curve), whose expression is reported in the image.}
	\label{calib ebt3 sallenoire}
\end{figure}\\
The film dose was corrected to obtain the dose delivered in the cells by measuring the attenuation caused by the \SI{800}{\micro\meter} cell plastic holder and the air gap. Precisely, the dose delivered to the cell is related to the dose absorbed by the film according to:
\begin{equation}
D_{cells}=D_{EBT3}\cdot R_{cells/EBT3}
\end{equation}
where \(R_{cells/EBT3}\) is the ratio between the dose absorbed by the cells and by the EBT3. This factor was measured with a double film arrangement, as shown in the insert of Fig.~\ref{double film dose}. A film was placed on the plastic holder surface, as during irradiations, and a second film was placed inside the plastic holder at the cell position. Since the ratio depends on the spectral shape of the beam, which varied slightly during the day and between days (see Fig.~\ref{spectra}), the double film measurement was carried out before each irradiation series. The results of the double film measurements, performed in five different days, are reported in Fig.~\ref{double film dose}. We measured a factor \(R_{cells/EBT3}\) of 0.71 with good reproducibility during the day and between days (5.8\% standard deviation). A mean dose rate of 1.1 Gy/s and 0.78 Gy/s was measured on the holder surface and at the cell position, respectively, with an instability of 8\% (standard deviation). As anticipated in the introduction, the mean dose rate achieved with such a kHz laser-driven electron source is higher than typical dose rates reported in literature with J-class lasers running at \SI{10}{Hz} (\(\sim\)Gy/min)~\cite{gizzi_laserplasma_2015,andreassi_radiobiological_2016,laschinsky_radiobiological_2012,oppelt_comparison_2015,beyreuther_establishment_2010,nicolai_realizing_2014}. A rough estimate of the peak dose rate in the pulse can be obtained by considering the temporal stretch of the beam energy components between \SI{0.5}{MeV} and \SI{1}{MeV} at the cell position. The calculation leads to a pulse duration of \SI{120}{ps} and a peak dose rate in the order of \SI{e7}{Gy/s}.
\begin{figure}[]
	\centering

	\includegraphics[width=1\columnwidth]{figures/Dose_doubleFilm2.jpg}

	\caption{Dose rate measured at the surface of the plastic holder and at the cell position with the double film arrangement shown in the insert.}
	\label{double film dose}
\end{figure}
\subsection{Dose homogeneity and uncertainty}
\sloppy
Since the position of the irradiation site inserted in the tube was fixed, the electron beam pointing needed to be carefully aligned by adjusting the incidence angle of the laser on the gas-jet before the beginning of the irradiations. Once the beam pointing was aligned with the cell holder, the beam showed a remarkably good stability in terms of dose distribution between irradiations, as shown in Fig.~\ref{film images}. \begin{figure*}[]
	\centering

	\includegraphics[width=0.8\textwidth]{figures/all_films.pdf}

	\caption{Dose distribution measured with the EBT3 film placed on the cell plastic support in six consecutive irradiations consisting in a dose escalation performed for the survival essay. The colormap interval is set between the maximum value and the 50\% percent level in all images. The red circles represent the ROI of \SI{1}{cm} diameter used to measure the dose as well as the surface where the cell response was studied. }
	\label{film images}
\end{figure*}
The figure shows typical dose distributions measured with the EBT3 during a dose escalation performed for the survival assay. We were able to precisely align the electron beam with the cell sample and keep the dose inhomogeneity below \(\pm7\)\% in most of the irradiations (85\%).\\The uncertainty on the dose delivered at the cell sample was evaluated by cross correlation between the dose measured by the EBT3 film and by the RNC placed behind the sample. In particular, the instability of the ratio between the two measured doses is an indicator of the spectrum stability between irradiations and of the uncertainty on the ratio \(R_{cells/EBT3}\) used to quantify the dose delivered at cell sample. Fig.~\ref{ratios} shows the typical behaviour of the ratio between the dose absorbed by the film and by the RNC during three dose escalation series performed for the survival assay.
\begin{figure}[]
	\centering

	\includegraphics[width=1\columnwidth]{figures/ratios_2407p53_2907WT_p53.pdf}

	\caption{Ratio between the dose meaured by the EBT3 film and by the Razor Nano Chamber as a function of the exposure time, i.e. the shutter opening time.}
	\label{ratios}
\end{figure}
As it is shown, we observed an increase of the ratio with the exposure time. A similar trend was also observed for the dose rate, which turned out to slightly decrease with the exposure time. It is possible that the observed patterns are due to the shutter used to block the laser and control the irradiation time. Opening the shutter, placed before the last four optics, after it is kept closed to block the laser during the sample replacement might generate a transient thermal phenomenon in the optics and affect the laser properties at the gas-jet. This effect will be verified and solved in future experiment by placing the laser shutter after the last focusing optics. The variation of the ratio between irradiations, typically below \(\pm10\%\), results in a comparable uncertainty on the \(R_{cells/EBT3}\) factor measured with the double film arrangement and, as a consequence, on the dose delivered to the cell sample.

\section{Radiation biology experiment}
\label{sec:2}
\subsection{Methods}
\sloppy
We performed irradiation of \textit{in-vitro} cancer cells with the \textit{Salle Noire} electron source both to investigate the biological response to kHz laser-accelerated electrons and to validate the dosimetry. To this aim, we irradiated colorectal cancer cells HCT116, both the wild type (HCT116 WT) and its radioresistant counterpart (HCT116 p53\(^{-/-}\)) mutated for the tumour suppressor gene p53. The use of both the wild type and the radioresistant counterpart allowed to verify the quality of the dosimetry protocol and of the radiobiology analysis. The colorectal cancer cells HCT116, wild type (WT) or mutated for the tumour suppressor gene p53 (p53\(^{-/-}\)), were cultured and grown as monolayers in plastic tissue culture disposable flasks (TPP) in McCoy's 5A (Modified) Medium with GlutaMAXTM (ThermoFisher Scientific), supplemented with 10\% fetal calf serum (PAA) and 1\% penicillin and streptomycin (ThermoFisher Scientific).\\All cell lineages were grown at 37\degree~C in a humidified atmosphere of 5\% CO\(_2\) in air. Based on the transverse dose homogeneity deducted from radiochromic films, a \SI{1}{cm} diameter circular region was delimited on the inside face of the culture dishes with Creamy Color long lasting lip pencil (Kiko) in which \SI{4e4}{cells} were seeded in \SI{100}{\micro\liter} of medium. The culture dish was positioned vertically behind the exit window. After exposure to the laser-driven electron beam, 1 mL of medium was added and the cell monolayers were incubated for 4~hours in standard conditions. Cells were harvested with Accutase (Merck), dispatched into 3 different wells of 12-well plate (TPP) in \SI{2.5}{\milli\liter} of medium and grown for five generations corresponding to five days (replication time below 20~h). Appropriate control samples were treated under the same conditions, including bringing the cell culture dish in a vertical position as for irradiation.\\ After this period, cells were harvested with \SI{250}{\micro\liter} of Accutase inactivated with \SI{250}{\micro\liter} of 1X PBS (ThermoFisher Scientific) supplemented with 10\% fetal calf serum. The final volume was adjusted to \SI{1}{\milli\liter} with 1X PBS and \SI{200}{\micro\liter} of each well were dispatched into a non-sterile Ubottom 96-well plate (TPP). In each well, \SI{2}{\micro\liter} of a propidium iodide solution (Sigma, \SI{100}{\micro g/\milli\liter} in 1X PBS) were added just before counting by flow cytometry technology. Cell acquisition was performed using Guava\textsuperscript{\tiny\textregistered} and the data analysis carried out with GuavaSoft (Merck Millipore).
\subsection{Survival assay results}
\sloppy
\begin{figure}[]
	\centering

	\includegraphics[width=1\columnwidth]{figures/survie.pdf}

	\caption{Survival curves for the HCT116 WT and HCT116 p53\(^{-/-}\) cell lines obtained in two independent experiments. In both graphs the points represent the mean of three radiobiological analysis carried out on the same irradiated sample and the vertical error bar corresponds to the standard deviation of the survival fraction obtained in the three analysis. The horizontal error bars represent the sum of the dose error due to the inhomogeneity within the analysed surface of the sample (circular surface of \SI{1}{cm} diameter) and of the uncertainty on the dosimetry.}
	\label{survival}
\end{figure}
The survival curves obtained for the HCT116 WT and HCT116 p53\(^{-/-}\) cell lines grown as monolayers are reported in Fig.\ref{survival}. Both the error due to the dose homogeneity and to the dosimetry uncertainty are included in the horizontal error bars. The experimental results were fitted with the Linear-Quadratic (LQ) model: \(\ln S= \alpha D-\beta D^2\), where D is the delivered dose and \(\alpha\) and \(\beta\) are adjustable parameters. The surviving fraction decreases with the delivered dose and the curves exhibit the typical shoulder shape of low LET radiations. Furthermore, the curve of the WT cell line lies well below the curve of the more radioresistant counterpart p53\(^{-/-}\), which confirms the quality of both the dosimetry and the biological analysis.
\begin{table}[]
	\centering
	\caption{D\(_{10}\) values for \textit{in-vitro} HCT116 cell lines, grown as monolayers, obtained with laser-accelerated electrons (from the survival curves of Fig.~\ref{survival}) and with \SI{662}{keV} photons of a \(^{137}\)Cs source at a conventional dose-rate, taken from Pommarel et al. (2017)\cite{pommarel_spectral_2017}.}
	\label{D10}
	\renewcommand{\arraystretch}{1.5}

	\begin{tabular}{c|cc}
		Cell line   &
		LAE & Photons \(^{137}\)Cs 662 keV \\ \hline
		HCCT116 WT  & \(4.34\pm0.73\) Gy                                & \(3.8\pm0.4\) Gy           \\
		HCCT116 p53 & \(5.05\pm0.66\) Gy                                  & \(5\pm0.4\) Gy             \\ \hline
	\end{tabular}
\end{table}\\It is relevant to compare the 10\% survival dose D\(_{10}\) obtained with the LAE beam with those previously obtained with low-LET \SI{662}{keV} photons at a conventional dose rate of \SI{1.4}{Gy/min} produced by a \(^{137}\)Cs source (integrated in an IBL 637 n\(^\circ\)9418 device manufactured by CIS Bio)\cite{pommarel_spectral_2017}.
The D\(_{10}\) values, determined from the LQ fitting curves (Kaleidagraph software) for the two cell types (WT and p53\(^{-/-}\)) are reported in Table~\ref{D10}. The D\(_{10}\) obtained with LAEs are close to the values
obtained with the \(^{137}\)Cs source, which is in line with the results reported in literature showing no remarkable difference in the RBE of LAEs in comparison with conventional sources of comparable LET~\cite{laschinsky_radiobiological_2012}.
\section{Conclusions and outlook}
	\sisetup{
		range-phrase= -}
We presented the first dosimetric characterisation and application to radiation biology of a kHz laser-driven electron beam. The use of both radiochromic films and an on-line Razor Nano Chamber allowed evaluation of the dose distribution delivered at the sample at each irradiation and to evaluate the uncertainty on the delivered dose. We showed good day-to-day reproducibility of the irradiation conditions in terms of dose rate and beam quality. A dose rate of \SI{1}{Gy/s} was achieved, which is higher than typical dose rates reported in literature with J-class lasers running at \SI{10}{Hz}. Besides the increase of the mean dose rate, the high repetition rate ensured a remarkably good stability of the dose distribution at the sample and indicates that such LAE source can be a promising alternative to electrons produced by \SI{10}{Hz} lasers for radiation biology applications. During the dose escalation sequence, we observed an impact of the exposure time on the beam spectrum, which introduced an uncertainty on dosimetry. The observed trends are likely to be due to the use of a beam shutter in the laser chain generating transient thermal phenomena in the optics. This phenomenon, which is not attributable to inherent source instability, can be significantly reduced in future studies by placing the shutter at the end of the laser chain. In addition, an ongoing upgrade of the \textit{Salle Noire} facility will lead to an increase of the laser energy up to \SIrange[range-units = single]{4}{5}{mJ/pulse}. This would allow to generate more energetic electrons (with a comparable charge per bunch) and to reduce the sensitivity of the dosimetry protocol to spectral instabilities thanks to the higher penetration depth.\\Regarding the radiobiological outcome, the survival curves obtained with HCT116 monolayer samples indicate no significant difference between the radiobiological effectiveness of kHz laser-driven electrons and conventional sources with similar LET, in line with previous findings reported in literature. This result further validates both the dosimetry and the biological analysis. Future experiments should explore the response of healthy cells and of spheroids samples as well as the impact of different fractionation modalities on the biological response~\cite{bayart_fast_2019} to provide a more complete indication of the therapeutic potential of such high peak dose rate sources. Moreover, boosting the energies up to few MeV without a significant decrease of the charge per bunch would allow to irradiate \textit{in vivo} biological targets, such as mice, and would increase the peak dose rate in the pulse by reducing the pulse temporal stretch to few picoseconds.

\bibliographystyle{./bibliography/unsrt-5aut}


text:
\section*
\thesection\protect\indent \parbox[t]{11.715cm} {#1}}
\addcontentsline{toc}{section}{Appendix\thesection\ \ \ #1} }

\newcommand{\C}{\mathbb C}
\newcommand{\complexs}{{\bbs C}}
\newcommand{\reals}{{\bbs R}}
\newcommand{\zeds}{{\bbs Z}}
\newcommand{\rat}{{\bibitem Q}}
\newcommand{\mat}{{\bibitem M}}
\newcommand{\mats}{{\bbs M}}
\newcommand{\NO}{\,\mbox{$\circ\atop\circ$}\,}
\newcommand{\id}{{1\!\!1}}
\newcommand{{\mathbb 1}}{{\mathbb 1}}

\def\semiprod{{\supset\!\!\!\!\!\!\!{\times}~}}
\def{\cal A}{{\cal A}}
\def{\cal H}{{\cal H}}

\def\bra#1{\left\langle #1\right|}
\def\ket#1{\left| #1\right\rangle}
\def\hat\sigma#1#2{\left\langle #1\right|\left. #2\right\rangle}
\def\vev#1{\left\langle #1\right\rangle}

\font\mathsm=cmmi9
\font\mybb=msbm10 at 12pt

\def\nonumber{\nonumber}
\newcommand{\rf}[1]{(\ref{#1})}
\newcommand{\eq}[1]{Eq.~(\ref{#1})}
\hyphenation{pre-print} \hyphenation{pre-prints}
\hyphenation{di-men-sion-al} \hyphenation{di-men-sion-al-ly}
\def\begin{eqnarray}{\begin{eqnarray}}
\def\end{eqnarray}{\end{eqnarray}}

\newcommand{{\textstyle{1\over 2}}}{{\textstyle{1\over 2}}}

\newcommand{\partial}{\partial}

\newcommand{{\mathrm i}}{{\mathrm i}}
\newcommand{{\mathrm d}}{{\mathrm d}}

\begin{document}

\begin{titlepage}

\begin{center}

\baselineskip=24pt

{\Large\bf Gribov horizon in Noncommutative QED}

\baselineskip=14pt

\vspace{1cm}

Ozório Holanda$^1$, Marcelo S. Guimaraes$^2$,  Luigi Rosa$^{3,5}$ and  Patrizia Vitale$^{4,5}$
\\ [6mm]
$^1${\it Centro de Matemática, Computação e Cognição - Universidade Federal do ABC, Santo André, SP, Brazil}\\
[4mm]
$^2$ {\it Departamento de F\'{\i }sica Te\'{o}rica, Instituto de F\'{\i }sica, UERJ - Universidade do Estado do Rio de Janeiro.\\ Rua S\~{a}o Francisco Xavier 524, 20550-013 Maracan\~{a}, Rio de Janeiro, Brasil.}\\
[4mm]
$^3${\it Dipartimento di Matematica e Applicazioni ``R.~Caccioppoli'', Universit\`{a} di Napoli
{\sl Federico II}}\\
{\it Monte S.~Angelo, Via Cintia, 80126 Napoli, Italy}
\\[4mm]
$^4${\it Dipartimento di Fisica ``E.~Pancini', Universit\`{a} di Napoli
{\sl Federico II}}\\
{\it Monte S.~Angelo, Via Cintia, 80126 Napoli, Italy}
\\[4mm]
$^5${\it INFN, Sezione di Napoli}\\
{\it Monte S.~Angelo, Via Cintia, 80126 Napoli, Italy}
\vskip 1truecm
{\small\tt
netoholanda91@gmail.com, msguimaraes@uerj.br,  luigi.rosa@unina.it, patrizia.vitale@unina.it}

\end{center}

\begin{abstract}
It is known that  Noncommutative QED (NCQED) exhibits Gribov ambiguities in the Landau gauge. These ambiguities are related to zero modes of the Faddeev-Popov operator and arise in the ghost propagator when it has a pole. In this work, we establish a positive Faddeev-Popov operator for NCQED and the condition for the ghost propagator  not to have  poles, the so-called Gribov no-pole condition. This condition is implemented in the path integral and allows for the calculation of the photon propagator in  momentum space, which is dependent on the squared non-commutativity parameter. In the commutative limit, the standard QED is recovered.
\end{abstract}

\end{titlepage}

\section{Introduction}
The occurrence of so-called  Gribov ambiguities for Yang-Mills theories was first discussed in \cite{Gribov1978}, where Gribov showed that for non-Abelian gauge theories on flat topologically trivial space-times,   gauge fixing is problematic. Namely, it is not possible to choose for the gauge potential one representative on each gauge orbit, that is, by considering the quotient of the space of gauge connections with respect to the gauge group. Gribov ambiguity amounts to the fact that there are in general different field configurations which obey the same gauge-fixing condition, but which are related by a gauge transformation, i.e., they are on the same gauge orbit. As first shown by Singer \cite{Singer1978} and independently by Narasimhan and Ramadas \cite{Narasimhan1979}, this occurrence can be given a precise mathematical characterisation in terms of topological non-triviality of $\mathcal{G}$, the pertinent gauge group involved. On this basis,  QED is singled out, with the relevant homotopy group $\Pi_1(\widehat{U}(1))$ being trivial. For Quantum Chromodynamics (QCD), however, Gribov copies pose a real threat to the accurate description of the low-energy regime of the theory,  the gauge group of SU(N) Yang-Mills theories being nontrivial for $N\ge 2$. Over the last few decades, several solutions have been proposed in order to restrict the functional integral to the region in the space of  connections where only one representative for each gauge orbit is present, the so-called  first Gribov region. The most notable efforts are the Gribov-Zwanziger approach \cite{Zwanziger:1989mf} and its refined version \cite{Dudal:2008sp}.

In \cite{Canfora2016, Kurkov2018} the problem has been   addressed  in the context of gauge theories on noncommutative space-time. It has been shown that, differently from the commutative case, noncommutative QED (NCQED) exhibits Gribov ambiguity and  gauge configurations have been found, with an infinite number of copies. The goal of the present work is to investigate the characterisation of the  first Gribov region.  Following the same approach as in non-Abelian gauge theories \cite{Dudal:2008sp}, the ensuing constraint on the physical space of connections is implemented  by requiring that  the ghost propagator does not develop poles. This implies a modification of the photon propagator.

The paper is organised as follows. In section \ref{gaugeth},  we shortly review the problem in the framework of standard gauge theory. We then recall in section \ref{NCQED} the formulation of QED in the noncommutative setting with Moyal type noncommutativity and review the derivation of the equation for Gribov copies introduced in \cite{Canfora2016, Kurkov2018}. In section \ref{Gribovregion}, we apply the techniques developed in \cite{Dudal:2008sp} to characterise the first Gribov region for NCQED. In section \ref{five}, we evaluate the no-pole condition in the Landau gauge and find the Gribov region function $V(\Omega)$. Finally, we calculate the photon propagator in noncommutative momentum space. We conclude with some remarks and future perspectives.

\section{Gribov ambiguity in Yang-Mills}\label{gaugeth}

Let $M=\mathbb R^4$ be the space-time  manifold and let us define the group of gauge transformations,   $\mathcal{G}=\widehat{U}(N)$, as the smooth maps $ {g:M\rightarrow U(N)}$,  with boundary condition $  g(x) \rightarrow 1$ as $ |{ x}| \rightarrow \infty$.

 A pure gauge  theory of fundamental interactions, without matter fields, is a theory where the dynamical fields are the gauge connections,  locally represented by Lie algebra valued one-forms $A$, with   $F= dA+ A\wedge A$ the local curvature two-form. When $M$ is the Euclidean space-time, the classical action describing the dynamics is
\begin{equation}
S=   \frac{1}{4}\int F^a_{\mu\nu} {F^{\mu\nu}}^a {\mathrm d}^4 x \label{act}
\end{equation}
where  $F= F^a_{\mu\nu} \tau_a dx^\mu\wedge dx^\nu$  and $\tau_a$ are the generators of the Lie algebra. The quadratic, free, sector of the theory, upon integrating by parts with suitable boundary conditions, can be written as
\begin{equation}
S= \frac{1}{2}\int {\mathrm d}^n x\, {\int {\mathrm d}^n y\,} A^a_\mu(x) \Delta_{ab}^{\mu\nu} (x,y) A^b_\nu(y), \label{action}
\end{equation}
with
\begin{equation}
\Delta_{ab}^{\mu\nu} (x,y)= (-\partial^2\delta^{\mu\nu} +\partial^\mu\partial^\nu) \delta^{(4)} (x-y)\delta_{ab} \label{Mop}
\end{equation}
and
\begin{equation}
Z[J]= \int_{\mathcal{A} }  [d \mu(\mathcal{A})] \exp\left(-\frac{1}{2}(S[A]+S_I[A,J])\right), \label{zeta}
 \end{equation}
 the generating functional and  $\mathcal{A}$ the space of gauge connections.  If  the operator $\Delta_{ab}^{\mu\nu}$ were invertible, as it is the case for scalar theories,  the Gaussian integral in \eqn{zeta} could be {formally} performed to give:
\begin{equation}
Z[J] = \left(\det \Delta\right)^{-\frac{1}{2}} \exp\left( \frac{1}{2} \int J \Delta^{-1} J\right)
\end{equation}
with $\Delta^{-1}$ the Euclidean propagator, but this is not the case for gauge theories unless we manage to integrate over equivalence classes of gauge connections. Indeed, because of gauge invariance of the free action under
$
A \rightarrow A^g= g A g^{-1} + dg g^{-1},
$
field configurations of the form  $dg g^{-1}$ (so-called pure gauge terms) are in the kernel of $\Delta_{ab}^{\mu\nu}$, such that
\begin{equation}
\Delta^{\mu\nu} \partial_\nu g g^{-1}= 0, \label{zeromodes}
\end{equation}
showing that the operator \eqn{Mop} has eigenvectors with zero eigenvalues (so-called zero modes).  To restrict the functional integration to gauge inequivalent potentials, one has to limit the functional integral \eqn{zeta} to the quotient space $\mathcal{B}= \mathcal{A}/\mathcal{G}.$
  Mathematically, this amounts to choosing a surface $\Sigma_{f} \subset \mathcal{A}$  which possibly intersects the gauge orbits only once.
The choice of $\Sigma_f$ is  physically rephrased as a {gauge fixing},  $f(A)= h$, for some chosen  functions $f, h$.  Within the standard quantization procedure of gauge theories, this is  achieved through the  introduction of the Faddeev-Popov determinant, which is only unambiguous when no topological obstructions arise.  But precisely the latter is at the origin of Gribov ambiguities. Therefore, we shortly summarise the standard procedure and the related topological issues.

The kinematical configuration space $\mathcal{A}$  is usually assumed to be  globally equivalent to the  product  $ \mathcal{B} \times \mathcal{G}$. In such a case, one has for the integration measure
\begin{equation} \label{measure}
 [d \mu(\mathcal{A})]= [d\mu(\mathcal{B})]  ~ [d\mu (\mathcal{G}) ]= [d\mu(\mathcal{B})]  ~[d\alpha ],
 \end{equation}
for gauge transformations close to the identity, $g(x) \simeq {\mathbf 1}+ \alpha^a (x) \tau_a$.  By    performing  a change of variables  $ [d\alpha ]\rightarrow  [d f(A) ]$, with the insertion of   the Jacobian
\begin{equation}
 \rm{Det} \frac{\delta f^{a}(x)}{\delta \alpha^{b}(y)}  \equiv \rm{Det}  \Delta (x,y) \label{FPdet},
\end{equation}
one arrives at
\begin{equation}
[d \mu(\mathcal{A})] \rm{Det} \Delta  =  [{ [d\mu(\mathcal{B})] }  [d\alpha ]  \rm{Det} \Delta   =  { [d\mu(\mathcal{B})] }\; [d f ].
\end{equation}
And, finally, by integrating over $[df]$ with the insertion of a delta function $\delta(f(A) -h(x))$ which implements the gauge choice, one obtains the measure on the quotient space:
\begin{equation}
  [d \mu(\mathcal{A})] \; \rm{Det} \Delta\;\; \delta(f(A) -h(x))={ [d\mu(\mathcal{B})] }. \label{measuregauged}
  \end{equation}
  The Jacobian in \eqn{FPdet} is the known Faddeev-Popov determinant.

However, the gauge fixing described above is not enough to remove unphysical degrees of freedom if the theory is non-Abelian.
Indeed, let us consider the gauge orbit
\begin{equation}
A^g =g  A g^{-1} + d g g^{-1}  \simeq  A+ D\alpha,
\end{equation}
with here  $ - (\partial_\mu D^\mu) \delta^{(4)}(x-y) \delta^{ab} $ is the FP operator for this choice of gauge fixing and $D\alpha=  d\alpha+ \alpha^a\wedge A^b [\tau_a, \tau_b]$.
 The gauge fixing condition  $\partial^\mu A^g _\mu= 0$  yields the so-called  {\it equation of copies} \cite{Gribov1978}
 \begin{equation}
 \partial_\mu D^\mu \alpha= 0 ,\label{copies}
 \end{equation}
 which may have nontrivial solutions, whenever the gauge group is non-Abelian, yielding to Gribov ambiguities.\footnote{\label{footnote3} In the Abelian case we only have trivial solutions, if we further assume that $\lim_{x\rightarrow\infty} \alpha(x)=0$.}
The  problem can be traced back to the topological non-equivalence of the kinematical space of connections ${\mathcal A}$ and the physical space ${\mathcal B}$, \cite{Singer1978}, \cite{Narasimhan1979}. Indeed \eqn{measuregauged} is only valid if $\mathcal{A}= \mathcal{B} \times \mathcal{G}$ globally, which, for $\mathcal{G}=\widehat{U}(N)$ is only true for $N=1$  namely for QED (see for example \cite{Nair2005} for a clear pedagogical discussion).

\section{The Gribov copies equation for  NC QED }\label{NCQED}
For each two functions $f, g$ defined on Euclidean space-time $\mathbb R^4$,  the
noncommutative Moyal star product $f\star g$ \cite{Moyal:1949sk} can be given the following asymptotic expansion
\begin{equation}
(f\star g)(x)=f(x)\exp \left\{ \frac{i}{2}\,\Theta ^{\rho \sigma }\overset{%
\leftarrow }{\partial _{\rho }}\overset{\rightarrow }{\partial _{\sigma }}%
\right\} g(x),  \label{star}
\end{equation}%
with $\rho ,\sigma =1,..,4$.
 The antisymmetric matrix $\Theta $, has the following nonzero components
\begin{equation}
\Theta _{1,2}=-\Theta _{2,1}=\Theta_1 ,\,\,\Theta _{3,4}=-\Theta _{4,3}=\Theta_2
 ,  \label{nctheta}
\end{equation}%
where $\Theta_i $ are real deformation parameters, in principle all different from each other,
characterizing noncommutativity. For simplicity, we perform a rescaling in order to make all parameters equal. When $\Theta_i \rightarrow 0$, the star
product tends to the standard commutative point-wise product of $f$
and $g$.
\subsection{Gauge transformations}
 Gauge theories  with gauge group ${\mathcal{G}}= {\rm Maps}(\mathbb R^4, U(N))$
are modified in the noncommutative setting by replacing the point-wise product with the non-local product~\eqn{star}.
The elements of the gauge group,  $U_\star(x)$, are defined as star exponentials
\begin{equation}
U_\star(x)=\exp_\star\left(i \alpha (x)^i T_i\right)= 1+ i \alpha^i
(x) T_i - \frac{1}{2}(\alpha^i\star \alpha^j)(x) T_i T_j + \ldots
\label{starexp}
\end{equation}
with $T_i$ the Lie algebra generators of the structure group, $T_i \in \mathfrak{u}(N)$.
The associated  matter fields transform under  deformed gauge transformations according to
\begin{equation}
\phi(x) \longrightarrow U_\star(x)\triangleright_\star \phi(x)=\exp_\star\left(i
\alpha^i (x) T_i\right) \triangleright_\star \phi (x). \label{gaugetransf}
\end{equation}
 with $\triangleright_\star$ indicating simultaneously  the appropriate representation of the Lie algebra generators  {\it and} the $\star$ multiplication in space-time.
 At the infinitesimal level, we have then
\begin{equation}
\phi(x) \longrightarrow  \phi(x) +i (\alpha\triangleright_\star \phi)(x),
\label{infgaugetransf}
\end{equation}
with
\begin{equation}
(\alpha\triangleright_\star\phi)(x)=i \left(\alpha^j(x)\star (T_j\triangleright \phi)
\right)(x).
\end{equation}
The gauge potential  transforms  as
\begin{equation}
A_\mu \rightarrow A^{\prime}_{\mu} = U_\star \star A_{\mu} \star U_\star^{\dagger}
+i\, U_\star \star \partial_{\mu} U_\star^{\dagger}.  \label{gtransgen}
\end{equation}
Specialising to NCQED, the infinitesimal transformation  reads therefore
\begin{equation}
A_\mu \rightarrow A^{\prime}_{\mu}[\alpha] = A_{\mu} + D_{\mu}\alpha + \mathcal{O%
}(\alpha),  \label{gtransinfin}
\end{equation}
where the covariant derivative $D_{\mu}$, now  only due to space-time {noncommutativity}, is given by

\begin{equation}
D_{\mu}\alpha = \partial_{\mu}\alpha + i \left(\alpha \star A_{\mu} - A_{\mu}\star \alpha\right)
\label{cd}
\end{equation}
and reduces to the standard Abelian form in the  commutative limit  $%
\Theta\rightarrow 0$.
 The field strength $F$  is given by
 \begin{equation}
 F_{\mu\nu}= D_{\mu} A_\nu -D_{\nu} A_\mu
 \end{equation}
 and transforms covariantly under star-gauge transformations.

It is precisely the occurrence of the covariant derivative which makes NCQED similar to non-Abelian gauge theories in commutative space-time, and, independently from the specific form of the gauge action $S[A]$ one deals with,\footnote{Besides the  natural candidate, which is obtained by replacing the point-wise product with the noncommutative one, i.e.
\begin{equation}
S[A]= \frac{1}{4}\int F_{\mu\nu}\star F^{\mu\nu}
\end{equation}
other proposals have been considered in the literature. See for example \cite{MVW14} and refs. therein.}  it is a meaningful question and a relevant issue for the quantization of the theory to investigate the existence of nontrivial solutions for the equation of Gribov copies. The latter is readily obtained by considering  $A_{\mu}$ and $A'_{\mu}$ satisfying  the same gauge condition, e.g. the Landau gauge,
\begin{eqnarray}
\partial^{\mu}A_{\mu}=\partial^{\mu}A'_{\mu}=0.
\end{eqnarray}
It follows that $\alpha$ yields a Gribov copy if
\begin{eqnarray}
-\partial^{\mu}D_{\mu}\alpha(x) &=& 0
\end{eqnarray}
The problem has been addressed in \cite{Canfora2016, Kurkov2018} where the equation has been solved for a class of particularly simple potentials, and an  infinite number of copies, which cannot be discarded by boundary conditions, i.e. they are non-trivial,  has been found.

 It is therefore clear that the standard  quantization procedure applied to NCQED does not provide a well-defined measure of integration over  gauge fields. In order to overcome the problem, we propose to apply the same techniques adopted for non-Abelian gauge theories, which have led to restricting the functional  integral  to an integration region free of copies: the Gribov region.

\section{Establishing the Gribov Region}\label{Gribovregion}
Once  the existence of copies has been established in terms of zero modes of the Faddeev-Popov operator, we need to characterise  the  region $\Omega$ in the space of connections where there are no ambiguities. The definition  which we adopt here for NCQED  was proposed by Gribov in 1978 for non-Abelian theories \cite{Gribov1978}. For the present case it reads:
\begin{eqnarray}
\Omega \equiv \{A_{\mu},\,\,\partial^{\mu} A_{\mu}=0,\,\, \Delta(x,y)>0\},
\end{eqnarray}
where the Faddeev-Popov operator is explicitly  given by
\begin{eqnarray}\label{FPstar}
\Delta(x,y)=-\partial_{\mu}D^{\mu}\delta(x-y)=-\partial^2\delta(x-y)+i\lambda[A_{\mu},\partial_{\mu}\delta(x-y)]_{\star}
\end{eqnarray}
and $ [f, g]_\star= f\star g -g\star f$.

The boundary $\partial \Omega$ represents  the first  Gribov horizon. For non-Abelian gauge theories it was shown \cite{DellAntonio:1989wae,DellAntonio:1991mms} that $\Omega$ is convex, bounded in all directions in the space of gauge connections  and that
each gauge orbit passes at least once in  $\Omega$. The inverse of the FP
operator, or equivalently the ghost propagator with external gauge field,
$G(k, A)= \Delta^{-1}(k, A)$,  with $k$ the momentum variable, can
be used to implement the restriction to $\Omega$ \cite{Gribov1978}. In the following we make the non-trivial assumption  that the results established
in \cite{DellAntonio:1989wae,DellAntonio:1991mms} may be extended to NCQED, namely that
\begin{itemize}
\item
The Gribov region is bounded in every direction  (in the functional space
of transverse gauge potentials).
\item  The Faddeev-Popov determinant changes sign at the Gribov horizon.
\item  Every gauge orbit passes inside the Gribov horizon.
\end{itemize}
Under these assumptions, we will adapt to NCQED   the procedure introduced in \cite{Capri:2012wx} to restrict  the
functional integration to the Gribov region. In the next subsection we provide some arguments to support these assumptions.

\subsection{ Properties of the NC Faddeev-Popov operator}

The FP operator is given by
\begin{align}
	\Delta(x,y, A) =-\partial_{\mu}D^{\mu}\delta(x-y) =-\partial^2\delta(x-y)+i\lambda[A_{\mu},\partial_{\mu}\delta(x-y)]_{\star}
\end{align}

\begin{itemize}
	\item \emph{\textbf{The Gribov region is limited in every direction}}

	First of all, note that  $\Delta(x,y, A)$ is hermitian in Euclidean space so that for real functions $f(x)$ and $g(x)$ one has
	\begin{align}
	\int d^Dx \int d^D y\; f(x)\star \Delta(x,y, A) \star g(y) \in \mathds{R}
	\end{align}
	Now, suppose that $\Delta(x,y, A)$ develops a zero mode for some value of $A$. This must come from some negtive contribution from the term ${\cal M}(x,y, A) = i\lambda[A_{\mu},\partial_{\mu}\delta(x-y)]_{\star}$ because the term $\partial^2\delta(x-y)$ is positive. Therefore one must have for some value of $A$ and some function $\psi(x)$ that
	\begin{align}
		\int d^Dx \int d^D y\; \psi(x)\star {\cal M}(x,y, A) \star \psi(y) = a < 0
	\end{align}
	Now introduce a real parameter $t$ multiplying $A$ so as to parameterize its amplitude isotropicaly and consider the expression
		\begin{align}
			\int d^Dx \int d^D y\; \psi(x)\star \Delta(x,y, t A) \star \psi(y)  = \int d^Dx \int d^D y\; \psi(x)(-\partial^2) \psi(x)  + t a
		\end{align}
	where we used that ${\cal M}(x,y, tA) = t {\cal M}(x,y, A)$. The first term is positive and since $a$ is negative, it is clear that this expression will become negative for some sufficiently high value of $t$. This means that the Gribov region, defined as the region in $A$ space such that $\Delta(x,y, A)$ is positive, will be eventually crossed for some finite value of $t$ and thus it must be limited in every direction.

	\item \emph{\textbf{The Faddeev-Popov determinant changes sign at the Gribov horizon.}} This follows from the definition of the Gribov horizon and can be seen from the above remarks. As the amplitude of $A$ becomes large the $FP$ reaches a zero mode (the Gribov horizon) and afterward it becomes negative as we saw by raising the parameter $t$ in the proof above.

	\item \emph{\textbf{Every gauge orbit passes inside the Gribov horizon.}}
	In order to prove this, one would need a functional $f(A)$ such that under an infinitesimal  gauge transformation $\delta A_{\mu} = D_{\mu} \alpha$ its extrema are determined by the gauge condition
	\begin{align}
	\delta f(A) = 0 \Rightarrow \partial_{\mu} A_{\mu} = 0
	\end{align}
	and its second variation is positive, thus making these extrema in fact the minima
	\begin{align}
		\delta^2 f(A) > 0 \Rightarrow \partial_{\mu}D_{\mu}  > 0
	\end{align}
	 The argument then go on to show that for every gauge orbit, as one moves along the orbit, the functional $f(A)$ will eventually attain its absolute minimum. But since the condition for the minimum is exactly the condition that defines the Gribov region, $\partial_{\mu}D_{\mu}  > 0$, it would follow that every gauge orbit will pass through the Gribov region. For the usual non-abelian case, the functional is $f(A) = \frac 12 \int d^D x A^a_{\mu}A^a_{\mu}$.

 We have verified that the same functional works for the NC case, by using the cyclicity of the product. Given
 \begin{eqnarray}
 	f(A)= \int d^Dx\, A_\mu \star A_\mu = \int d^Dx\, A_\mu  A_\mu .
 \end{eqnarray}
The first variation is
\begin{eqnarray}
	0= \delta f(A) &=&  \int d^Dx\, A_\mu \star \delta A_\mu =  \int d^Dx\, A_\mu \star (\partial_\mu \alpha + i [ \alpha, A_\mu]_\star)= \int d^Dx\, A_\mu \star \partial_\mu \alpha \nonumber\\
	&=& \int d^Dx\, A_\mu\cdot  \partial_\mu \alpha    =\int d^Dx\, \partial _\mu A_\mu \cdot \alpha=\int d^Dx\, \partial _\mu A_\mu \star \alpha,
\end{eqnarray}
which implies $\partial _\mu A_\mu= 0$.  Here  we have used the cyclicity of Moyal  product under the integral  ($\int a\star b\star c$ is cyclic) and $\int a\star b= \int b\star a$. Notice that the last property is in general not true for other star products. Then the second variation gives
\begin{eqnarray}
	0<\delta^2 f(A)= \int dx \delta \partial_\mu  A_\mu \star \alpha= \int \partial_\mu D_\mu \alpha \star \alpha = \int \partial_\mu D_\mu \alpha \cdot \alpha
\end{eqnarray}
which, for arbitrary $\alpha$ implies the positivity of  $\partial_\mu D_\mu$, namely the condition defining the Gribov region.  As mentioned above, to provide a complete proof one would need the show that the functional $f(A)$ attains its absolute minimum along every orbit in the NC case as well. This is a reasonable expectation but to pursue such a formal proof is beyond the scope of this work.
\end{itemize}

\section{NCQED no-pole condition}\label{five}

 The requirement that the FP operator has no zero eigenvalues is implemented in the standard approach \cite{Gribov1978}  in terms of the request that the
inverse of the operator doesn't have poles: this is the  so-called {\it no pole condition}. The latter is seen to restrict
the region of integration to  the Gribov region  by introducing   a factor $V(\Omega)$ in the generating functional,
\begin{equation}\label{fungen}
Z(A,c,\bar c)=\int_{\Omega}\,\, [dA][d\bar{c}][dc]\,V(\Omega)\, \delta(\partial\cdot A)
\exp\left[-S[A]-\int\,dx\,dy \,\bar{c}(x)\star \Delta(x,y) \star c(y) \right]
\end{equation}
with $\Delta(x,y)$ given in \eqref{FPstar}. We recall that we are in the Euclidean setting.
We work in the Landau gauge, $\delta(\partial\cdot A)$. To determine $V(\Omega)$ we use the relationship between the ghost sector of the theory and the FP operator, which emerges from calculating the exact ghost propagator
\begin{equation}
\langle\bar c(p) \,c (-p)\rangle= \frac{\delta}{\delta J_c(y)}\frac{\delta}{\delta J_{\bar c}(x)} Z[J]
\end{equation}
with
\begin{eqnarray}
Z[J]&=&\int_{\Omega}\,\, [dA][d\bar{c}][dc]\,V(\Omega)\, \delta(\partial_{\mu}A_{\mu}) \\
&\times & \exp\left[-S[A]-\int\,dx\,dy \, \bar{c}(x)\star \Delta(x,y) \star c(y) +\int \, dx\, \left(J_c(x) \star c(x)+ \bar c(x) \star J(x)\right) \nonumber\right].
\end{eqnarray}
Thanks to the  cyclicity and closure properties of the Moyal product\footnote{They amount respectively  to the fact that $\int f_1\star ...\star f_n= \int f_n\star f_1\star...$ and $\int f\star g= f \cdot g$. }\label{cycl} we may perform the Gaussian integral over ghost fields in standard fashion, obtaining
\begin{eqnarray}
Z(J)=C \det (\Delta) \exp\left[-\int d^d x d^d y J_c(x)\star \Delta^{-1}(x,y)\star J_{\bar{c}}(y)\right]
\end{eqnarray}
so that
\begin{equation}
\langle\bar c(p) \,c (-p)\rangle= \int [d A] V(\Omega) \delta(\partial^\mu A_\mu) \det(-\partial^\mu D_\mu)  \Delta^{-1}(x,y) \exp(-S[A])
\end{equation}
with  $\Delta^{-1}(x,y)$ the ghost propagator at first order in $\hbar$.
According to \cite{Capri:2012wx} it is convenient to interpret $\Delta(x,y)$ as the space representation of an abstract operator $\Delta$
\begin{equation}
\Delta(x,y)= \langle x|\Delta|y \rangle.
\end{equation}
Inserting completeness relations  $\int d^dx\, |x\rangle\langle x|=1$, we may obtain   $\Delta(x,y)$  in momentum space
\begin{equation}
\Delta(p,q)= \int d^dx\int d^dy \langle p|x\rangle\star\langle x|\Delta |y\rangle\star \langle y|q\rangle=\int d^d x \int d^d y \,e^{-ipx}\star\Delta(x,y)\star e^{iqy}.
\end{equation}
The explicit calculation of star products, performed in  App. \ref{appA} yields
\begin{equation}\label{mpq}
\Delta(p,q)=
q^2\delta(p-q)+2i\lambda \sin\left(\frac{1}{2}\Theta_{\rho\sigma}p_{\rho}q_{\sigma}\right)\tilde A_{\mu}(p-q) q_{\mu}.
\end{equation}
where
\begin{eqnarray}
\langle x|p\rangle&=& \frac{e^{ipx}}{(2\pi)^{d/2}}\nonumber\\
\tilde A_{\mu}(p-q)&=& \int \frac{d^d x }{(2\pi)^d} e^{-i(p-q)x}\star A_\mu(x)= \int \frac{d^d x }{(2\pi)^d} \,e^{-i(p-q)x} A_\mu(x)
\end{eqnarray}

Upon introducing the matrix notation
\begin{eqnarray}
\mathbbm{1}_{pq}&=&\braket{p\vert \mathbbm{1}\vert q}=\delta(p-q),\\
\mathbbm{A}_{pq}&=&\braket{p\vert \mathbbm{A}\vert q}=-2i\sin\left(\frac{1}{2}\Theta_{\rho\sigma}p_{\rho}q_{\sigma}\right)\tilde A_{\mu}(p-q)\frac{q_{\mu}}{q^2},
\end{eqnarray}
we may rewrite Eq. \eqn{mpq} as
\begin{eqnarray}
\Delta_{pq}=q^2(\mathbbm{1}_{pq}-\lambda\mathbbm{A}_{pq}).
\end{eqnarray}

The no-pole condition works as a bound on the allowed amplitudes of the gauge field configurations and  implies that inverse of this expression makes sense. Then, the inverse of the operator $\Delta$ can be written as
\begin{eqnarray}
\Delta_{pq}^{-1}=\frac{1}{q^2(\mathbbm{1}_{pq}-\lambda\mathbbm{A}_{pq})}=\frac{1}{q^2}\sum_{n=0}^{\infty}[(\lambda\mathbbm{A})^n]_{pq}.\label{eq8}
\end{eqnarray}
We are interested in the poles of this expression, it is then convenient to study the corresponding normalized trace
\begin{eqnarray}
\frac{1}{V}\left.\Delta_{pq}^{-1}\right\vert_{p=q=k}\equiv \frac{1}{k^2}\left(1 +  \sigma(k,A)\right).\label{formfactor}
\end{eqnarray}
with $V=\int \frac{d^dx}{(2\pi)^d}$ the infinite volume factor. This equation defines the so-called form factor $\sigma(k,A)$, which encapsulates all non-trivial information about the pole structure of the full ghost propagator, in fact, if it vanishes we recover the free propagator expression. In order to implement the no-pole condition one has to integrate over all gauge field configurations leading to
\begin{eqnarray}\label{1PI}
\Delta^{-1}_{kk}(k) =\langle \Delta^{-1}_{kk}(k, A) \rangle_{conn} = \frac{V}{k^2} \left(1 + \langle \sigma(k,A) \rangle_{conn} \right) = \frac{V}{k^2} \frac{1}{1- \langle \sigma(k, A)\rangle_{1PI} }
\end{eqnarray}
where $\langle \cdots \rangle_{conn}$ is the sum of  all connected diagrams and $\langle \cdots \rangle_{1PI}$ is the sum of all one-particle-irreducible diagrams (see \cite{Capri:2012wx} for details). The no-pole condition is the requirement that the ghost propagator does not develop a pole for any value of $k\neq 0$. The evaluation of the exact no-pole condition thus amounts to the computation of the exact value $\langle \sigma(k, A)\rangle_{1PI}$. This is of course a very difficult problem ( see the discussion in Appendix (\ref{appALLORDERS})) and we will study here an approximate solution by considering only the first non-trivial contribution.
From expression \eqref{eq8}, we get
\begin{eqnarray}\label{M-1}
k^2\Delta^{-1}_{kk}(k, A)=\mathbbm{1}_{kk}+g\mathbbm{A}_{kk}+\lambda^2\mathbbm{A}^2_{kk}+\mathcal{O}(\lambda^3).
\end{eqnarray}
Then, we can write $\sigma(k,A)$ as
\begin{eqnarray}
\sigma(k,A)&=&\frac{1}{V }(\mathbbm{1}_{kk}+\lambda \mathbbm{A}_{kk}+\lambda^2\mathbbm{A}^2_{kk}+ \mathcal{O}(\lambda^3))-1.\\
&=& \frac{\lambda^2}{V} \int \frac{d^d p}{(2\pi)^{d}} \int \frac{d^d q}{(2\pi)^{d}}\mathbbm{A}_{kp}\left(  q^2 \Delta_{pq}^{-1}  \right)\mathbbm{A}_{qk}
\label{eq540}
\end{eqnarray}
where we used \eqref{eq8} and the relations
\begin{eqnarray}
\mathbbm{1}_{kk}&=&V,\\
\mathbbm{A}_{kk}&=&-2\sin(0)A_{\mu}(0)\frac{ik_{\mu}}{k^2}=0,
\end{eqnarray}
The first nontrivial contribution to $\sigma(k,A)$ is thus obtained approximating the  ghost propagator inside the integral \eqref{eq8} by its tree-level expression, that is  $q^2 \Delta_{pq}^{-1} \approx \mathbbm{1}_{pq}$. Therefore, noting that
\begin{eqnarray}
\mathbbm{A}^{2}_{kk}&=&\int \frac{d^d q}{(2\pi)^{d}}\, \braket{k\vert \mathbbm{A}\vert q}\braket{q\vert \mathbbm{A}\vert k}\nonumber\\
&=&-4\int \frac{d^d q}{(2\pi)^{d}}\sin\left(\frac{1}{2}\Theta_{\rho\sigma}k_{\rho}q_{\sigma}\right)\tilde A_{\mu}(k-q)\frac{q_{\mu}}{q^2}\sin\left(\frac{1}{2}\Theta_{\rho\sigma}q_{\rho}k_{\sigma}\right)\tilde A_{\nu}(q-k)\frac{k_{\nu}}{k^2}\label{eq654}
\end{eqnarray}
equation \eqref{eq540} becomes
\begin{eqnarray}
\sigma(k,A,\Theta)=-\frac{4\lambda^2}{V}\int \frac{d^d q}{(2\pi)^{d}}\sin\left(\frac{1}{2}\Theta_{\rho\sigma}k_{\rho}q_{\sigma}\right)\tilde A_{\mu}(k-q)\frac{q_{\mu}}{q^2}\sin\left(\frac{1}{2}\Theta_{\rho\sigma}q_{\rho}k_{\sigma}\right)\tilde A_{\nu}(q-k)\frac{k_{\nu}}{k^2}
\end{eqnarray}
where the $\Theta$ dependence was highlighted in the argument of $\sigma$. Considering the Landau gauge $q_{\mu}A_{\mu}(k-q)=k_{\mu}A_{\mu}(k-q)$ and changing $q\rightarrow q+k$, we obtain,
\begin{eqnarray}
\sigma(k,A,\Theta)&=&\frac{4\lambda^2}{V}\frac{k_{\mu}k_{\nu}}{k^2} \int \frac{d^d q}{(2\pi)^{d}} \sin^2\left(\frac{1}{2}\Theta_{\rho\sigma}q_{\rho}k_{\sigma}\right)\tilde A_{\mu}(-q)\frac{1}{(k+q)^2}\tilde A_{\nu}(q).
\end{eqnarray}
The latter may be further simplified by observing that, in the Landau gauge, $\tilde A_\mu(-q) \tilde A_\nu(q)$ is transversal,
\begin{equation}
\tilde A_\mu(-q) \tilde A_\nu(q)=\omega(A)\left(\delta_{\mu\nu}-\frac{q_\mu q_\nu}{q^2}\right).
\end{equation}
Multiplying by $\delta_{\mu\nu}$ the factor $\omega(A)$ is found to be    $\omega(A)= \frac{1}{d-1}|\tilde A|^2$.
Thus we obtain
\begin{eqnarray}
\sigma(k, A,\Theta)=\frac{4\lambda^2}{V}\frac{1}{d-1}\frac{k_{\mu}k_{\nu}}{k^2} \int \frac{d^d q}{(2\pi)^{d}} \sin^2\left(\frac{1}{2}\Theta_{\rho\sigma}q_{\rho}k_{\sigma}\right)\frac{\vert A_{\lambda}(q)\vert^2}{(k+q)^2}\left(1-\frac{q_{\mu}q_{\nu}}{q^2}\right)\label{nopole}.
\end{eqnarray}

Observing that $\dfrac{\sin^2(x)}{x^2}\leq 1$ we have $\sigma(k,\Theta)\leq I(k,\Theta)$. The function $I(k,\Theta)$ is defined as
\begin{eqnarray}
I(k, A,\Theta) &=&\frac{\lambda^2}{V}\frac{\Theta_{\rho\sigma}\Theta_{\omega\lambda}}{d-1}\frac{k_{\mu}k_{\nu}k_{\sigma}k_{\lambda}}{k^4} \int \frac{d^d q}{(2\pi)^{d}} \vert \tilde A_{\lambda}(q)\vert^2\frac{k^2 q^2}{(k+q)^2}\left(\delta_{\mu\nu}-\frac{q_{\mu}q_{\nu}}{q^2}\right)\frac{q_{\rho}q_{\omega}}{q^2}\\
&=&\frac{k_{\mu}k_{\nu}k_{\lambda}k_{\sigma}}{k^4} \mathcal{I}_{\mu\nu\lambda\sigma}(k,\Theta),\label{nopolesmall}
\end{eqnarray}
with
\begin{eqnarray}
\mathcal{I}_{\mu\nu\lambda\sigma}(k,A,\Theta )&=&\frac{\lambda^2}{V}\frac{\Theta_{\rho\sigma}\Theta_{\omega\lambda}}{d-1} \int \frac{d^d q}{(2\pi)^{d}} \vert \tilde A_{\lambda}(q)\vert^2\frac{k^2 q^2}{(k+q)^2}\left(\delta_{\mu\nu}-\frac{q_{\mu}q_{\nu}}{q^2}\right)\frac{q_{\rho}q_{\omega}}{q^2}\label{functionnopolesmall}.
\end{eqnarray}
We are interested in finding an upper-bound for the latter, so to have a refined no-pole condition $I(k,\Theta)<1$.
In Appendix \eqref{apa} and \eqref{apb}  we perform an explicit analysis in $d=2$ and $d=4$ and show that $I(k,A,\Theta)$ is an increasing function of $k$, with the maximum value  attained at $k\rightarrow \infty$. In what follows we shall therefore assume $d=2,4$, although  the analysis could  be probably extended to other dimensions.
Because of that, the  no-pole condition is certainly satisfied if we impose
 \begin{eqnarray}
 I_{max}(A,\Theta )<1,
 \end{eqnarray}
 with
 \begin{eqnarray}
 I_{max}(A,\Theta )=\frac{k_{\mu}k_{\nu}k_{\lambda}k_{\sigma}}{k^4} \mathcal{I}_{\mu\nu\lambda\sigma}(\infty,A,\Theta ).\label{lema}
 \end{eqnarray}
Taking the limit $k^2\rightarrow \infty$ in \eqref{functionnopolesmall}, we obtain
\begin{eqnarray}
\mathcal{I}_{\mu\nu\lambda\sigma}(\infty, A,\Theta)&=&\frac{\lambda^2}{V}\frac{\Theta_{\rho\sigma}\Theta_{\omega\lambda}}{d-1} \int \frac{d^d q}{(2\pi)^{d}} \vert \tilde A_{\lambda}(q)\vert^2\left(\delta_{\mu\nu}q_{\rho}q_{\omega}-\frac{q_{\mu}q_{\nu}q_{\rho}q_{\omega}}{q^2}\right). \label{fff}
\end{eqnarray}
Using the following properties of regularized momentum integrals \cite{Anacleto:2014aha}
\begin{eqnarray}
\int \frac{d^d q}{(2\pi)^{d}} f(q^2) q_{\rho}q_{\omega}&=&\frac{\delta_{\rho\omega}}{d}\int \frac{d^d q}{(2\pi)^{d}} \, q^2 f(q^2),\label{id1}\\
\int \frac{d^d q}{(2\pi)^{d}} f(q^2) q_{\omega}q_{\rho}q_{\mu}q_{\nu}&=&\frac{\delta_{\mu\nu}\delta_{\omega\rho}+\delta_{\mu\rho}\delta_{\nu\omega}+\delta_{\mu\omega}\delta_{\rho\nu}}{d(d+2)}\int \frac{d^d q}{(2\pi)^{d}} \, q^4 f(q^2),
\end{eqnarray}
we rewrite \eqref{fff} as
\begin{eqnarray}
\mathcal{I}_{\mu\nu\lambda\sigma}(\infty, A,\Theta )&=&\frac{\lambda^2}{V}\Theta_{\rho\sigma}\Theta_{\omega\lambda} \left(\frac{(d+1)\delta_{\mu\nu}\delta_{\rho\omega}-\delta_{\mu\rho}\delta_{\nu\omega}-\delta_{\mu\omega}\delta_{\rho\nu}}{d(d+2)(d-1)}\right)\int \frac{d^d q}{(2\pi)^{d}}\, q^2\vert \tilde A_{\lambda}(q)\vert^2\nonumber\\
&=&\frac{\lambda^2}{V}\left(\frac{(d+1)\delta_{\mu\nu}\Theta^2_{\sigma\lambda}-\Theta_{\mu\sigma}\Theta_{\nu\lambda}-\Theta_{\mu\lambda}\Theta_{\nu\sigma}}{d(d+2)(d-1)}\right) \int \frac{d^d q}{(2\pi)^{d}}\, q^2\vert \tilde A_{\lambda}(q)\vert^2. \label{eq764}
\end{eqnarray}
 Substituing \eqref{eq764} in \eqref{lema}, we conclude that
\begin{eqnarray}\label{Imax}
I_{max}( A,\Theta )&=&\frac{\lambda^2}{V}\frac{k_{\mu}k_{\nu}k_{\lambda}k_{\sigma}}{k^4} \left(\frac{(d+1)\delta_{\mu\nu}\Theta^2_{\sigma\lambda}-\Theta_{\mu\sigma}\Theta_{\nu\lambda}-\Theta_{\mu\lambda}\Theta_{\nu\sigma}}{d(d+2)(d-1)}\right) \int \frac{d^d q}{(2\pi)^{d}} \, q^2 \vert \tilde A_{\lambda}(q)\vert^2\nonumber\\
&=&\frac{\lambda^2}{V} \left(\frac{(d+1)\Theta^2}{d(d+2)(d-1)}\right) \int \frac{d^d q}{(2\pi)^{d}} \, q^2 \vert \tilde A_{\lambda}(q)\vert^2 ,
\end{eqnarray}
where we have used  $k_{\sigma}k_{\lambda}\Theta^2_{\sigma\lambda}=k^2\Theta^2$ and $k_{\mu}k_{\sigma}\Theta_{\mu\sigma}=0$.

Thanks to the latter, we may  estimate  the first Gribov region and exhibit its dependence on $\Theta$. Indeed, since $\sigma(k,A,\Theta)$ is always smaller or equal than $I(k,A,\Theta)$, which we have shown to be an increasing function of $k$, we may choose in Eq. \eqn{fungen}
\begin{equation}
V(\Omega)=\vartheta(1-I_{max}(A, \Theta)).
\end{equation}
or using the Heaviside function integral form,
\begin{eqnarray}
V(\Omega)= \int_{-\infty +\epsilon}^{+\infty +\epsilon}\frac{d\tau}{2\pi i \tau}e^{\tau(1-I_{max}(\Theta))}.
\end{eqnarray}
We can insert this into the path integral \eqref{action} which takes the form
\begin{eqnarray}\label{Zfunctional}
Z(J)=C\int \frac{d\tau}{2\pi i \tau}e^{\tau}\int[dA]e^{-\tau I_{max}(A,\Theta)-\int d^dx\,\,\frac{1}{4}F_{\mu\nu}\star F_{\mu\nu} + \frac{1}{2\alpha}b \star\partial_{\mu}A_{\mu} +A_{\mu}\star J_{\mu}}.
\end{eqnarray}
The tree-level photon propagator can be read from the quadractic part and is given by
\begin{eqnarray}
\Delta_{\mu\nu}(k^2)&=&\left((1+\gamma)k^2 \delta^{\mu\nu}+\left(\frac{1}{\alpha}-1\right)k_{\mu}k_{\nu}\right),\,\,\gamma \equiv\frac{2(d+1)}{d(d+2)(d-1)}\frac{\tau\,\Theta^2 \lambda^2}{V}.
\end{eqnarray}
Therefore, the transverse photon propagator is given by
\begin{eqnarray}
\langle \tilde A_{\mu}(-k)\tilde A_{\nu}(k)\rangle &=&\frac{1}{(1+\gamma)k^2}\left(\delta_{\mu\nu}-{k_{\mu}k_{\nu}\over k^2}\right).
\end{eqnarray}
We see that, {\it at least within Moyal type non-commutativity}\footnote{Indeed, if the star product is not closed, already the quadratic part of the action, consisting of $\int (\partial_\mu A_\nu-\partial_\nu A_\mu) \star  \partial_\mu A_\nu-\partial_\nu A_\mu)$ is deformed. Namely, the $\star$ product cannot be removed, resulting in a modification of $\Delta_{\mu\nu}(k^2)$.} there is no qualitative modification of the photon propagator, which is changed only by a scale factor that can be absorbed by a $\Theta$ dependent renormalization of the gauge field. So the effect of Gribov copies is not felt and, in particular, there is no confining phase due to Gribov effects. It is a situation very similar to the case of the ${\cal N} = 4$ Super Yang-Mills, where Gribov effects are suppressed and no scale is generated \cite{Capri:2014tta}, thus maintaining the conformal invariance of the theory even when Gribov copies are taken into account. In the present case, the scale associated with the Gribov region is the noncommutative parameter $\Theta$, which is not dynamically generated but is already present as a fundamental scale of the theory from the start.

 This is to be contrasted with the original (commutative space) non-abelian Gribov analysis \cite{Gribov1978}, where the modification of the quadratic part of the action coming from the form factor is $\sim \frac{\vert \tilde A_{\lambda}(q)\vert^2}{q^2}$ and it leads to an effective gluon propagator of the form  $\sim \frac{q^2}{q^4 + \kappa^4}$, where $\kappa$ is the Gribov parameter determining the scale associated with the Gribov horizon. This scale is dynamically generated and can be computed by a gap equation,
 \begin{eqnarray}
 \int \frac{d^d q}{(2\pi)^{d}} \frac{1}{q^4 + \kappa^4} = 1
 \end{eqnarray}
 whereas in the present case, the modification of the quadratic part induced by $I_{max}( A,\Theta )$, \eqref{Imax},  is of the form $\sim \gamma q^2\vert \tilde A_{\lambda}(q)\vert^2$, with $\gamma$ proportional to $\Theta^2$.  This produces an ill-defined gap equation. In fact, the gap equation is obtained by extremising the vacuum energy functional, ${\cal E}(\tau) = - \ln Z$, with respect to $\tau$. So, after integrating over the gauge fields in \eqref{Zfunctional} considering only the quadratic part, we would need to find the extremum of
 \begin{eqnarray}
 f(\tau) = -\ln \tau + \tau + \int \frac{d^d q}{(2\pi)^{d}} \ln \left((1+\gamma )q^2\right)
 \end{eqnarray}
 The integral in this expression is ill-defined and in fact it is zero under dimensional regularization. Considering further that we eventually want to take the thermodynamic limit $V\rightarrow \infty$, keeping $\gamma$ finite (which means take $\tau \rightarrow \infty$), we see that this equation provides no consistent solution in this limit. Therefore, within the approximations considered, there is no dynamically generated scale  associated with the Gribov region in  NCQED. We conclude that the scale associated with the Gribov parameter in the non-abelian setting  is replaced by the nondynamical noncommutative scale here. We can observe that in the commutative limit,  when $\Theta \to 0$, all Gribov features are removed and the theory returns to  standard QED.

As a final, important, remark, notice that our analysis is not exact. Besides the approximation of working only to leading nontrivial order in the form factor $\sigma(k, A, \Theta)$, we also consider a stronger no-pole condition than the needed one, since $I(k, A, \Theta)$ is an upper bound for $\sigma(k, A, \Theta)$. The result could be made more precise by maximising directly \eqref{nopole}. There are however technical problems which we haven't succeeded to solve for the moment.

\section{Conclusion}

We have analysed  the no-pole condition for NCQED with Moyal type non-commutativity.  By exploiting  the formal analogy of the latter with non-Abelian gauge theories, we have restricted  the effective action to the Gribov region by introducing a constraining factor, $V(\Omega)$, which is obtained through the  ghost propagator in a fashion similar to $SU(N)$ gauge theories. As already stressed, our result is to be intended as a first  estimate of the Gribov region, not only because we work to leading nontrivial order in the form factor $\sigma(k,A,\Theta)$, but also because the request that   $I(k,\Theta)<1$ is a stronger no-pole condition than the needed one, $I(k,\Theta)$  representing  an upper bound for $\sigma(k, A, \Theta)$. This is due  to technical difficulties in maximising $\sigma(k,\Theta)$ in a $k$-invariant way. Still, this upper bound  guarantees the absence of poles for the ghost propagation. Also, it permits the calculation of a non-commutative transversal photon propagator proportional to $\dfrac{1}{(1+\gamma)k^2}$, which results to be  dependent on $\Theta^2$, the squared non-commutative parameter.  We therefore conclude that, in the approximation chosen,  there is no qualitative modification of the photon propagator, it being only  changed  by a scale factor that can be absorbed by a $\Theta$ dependent renormalisation of the gauge field. Let us stress, however, that this result has only been proven for the Moyal star product. Already mild modifications of the latter, such as the Wick-Voros product, would introduce a momentum dependent weight which could modify the behaviour of the propagator (see  for example \cite{Galluccio:2008wk} for an application to scalar field theory). It would be interesting to analyse this aspect more in detail.
Some immediate research questions can be investigated following these results. First, naturally suggested by the present analysis we would like to find a more precise solution of $\sigma(k,\Theta)<1$  dependent on the non-commutative parameter. Another possible path is to extend the analysis beyond the original Gribov semiclassical treatment and construct the analogous of the full Gribov-Zwanziger action \cite{zwanziger1989action, Zwanziger:1989mf, zwanziger1993renormalizability} for the NCQED. Finally, for the reasons outlined above, it would be interesting to repeat the analysis for a different kind of non-commutativity.

\section*{Acknowledgments}
O.H. acknowledges support  by the São Paulo Research Foundation (FAPESP), by the grant 2019/26291-8.  M. S. Guimaraes is a level 2 CNPq researcher under the contract 310049/2020-2.
L. Rosa and P. Vitale are partially supported by INFN.

\begin{appendices}
\section{Faddeev-Popov operator in   momentum space}\label{appA}
We start writing $\Delta (x,y)$ in the momentum space as
\begin{eqnarray}
\Delta(p,q)&=&\int d^d x \int d^d y \,e^{-ipx}\star\Delta(x,y)\star e^{iqy}.
\end{eqnarray}
We can write that as
\begin{eqnarray}
\Delta(p,q)&=&\int d^d x \int d^d
y \,e^{-ipx}\star(-\partial^2\delta(x-y)+i\lambda[A_{\mu},\partial_{\mu}\delta(x-y)]_{\star})\star e^{iqy},\nonumber\\
&=&\int d^d x \int d^d
y \,e^{-ipx}\star(-\partial^2\delta(x-y)+i\lambda A_{\mu}(x)\partial_{\mu}\star\delta(x-y)\nonumber\\&-&i\lambda\partial_{\mu}\star\delta(x-y)\star A_{\mu}(x))\star e^{iqy}.\label{faddeev}
\end{eqnarray}
By substituting the delta function in momentum space
\begin{eqnarray}
\delta(x-y)=\int \frac{dk}{2\pi ^{d}} e^{i k (x-y)}\label{deltafunction}
\end{eqnarray}
we arrive at
\begin{eqnarray}
\Delta(p,q) &=&\int d^d x \int d^d
y\int \frac{dk}{(2\pi)^{d}} \,e^{-ipx}\star\Bigl(k^2e^{i k (x-y)}-\lambda A_{\mu}(x)\star k_{\mu}e^{i k (x-y)}\nonumber\\&+&\lambda k_{\mu}e^{i k (x-y)}\star A_{\mu}(x)\Bigr)\star e^{iqy}.
\end{eqnarray}
On using  Eq. \eqn{star} to perform the  star   product in the $x$ and $y$ variables we get as an intermediate step
\begin{eqnarray}
\Delta(p,q) =\int d^d x \int d^d
y\int \frac{d^dk}{(2\pi)^{d}} \,e^{-ipx}\star\Bigl(k^2e^{i k (x-y)+iqy}e^{\frac{i}{2}\Theta_{\rho\sigma}k_{\rho}q_{\sigma}}\nonumber\\-\lambda A_{\mu}(x)\star k_{\mu}e^{i k (x-y)+iqy}e^{\frac{i}{2}\Theta_{\rho\sigma}k_{\rho}q_{\sigma}}\Bigr)+\lambda k_{\mu}e^{i k (x-y)-ipx}e^{\frac{i}{2}\Theta_{\rho\sigma}k_{\rho}p_{\sigma}}\star A_{\mu}(x)\star e^{iqy}\label{moyalfp}
\end{eqnarray}
By computing the  first term in \eqref{moyalfp} and using the  closure of the star product (see footnote \ref{cycl}), we obtain
\begin{eqnarray}
&&\int  d^d x \int d^d
y\int \frac{d^dk}{(2\pi)^{d}} \,e^{-ipx}\star k^2e^{i k (x-y)+iqy}e^{\frac{i}{2}\Theta_{\rho\sigma}k_{\rho}q_{\sigma}}\nonumber\\
&=& \int  d^d x\int \frac{d^dk}{(2\pi)^{d}} e^{-i p x} k^2 e^{ikx}\delta(q-k)e^{\frac{i}{2}\Theta_{\rho\sigma}k_{\rho}q_{\sigma}}\nonumber\\
&=& \int  d^d x q^2  e^{-i p x} e^{iqx}=
q^2 \delta(q-p).\label{moyalfp1}
\end{eqnarray}
From the second term in \eqref{moyalfp}, using the cyclicity and closure of the product
\begin{eqnarray}
&-&\lambda\int d^d x \int d^d y\int \frac{d^dk}{(2\pi)^{d}} \,k_{\mu} e^{\frac{i}{2}\Theta_{\rho\sigma}k_{\rho}q_{\sigma}}\, e^{-ipx}\star A_{\mu}(x)\star e^{i k (x-y)+iqy}\nonumber\\
&=& -\lambda\int d^d x \int \frac{d^dk}{(2\pi)^{d}} \,k_{\mu}\delta(q-k) e^{\frac{i}{2}\Theta_{\rho\sigma}k_{\rho}q_{\sigma}}\,e^{-ipx}\star A_{\mu}(x)\star e^{ikx}\nonumber\\
&=& -\lambda\int d^d x \, q_{\mu}\,e^{-ipx}\star A_{\mu}(x)\star e^{iqx}\nonumber\\
&=& -\lambda\int d^d x \, q_{\mu}\,e^{iqx}\star e^{-ipx}\star A_{\mu}(x)\nonumber\\
&=& -\lambda\int d^d x \, q_{\mu}\, e^{\frac{i}{2}\Theta_{\rho \sigma}q_\rho p_\sigma} e^{i(q-p)x}\star A_{\mu}(x)=-\lambda\int d^d x \, q_{\mu}\, e^{\frac{i}{2}\Theta_{\rho \sigma}q_\rho p_\sigma} e^{i(q-p)x} A_{\mu}(x)\nonumber\\
&=& -\lambda q_{\mu}A_{\mu}(p-q) e^{\frac{i}{2}\Theta_{\rho \sigma}q_\rho p_\sigma}.\label{moyalfp2}
\end{eqnarray}
Following the same procedure in the third term in \eqref{moyalfp},
\begin{eqnarray}
&&\lambda\int d^d x \int d^d y\int \frac{dk}{(2\pi)^{d}} \, k_{\mu}e^{\frac{i}{2}\Theta_{\rho\sigma}k_{\rho}p_{\sigma}}\, e^{i k (x-y)-ipx}\star A_{\mu}(x)\star e^{iqy}\nonumber\\
&=&\lambda q_{\mu}A_{\mu}(p-q) e^{\frac{i}{2}\Theta_{\rho \sigma}p_\rho q_\sigma}\label{moyalfp3}
\end{eqnarray}
Inserting \eqref{moyalfp1}, \eqref{moyalfp2} and \eqref{moyalfp3} in \eqref{moyalfp}, we prove that
\begin{eqnarray}
\Delta(p,q) &=&q^2 \delta(p-q)+\lambda q_{\mu}A_{\mu}(p-q)\left(e^{\frac{i}{2}\Theta_{\rho \sigma}p_\rho q_\sigma}-e^{-\frac{i}{2}\Theta_{\rho \sigma}p_\rho q_\sigma}\right)\nonumber\\
&=&q^2 \delta(p-q)+2i\lambda q_{\mu}A_{\mu}(p-q)\sin\left(\frac{1}{2}\Theta_{\rho\sigma}p_{\rho}q_{\sigma}\right).\label{deltapq}
\end{eqnarray}
\section{Analysis of $\sigma(k, \Theta)$ in two dimensions}\label{apa}
In this Section we consider the function $\sigma(k, \Theta)$ in $d=2$ dimensions.  In polar coordinates of the plane, $(q, \alpha)$, $q^2= q^\mu q_\mu, \alpha\in [0, 2\pi]$,  Eq. \eqn{nopole} becomes
\begin{equation}\label{sigma2}
\sigma(k, \Theta)=\frac{4 \lambda^2}{V}\int_0^\infty \frac{dq}{4\pi^2}\;q|A(q)|^2\int_0^{2\pi} d\alpha \frac{\sin^2{\left(\frac{1}{2}\Theta q k \sin{\alpha}
\right)  } \sin^2{\alpha} }{(k-q)^2},
\end{equation}
where we have chosen $k_\mu = (0,k)$.
Observing  that $\frac{\sin^2x}{x^2}\leq1$, we have
\begin{eqnarray}
\sigma(k, \Theta) &\leq& I(k, \Theta)  =\frac{4 \lambda^2}{V}\int_0^\infty \frac{dq}{4\pi^2}\;q|A(q)|^2\int_0^{2\pi} d\alpha \frac{\left( {\frac{1}{2}\Theta q k \sin{\alpha}   } \right)^2 \sin^2{\alpha} }{(k-q)^2} \\
&& \nonumber\\
&=&\frac{\lambda^2}{4\pi V}\int_0^\infty dq\;q|A(q)|^2\frac{1}{4}\Theta^2 q^2 k^2 \left(
\frac{\left(3 k^2-q^2\right)}{4 k^4} \vartheta (k-q)-\frac{\left(k^2-3 q^2\right)}{4 q^4} \vartheta (q-k)
\right) \nonumber
\end{eqnarray}
where $\vartheta(x)$ is the Heaviside step function. We are interested in finding an upper-bound for the latter. Therefore we compute
\begin{eqnarray}
\frac{\partial I(k,\Theta)}{\partial k}&=&\frac{ \lambda^2}{V}\Theta^2 \int_0^\infty dq\;q^3|A(q)|^2 \left(
\frac{k^2 \left(k^2-3 q^2\right) \delta (q-k)}{16 q^2}+\frac{q^2 \left(3 k^2-q^2\right)
   \delta (k-q)}{16 k^2} \right. \\
   &&\nonumber \\
   &&\left.+\frac{k^4 \left(3 q^2-2 k^2\right) \vartheta (q-k)+q^6 \theta
   (k-q)}{8 k^3 q^2}\right)  \nonumber \\
      &&\nonumber \\
   &=&\frac{ \lambda^2}{V}\Theta^2 \int_0^\infty dq\;q^3|A(q)|^2 \left(\frac{k^4 \left(3 q^2-2 k^2\right) \theta (q-k)+q^6 \vartheta
   (k-q)}{8 k^3 q^2}\right) >0 (\mbox{   unless  } |A(q)|=0 ) \nonumber
\end{eqnarray}
This is always positive, therefore $I(k, \Theta)$ is an increasing function of $k$ and it reaches its maximum at $k\rightarrow \infty$. We have
\begin{equation}
I_{max}(\Theta)=\lim_{k\rightarrow\infty}I(k,  \Theta)=\frac{ \lambda^2}{V}\Theta^2\frac{3}{16\pi} \int_0^\infty dq\;q^3|A(q)|^2
\end{equation}
and the no-pole condition is certainly verified if we impose
\begin{equation}
I_{max}( \Theta)<1.
\end{equation}
Thanks to the latter, we may  estimate  the first Gribov region and exhibit its dependence on $\Theta$.

\section{Analysis of $\sigma(k,  A, \Theta)$ in four Dimensions}\label{apb}
Let us now address the problem  in $d=4$ dimensions, in the Euclidean setting as previously done. We  work in spherical coordinates $(q, \psi, \alpha, \phi)$ with $q^2= q^\mu q_\mu, \psi,\alpha \in[0,\pi], \phi\in[0, 2\pi]$ and choose $k_\mu=(k,0,0,0)$. Since the computation is Euclidean we are free to choose the reference frame as we like, so we choose the easiest one  and the result will not the depend on the choice. As for the noncommutativity matrix $\Theta$, according to the choice we have made with   \eqn{nctheta}, the only contributing factor will be $\Theta^{12}$, it being
\begin{equation}
\Theta^{\mu\nu}q_\mu k_\nu=  - \Theta^{12}q_2 k_1= \Theta k q \sin \psi \cos \alpha
\end{equation}
with $\Theta$ a real parameter.  Eq. \eqn{nopole} becomes then
\begin{eqnarray}\label{4dsigm}
\sigma(k, A, \Theta)&=&\frac{4 \lambda^2}{3V}\int_0^\infty \frac{dq}{16\pi^4}\;q^3|A(q)|^2\int_0^{2\pi} d\phi\int_0^\pi \sin{\alpha}\,d\alpha\int_0^\pi \sin^4{\psi}\,d\psi \frac{\sin^2{\left(\frac{1}{2}\Theta q k \sin{\psi}\cos{\alpha} \right)  } }{(k^2+q^2-2k q\cos{\psi})} \nonumber\\
&=&\frac{8\pi \lambda^2}{3V}\int_0^\infty \frac{dq}{16\pi^4}\;q^3|A(q)|^2\int_0^\pi \sin{\alpha}\, d\alpha\int_0^\pi \sin^4{\psi}\, d\psi \frac{\sin^2{\left(\frac{1}{2}\Theta q k \sin{\psi}\cos{\alpha} \right)  } }{(k^2+q^2-2k q\cos{\psi})}\nonumber
\end{eqnarray}
As for the two-dimensional case,  since   $\frac{\sin^2(x)}{x^2}\leq1$ we
have
{\small{\begin{eqnarray}
\sigma(k, A, \Theta) &\leq& I(k,  A, \Theta)  =
\frac{8\pi \lambda^2}{3V}\int_0^\infty \frac{dq}{16\pi^4}\;q^3|A(q)|^2\int_0^\pi \sin{\alpha}\,d\alpha\int_0^\pi \sin^4{\psi}d\psi \frac{{\left(\frac{1}{4}\Theta^2 q^2 k^2 \sin^2{\psi}\cos^2{\alpha} \right)  } }{(k^2+q^2-2k q\cos{\psi})} \nonumber\\
&=&\frac{2\pi \lambda^2}{V}\Theta^2 \int_0^\infty \frac{dq}{16\pi^4}\;q^5|A(q)|^2\int_0^\pi \cos^2{\alpha} \sin{\alpha}\,d\alpha\int_0^\pi d\psi \frac{{\left( k^2 \sin^6{\psi}\right)  } }{(k^2+q^2-2k q\cos{\psi})} \nonumber\\
&=&\frac{4\pi \lambda^2\Theta^2}{9V}\int_0^\infty \frac{dq}{16\pi^4}\;q^5|A(q)|^2\int_0^\pi d\psi \frac{{\left( k^2 \sin^6{\psi}\right)  } }{(k^2+q^2-2k q\cos{\psi})} \\
&& \nonumber\\
&=&\frac{ \lambda^2\Theta^2}{36\pi^2 V}\int_0^\infty dq\;q^5|A(q)|^2 \left(
\frac{\left(q^4-5 k^2 q^2+10 k^4\right) \vartheta (k-q)}{16 k^4}+
\frac{k^2 \left(10 q^4-5 k^2 q^2+k^4\right) \vartheta (q-k)}{32 q^6}
\right) \nonumber
\end{eqnarray}
}}
where $\vartheta(x)$ is the Heaviside step function. We now proceed as in
the previous section, by  analysing the behaviour of the $I(k,A,\Theta)$ as a function of $k$. The derivative of $I(k, A,  \Theta)$ with respect to
$k$ is given by:
\begin{eqnarray}
\frac{\partial I}{\partial k}&=&\frac{4 \pi^2 \lambda^2\Theta^2}{9V}\int_0^\infty dq\,q^5|A(q)|^2
\left(\frac{10 k^4+q^4 -5 k^2 q^2 }{16 k^4}\delta (k-q)-\frac{10
   k^2 q^4+k^6-5 k^4 q^2}{32 q^6} \delta (q-k) \right. \nonumber\\
   &&\left.+\frac{2 q^8 \left(5 k^2-2 q^2\right) \vartheta
   (k-q)+k^6 \left(-10 k^2 q^2+3 k^4+10 q^4\right) \vartheta (q-k)}{16 k^5 q^6}\right)\nonumber\\
   &=&\frac{4 \pi^2 \lambda^2\Theta^2}{9V}\int_0^\infty dq\,q^5|A(q)|^2
   \left(\frac{2 q^8 \left(5 k^2-2 q^2\right) \vartheta
   (k-q)+k^6 \left(3 k^4+10 q^4-10 k^2 q^2\right) \vartheta (q-k)}{16 k^5
q^6}\right)\nonumber \\
   &+&\frac{ \pi \lambda^2\Theta^2}{12V} k^5 \vartheta (k) |A(k)|^2.
\end{eqnarray}
The latter is always strictly positive for $|A(q)|\neq 0$, therefore
 $I(k; A; \Theta)$ is an increasing function of $k$ and its maximum is obtained by performing  the limit $k\rightarrow\infty$:
\begin{equation}
I_{max}( A, \Theta)=\frac{5  \lambda^2\Theta^2}{288 \pi^2 V}\int_0^\infty dq\;q^5|A(q)|^2.
\end{equation}
Therefore the  no-pole condition is attained by imposing
\begin{equation}
I_{max}(A, \Theta)<1
\end{equation}
and the first Gribov region is taken into account in the functional integral  \eqn{fungen} by posing
\begin{equation}
V(\Omega)=\vartheta(1-I_{max}(A, \Theta)),
\end{equation}
with $\vartheta$ the Heaviside function. The same conclusions as for the two dimensional case apply: the horizon depends on the square of the non-commutativity parameter and is removed in the commutative limit. The estimate can be probably made more precise by maximising directly the function $\sigma(k; A, \Theta)$, provided one is able to overcome the even more challenging technical difficulties in dealing with the higher-dimensional
integral
\eqn{4dsigm}.

\section{ What to expect of an all order computation?}\label{appALLORDERS}
 Here we discuss perspectives on the extension of analysis to include  higher order terms. We start with equation \eqref{eq654}
\begin{eqnarray}
\mathbb{A}^{2}_{kk}&=&\int \frac{d^d q}{(2\pi)^{d}}\, \langle k\vert \mathbb{A}\vert q\rangle\langle q\vert \mathbb{A}\vert k \rangle \nonumber\\
&=&-4\int \frac{d^d q}{(2\pi)^{d}}\sin\left(\frac{1}{2}\Theta_{\rho\sigma}k^{\rho}q^{\sigma}\right)\tilde A_{\mu}(k-q)\frac{q_{\mu}}{q^2}\sin\left(\frac{1}{2}\Theta_{\rho\sigma}q^{\rho}k^{\sigma}\right)\tilde A_{\nu}(q-k)\frac{k_{\nu}}{k^2}.
\end{eqnarray}
One can easily proceed and write a general formal expression for the term of order $n$
\begin{eqnarray}
	\mathbb{A}^{n}_{kk}&=&\int \frac{d^d q_1}{(2\pi)^{d}}\,\cdots  \int \frac{d^d q_{n-1}}{(2\pi)^{d}}\, \langle k\vert \mathbb{A}\vert q_1\rangle\langle q_1\vert \mathbb{A}\vert q_2 \rangle \cdots \langle q_{n-2}\vert \mathbb{A}\vert q_{n-1}\rangle\langle q_{n-1}\vert \mathbb{A}\vert k \rangle \nonumber\\
	&=&(-2i)^n\int \frac{d^d q_1}{(2\pi)^{d}}\,\cdots  \int \frac{d^d q_{n-1}}{(2\pi)^{d}}\,\sin\left(\frac{1}{2}\Theta_{\rho\sigma}k^{\rho}q_1^{\sigma}\right)\sin\left(\frac{1}{2}\Theta_{\rho\sigma}q_1^{\rho}q_2^{\sigma}\right)\cdots \sin\left(\frac{1}{2}\Theta_{\rho\sigma}q_{n-1}^{\rho}k^{\sigma}\right) \nonumber\\
	&&\tilde A_{\mu}(k-q_1)\tilde A_{\sigma_1}(q_1-q_2)\cdots \tilde A_{\sigma_{\nu}}(q_{n-1}-k)
	\frac{q_1^{\mu}q_2^{\sigma_1}\cdots q_{n-1}^{\sigma_{n-1}}}{q_1^2 q_2^2\cdots q_{n-1}^2}\frac{k_{\nu}}{k^2}\nonumber\\
	&=&(-2i)^n \frac{k_{\mu}k_{\nu}}{k^2}\int \frac{d^d q_1}{(2\pi)^{d}}\,\cdots  \int \frac{d^d q_{n-1}}{(2\pi)^{d}}\,\sin\left(\frac{1}{2}\Theta_{\rho\sigma}k^{\rho}q_1^{\sigma}\right)\sin\left(\frac{1}{2}\Theta_{\rho\sigma}q_1^{\rho}q_2^{\sigma}\right)\nonumber\\
	&&\cdots \sin\left(\frac{1}{2}\Theta_{\rho\sigma}q_{n-1}^{\rho}k^{\sigma}\right) \tilde A_{\mu}(k-q_1)\tilde A_{\sigma_1}(q_1-q_2)\cdots \tilde A_{\nu}(q_{n-1}-k)
	\frac{q_2^{\sigma_1}\cdots q_{n-1}^{\sigma_{n-1}}}{q_1^2 q_2^2\cdots q_{n-1}^2}.
\end{eqnarray}
where the Landau gauge property $q_{\mu}A_{\mu}(k-q)=k_{\mu}A_{\mu}(k-q)$  was used. It then follows that the exact ghost form factor is given by
\begin{eqnarray}
	\sigma(k,A) = \sum_n \sigma^{(n)}(k,A),
\end{eqnarray}
where
\begin{eqnarray}
	\sigma^{(n)}(k,A) &=&\frac{(-2i\lambda)^n}{V} \frac{k_{\mu}k_{\nu}}{k^2}\int \frac{d^d q_1}{(2\pi)^{d}}\,\cdots  \int \frac{d^d q_{n-1}}{(2\pi)^{d}}\,\sin\left(\frac{1}{2}\Theta_{\rho\sigma}k^{\rho}q_1^{\sigma}\right)\sin\left(\frac{1}{2}\Theta_{\rho\sigma}q_1^{\rho}q_2^{\sigma}\right)\nonumber\\
	&&\cdots \sin\left(\frac{1}{2}\Theta_{\rho\sigma}q_{n-1}^{\rho}k^{\sigma}\right) \tilde A_{\mu}(k-q_1)\tilde A_{\sigma_1}(q_1-q_2)\cdots \tilde A_{\nu}(q_{n-1}-k)
	\frac{q_2^{\sigma_1}\cdots q_{n-1}^{\sigma_{n-1}}}{q_1^2 q_2^2\cdots q_{n-1}^2}.\nonumber\\
\end{eqnarray}

Now, we worked semiclassically because we stopped at the lowest nontrivial order, which is the quadratic order $\sigma^{(2)}(k,A)$. But in order to go to higher orders one cannot ignore interactions and the problem must be addressed as discussed around \eqref{1PI}, that is, one must deal with $\sigma(k) = \langle \sigma(k,A) \rangle_{1PI}$. This amounts to compute at each order $n$ the following expression
\begin{eqnarray}
	\sigma^{(n)}(k)&=&\frac{(-2i\lambda)^n}{V} \frac{k_{\mu}k_{\nu}}{k^2}\int \frac{d^d q_1}{(2\pi)^{d}}\,\cdots  \int \frac{d^d q_{n-1}}{(2\pi)^{d}}\,\sin\left(\frac{1}{2}\Theta_{\rho\sigma}k^{\rho}q_1^{\sigma}\right)\sin\left(\frac{1}{2}\Theta_{\rho\sigma}q_1^{\rho}q_2^{\sigma}\right)\nonumber\\
	&&\cdots \sin\left(\frac{1}{2}\Theta_{\rho\sigma}q_{n-1}^{\rho}k^{\sigma}\right) \langle \tilde A_{\mu}(k-q_1)\tilde A_{\sigma_1}(q_1-q_2)\cdots \tilde A_{\nu}(q_{n-1}-k)\rangle_{1PI}
	\frac{q_2^{\sigma_1}\cdots q_{n-1}^{\sigma_{n-1}}}{q_1^2 q_2^2\cdots q_{n-1}^2}.\nonumber\\
\end{eqnarray}
Therefore, one must be able to compute the $n$-point function $\langle \tilde A_{\mu}(k-q_1)\tilde A_{\sigma_1}(q_1-q_2)\cdots \tilde A_{\nu}(q_{n-1}-k)\rangle_{1PI}$, a not so easy task indeed.

The desired results rests on a proper understanding of the behavior of  $\sigma(k)$ as a function of $k$. In the semiclassical quadractic case, analyzed in the paper, the result indicating no qualitative change in the propagator can be traced to the presence of sines in the expression for  $\sigma^{(2)}(k)$, that leads to a stronger bound of the integral suppressing the Gribov effects. The expression for $\sigma^{(n)}(k)$  also displays sines and these should similarly provide stronger bounds in comparison to the non-abelian commutative case. So, one is led to speculate that Gribov effects could be suppressed at higher order as well. But, on the other hand, this will depend on the behavior of the gauge field $n$-point function in a complicated self-consistent way, thus making it difficult to settle for a definite answer as to what happens at higher orders.
\end{appendices}

redpajama_set_name: RedPajamaArXiv