# Laser-induced electron emission from Au nanowires: a probe for orthogonal
polarizations

Eric R. Jones<EMAIL_ADDRESS>Department of Physics and
Astronomy, University of Nebraska-Lincoln, Lincoln, Nebraska 68588, USA  Wayne
C. Huang  Department of Physics and Astronomy, University of Nebraska-Lincoln,
Lincoln, Nebraska 68588, USA  Gobind Basnet  Department of Physics, Kansas
State University, Manhattan, Kansas 66506, USA  Bret N. Flanders  Department
of Physics, Kansas State University, Manhattan, Kansas 66506, USA  Herman
Batelaan<EMAIL_ADDRESS>Department of Physics and Astronomy, University of
Nebraska-Lincoln, Lincoln, Nebraska 68588, USA

###### Abstract

Photoelectron field emission, induced by femtosecond laser pulses focused on
metallic nanotips, provides spatially coherent and temporally short electron
pulses. Properties of the photoelectron yield give insight into both the
material properties of the nanostructure and the exciting laser focus.
Ultralong nanoribbons, grown as a single crystal attached to a metallic taper,
are sources of electron field emission that have not yet been characterized.
In this report, photoemission from gold nanoribbon samples is studied and
compared to emission from tungsten and gold tips. We observe that the emission
from sharp tips generally depends on one transverse component of the exciting
laser field, while the emission of a blunted nanoribbon is found to be
sensitive to both components. We propose that this property makes
photoemission from nanoribbons a candidate for position-sensitive detection of
the longitudinal field component in a tightly focused beam.

A consequence of tightly focusing a beam of light is that the beam will become
longitudinally polarized near the focus.Richards and Wolf (1959); Cicchitelli,
Hora, and Postle (1990); Dorn, Quabis, and Leuchs (2003) Longitudinally
polarized beams are desirable because their focal widths can be below the
diffraction limit, and the longitudinal component of the electric field does
not contribute to the energy flow along the direction of beam propagation.Wang
_et al._ (2008); Yu  _et al._ (2016) These features find applications in high-
resolution optical microscopy,Yoshiki  _et al._ (2007); Terakado, Watanabe,
and Kano (2009) optical data storage,Kim  _et al._ (2007) particle
trapping,Zhan (2004, 2012) charged particle acceleration,Gupta  _et al._
(2007) material ablation,Hnatovsky  _et al._ (2011) and pushing the high-
intensity frontier. The longitudinal field component of a laser focus has been
characterized in situ by atomic fluorescenceNovotny  _et al._ (2001) and near-
field microscopy,Yu  _et al._ (2016) and ex situ via imaging of material
damageHnatovsky  _et al._ (2011) and atomic force microscopy of thin film
deformation.Gilbert  _et al._ (2006) As the in situ methods of characterizing
the longitudinal field are limited by intensity or to a resonant wavelength, a
flexible alternative would be preferable.

Photoelectron field emission, induced by focusing femtosecond laser pulses
onto sharp metallic tapers with nanometric radii of curvature, Hommelhoff  _et
al._ (2006); Ropers  _et al._ (2007); Barwick  _et al._ (2007) has a broad
range of applications. Temporally short electron wave packetsKrüger, Schenk,
and Hommelhoff (2011) with high spatial coherenceHommelhoff  _et al._ (2009);
Ehberger  _et al._ (2015) can be achieved with moderate intensities. Tip
sources have thus been integrated into electron microscopes to obtain sub-
micron spatial resolution with femtosecond timing.Barwick, Flannigan, and
Zewail (2009); Gulde  _et al._ (2014); Piazza  _et al._ (2015) Femtosecond
electron pulses have been used to study fundamental quantum mechanics, as in
testing the existence of forces in the Aharonov-Bohm effect,Caprez, Barwick,
and Batelaan (2007) and observing diffraction in time.Krüger, Schenk, and
Hommelhoff (2011); Krüger  _et al._ (2012) Electron emission from nanotips is
obtained for a range of laser intensities and wavelengths.Herink, Wimmer, and
Ropers (2014); Förster  _et al._ (2016); Huang  _et al._ (2017) Nanotip
emission might then provide an alternative for characterizing the longitudinal
component of a laser focus. However, as nanotip emission is dominated by a
single transverse component of the polarization of the exciting field, it has
not yet been utilized as a detector of longitudinal fields.

Ultralong gold nanoribbons,Basnet  _et al._ (2017) grown by the method of
directed electrochemical nanowire assembly (DENA),Ozturk  _et al._ (2007);
Flanders (2012)are as of yet unstudied sources for ultrafast electron
photoemission. The DENA methodology results in single-crystalline samples, as
confirmed by electron diffraction. Previous studies into the optical damage
threshold of similar nanowire samples indicated that single-crystalline gold
nanowires could tolerate high laser pulse peak intensities before melting, but
that they exhibited long cooling times.Summers  _et al._ (2014) From these
studies, the exchange and dissipation of heat between the nanowire electrons
and the crystalline lattice could be modelled. A logical next step would then
be to characterize electrons emitted from nanoribbons in response to
ultrashort pulse illumination. Photoelectron emission is known to carry
information on material effects, such as plasmonic dynamics and laser
heating.Vogelsang  _et al._ (2015); Kealhofer  _et al._ (2012) Temporally
short electron pulses have been observed from plasmonic nanostructures,Irvine,
Dechant, and Elezzabi (2004); Dombi  _et al._ (2010) and the photoemission
yield has been used as a sensitive probe of the plasmonic field enhancement
from nanostructures.Rácz  _et al._ (2017) We measured photoelectron emission
from nanoribbons in an attempt to determine their plasmonic or thermal
properties.Vogelsang  _et al._ (2015); Rácz  _et al._ (2017); Kealhofer  _et
al._ (2012) In the following, nanoribbon samples are characterized by
photoelectron emission and compared to standard single-crystalline gold and
tungsten samples. The resulting photoelectron emission spectra reveal that
nanoribbons can be employed as position- and polarization-sensitive detectors
within a laser focus, providing a potential in situ sub-wavelength probe for
longitudinal polarization.

A schematic for the system used to characterize tip samples is given in Fig.
1.

Figure 1:  Experimental schematic showing pump/probe configuration for tip and
nanoribbon photoelectron emission.

The intensity of the output from a Ti:Saph oscillator (Spectra Physics
Tsunami, $80\text{\,}\mathrm{MHz}$ repetition rate, $800\text{\,}\mathrm{nm}$
central wavelength, $100\text{\,}\mathrm{fs}$ pulse width) is controlled by a
variable attenuator (VA), which consists of a $\lambda/2$ plate and a Brewster
window. The attenuated beam is split into a pump and a probe beam by a
balanced Mach-Zehnder interferometer (IFM). Photoelectron emission is observed
in pump/probe and single beam experiments. The time delay between pump and
probe pulses, $\uptau$, can be varied from $-4\text{\,}\mathrm{ps}$ to
$4\text{\,}\mathrm{ps}$ manually by a micrometer. The additive ratio, measured
as a function of the delay $\uptau$, is obtained by dividing the emission rate
measured with both pump and probe pulses delivered to the samples by the sum
of the individual rates due to the pump and probe pulses. A second $\lambda/2$
plate (HWP) rotates the polarization of both beams prior to delivery to the
experimental chamber. A rotational stepper motor is used to scan the beam
power and the polarization.

The experimental chamber, which is detailed in Ref. Barwick  _et al._ , 2007,
is maintained at $2\text{\times}{10}^{-7}$ Torr. The beams are focused within
the chamber by an off-axis parabolic mirror (OAPM) to a full-width half
maximum of $3.6\text{\,}\mathrm{\SIUnitSymbolMicro m}$. A 3-axis stage,
coupled to the chamber by flexible bellows, positions tip samples into the
focus. Mounted tip samples were biased at $-100\text{\,}\mathrm{V}$, as this
was lower than the threshold for Fowler-Nordheim field emission. Electrons
were collimated through two $4\text{\,}\mathrm{mm}$ apertures before being
detected by a microchannel plate (MCP).

Electron pulses from the MCP were amplified and discriminated. Discriminator
pulses were counted by a multichannel scaler, and used as the start trigger
for a time-to-amplitude converter (TAC). The output reference signal from the
Ti:Saph oscillator was used as the TAC stop trigger to measure the arrival
time of electrons. Timing spectra were obtained by sending the TAC output
pulses to a multichannel analyzer (MCA).

Nanoribbon samples were prepared using the DENA methodology. Ozturk  _et al._
(2007); Flanders (2012) Nanoribbon samples are reported to have a thickness of
$40\text{\,}\mathrm{nm}$, and widths ranging from
$130\text{\,}\mathrm{nm}$–$360\text{\,}\mathrm{nm}$ along the length of the
wire. The tip can have a radius of curvature of
$10\text{\,}\mathrm{nm}$.Basnet  _et al._ (2017) These dimensions can be
tailored during the growth process to make nanoribbons that are well-suited
for photoemission.

In order to distinguish which photoemission properties arise due to the
material or geometry of the nanoribbon samples, single-crystal tungsten (W)
and gold (Au) tips were prepared for comparison. Samples of W wire
($200\text{\,}\mathrm{\SIUnitSymbolMicro m}$ diameter) were annealedGreiner
and Kruse (2007), and then etched via the lamella drop-off method.Müller  _et
al._ (1999) Au wire samples ($99.95$% purity, Ted Pella,
$200\text{\,}\mathrm{\SIUnitSymbolMicro m}$ diameter) were annealedRoy,
Williams, and Mingard (2010), and then etched as according to Refs. Gingery
and Bühlmann, 2007; Eligal  _et al._ , 2009.

The tip and nanoribbon samples were mounted to SEM pin stubs with silver
paste. SEM images of the samples are given in Fig. 2(a) along with plots of
the beam focus (red) and intensity profile (white) as measured by
photoemission. From left to right is shown W (I), Au (II), an undamaged
$23\text{\,}\mathrm{\SIUnitSymbolMicro m}$ Au nanoribbon (III), and an
$11\text{\,}\mathrm{\SIUnitSymbolMicro m}$ Au nanoribbon (IV) obtained after
the $23\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon was blunted during
pump/probe experiments. Images were taken before and after experimental
characterization to determine the extent of damage due to laser illumination.

Photoemission data from single-beam experiments are shown in Fig. 2(b)–(d). In
Fig. 2(b), the emission rate is shown as samples were translated through the
laser focus. The W (I) and Au (II) tip samples show emission localized at the
tip apex only, while the nanoribbon samples (III and IV) can emit from
multiple locations along their length. This feature confirmed that a
nanoribbon remained attached to the Au substrate after imaging and transfer to
the experimental chamber. Thin lines between data points serve as a guide to
the eye. Fig. 2(c) shows the dependence of electron yield on the average power
of the beam, plotted on a log-log scale. The value, $n$, of the power
dependence ($\propto I^{n}$), is often used to identify the emission process
of a tip. The W tip and $23\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon
have slopes of $n=3$, while the Au tip is found to have a slope of $n=3$ for
low power, and $n=5$ for higher power. Such behavior, that is, the increase in
power law slope with increasing laser power, has been observed in W tips and
studies of above threshold photoemission.Barwick  _et al._ (2007); Schenk,
Krüger, and Hommelhoff (2010) The $11\text{\,}\mathrm{\SIUnitSymbolMicro m}$
nanoribbon has a slope of $n=5$. With these values of $n$, the position
dependence of the samples in Fig. 2 (b) can be fit with a gaussian function to
determine the size of the focal waist (bold lines). The focal waist has a
fitted full width at half maximum of $3.6\text{\,}\mathrm{\SIUnitSymbolMicro
m}$ from the W data. Fig. 2(d) shows the variation of emission rate as the
polarization of the beam is rotated by a $\lambda/2$ plate. The high contrast
$None$ spaced peaks in the tip samples and
$23\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon support that the sample
geometry is well-defined with respect to the laser polarization in the focus,
and that the emission process is dominated by a preferred laser polarization.
The broadened

Figure 2:  Tip sample characteristics. (a) Scanning electron microscope (SEM)
images of the tip samples: from left to right is shown annealed W (I),
annealed Au (II), a $23\text{\,}\mathrm{\SIUnitSymbolMicro m}$ Au nanoribbon
(III), and an $11\text{\,}\mathrm{\SIUnitSymbolMicro m}$ Au nanoribbon (IV).
The laser spot size (red) and intensity profile (white) have a full width half
maximum of $3.6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ as fit from the W data
in (b) and (c). (b) Scaled electron counts as function of tip position in
focus. Gaussian fits to the data are shown in bold lines. Thin lines are
guides to the eye. (c) Power dependence of electron emission. (d) Polarization
dependence of electron emission.

peaks and reduced contrast of the $11\text{\,}\mathrm{\SIUnitSymbolMicro m}$
nanoribbon electron yield (IV) indicate that the emission process depends on
both transverse components of the exciting field. A feature consistent with
multiphoton emission is that the power law slope, $n$, will agree with the
polarization dependence on the emission rate. The emission rate is $\propto
cos^{2n}(\theta)$ in a multiphoton emission model, where $\theta$ is the
polarization of the laser relative to the tip direction. The W and Au
polarization data agree with a fit with $n=3$, while both nanoribbon samples
require a combination of $n=1$ and $n=3$. This peculiarity of both nanoribbon
samples indicates a deviation from typical multiphoton emission.

Pump/probe and single beam experiments with the TAC/MCA configuration revealed
the timing features of electron emission. The additive ratio of emission from
the samples is plotted in Fig. 3(a) as a function of $\uptau$. The
polarization of the focus was chosen for the optimum electron yield from each
sample. A ratio of 1 indicates that the emission yield from the probe pulses
are independent from the pump pulses. A ratio significantly greater than 1
indicates emission processes that are slower than the time delay between
pulses.Barwick  _et al._ (2007) When the delay is shorter than the pulse
duration, the additive ratio can vary due to interference between the pulses.
The W tip (blue triangles), Au tip (green triangles), and
$23\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon all have additive ratios
that are close to 1 when the pulse delay is outside of the
$\pm$200\text{\,}\mathrm{fs}$$ interference window, so the emission processes
are as fast as the $100\text{\,}\mathrm{f}\mathrm{s}$ laser pulse duration and
thus prompt. The $11\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon (black
squares) has an average additive ratio of 14.9 for delays longer than the
pump/probe interference window, therefore the process is not prompt.
Measurement of the ratio for delays with high constructive interference were
avoided to prevent damage to the Au tip and
$11\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon samples.

Normalized time spectra of electron emission from tip samples are plotted in
Fig. 3(b). Shown, grouped from left to right, are the spectra from the W tip
(blue line), Au tip (green line), $23\text{\,}\mathrm{\SIUnitSymbolMicro m}$
nanoribbon (red line), and the $11\text{\,}\mathrm{\SIUnitSymbolMicro m}$
nanoribbon (black line). The peak separation for each sample shows the
$13\text{\,}\mathrm{ns}$ pulse separation of the oscillator. Sharp peaks
indicate pulsed electron emission, while sustained signal after the laser
pulse indicates background emission. The
$11\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon has a significant
background as compared to the other samples, indicating that electron emission
is continuing after the exciting laser pulse is gone. The emission process of
the $11\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon is therefore ruled
out as purely multiphoton, and is likely due to both multiphoton and laser
heating of the nanoribbon structure.

To further investigate this feature, the pulsed and background contributions
to emission from the $11\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon are
plotted as a function of $\lambda/2$ angle in Fig. 3(c).

Figure 3:  Pump/probe, timing, and polarization control of emission processes.
(a) Additive ratio of electron emission vs. pulse delay. Color and marker
conventions follow Fig. 2. Ratio values of 1 and 14.9 are highlighted by
dashed black lines. (b) Normalized time spectra of electron emission from tip
samples. Grouped from left to right are the W tip (blue), Au tip (green),
$23\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon (red), and
$11\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon (black). (c)
Polarization dependence of summed counts of time spectra background (blue
hatched) and pulsed (red hatching) emission vs. $\lambda/2$ angle for the
$11\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon. Colored lines are
guides to the eye. The inset shows the background and pulsed signal for a
single timing peak taken at $\lambda/2$ angle $=124^{\circ}$.

Time spectra were recorded for each $\lambda/2$ angle. The process for
dividing each spectrum into pulsed and background contributions is illustrated
in the inset of Fig. 3(c), which shows a portion of the time spectrum taken at
$\lambda/2$ angle of $None$. The red hatched region of the inset indicates the
pulsed contribution, and the blue hatched region indicates the background. The
background regions are defined by taking linear fits to the tails of the
timing spectra, and extending those fits to the rising edges of the timing
peaks. This procedure is performed for two $13\text{\,}\mathrm{ns}$ oscillator
periods. The counts in the background regions are summed, giving the data
points marked by the blue hatched squares. The background contributions are
then subtracted from the total counts in each spectrum. This results in the
data points marked by the red hatched squares. The maximum emission for the
background process occurs at a $\lambda/2$ angle that is shifted relative to
the pulsed process. For example, for a $\lambda/2$ angle of $None$, the
photoemission signal is dominated by the pulsed process, as the number of
counts in the tails of the electron time of flight spectra is low. For a
$\lambda/2$ angle of $None$, the photoemission signal is dominated by the
background process, which is marked by a comparatively high number of counts
in the delayed tail of the electron time of flight spectrum. The explanation
for the variation in electron signal is likely due to the nanoribbon absorbing
more of the incident laser pulse energy at $None$, and less at $None$. The
higher degree of energy absorption in the
$11\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon requires a higher degree
of energy dissipation, which occurs by an additional process–thermal
dissipation–that is much slower than multiphoton-driven processes. This
interpretation is consistent with the results of Ref. Kealhofer  _et al._ ,
2012, which demonstrated polarization control of thermally enhanced
photoemission from nanotips. The polarization control of these multiphoton
(pulsed) and thermal (background) processes indicates that they respond to
different components of the laser field in the focus.

Ultralong Au nanoribbons are unique nanostructures for the study of electron,
thermal, and plasmonic transport by laser-induced electron emission. Previous
work indicated that nanoribbons are resilient to damage by laser intensities
on the order of $\text{\,}\mathrm{TW}\text{/}{\mathrm{cm}}^{2}$,Summers  _et
al._ (2014) but we found that damage can occur with lower intensities. This
observation resulted in the $11\text{\,}\mathrm{\SIUnitSymbolMicro m}$
nanoribbon, which differed from the original
$23\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon by having a less defined
apex. An immediate consequence of the laser damage to the
$23\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon was that the shorter
$11\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon required higher incident
laser power to emit, which is evident in Fig. 2(c). The change in the
nanoribbon’s apex geometry also coincided with emission from the
$11\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon being superadditive for
emission rates similar to the other samples. The superadditive emission in
pump/probe experiments was accompanied by delayed tails in the photoelectron
time of flight spectra in the $11\text{\,}\mathrm{\SIUnitSymbolMicro m}$
nanoribbon. Superadditive and delayed emission are not consistent with
plasmon-induced field emission, as plasmonic emission is reported to have a
standard pump/probe cross-correlation and thus is as fast as the exciting
laser pulses.Dombi  _et al._ (2010); Vogelsang  _et al._ (2015) Such tails
were not observed at the $23\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon
apex before it was damaged. These features, the superadditive emission and the
delayed arrival times, are consistent with the
$11\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon being more susceptible
to laser heating than the $23\text{\,}\mathrm{\SIUnitSymbolMicro m}$
nanoribbon. The higher susceptibility to heating made the
$11\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon sensitive to both
transverse components of the focused laser field, as shown in Fig. 3(c). The
sensitivity to both transverse components of the focused field is a unique
feature of the thermal emission of the nanoribbon that was not observed in W
or Au tips. The cone structure of the nanotip samples leads to much faster
cooling times than the nanoribbon samples. This suggests that a tailored
nanoribbon could be oriented to probe other polarization components in a focus
as well. A schematic for an oriented nanoribbon as a probe of the longitudinal
component of a focused non-paraxial beam is given in Fig. 4.

Figure 4:  (a) Micrographs of a long nanoribbon and sample (IV) (inset). The
nanoribbon is twisted to highlight that it is ribbon shaped. (b) Schematic of
nanoribbon as a detector of longitudinal ($E_{L}$) and transverse ($E_{T}$)
electric fields.

Shown in Fig. 4(a) is a high-resolution SEM image of a nanoribbon and our
$11\text{\,}\mathrm{\SIUnitSymbolMicro m}$ nanoribbon as the inset. The
nanoribbon is twisted in the high-resolution image to highlight its shape.
With the nanoribbon oriented as shown in Fig. 4(b), the broad side of the
nanoribbon would heat depending on the strength of the longitudinal component
of the exciting field, $E_{L}$, and the delayed electron yield would therefore
depend on $E_{L}$, while the peaked electron yield could depend on a
combination of the transverse component, $E_{T}$, and also $E_{L}$. The
capability to distinguish the slow thermal electron yield at the nanoribbon
apex makes electron emission from a nanoribbon a subwavelength probe of
orthogonal polarizations.

###### Acknowledgements.

We gratefully acknowledge funding by NSF EPSCoR NE-KS Track-II, Award No. EPS
1430519, and NSF grant No. 1602755.

## References

  * Richards and Wolf (1959) B. Richards and E. Wolf, Proc. R. Soc. London, Ser. A 253, 358 (1959).
  * Cicchitelli, Hora, and Postle (1990) L. Cicchitelli, H. Hora, and R. Postle, Phys. Rev. A 41, 3727 (1990).
  * Dorn, Quabis, and Leuchs (2003) R. Dorn, S. Quabis, and G. Leuchs, J. Mod. Opt. 50, 1917 (2003).
  * Wang  _et al._ (2008) H. Wang, L. Shi, B. Lukyanchuk, C. Sheppard, and C. T. Chong, Nat. Photonics 2, 501 (2008).
  * Yu  _et al._ (2016) A.-p. Yu, G. Chen, Z.-h. Zhang, Z.-q. Wen, L.-r. Dai, K. Zhang, S.-l. Jiang, Z.-x. Wu, Y.-y. Li, C.-t. Wang, and X.-g. Luo, Sci. Rep. 6, 38859 (2016).
  * Yoshiki  _et al._ (2007) K. Yoshiki, K. Ryosuke, M. Hashimoto, N. Hashimoto, and T. Araki, Opt. Lett. 32, 1680 (2007).
  * Terakado, Watanabe, and Kano (2009) G. Terakado, K. Watanabe, and H. Kano, Appl. Opt. 48, 1114 (2009).
  * Kim  _et al._ (2007) W.-C. Kim, N.-C. Park, Y.-J. Yoon, H. Choi, and Y.-P. Park, Opt. Rev. 14, 236 (2007).
  * Zhan (2004) Q. Zhan, Opt. Express 12, 3377 (2004).
  * Zhan (2012) Q. Zhan, Opt. Express 20, 6058 (2012).
  * Gupta  _et al._ (2007) D. N. Gupta, N. Kant, D. E. Kim, and H. Suk, Phys. Lett. A 368, 402 (2007).
  * Hnatovsky  _et al._ (2011) C. Hnatovsky, V. Shvedov, W. Krolikowski, and A. Rode, Phys. Rev. Lett. 106, 123901 (2011).
  * Novotny  _et al._ (2001) L. Novotny, M. R. Beversluis, K. S. Youngworth, and T. G. Brown, Phys. Rev. Lett. 86, 5251 (2001).
  * Gilbert  _et al._ (2006) Y. Gilbert, R. Bachelot, P. Royer, A. Bouhelier, G. P. Wiederrecht, and L. Novotny, Opt. Lett. 31, 613 (2006).
  * Hommelhoff  _et al._ (2006) P. Hommelhoff, Y. Sortais, A. Aghajani-Talesh, and M. A. Kasevich, Phys. Rev. Lett. 96, 077401 (2006).
  * Ropers  _et al._ (2007) C. Ropers, D. R. Solli, C. P. Schulz, C. Lienau, and T. Elsaesser, Phys. Rev. Lett. 98, 043907 (2007).
  * Barwick  _et al._ (2007) B. Barwick, C. Corder, J. Strohaber, N. Chandler-Smith, C. Uiterwaal, and H. Batelaan, New J. Phys. 9, 142 (2007).
  * Krüger, Schenk, and Hommelhoff (2011) M. Krüger, M. Schenk, and P. Hommelhoff, Nature 475, 78 (2011).
  * Hommelhoff  _et al._ (2009) P. Hommelhoff, C. Kealhofer, A. Aghajani-Talesh, Y. R. Sortais, S. M. Foreman, and M. A. Kasevich, Ultramicroscopy 109, 423 (2009), iFES 2008\.
  * Ehberger  _et al._ (2015) D. Ehberger, J. Hammer, M. Eisele, M. Krüger, J. Noe, A. Högele, and P. Hommelhoff, Phys. Rev. Lett. 114, 227601 (2015).
  * Barwick, Flannigan, and Zewail (2009) B. Barwick, D. J. Flannigan, and A. H. Zewail, Nature 462, 902 (2009).
  * Gulde  _et al._ (2014) M. Gulde, S. Schweda, G. Storeck, M. Maiti, H. K. Yu, A. M. Wodtke, S. Schäfer, and C. Ropers, Science 345, 200 (2014).
  * Piazza  _et al._ (2015) L. Piazza, T. T. A. Lummen, E. Quiñonez, Y. Murooka, B. W. Reed, B. Barwick, and F. Carbone, Nat. Commun. 6, 6407 (2015).
  * Caprez, Barwick, and Batelaan (2007) A. Caprez, B. Barwick, and H. Batelaan, Phys. Rev. Lett. 99, 210401 (2007).
  * Krüger  _et al._ (2012) M. Krüger, M. Schenk, M. Förster, and P. Hommelhoff, J. Phys. B 45, 074006 (2012).
  * Herink, Wimmer, and Ropers (2014) G. Herink, L. Wimmer, and C. Ropers, New J. Phys. 16, 123005 (2014).
  * Förster  _et al._ (2016) M. Förster, T. Paschen, M. Krüger, C. Lemell, G. Wachter, F. Libisch, T. Madlener, J. Burgdörfer, and P. Hommelhoff, Phys. Rev. Lett. 117, 217601 (2016).
  * Huang  _et al._ (2017) W. C.-W. Huang, M. Becker, J. Beck, and H. Batelaan, New J. Phys. 19, 023011 (2017).
  * Basnet  _et al._ (2017) G. Basnet, K. R. Panta, P. S. Thapa, and B. N. Flanders, Appl. Phys. Lett. 110, 073106 (2017).
  * Ozturk  _et al._ (2007) B. Ozturk, B. N. Flanders, D. R. Grischkowsky, and T. D. Mishima, Nanotechnology 18, 175707 (2007).
  * Flanders (2012) B. N. Flanders, Mod. Phys. Lett. B 26, 1130001 (2012).
  * Summers  _et al._ (2014) A. M. Summers, A. S. Ramm, G. Paneru, M. F. Kling, B. N. Flanders, and C. A. Trallero-Herrero, Opt. Express 22, 4235 (2014).
  * Vogelsang  _et al._ (2015) J. Vogelsang, J. Robin, B. J. Nagy, P. Dombi, D. Rosenkranz, M. Schiek, P. Groß, and C. Lienau, Nano Lett. 15, 4685 (2015).
  * Kealhofer  _et al._ (2012) C. Kealhofer, S. M. Foreman, S. Gerlich, and M. A. Kasevich, Phys. Rev. B 86, 035405 (2012).
  * Irvine, Dechant, and Elezzabi (2004) S. E. Irvine, A. Dechant, and A. Y. Elezzabi, Phys. Rev. Lett. 93, 184801 (2004).
  * Dombi  _et al._ (2010) P. Dombi, S. E. Irvine, P. Rácz, M. Lenner, N. Kroó, G. Farkas, A. Mitrofanov, A. Baltuška, T. Fuji, F. Krausz, and A. Y. Elezzabi, Opt. Express 18, 24206 (2010).
  * Rácz  _et al._ (2017) P. Rácz, Z. Pápa, I. Márton, J. Budai, P. Wróbel, T. Stefaniuk, C. Prietl, J. R. Krenn, and P. Dombi, Nano Lett. 17, 1181 (2017).
  * Greiner and Kruse (2007) M. Greiner and P. Kruse, Rev. Sci. Instrum. 78, 026104 (2007).
  * Müller  _et al._ (1999) A. Müller, F. Müller, M. Hietschold, F. Demming, J. Jersch, and K. Dickmann, Rev. Sci. Instrum. 70, 3970 (1999).
  * Roy, Williams, and Mingard (2010) D. Roy, C. M. Williams, and K. Mingard, J. Vac. Sci. Technol. B Nanotechnol. Microelectron. 28, 631 (2010).
  * Gingery and Bühlmann (2007) D. Gingery and P. Bühlmann, Rev. Sci. Instrum. 78, 113703 (2007).
  * Eligal  _et al._ (2009) L. Eligal, F. Culfaz, V. McCaughan, N. I. Cade, and D. Richards, Rev. Sci. Instrum. 80, 033701 (2009).
  * Schenk, Krüger, and Hommelhoff (2010) M. Schenk, M. Krüger, and P. Hommelhoff, Phys. Rev. Lett. 105, 257601 (2010).


e1e-mail<EMAIL_ADDRESS>e2e-mail:
<EMAIL_ADDRESS>e3e-mail: helios.sanchis-
<EMAIL_ADDRESS>
11institutetext: Institut für Theoretische Physik, Justus-Liebig-Universität
Giessen, 35392 Giessen, Germany 22institutetext: Helmholtz Forschungsakademie
Hessen für FAIR (HFHF), GSI Helmholtzzentrum für Schwerionenforschung, Campus
Gießen, 35392 Gießen, Germany 33institutetext: Institute of Physics,
University of Graz, NAWI Graz, Universitätsplatz 5, 8010 Graz, Austria
44institutetext: Silicon Austria Labs GmbH, Inffeldgasse 33, 8010 Graz,
Austria

# Spectrum of scalar and pseudoscalar glueballs from functional methods

Markus Q. Hubere1,addr1  Christian S. Fischere2,addr1,addr2  Hèlios Sanchis-
Alepuze3,addr3,addr4

###### Abstract

We provide results for the spectrum of scalar and pseudoscalar glueballs in
pure Yang-Mills theory using a parameter-free fully self-contained truncation
of Dyson-Schwinger and Bethe-Salpeter equations. The only input, the scale, is
fixed by comparison with lattice calculations. We obtain ground state masses
of $1.9\,\text{GeV}$ and $2.6\,\text{GeV}$ for the scalar and pseudoscalar
glueballs, respectively, and $2.6\,\text{GeV}$ and $3.9\,\text{GeV}$ for the
corresponding first excited states. This is in very good quantitative
agreement with available lattice results. Furthermore, we predict masses for
the second excited states at $3.7\,\text{GeV}$ and $4.3\,\text{GeV}$. The
quality of the results hinges crucially on the self-consistency of the
employed input. The masses are independent of a specific choice for the
infrared behavior of the ghost propagator providing further evidence that this
only reflects a nonperturbative gauge completion.

###### Keywords:

glueballs, bound states, correlation functions, Dyson-Schwinger equations,
Yang-Mills theory, 3PI effective action

###### pacs:

12.38.Aw, 14.70.Dj, 12.38.Lg

††journal: Eur. Phys. J. C

##  1 Introduction

Glueballs, i.e. hadrons that consist of gluons only, are extremely fascinating
objects to study. They arise due to the non-Abelian nature of Yang-Mills
theory which allows for the formation of gauge invariant states of gluons that
interact strongly amongst each other. The properties of glueballs have been
studied in many models since their prediction in the 1970s Fritzsch:1972jv ;
Fritzsch:1975tx . Today, many glueball masses in pure Yang-Mills theory are
known rather accurately owing to high statistics quenched lattice calculations
Bali:1993fb ; Morningstar:1999rf ; Chen:2005mg ; Athenodorou:2020ani .
Unquenched lattice calculations of glueball masses are still on the
exploratory level with considerable uncertainties due to severe problems with
the signal to noise ratio, see, e.g. Gregory:2012hu  and references therein.
Alternative theoretical frameworks, such as Hamiltonian many body methods
Szczepaniak:1995cw ; Szczepaniak:2003mr  or chiral Lagrangians Janowski:2011gt
; Eshraim:2012jv , have shed some light on potential mass patterns and
identifications of experimental states dominated by their glueball content,
see Ochs:2013gi  for a comprehensive review. However, it seems fair to state
that our detailed understanding of glueball formation from the underlying
dynamics of Yang-Mills theory is still far from complete. In this work, we
provide an additional, complementary perspective from functional methods.

While the calculation of mesons and baryons from functional bound state
equations is a very active field, see, e.g. Cloet:2013jya ; Eichmann:2016yit
and references therein, studies of exotic states in this framework are less
abundant. This is particularly true for glueballs due to the inherent
complexity of gauge-fixed Yang-Mills theories. Some results have been reported
in Meyers:2012ka ; Sanchis-Alepuz:2015hma ; Souza:2019ylx ; Kaptari:2020qlt ,
but they remain on an exploratory level due to the ansaetze used for the
input. An alternative approach to extract glueball masses from Landau gauge
correlation functions was followed in Dudal:2010cd ; Dudal:2013wja .

Figure 1:  The coupled set of BSEs for a glueball made from two gluons and a
pair of Faddeev-Popov (anti-)ghosts. Wiggly lines denote dressed gluon
propagators, dashed lines denote dressed ghost propagators. The gray boxes
represent interaction kernels given in Fig. 2. The Bethe-Salpeter amplitudes
of the glueball are denoted by gray disks.

Functional methods allow direct studies of the internal structures of bound
states as determined by the dynamics of QCD. As an example, let us mention
tetraquarks where the role of different two-body clusterings can be studied,
see Ref. Eichmann:2020oqt  for a summary. Another promising aspect of
functional methods is that they provide direct access to analytic properties.
The necessary numerical techniques are not fully developed yet, but they have
already been applied to some interesting questions like dynamical resonances
or spectral functions, see, for instance, Maris:1995ns ; Alkofer:2003jj ;
Windisch:2012sz ; Strauss:2012dg ; Windisch:2013dxa ; Pawlowski:2015mia ;
Strodthoff:2016pxx ; Eichmann:2019dts ; Fischer:2020xnb . As a third example
where functional methods can provide useful insights into QCD we mention the
study of its phases and the transitions between them, see, e.g.,
Pawlowski:2010ht ; Fischer:2018sdj . In the present work we add another
application of functional equations to the list. It is set apart from the
previous examples by delivering quantitative results without the need for any
model parameters which typically substitute for missing information in the
considered truncated systems of equations. This is to our knowledge the first
such calculation. If such setups can be devised and solved also for other
quantities of interest, this would constitute another major step forward.

What enabled us to perform a parameter-free calculation is the solid
understanding of the properties of Yang-Mills correlation functions obtained
in the last decade and a refinement of the methods to compute them from
functional equations, e.g. Huber:2012kd ; Blum:2014gna ; Eichmann:2014xya ;
Huber:2016tvc ; Cyrol:2016tym ; Huber:2017txg ; Aguilar:2018csq ;
Huber:2018ned ; Aguilar:2019jsj ; Aguilar:2019kxz ; Huber:2020keu . In
particular, we employ recent results from a fully self-contained truncation of
Dyson-Schwinger equations (DSEs) for propagators and vertices that leads to
good agreement with corresponding gauge fixed lattice results without any
tuning Huber:2020keu . As we will see in the following, the glueball spectrum
extracted from the corresponding set of bound state Bethe-Salpeter equations
(BSEs) agrees quantitatively with corresponding lattice results.

The remainder of the article is organized as follows. In Sec. 2 we introduce
the bound state equations we solve and the employed input. The extraction of
the spectrum is explained in Sec. 3 and the results are presented and
discussed in Sec. 4. We summarize and provide an outlook in Sec. 5. The
employed methods are illustrated in the appendix for a meson system.

##  2 Framework and input

Figure 2:  Interaction kernels from the three-loop 3PI effective action. All
propagators are dressed; black disks represent dressed vertices. In our
calculation, we include the diagrams inside the red rectangles. Details are
discussed in the text.

We consider Landau gauge Yang-Mills theory with gauge fixing via the Faddeev-
Popov procedure. The corresponding bound state equation for a two-gluon
glueball is depicted in Fig. 1. Since for some quantum numbers the two-gluon
state also couples to a ghost–anti-ghost state, we also need to consider the
corresponding BSE. For simplicity, we refer to them as glueball-part and
ghostball-part. The forms of the interaction kernels can be derived from an
effective action and its truncation. We use the kernels obtained from a 3PI
effective action truncated to three loops Berges:2004pu ; Carrington:2010qq .
This choice is motivated by excellent results that have been obtained in this
truncation scheme for the coupled set of DSEs for the dressed gluon
propagator, the dressed ghost propagator and the dressed ghost-gluon as well
as three-gluon vertices Huber:2020keu  discussed below. A truncation on a
similar level has also been used in the quark sector in Ref. Williams:2015cvx
. The corresponding kernels, shown in Fig. 2, are derived by performing
appropriate functional derivatives of the action Fukuda:1987su ; McKay:1989rk
; Sanchis-Alepuz:2015hma ; Sanchis-Alepuz:2015tha .

There are some noteworthy and important differences between a setup based on
the 1PI effective action and our 3PI setup. The 3PI effective action depends
on dressed propagators and three-point functions which can be calculated from
their equations of motion corresponding to stationarity conditions of the
action. Typically, one employs a (dressed) loop expansion for the 3PI
effective action with a finite number of terms Berges:2004pu ;
Carrington:2010qq . The truncated DSE system derived therefrom is then a
finite and closed system of equations for propagators and three-point
functions with all higher-order vertices bare. Thus, it contains non-
perturbative information for the three-point functions in a significantly
different fashion than a truncation derived from the 1PI effective action. In
the latter, all propagators and vertices are treated on the same footing, and
for any truncation of the 1PI effective action, there are infinitely many
equations to be solved. Thus, one normally performs the truncation at the
level of the system of equations (derived from the 1PI effective action)
instead of the action itself. If the 3PI effective action is truncated to
three loops, the equations of motions for the three-point functions have a
similar structure as the equations from the 1PI effective action without two-
loop terms Berges:2004pu . However, there are no bare three-point functions in
the former case. This seems like a small difference, but it turned out that
the equations of motion of the three-loop expanded 3PI effective action, which
are one-loop, outperform the results from equations of 1PI equations truncated
to one loop. An explicit comparison for the three-gluon vertex showed that
including the two-loop diagrams for the latter leads to good agreement.
However, such two-loop calculations are much harder to realize Huber:2020keu .

The similar structure of the equations of motion from 1PI and 3PI effective
actions persists also for the kernels of the bound state equations. A typical
kernel for a 1PI calculation contains only one-particle exchanges with one
bare vertex, but some studies beyond that exist as well, e.g., Bender:1996bb ;
Fischer:2009jm . The lowest order of 3PI kernels contain only one-particle
exchanges but in contrast to the 1PI case all vertices are dressed. The next
order also contains two-particle exchanges. For now, these are not taken into
account due to their technical complexity.

Of course, a priori there is no known way to guarantee that higher order
corrections are small in a dressed loop expansion. The results of this work,
in particular the excellent agreement of our glueball masses with
corresponding lattice results, provides some circumstantial evidence that this
may be the case. Further evidence has been collected in Eichmann:2014xya ;
Huber:2016tvc ; Huber:2017txg ; Huber:2018ned , where selected extensions of
the present truncation scheme have already been calculated and found to be
small indeed. Nevertheless, more systematic tests in dedicated calculations
need to be done in the future to further elaborate on this issue.

Having specified the underlying equations, we can now turn to solving them. A
solution to the coupled set of BSEs for glueballs can be found by treating
them as an eigenvalue equation for the matrix $\mathcal{K}$ of interaction
kernels with the eigenvector $\Gamma=(\Gamma_{\mu\nu},\Gamma_{gh})$ combining
the glueball-part $\Gamma_{\mu\nu}$ and the ghostball-part $\Gamma_{gh}$ of
the Bethe-Salpeter amplitude. This leads to

$\displaystyle\mathcal{K}\cdot\Gamma(P,p)=\lambda(P)\,\Gamma(P,p),$ (1)

where $P$ and $p$ are the total and relative momenta. For a solution, the
eigenvalue $\lambda(P)$ must equal $1$. The mass is determined from the
corresponding value of the total momentum, $M^{2}=-P^{2}$.

This structure is very general and valid for all allowed quantum numbers. Here
we treat the simplest glueballs, the scalar and pseudoscalar ones. To this
end, we need to specify the forms for the amplitudes. The glueball-part
$\Gamma_{\mu\nu}$ has two open Lorentz indices and is transverse. For the
scalar glueball ($J^{PC}=0^{++}$) one can find two independent tensors with
these properties given by Meyers:2012ka :

$\displaystyle\Gamma^{++}_{\mu\nu}(p,P)=h^{++}_{2}(p,P)\left(g_{\mu\nu}-\frac{p_{2\mu}p_{1\nu}}{p_{1}\cdot
p_{2}}\right)$
$\displaystyle\quad+h^{++}_{1}(p,P)\frac{(p_{1}^{2}\,p_{2\mu}-p_{1}\cdot
p_{2}\,p_{1\mu})(p_{2}^{2}\,p_{2\nu}-p_{1}\cdot p_{2}\,p_{1\nu})}{(p_{1}\cdot
p_{2})^{3}-p_{1}\cdot p_{2}\,p_{1}^{2}\,p_{2}^{2}}$ (2)

with $p_{1/2}=p\pm P/2$. For the pseudoscalar glueball ($J^{PC}=0^{-+}$), only
one tensor exists which we choose as

$\displaystyle\Gamma^{-+}_{\mu\nu}(p,P)=h^{-+}(p,P)\,\epsilon_{\mu\nu\rho\sigma}\hat{p}^{T}_{\rho}\hat{P}_{\sigma}.$
(3)

The hat indicates normalization and the superscript $T$ that the vector is
made transverse with respect to $P$. For the ghostball-part, which is a scalar
in Lorentz space, the amplitude is simply given by

$\displaystyle\Gamma^{++}_{gh}(p,P)=h^{++}_{3}(p,P).$ (4)

There is no corresponding amplitude with negative parity. This simplifies the
BSE for the pseudoscalar glueball where the ghostball-part of the amplitude
drops out.

The input required to solve the BSEs are the dressed gluon and ghost
propagators, $D_{\mu\nu}$ and $D_{G}$, respectively, given by

$D_{\mu\nu}(p)=\left(\delta_{\mu\nu}-\frac{p_{\mu}p_{\nu}}{p^{2}}\right)\frac{Z(p^{2})}{p^{2}},\,\,\,\,\,\,\,D_{G}(p)=\frac{-G(p^{2})}{p^{2}},$
(5)

as well as the dressed three-gluon and ghost-gluon vertices. For these we use
numerical results from a DSE system also derived from the three-loop 3PI
effective action. A graphical representation of this truncation together with
a thorough discussion of all technical details and merits can be found in
Huber:2020keu . Here, we only wish to mention that the scheme is self-
contained, i.e., it can be solved without any ad-hoc ansaetze and parameters.
Thus, either correlation functions are taken into account and solved for self-
consistently, or they are consistently neglected.

Figure 3: Gluon dressing function $Z(p^{2})$ (left) and gluon propagator
$D(p^{2})$ (right) in comparison to lattice data Sternbeck:2006rd . For the
sake of comparison, the functional results were renormalized to agree with the
lattice results at $6\,\text{GeV}$. Different lines correspond to different
decoupling/scaling solutions as explained in the text.

Figure 4: Left: Ghost dressing function $G(p^{2})$ in comparison to lattice
data Sternbeck:2006rd . For the sake of comparison, the functional results
were renormalized to agree with the lattice results at $8\,\text{GeV}$. Right:
Ghost-gluon vertex dressing function at the symmetric point in comparison to
$SU(2)$ lattice data Maas:2019ggf . Different lines correspond to different
decoupling/scaling solutions as explained in the text. Figure 5: Solutions for
the three-gluon vertex dressing function at the symmetric point in comparison
to lattice data Cucchieri:2008qm ; Sternbeck:2017ntv ). For the sake of
comparison, all data were renormalized to $1$ at $5\,\text{GeV}$. Different
lines correspond to different decoupling/scaling solutions as explained in the
text.

A truncation independent property of the Yang-Mills correlation functions in
the continuum is that they appear as a one-parameter family of so-called
decoupling/scaling solutions with the scaling limit as an end-point
Boucaud:2008ji ; Aguilar:2008xm ; Fischer:2008uz ; Alkofer:2008jy . For the
ghost and gluon propagators, these solutions agree at large momenta but start
to differ around the scale of $1\,\text{GeV}$ and below. The appearance of
this one-parameter family has been discussed to be caused by incomplete
nonperturbative gauge fixing in the Landau gauge/Faddeev-Popov setup due the
Gribov problem Fischer:2008uz ; Maas:2009se . Indeed, on the lattice it is
well studied that different samplings of the Gribov region influence the low
momentum behavior of correlation functions, see, e.g. Maas:2009se ;
Maas:2011se ; Sternbeck:2012mf  and references therein. If this picture is
correct, all solutions of the family should lead to the same results for
physical quantities. The glueballs studied in this work provide an ideal
testing ground for this hypothesis.

We thus used a number of different solutions within this one-parameter family
including the scaling solution taken from Ref. Huber:2020keu  for our
calculation of glueball masses. The corresponding gluon propagators and
dressing functions are shown in Fig. 3. Since the precise correspondence of
continuum and lattice gauge fixing methods is an open issue Maas:2019ggf , we
show a selection of functional results and lattice results for one particular
prescription to deal with Gribov copies. All employed gluon propagators
feature a maximum at $p^{2}>0$, but they are so shallow for two cases that
they are hardly visible in the plots. The ghost dressing function is shown in
Fig. 4. For the plots, the propagators were renormalized such that they agree
with the lattice data of Ref. Sternbeck:2006rd  at some high momentum. This
rescaling is unrelated to the renormalization employed in the actual
calculation, where different renormalization conditions were used. Especially
for the vertices the consistency of the renormalization procedure is an
important issue as discussed in Huber:2020keu . Here, we can directly use
these renormalized results.

The ghost-gluon vertex dressing function is shown in Fig. 4. The plot only
shows the symmetric point at which all momentum squares are equal, but its
full kinematic dependence was taken into account in the calculation. The
ghost-gluon vertex is finite in the Landau gauge Taylor:1971ff  and thus no
renormalization is necessary. The three-gluon vertex is shown in Fig. 5 where
it is compared to lattice data from Refs. Cucchieri:2008qm ; Sternbeck:2017ntv
; see Refs. Athenodorou:2016oyh ; Boucaud:2017obn  for similar results. Also
here the full kinematic dependence was taken into account, but its angular
dependence is very weak. For the sake of comparison, all data were
renormalized to $1$ at $5\,\text{GeV}$.

With the functional renormalization group, a truncation very similar to the
one used here leads to equally good results for the propagators and vertices
Cyrol:2016tym  which further supports the trustworthiness of the present
truncation. While improvements of this truncation were partially already
tested and found to be small, e.g., the impact of neglected contributions like
four-point functions Huber:2016tvc ; Huber:2017txg ; Huber:2018ned  and a full
basis for the three-gluon vertex Eichmann:2014xya , more tests should be
performed in the future.

For quark anti-quark bound states, the consistency of the truncations for the
DSE and BSE system is known to be crucial. This is true in particular for the
chiral properties, where the axial-vector Ward-Takahashi identity relates the
self-energies and integration kernels, see, e.g. Eichmann:2016yit . In the
present case of pure Yang-Mills glueballs, similar relations do not exist to
our knowledge, but it is certainly important to consider the impact of
consistency between the truncations of the DSEs and the BSEs.

Figure 6: Eigenvalues for off-shell values of $P^{2}=-M^{2}$ for scalar (left)
and pseudoscalar (right) glueballs. Shown are the eigenvalues for the ground
state and the first two excited states.

##  3 Extraction of the spectrum

Bethe-Salpeter equations describe bound states with time-like total momenta.
As a consequence, propagators and vertices inside Bethe-Salpeter equations are
tested at complex squared momenta. While calculations in the complex plane are
in principle possible with functional methods, see, e.g. Maris:1995ns ;
Fischer:2005en ; Fischer:2008sp ; Krassnigg:2009gd ; Strauss:2012dg ;
Kamikado:2013sia ; Pawlowski:2015mia ; Fischer:2020xnb , such calculations are
realized only with comparatively simple truncations. More advanced truncation
schemes, like the one we use in this work, have been solved for space-like
Euclidean momenta only. Analytic continuations of correlation functions have
been explored Dudal:2013yva ; Cyrol:2018xeq ; Binosi:2019ecz ; Dudal:2019gvn ,
but definite conclusions are still lacking.

Thus, we follow here an alternative approach and instead analytically continue
the eigenvalues calculated for Euclidean momenta $P^{2}>0$ into the time-like
momentum domain $P^{2}<0$. For every bound state we have calculated the
eigenvalue $\lambda(P^{2})$ of the BSE for $100$ space-like and real momenta
$P^{2}\in[10^{-4},0.25]\,\text{GeV}^{2}$ and extrapolate to time-like momenta
using Schlessinger’s method based on continued fractions Schlessinger:1968spm
; Tripolt:2018xeo . The extrapolation error is dominated by (small) numerical
inaccuracies in the input data and the limitations of the chosen extrapolation
function itself. We partially quantify the extrapolation error by a bootstrap-
like procedure: We take only a random subset of 80 points, which is enough to
get a good extrapolation function, and calculate the mass. Excluding
exceptional extrapolations that do not give a bound state, we repeat that
process 100 times and average the results. The one standard deviation errors
indicated in our results for the bound state masses stem from this procedure.
Of course, this does not include additional errors due to the truncation and
other sources.

We tested the extrapolation procedure in a case where we are able to compare
directly with a solution. To this end, we used a quark–anti-quark BSE in
rainbow-ladder approximation that is reviewed in Ref. Eichmann:2016yit . The
comparison between the extrapolated and the calculated eigenvalues is
described in A. Within errors both results agree very well.

One complication which needs to be taken into account is that not every
eigenvalue curve in the space-like domain can be extrapolated to a physical
state with a positive real mass. We encountered (and discarded) some complex
eigenvalue curves and one example for a (spurious) tachyonic state in case of
the pseudoscalar glueball. Furthermore, eigenvalue curves may cross. Hence,
the hierarchy of eigenvalues at space-like momenta does not need to correspond
to the hierarchy of the solutions.

The extrapolation is illustrated in Fig. 6 where the solutions for off-shell
values of $P^{2}=-M^{2}$ are shown. We chose deliberately to print the
eigenvalues over $M^{2}$ and not $M$, as the extrapolation is done in $P^{2}$.
For this plot, Eq. (1) was solved for various values of $\lambda\neq 1$. As
explained above, the errors are calculated by averaging over 100 random
subsets of the input points. The Schlessinger extrapolation sometimes contains
artifacts, e.g., very narrow poles. These can affect the average locally which
can lead to the small roughness observed in the errors. It might be possible
to refine the averaged extrapolations by applying further statistical methods
or other optimization techniques as, for example, in Ref. Binosi:2019ecz . For
now, the obtained precision is sufficient.

| Morningstar:1999rf  | Chen:2005mg  | Athenodorou:2020ani  | This work
---|---|---|---|---
State | $M\,[\text{MeV}]$ | $M/M_{0^{++}}$ | $M\,[\text{MeV}]$ | $M/M_{0^{++}}$ | $M\,[\text{MeV}]$ | $M/M_{0^{++}}$ | $M\,[\text{MeV}]$ | $M/M_{0^{++}}$
$0^{++}$ | $1760(50)$ | $1(0.04)$ | $1740(50)$ | $1(0.04)$ | $1651(23)$ | $1(0.02)$ | $1850(130)$ | $1(0.1)$
$0^{{}^{*}++}$ | $2720(180)$ | $1.54(0.11)$ | – | – | $2840(40)$ | $1.72(0.034)$ | $2570(210)$ | $1.39(0.15)$
$0^{{}^{**}++}$ | – | – | – | – | $3650(60)^{\dagger}$ | $2.21(0.05)^{\dagger}$ | $3720(160)$ | $2.01(0.16)$
$3580(150)^{\dagger}$ | $2.17(0.1)^{\dagger}$
$0^{-+}$ | $2640(40)$ | $1.5(0.05)$ | $2610(40)$ | $1.5(0.05)$ | $2600(40)$ | $1.574(0.032)$ | $2580(180)$ | $1.39(0.14)$
$0^{{}^{*}-+}$ | $3710(60)$ | $2.1(0.07)$ | – | – | $3540(80)$ | $2.14(0.06)$ | $3870(120)$ | $2.09(0.16)$
$0^{{}^{**}-+}$ | – | – | – | – | $4450(140)^{\dagger}$ | $2.7(0.09)^{\dagger}$ | $4340(200)$ | $2.34(0.19)$
$4540(120)^{\dagger}$ | $2.75(0.08)^{\dagger}$
Table 1: Ground and excited state masses $M$ of scalar and pseudoscalar
glueballs. Compared are lattice results from Morningstar:1999rf ; Chen:2005mg
; Athenodorou:2020ani  with the results of this work. For Morningstar:1999rf ;
Chen:2005mg , the errors are the combined errors from statistics and the use
of an anisotropic lattices. For Athenodorou:2020ani , the error is statistical
only. In our results, the error comes from the extrapolation method. All
results use the same value for $r_{0}=1/(418(5)\,\text{MeV})$, see text for
details. The related error is not included in the table. Masses with † are
conjectured to be the second excited states.

Figure 7:  Results for scalar and pseudoscalar glueball ground states and
excited states from lattice simulations Morningstar:1999rf ;
Athenodorou:2020ani  and this work. In the left diagram we display the
glueball masses on an absolute scale set by $r_{0}=1/418(5)\,\text{MeV})$. In
the right diagram we display the spectrum relative to the ground state.

Finally, the scalar glueball BSE leads to some additional insight with regard
to the ultraviolet behavior of the Bethe-Salpeter amplitudes. The integral in
the BSE is convergent only if the amplitude falls off polynomially at large
momenta. This is a standard behavior of all Bethe-Salpeter amplitudes in the
quark sector Eichmann:2016yit  and also expected for the glueball one.
However, we do find eigenvalue curves related to states with amplitudes that
run logarithmically and therefore do not satisfy this constraint111One
possible reason for the appearance of such solutions is the technical
implementation of the BSE problem. In practice, a cutoff for the radial
integral in the BSE and discrete quadrature rules are used. Thus, the
eigenproblem solved is that of a large but finite matrix which may contain
additional solutions not corresponding to a solution of the original BSE and
not constrained by the convergence properties stated above.. As a consequence,
these eigenvalue curves are cut-off dependent. In the gluon DSE, a similar
problem is associated with terms that break gauge covariance and lead to
quadratic divergences. A range of methods have been applied to remove these
terms, see e.g. Huber:2014tva  and references therein. In the BSE, the
simplest way to deal with this problem is to identify the solutions with the
wrong asymptotics as artifacts and discard them. We adopted this strategy in
this work for both scalar and pseudoscalar glueballs. It remains to be seen
whether the problem disappears entirely in even more advanced truncation
schemes.

##  4 Results and discussion

We show our results for the ground and excited states of scalar and
pseudoscalar glueball masses in Tab. 1 and Fig. 7 together with results from
lattice calculations. After the first version of this article appeared on
arxiv.org, new lattice data became available Athenodorou:2020ani  which we
also include now in the comparison. For the sake of comparability, we adjusted
all scales to the same value of the Sommer parameter $r_{0}$. Originally, we
set the scale of our results from the gluon propagator of Ref.
Sternbeck:2006rd  where $r_{0}=0.5\,\text{fm}$ was used, while for the lattice
results of Refs. Morningstar:1999rf ; Chen:2005mg  it was
$r_{0}=0.481(23)\,\text{fm}$. In Ref. Athenodorou:2020ani , $r_{0}$ was set to
$0.472(5)\,\text{fm}$ to express the obtained results in physical units. Here,
we adopt this value for the comparison of all results and rescaled the masses
accordingly. For our results this means that we multiply the masses by the
factor $1.059$. In addition, we show the results in terms of the corresponding
scalar ground state masses.

For the scalar glueball we find a ground state mass of
$1.850(130)\,\text{GeV}$. The extrapolation is very stable in this case. The
first excited state is at $2.570(210)\,\text{GeV}$, and we also find a
candidate for the second excited state at $3.720(160)\,\text{GeV}$. The
pseudoscalar glueball ground state is at $2.580(180)\,\text{GeV}$. The first
two excited states are at $3.870(120)\,\text{GeV}$ and
$4.340(200)\,\text{GeV}$. The masses for the scalar and pseudoscalar second
excited states were, to our knowledge, predicted here for the first time.

Ref. Athenodorou:2020ani  contains lists of states for the different
irreducible representations of the octahedral subgroup to which the rotational
symmetry is broken on the lattice. The spin of a glueball state is identified
by searching for nearly degenerate states in corresponding representations.
The scalar and pseudoscalar glueballs only appear in the representation
$A_{1}$. Four states with the corresponding quantum numbers $P=1$, $C=\pm 1$
are listed. The lower two are identified as the (pseudo)scalar glueball and
its first excitation. The other two are very close to each other. Since their
values are similar to the ones found here for the second excited states, we
conjecture that one of each pair indeed is such a state.

The comparison with the lattice results of Refs. Morningstar:1999rf ;
Chen:2005mg ; Athenodorou:2020ani  shows that we do not only recover the same
hierarchy but that our results also agree quantitatively very well with the
lattice results even including second excited states.

The amplitudes for the ground states and the two excited states at the lowest
calculated value for $P^{2}>0$ are shown in Fig. 8. Although they do not
correspond to the physical amplitudes, for which $P^{2}=-M^{2}$, they provide
a lot of insight. First of all, we see for the meson example described in
detail in A that the amplitudes show only a small dependence on the momentum
variable $P^{2}$. In that case, we can even use the results at space-like
momenta to extrapolate the results to time-like momenta. Second, off-shell
amplitudes are required for calculations where the corresponding bound states
are intermediate states over which one integrates.

Figure 8: Leading amplitudes for the scalar (top) and pseudoscalar (bottom)
ground states (left), first excited states (middle) and second excited states
(right) for the lowest calculated $P^{2}>0$.

All amplitudes in the ground states show a distinct peak around
$1\,\text{GeV}$. The first and second radial excitations contain one and two
nodes, respectively. For the scalar glueball we observe an interesting
interplay of the three different components given in Eqs. (2) and (4). Whereas
the ghostball-part plays an important role in the ground state, it is
negligible for the two excited states. Thus, neglecting the ghost contribution
would distort the results quantitatively (but not qualitatively) as we tested
explicitly. Furthermore, for ground and excited states different gluon
amplitudes are largest.

We explicitly verified that the masses extracted for the glueballs are
invariant (within error bars) of our choice of input within the one-parameter
family shown in Figs. 3-5. The endpoint of this family is the scaling solution
for which we obtain also the same results. This supports the view that the
family corresponds to different gauge choices within Landau gauge Maas:2009se
. At the same time, it constitutes a nontrivial test of the reliability of the
employed truncation, as any inconsistency could destroy the gauge independence
easily. We explicitly verified that the spectrum gets distorted if
inconsistent input (e.g. propagators from one member of the family and
vertices from another) is used. This sensitivity on the employed input
indicates that it would be challenging to replace part or all of the input by
models.

In this context, it is also interesting to mention that in the process of
devising the setup employed in this work, we also tried various models for
propagators and/or vertices for testing purposes. However, not only were we
never able to produce results somewhat close to the present ones, also
qualitative features of the eigenvalue spectrum were different. As an example,
let us mention Ref. Sanchis-Alepuz:2015hma , in which, using model vertices, a
reasonable result for the scalar glueball was obtained, but the pseudoscalar
one was above $4\,\text{GeV}$. We thus conclude that the quality of the input
is very important for a good overall picture.

##  5 Summary and outlook

We calculated the masses of the three lowest states for scalar and
pseudoscalar glueballs in pure Yang-Mills theory from propagators and vertices
in the Landau gauge obtained from Dyson-Schwinger equations. We gave special
emphasis to the consistency between the DSE and BSE setups. The truncation we
employed is completely self-contained and does not depend on any parameter or
ansatz. The scale is inherited from the comparison with lattice results for
the gluon propagator. The gauge dependent input for the propagators and
vertices translates into invariant results for the masses within a one-
parameter family of nonperturbative completions of Landau gauge. Our results
agree quantitatively very well with the lattice results of Refs.
Morningstar:1999rf ; Chen:2005mg  and we add two more states to the known
spectrum of pure Yang-Mills theory. It is encouraging that these new states
can be matched with subsequently published lattice results Athenodorou:2020ani
.

To our mind, the present results constitute a considerable step forward to
describe QCD bound states in a continuum approach from first principles. There
are many applications to which this setup can and should be extended. In
particular, the quark sector should be included to get access to real-world
glueballs. The calculation of higher spin glueballs or even other exotics like
hybrids may be possible within the same framework as well. At the same time,
these results motivate further studies of functional equations for complex
momenta to eliminate the need of the extrapolation of the eigenvalue curves
from space-like Euclidean data. A step in this direction was taken in Refs.
Strauss:2012dg ; Fischer:2020xnb .

## Acknowledgments

We thank Richard Williams for useful discussions. This work was supported by
the DFG (German Research Foundation) grant FI 970/11-1, by the BMBF under
contract No. 05P18RGFP1, and by the FWF (Austrian Science Fund) under Contract
No. P29216-N36. The HPC Cluster at the University of Graz was used for the
numerical computations.

##  Appendix A Comparison between calculation and extrapolation at time-like
momenta for a meson example

We illustrate the extrapolation procedure to time-like momenta via
Schlessinger’s method with the example of a meson. For the interaction we use
the Maris-Tandy model Maris:1999nt  that allows also a direct calculation so
that we can compare the numerically exact results with the extrapolation.
Since we are only interested in the reliability of the method, we do not
consider a specific physical example but simply use a system of a quark and an
anti-quark such that we obtain a meson mass of $M=2.62$ GeV (with $\eta=1.8$
and $\Lambda=0.72\,\text{GeV}$ in the Maris-Tandy model). This is in the same
range as the glueball masses. Our extrapolation procedure leads to $M=2.67(4)$
GeV.

Figure 9: Extrapolation exemplified with a meson of mass $M=2.62\,\text{GeV}$.
The red dots represent the exact solutions for the eigenvalues, the star the
physical one. The orange line is the averaged extrapolation with errors
indicated by the band. The horizontal error bars represent the errors for
specific mass values. Below $2\,\text{GeV}$, the agreement is so good that the
points lie on top of each other.

Figure 10: Top: Comparison of calculated (line) and extrapolated (dots)
leading meson amplitudes at selected values of $P^{2}$. Bottom: Leading meson
amplitude at time-like momenta $P^{2}$.

We illustrate the extrapolation in Fig. 9 which shows the eigenvalues
calculated for real masses (red points) and the extrapolation from the
calculation on the space-like side (orange line with error band). One has to
distinguish two ways of measuring errors. The orange band in Fig. 9 is
calculated from averaging over extrapolations using random subsets of input
points. The mass is determined from the position where the eigenvalue curve
crosses one. To calculate an error for the mass, we average over corresponding
solutions. Since the distribution of extrapolation curves is not necessarily
uniformly distributed, this error does not have to agree with the orange band.
In the plot, we illustrate this by solving for would-be masses, viz., we solve
for a range of $\lambda$’s. This gives us the blue dots with the horizontal
error bars which correspond to the errors given for the glueballs results. All
indicated errors correspond to one standard deviation.

In the main text, we showed results for the glueball amplitudes at an off-
shell value for $P^{2}$. For the meson example discussed here, we can
explicitly illustrate that the $P^{2}$-dependence of the amplitude is small.
This is shown in Fig. 10. In addition, we also extrapolated the amplitude from
the space-like to the time-like side and found only small quantitative
differences. Since this is sufficient for illustration purposes, we do not
average over several extrapolations as for the eigenvalue curve and simply use
all input points for the extrapolation.

## References

  * (1) H. Fritzsch and M. Gell-Mann, “Current algebra: Quarks and what else?”, eConf C720906V2 (1972) 135–165, arXiv:hep-ph/0208010 [hep-ph].
  * (2) H. Fritzsch and P. Minkowski, “Psi Resonances, Gluons and the Zweig Rule”, Nuovo Cim. A30 (1975) 393.
  * (3) UKQCD Collaboration, G. S. Bali, K. Schilling, A. Hulsebos, A. C. Irving, C. Michael, and P. W. Stephenson, “A Comprehensive lattice study of SU(3) glueballs”, Phys. Lett. B309 (1993) 378–384, arXiv:hep-lat/9304012 [hep-lat].
  * (4) C. J. Morningstar and M. J. Peardon, “The Glueball spectrum from an anisotropic lattice study”, Phys. Rev. D60 (1999) 034509, arXiv:hep-lat/9901004 [hep-lat].
  * (5) Y. Chen et al., “Glueball spectrum and matrix elements on anisotropic lattices”, Phys. Rev. D73 (2006) 014516, arXiv:hep-lat/0510074 [hep-lat].
  * (6) A. Athenodorou and M. Teper, “The glueball spectrum of SU(3) gauge theory in 3+1 dimension”, arXiv:2007.06422 [hep-lat].
  * (7) E. Gregory, A. Irving, B. Lucini, C. McNeile, A. Rago, C. Richards, and E. Rinaldi, “Towards the glueball spectrum from unquenched lattice QCD”, JHEP 10 (2012) 170, arXiv:1208.1858 [hep-lat].
  * (8) A. Szczepaniak, E. S. Swanson, C.-R. Ji, and S. R. Cotanch, “Glueball spectroscopy in a relativistic many body approach to hadron structure”, Phys. Rev. Lett. 76 (1996) 2011–2014, arXiv:hep-ph/9511422 [hep-ph].
  * (9) A. P. Szczepaniak and E. S. Swanson, “The Low lying glueball spectrum”, Phys. Lett. B577 (2003) 61–66, arXiv:hep-ph/0308268 [hep-ph].
  * (10) S. Janowski, D. Parganlija, F. Giacosa, and D. H. Rischke, “The Glueball in a Chiral Linear Sigma Model with Vector Mesons”, Phys. Rev. D84 (2011) 054007, arXiv:1103.3238 [hep-ph].
  * (11) W. I. Eshraim, S. Janowski, F. Giacosa, and D. H. Rischke, “Decay of the pseudoscalar glueball into scalar and pseudoscalar mesons”, Phys. Rev. D87 no. 5, (2013) 054036, arXiv:1208.6474 [hep-ph].
  * (12) W. Ochs, “The Status of Glueballs”, J. Phys. G40 (2013) 043001, arXiv:1301.5183 [hep-ph].
  * (13) I. C. Cloet and C. D. Roberts, “Explanation and Prediction of Observables using Continuum Strong QCD”, Prog. Part. Nucl. Phys. 77 (2014) 1–69, arXiv:1310.2651 [nucl-th].
  * (14) G. Eichmann, H. Sanchis-Alepuz, R. Williams, R. Alkofer, and C. S. Fischer, “Baryons as relativistic three-quark bound states”, Prog. Part. Nucl. Phys. 91 (2016) 1–100, arXiv:1606.09602 [hep-ph].
  * (15) J. Meyers and E. S. Swanson, “Spin Zero Glueballs in the Bethe-Salpeter Formalism”, Phys. Rev. D87 no. 3, (2013) 036009, arXiv:1211.4648 [hep-ph].
  * (16) H. Sanchis-Alepuz, C. S. Fischer, C. Kellermann, and L. von Smekal, “Glueballs from the Bethe-Salpeter equation”, Phys. Rev. D92 (2015) 034001, arXiv:1503.06051 [hep-ph].
  * (17) E. V. Souza, M. N. Ferreira, A. C. Aguilar, J. Papavassiliou, C. D. Roberts, and S.-S. Xu, “Pseudoscalar glueball mass: a window on three-gluon interactions”, Eur. Phys. J. A56 no. 1, (2020) 25, arXiv:1909.05875 [nucl-th].
  * (18) L. Kaptari and B. Kämpfer, “Mass spectrum of pseudo-scalar glueballs from a Bethe-Salpeter approach with the rainbow-ladder truncation”, Few Body Syst. 61 no. 3, (2020) 28, arXiv:2004.06523 [hep-ph].
  * (19) D. Dudal, M. S. Guimaraes, and S. P. Sorella, “Glueball masses from an infrared moment problem and nonperturbative Landau gauge”, Phys. Rev. Lett. 106 (2011) 062003, arXiv:1010.3638 [hep-th].
  * (20) D. Dudal, M. S. Guimaraes, and S. P. Sorella, “Pade approximation and glueball mass estimates in $3d$ and $4d$ with $N_{c}=2,3$ colors”, Phys. Lett. B732 (2014) 247–254, arXiv:1310.2016 [hep-ph].
  * (21) G. Eichmann, C. S. Fischer, W. Heupel, N. Santowsky, and P. C. Wallbott, “Four-quark states from functional methods”, arXiv:2008.10240 [hep-ph].
  * (22) P. Maris, “Confinement and complex singularities in QED in three-dimensions”, Phys. Rev. D52 (1995) 6087–6097, arXiv:hep-ph/9508323 [hep-ph].
  * (23) R. Alkofer, W. Detmold, C. S. Fischer, and P. Maris, “Analytic properties of the Landau gauge gluon and quark propagators”, Phys. Rev. D70 (2004) 014014, arXiv:hep-ph/0309077.
  * (24) A. Windisch, M. Q. Huber, and R. Alkofer, “On the analytic structure of scalar glueball operators at the Born level”, Phys. Rev. D87 no. 6, (2013) 065005, arXiv:1212.2175 [hep-ph].
  * (25) S. Strauss, C. S. Fischer, and C. Kellermann, “Analytic structure of the Landau gauge gluon propagator”, Phys.Rev.Lett. 109 (2012) 252001, arXiv:1208.6239 [hep-ph].
  * (26) A. Windisch, M. Q. Huber, and R. Alkofer, “How to determine the branch points of correlation functions in Euclidean space”, Acta Phys. Polon. Supp. 6 no. 3, (2013) 887–892, arXiv:1304.3642 [hep-ph].
  * (27) J. M. Pawlowski and N. Strodthoff, “Real time correlation functions and the functional renormalization group”, Phys. Rev. D92 no. 9, (2015) 094009, arXiv:1508.01160 [hep-ph].
  * (28) N. Strodthoff, “Self-consistent spectral functions in the O(N) model from the functional renormalization group”, Phys. Rev. D95 no. 7, (2017) 076002, arXiv:1611.05036 [hep-th].
  * (29) G. Eichmann, P. Duarte, M. Peña, and A. Stadler, “Scattering amplitudes and contour deformations”, Phys. Rev. D 100 no. 9, (2019) 094001, arXiv:1907.05402 [hep-ph].
  * (30) C. S. Fischer and M. Q. Huber, “Landau gauge Yang-Mills propagators in the complex momentum plane”, Phys. Rev. D 102 no. 9, (2020) 094005, arXiv:2007.11505 [hep-ph].
  * (31) J. M. Pawlowski, “The QCD phase diagram: Results and challenges”, AIP Conf.Proc. 1343 (2011) 75–80, arXiv:1012.5075 [hep-ph].
  * (32) C. S. Fischer, “QCD at finite temperature and chemical potential from Dyson–Schwinger equations”, Prog. Part. Nucl. Phys. 105 (2019) 1–60, arXiv:1810.12938 [hep-ph].
  * (33) M. Q. Huber and L. von Smekal, “On the influence of three-point functions on the propagators of Landau gauge Yang-Mills theory”, JHEP 1304 (2013) 149, arXiv:1211.6092 [hep-th].
  * (34) A. Blum, M. Q. Huber, M. Mitter, and L. von Smekal, “Gluonic three-point correlations in pure Landau gauge QCD”, Phys. Rev. D 89 (2014) 061703(R), arXiv:1401.0713 [hep-ph].
  * (35) G. Eichmann, R. Williams, R. Alkofer, and M. Vujinovic, “The three-gluon vertex in Landau gauge”, Phys.Rev. D89 (2014) 105014, arXiv:1402.1365 [hep-ph].
  * (36) M. Q. Huber, “Correlation functions of three-dimensional Yang-Mills theory from Dyson-Schwinger equations”, Phys. Rev. D93 no. 8, (2016) 085033, arXiv:1602.02038 [hep-th].
  * (37) A. K. Cyrol, L. Fister, M. Mitter, J. M. Pawlowski, and N. Strodthoff, “Landau gauge Yang-Mills correlation functions”, Phys. Rev. D94 no. 5, (2016) 054005, arXiv:1605.01856 [hep-ph].
  * (38) M. Q. Huber, “On non-primitively divergent vertices of Yang–Mills theory”, Eur. Phys. J. C77 no. 11, (2017) 733, arXiv:1709.05848 [hep-ph].
  * (39) A. C. Aguilar, M. N. Ferreira, C. T. Figueiredo, and J. Papavassiliou, “Nonperturbative structure of the ghost-gluon kernel”, Phys. Rev. D99 no. 3, (2019) 034026, arXiv:1811.08961 [hep-ph].
  * (40) M. Q. Huber, “Nonperturbative properties of Yang–Mills theories”, Phys. Rept. 879 (2020) 1–92, arXiv:1808.05227 [hep-ph].
  * (41) A. C. Aguilar, M. N. Ferreira, C. T. Figueiredo, and J. Papavassiliou, “Nonperturbative Ball-Chiu construction of the three-gluon vertex”, Phys. Rev. D99 no. 9, (2019) 094010, arXiv:1903.01184 [hep-ph].
  * (42) A. C. Aguilar, M. N. Ferreira, C. T. Figueiredo, and J. Papavassiliou, “Gluon mass scale through nonlinearities and vertex interplay”, Phys. Rev. D100 no. 9, (2019) 094039, arXiv:1909.09826 [hep-ph].
  * (43) M. Q. Huber, “Correlation functions of Landau gauge Yang-Mills theory”, Phys. Rev. D 101 no. 11, (2020) 11, arXiv:2003.13703 [hep-ph].
  * (44) J. Berges, “n-PI effective action techniques for gauge theories”, Phys. Rev. D70 (2004) 105010, arXiv:hep-ph/0401172.
  * (45) M. Carrington and Y. Guo, “Techniques for n-Particle Irreducible Effective Theories”, Phys.Rev. D83 (2011) 016006, arXiv:1010.2978 [hep-ph].
  * (46) R. Williams, C. S. Fischer, and W. Heupel, “Light mesons in QCD and unquenching effects from the 3PI effective action”, Phys. Rev. D93 no. 3, (2016) 034026, arXiv:1512.00455 [hep-ph].
  * (47) R. Fukuda, “Stability Conditions in Quantum System. A General Formalism”, Prog. Theor. Phys. 78 (1987) 1487–1507.
  * (48) D. W. McKay and H. J. Munczek, “Composite Operator Effective Action Considerations on Bound States and Corresponding S Matrix Elements”, Phys. Rev. D40 (1989) 4151.
  * (49) H. Sanchis-Alepuz and R. Williams, “Hadronic Observables from Dyson-Schwinger and Bethe-Salpeter equations”, J. Phys. Conf. Ser. 631 no. 1, (2015) 012064, arXiv:1503.05896 [hep-ph].
  * (50) A. Bender, C. D. Roberts, and L. Von Smekal, “Goldstone theorem and diquark confinement beyond rainbow ladder approximation”, Phys. Lett. B 380 (1996) 7–12, arXiv:nucl-th/9602012.
  * (51) C. S. Fischer and R. Williams, “Probing the gluon self-interaction in light mesons”, Phys. Rev. Lett. 103 (2009) 122001, arXiv:0905.2291 [hep-ph].
  * (52) A. Sternbeck, arXiv:hep-lat/0609016, PhD thesis, Humboldt-Universität zu Berlin, 2006\.
  * (53) A. Maas, “Constraining the gauge-fixed Lagrangian in minimal Landau gauge”, SciPost Phys. 8 no. 5, (2020) 071, arXiv:1907.10435 [hep-lat].
  * (54) A. Cucchieri, A. Maas, and T. Mendes, “Three-point vertices in Landau-gauge Yang-Mills theory”, Phys. Rev. D77 (2008) 094510, arXiv:0803.1798 [hep-lat].
  * (55) A. Sternbeck, P.-H. Balduf, A. Kızılersu, O. Oliveira, P. J. Silva, J.-I. Skullerud, and A. G. Williams, “Triple-gluon and quark-gluon vertex from lattice QCD in Landau gauge”, PoS LATTICE2016 (2017) 349, arXiv:1702.00612 [hep-lat].
  * (56) P. Boucaud et al., “IR finiteness of the ghost dressing function from numerical resolution of the ghost SD equation”, JHEP 06 (2008) 012, arXiv:0801.2721 [hep-ph].
  * (57) A. Aguilar, D. Binosi, and J. Papavassiliou, “Gluon and ghost propagators in the Landau gauge: Deriving lattice results from Schwinger-Dyson equations”, Phys.Rev. D78 (2008) 025010, arXiv:0802.1870 [hep-ph].
  * (58) C. S. Fischer, A. Maas, and J. M. Pawlowski, “On the infrared behavior of Landau gauge Yang-Mills theory”, Annals Phys. 324 (2009) 2408–2437, arXiv:0810.1987 [hep-ph].
  * (59) R. Alkofer, M. Q. Huber, and K. Schwenzer, “Infrared singularities in landau gauge yang-mills theory”, Phys. Rev. D81 (2010) 105010, arXiv:0801.2762 [hep-th].
  * (60) A. Maas, “Constructing non-perturbative gauges using correlation functions”, Phys. Lett. B689 (2010) 107–111, arXiv:0907.5185 [hep-lat].
  * (61) A. Maas, “Describing gauge bosons at zero and finite temperature”, Phys.Rept. 524 (2013) 203–300, arXiv:1106.3942 [hep-ph].
  * (62) A. Sternbeck and M. Müller-Preussker, “Lattice evidence for the family of decoupling solutions of Landau gauge Yang-Mills theory”, Phys.Lett. B726 (2013) 396–403, arXiv:1211.3057 [hep-lat].
  * (63) J. C. Taylor, “Ward identities and charge renormalization of the yang- mills field”, Nucl. Phys. B33 (1971) 436–444.
  * (64) A. Athenodorou, D. Binosi, P. Boucaud, F. De Soto, J. Papavassiliou, J. Rodriguez-Quintero, and S. Zafeiropoulos, “On the zero crossing of the three-gluon vertex”, Phys. Lett. B761 (2016) 444–449, arXiv:1607.01278 [hep-ph].
  * (65) P. Boucaud, F. De Soto, J. Rodríguez-Quintero, and S. Zafeiropoulos, “Refining the detection of the zero crossing for the three-gluon vertex in symmetric and asymmetric momentum subtraction schemes”, Phys. Rev. D95 no. 11, (2017) 114503, arXiv:1701.07390 [hep-lat].
  * (66) C. S. Fischer, P. Watson, and W. Cassing, “Probing unquenching effects in the gluon polarisation in light mesons”, Phys. Rev. D72 (2005) 094025, arXiv:hep-ph/0509213 [hep-ph].
  * (67) C. S. Fischer, D. Nickel, and R. Williams, “On Gribov’s supercriticality picture of quark confinement”, Eur. Phys. J. C60 (2009) 47–61, arXiv:0807.3486 [hep-ph].
  * (68) A. Krassnigg, “Excited mesons in a Bethe-Salpeter approach”, PoS CONFINEMENT8 (2008) 075, arXiv:0812.3073 [nucl-th].
  * (69) K. Kamikado, N. Strodthoff, L. von Smekal, and J. Wambach, “Real-time correlation functions in the $O(N)$ model from the functional renormalization group”, Eur. Phys. J. C74 no. 3, (2014) 2806, arXiv:1302.6199 [hep-ph].
  * (70) D. Dudal, O. Oliveira, and P. J. Silva, “Källén-Lehmann spectroscopy for (un)physical degrees of freedom”, Phys. Rev. D89 no. 1, (2014) 014010, arXiv:1310.4069 [hep-lat].
  * (71) A. K. Cyrol, J. M. Pawlowski, A. Rothkopf, and N. Wink, “Reconstructing the gluon”, SciPost Phys. 5 no. 6, (2018) 065, arXiv:1804.00945 [hep-ph].
  * (72) D. Binosi and R.-A. Tripolt, “Spectral functions of confined particles”, Phys. Lett. B 801 (2020) 135171, arXiv:1904.08172 [hep-ph].
  * (73) D. Dudal, O. Oliveira, M. Roelfs, and P. Silva, “Spectral representation of lattice gluon and ghost propagators at zero temperature”, Nucl. Phys. B952 (2020) 114912, arXiv:1901.05348 [hep-lat].
  * (74) L. Schlessinger, “Use of Analyticity in the Calculation of Nonrelativistic Scattering Amplitudes”, Phys. Rev. 167 no. 3, (1968) 1411\.
  * (75) R.-A. Tripolt, P. Gubler, M. Ulybyshev, and L. Von Smekal, “Numerical analytic continuation of Euclidean data”, Comput. Phys. Commun. 237 (2019) 129–142, arXiv:1801.10348 [hep-ph].
  * (76) M. Q. Huber and L. von Smekal, “Spurious divergences in Dyson-Schwinger equations”, JHEP 1406 (2014) 015, arXiv:1404.3642 [hep-ph].
  * (77) P. Maris and P. C. Tandy, “Bethe-Salpeter study of vector meson masses and decay constants”, Phys. Rev. C60 (1999) 055214, arXiv:nucl-th/9905056 [nucl-th].