Reference : A. Bazavov, N. Brambilla, P. Petreczky, A. Vairo, J. H. Weber (TUMQCD Collaboration),
Color screening in (2+1)-flavor QCD, Phys. Rev. D 98, 054511 (2018), arXiv:1804.10600.

Companion page (equation by equation): color screening in 2+1 flavor QCD at finite temperature
| Research Papers | Debye mass | EQCD | quark-gluon plasma

This page answers a different question from the companion page. Not what each equation says, but
*why the authors performed each step of the analysis, and why the paper needs fifteen main figures
plus a small forest of appendix figures.* Read it as: every plot in this paper exists to close one
specific loophole.

• #1. The paper in one paragraph
• #2. The physics question, and why one plot cannot answer it
• #3. The chain of reasoning (map of the paper)
• #4. What is measured, and why three different operators
• #5. Why subtracted free energies, and why 2F_Q
• #6. Why so much machinery before a single physics plot
• #7. Figure by figure
• #8. The appendix figures, grouped by the doubt they kill
• #9. Why so many plots -- the short answer
• #10. Results worth remembering

Two infinitely heavy (static) colour charges are placed a distance \(r\) apart inside a thermal bath of
gluons and \(2+1\) dynamical quarks. On the lattice one measures the correlator of their worldlines and
converts it into a free energy. The paper asks two questions: (i) at what distance does the medium
begin to screen the colour force, and (ii) once it does, can weak-coupling effective field theories
(pNRQCD, EQCD) describe what is seen? The answers, obtained over the very wide range
\(140 \le T \le 5814\) MeV with the HISQ action and six values of \(N_\tau\), are: screening switches on at
\(rT \approx 0.3\), and in the window \(0.3 \lesssim rT \lesssim 0.6\) weak coupling works.

Everything else in the 45 pages is the price of making those two sentences defensible.

At high temperature QCD has three scales, not one:

scale	name	physics living there	can perturbation theory reach it?
\(1/r\) (short distances)	hard/vacuum	ordinary Coulomb attraction, asymptotic freedom	yes
\(\pi T\)	hard thermal	lowest Matsubara mode	yes
\(gT\)	soft / electric	Debye screening, mass \(m_D\)	yes, after resummation (EQCD)
\(g^2T\)	ultrasoft / magnetic	confining magnetic sector	no (Linde problem)

Each scale owns a different effective theory, and each effective theory is valid only in its own
window of \(r\). Therefore:

• A single plot of "free energy versus \(r\)" cannot decide anything, because at different \(r\) one is in different theories.
• The natural variable is not \(r\) but the dimensionless \(rT\): it says which scale you are probing. This is why almost every physics figure in the paper is plotted against \(rT\) and not against \(r\) in fm.
• "Where does screening start?" is really the question *where does the \(T=0\) static energy \(V_S(r)\) stop describing the finite-\(T\) free energy?* So the \(T=0\) curve must appear in the same plot -- and it does, as the black band in Figs. 1, 2 and 4.

A second complication: what one can measure on a lattice is not unique. The colour-averaged free
energy \(F_{Q\bar Q}\), the singlet free energy \(F_S\) and the cyclic Wilson loop \(F_W\) are three
different objects that agree only in limits. Each of them therefore needs its own comparison figure.

The whole paper is one pipeline. Each arrow below is a place where a plot is required, either to
justify the step or to show its result.

1. Measure bare correlators on many ensembles. (Sec. II)
2. Renormalize them using the Polyakov loop. (Sec. II, Eqs. 10-14)
3. Remove short-distance lattice artifacts using \(T=0\) data. (App. B)
4. Interpolate in \(r\) and \(T\), then extrapolate \(a\to0\). (App. C)
5. Check finite volume and sea-quark-mass dependence. (App. A, App. D)
6. Result plots: \(F_S\), \(F_{Q\bar Q}\) vs \(V_S\); effective coupling. (Figs. 1-4)
7. Cross-check gauge dependence with cyclic Wilson loops. (Fig. 5)
8. Compare small \(r\) with pNRQCD. (Figs. 6, 7, 9, 10, 11)
9. Compare screening regime with EQCD. (Figs. 8, 12)
10. Fit the asymptotic exponential decay, extract screening masses. (Figs. 13-15)

Steps 3-5 produce no physics on their own. They exist because a systematic error at the \(5\)--\(10\%\)
level would be exactly the size of the effect being claimed, and would therefore be
indistinguishable from screening. That is the single most important reason the paper has so many plots.

\[ C_P(r,T) = \left\langle \mathrm{Tr}\,P(\mathbf{0})\;\mathrm{Tr}\,P^\dagger(\mathbf{r}) \right\rangle = \exp\left(-\frac{F_{Q\bar Q}(r,T)}{T}\right) \]

Why this one. It is gauge invariant and it is the object that appears in the partition function with
one static quark and one static antiquark inserted, so its free energy is unambiguous. Its drawback is
that it is a colour average: it contains the attractive singlet channel and the repulsive octet
channel together, weighted \(1/9\) and \(8/9\). Hence the normalization offset \(T\ln 9\) that appears in
the right panel of Fig. 1.

The colour average is also why \(F_{Q\bar Q}^{\rm sub}\) falls like \(1/r^2\) and not like \(1/r\) at short
distance: single-gluon exchange cancels between the two channels, and the leading survivor is
two-gluon exchange, \(\sim \alpha_s^2/(r^2T)\). That single fact dictates the axis choice of Figs. 3, 10,
11 and 12, where the plotted quantity is \(-r^2T F^{\rm sub}\) or \(-9r^2T F^{\rm sub}\): in that
combination the leading order is a pure number, so any \(r\)-dependence left in the plot is either the
running of \(\alpha_s\) or genuine screening. A flat curve means "still \(1/r^2\)"; a curve bending upward
means the singlet Coulomb \(1/r\) has taken over; a curve bending down and dying means screening.

Why this one too. Quarkonium in a medium is a singlet object; the potential that enters a
Schroedinger equation is the singlet one, not the colour average. \(F_S\) is also the quantity that has a
clean \(T\to0\) limit, namely the static energy \(V_S(r)\), which is known very precisely from \(T=0\)
lattice work. So \(F_S\) is the quantity in which "medium effect" can be defined as
\(V_S(r) - F_S(r,T)\) -- the object of Figs. 6 and 7.

The price. \(C_S\) is defined in Coulomb gauge, so it is not gauge invariant beyond leading order.
This is a genuine loophole, and it is closed by the third operator.

Why a third one. It is gauge invariant and reduces to the same singlet quantity, so if the
Coulomb-gauge \(F_S\) were an artifact of the gauge choice, the two would disagree. Fig. 5 is the plot
that tests this, and it is also why HYP smearing levels are scanned there: the Wilson loop needs
smearing to have any signal at all, so one must show that the agreement survives the smearing.

A bare Polyakov loop contains a self-energy divergence \(\propto 1/a\) for each static quark. Two
strategies exist; the paper uses the cheap and exact one: divide the correlator by the square of the
bare Polyakov loop,

\[ C^{\rm sub}(r,T) = \frac{C^{\rm bare}(r,T)}{\left(L^{\rm bare}(T)\right)^2},\qquad F^{\rm sub} = -T\ln C^{\rm sub}, \]

so that the divergence cancels identically, configuration by configuration. All renormalization is
then pushed into one single number per temperature, \(F_Q(T)\), defined by \(L = e^{-F_Q/T}\), which is
fixed by matching to the \(T=0\) static energy in the \(r_1\) scheme. The physical free energies are then

\[ F_S = F_S^{\rm sub} + 2F_Q, \qquad F_{Q\bar Q} = F_{Q\bar Q}^{\rm sub} + 2F_Q . \]

Consequences you can see in the figures.

• \(F^{\rm sub}\to0\) as \(r\to\infty\) by construction, so \(F\to 2F_Q\): that is the horizontal band in Fig. 1.
• Because \(F^{\rm sub}\) is a ratio, it is the statistically cleanest quantity, which is why every screening study (Figs. 3, 8, 10-15) is done with \(F^{\rm sub}\) and not with \(F\).
• Because the renormalization is one constant, a mistake in it shifts curves rigidly but cannot create or destroy screening. This separation of "additive constant" from "shape" is a deliberate design choice.

The claim of the paper is that medium effects on \(F_S\) are invisible below \(rT\approx0.3\) and visible
above it. At \(T\sim1\) GeV, \(rT=0.3\) means \(r\approx0.06\) fm -- a handful of lattice spacings. At such
short distances the dominant error is not statistics, it is discretization. So four separate
systematic studies must be done, and each one generates its own set of plots:

systematic	why it could fake the signal	where it is dealt with
cutoff effects at short \(r\)	a lattice Coulomb potential deviates from \(1/r\) by tens of percent at \(r\sim a\); that deviation is \(r\)-dependent and would look exactly like a medium effect	App. B, non-perturbative correction factors \(\langle K\rangle\) built from \(T=0\) data
continuum limit	different \(N_\tau\) at the same \(T\) means different \(a\); without \(a\to0\) one cannot separate \(a\)-dependence from \(T\)-dependence	App. C, fits in \(1/N_\tau^2,\,1/N_\tau^4\) (HISQ has only even powers)
finite volume	screening at large \(r\) is precisely where the box boundary is felt; a too-small box would mimic faster screening	aspect ratios \(N_\sigma/N_\tau = 4\) and \(6\) compared, App. C and D
sea quark mass	ensembles with \(m_l=m_s/20\) and \(m_s/5\) exist; if the free energy depended on it, high-\(T\) conclusions would be unsafe	App. A, Fig. 16

Plus one purely technical obstacle: different \(N_\tau\) give data on different grids of \((r,T)\), and
nothing can be extrapolated point by point unless the points line up. Hence the interpolation
machinery of App. C (local versus global fits), and hence the plots comparing the two interpolation
strategies -- if local and global fits disagreed, the continuum limit would be model-dependent.

For each figure: plotted = what is on the axes, why = the question it was made to answer,
result = what one is supposed to read off, feeds = what later step depends on it.

Plotted. Left: continuum-extrapolated singlet free energy \(F_S(r,T)\) in MeV, one curve per
temperature. Right: the colour-averaged \(F_{Q\bar Q}(r,T)\), with the \(T\ln9\) convention subtracted.
Both panels carry the \(T=0\) static energy \(V_S(r)\) as a black band, and a horizontal band at \(2F_Q(T)\).

Why. This is the defining plot of the paper. "Screening" is defined operationally as the distance
at which the finite-\(T\) curve peels away from the \(T=0\) curve, so the two must be drawn on the same
axes with the same normalization. That is also why so much work went into fixing \(F_Q\): without a
correct additive constant the curves would peel away trivially.

Result. At short \(r\) all temperatures lie on \(V_S(r)\): the medium is blind to a small colour-neutral
pair. Each curve then bends over and flattens at \(2F_Q(T)\), the free energy of two independent screened
quarks. The bend happens at smaller \(r\) for higher \(T\). \(F_S\) follows \(V_S\) noticeably longer than
\(F_{Q\bar Q}\) does, because \(F_{Q\bar Q}\) contains the octet channel, which is sensitive to the medium
already at short distance.

Feeds. Everything. Figs. 2 and 4 differentiate these curves; Figs. 6-7 subtract them.

Plotted. \(\alpha_{Q\bar Q}(r,T) = \frac{1}{C_F}r^2\,\partial_r F_S(r,T)\) versus \(r\), one curve per
temperature, against the vacuum \(\alpha_{Q\bar Q}(r)\) from \(V_S\) (black band).

Why. Two reasons. First, a derivative removes the additive constant \(2F_Q\) entirely, so this plot is
free of the renormalization scheme -- it is an independent way of asking "when does the finite-\(T\)
curve leave the vacuum curve?", one that cannot be faked by a wrong constant. Second, the value of
the coupling at its maximum answers a physically different question: is weak coupling legitimate at all?

Result. Each finite-\(T\) curve follows the vacuum band, peaks near \(r_{\max}\approx0.4/T\), and falls.
The peak position is the screening scale seen from another angle; the peak value is below \(0.5\) for
\(T\gtrsim320\) MeV, which is the paper's quantitative licence to use weak-coupling EFTs above roughly
\(300\) MeV -- and the basis for the statement that weak-coupling pNRQCD may be applied to quarkonium
there.

Plotted. The subtracted colour-averaged free energy times \(-r^2T\), on a log-log scale, at small \(r\).

Why. To test the functional form rather than the size. As explained in Sec. 4 above, the leading
weak-coupling prediction for the colour-averaged correlator at \(T\gg\alpha_s/r\) is \(\propto1/r^2\),
because single-gluon exchange cancels in the colour average. The chosen combination turns that
prediction into a horizontal line, so the eye can distinguish two regimes at a glance:
flat \(\Rightarrow\) the \(1/r^2\) two-gluon regime; rising linearly \(\Rightarrow\) the singlet-dominated
\(1/r\) Coulomb regime.

Result. Both regimes are visible, and the changeover happens where expected. This validates the
short-distance picture before any EFT curve is drawn on top of the data.

Plotted. Same construction as Fig. 2, but built from the colour-averaged free energy, again against
the vacuum band.

Why. Fig. 2 alone could be accused of being a Coulomb-gauge artifact, since \(F_S\) is gauge dependent.
Repeating the analysis with a manifestly gauge-invariant quantity checks that the qualitative picture
(follow the vacuum, peak, fall) is not an artifact.

Result. The colour-averaged coupling departs from the vacuum band earlier and peaks lower, exactly as
the octet contamination predicts. Same physics, quantitatively shifted -- which is itself informative:
it tells you how much of "the medium effect" in \(F_{Q\bar Q}\) is really just colour averaging.

Plotted. \(\ln(-rF_W^{\rm sub})\) from smeared and unsmeared cyclic Wilson loops, with the
Coulomb-gauge singlet result overlaid as bursts; four panels: \(T=200\) and \(480\) MeV, \(N_\tau=8\) and \(12\).

Why. This is the gauge-dependence audit announced in Sec. 4. The cyclic Wilson loop is gauge
invariant but noisy and needs HYP smearing; smearing itself distorts short distances. So the figure has
to show both that smeared Wilson loops converge to the Coulomb-gauge singlet, and that the
convergence is not an accident of one smearing level, one temperature or one lattice spacing. That is
why it is four panels and several smearing levels rather than one curve.

Result. With enough smearing the two agree. The singlet free energy used everywhere else in the paper
is therefore not a gauge artifact.

Plotted. The difference between the \(T=0\) static energy and the finite-\(T\) singlet free energy, in
MeV (left) and in units of \(T\) (right), versus \(r\).

Why. Fig. 1 shows two curves lying on top of each other at small \(r\); a difference plot is the only
honest way to see how well they really coincide, because a \(10\) MeV effect is invisible on a plot whose
range is hundreds of MeV. The two panels exist because the comparison to be made next (pNRQCD)
naturally predicts the difference in units of \(T\), while the experimental intuition ("how many MeV?")
lives in the left panel.

Result. The difference is tiny and nearly \(r\)-independent below \(rT\approx0.3\), then grows quickly.
This is the sharpest form of the paper's central claim.

Plotted. The same difference, now with the weak-coupling pNRQCD prediction overlaid; dotted line for
renormalization scale \(\mu=2\pi T\), band from \(\mu=\pi T\) to \(4\pi T\). The perturbative curve is shifted
by a small constant.

Why. Having established the size of the medium effect, one asks whether perturbation theory
reproduces it. Three details of the figure are themselves arguments:

• The band is not an error bar, it is a scale variation: in a truncated perturbative series the choice of \(\mu\) is arbitrary, and the spread it produces is the honest estimate of the missing orders. A narrow band at high \(T\) is evidence that the expansion is converging.
• The constant shift is not fudging: the pNRQCD expression is missing an \(r\)-independent \(\mathcal{O}(g^6)T\) matching term, so only the \(r\)-dependence is predicted. Shifting makes the comparison test what is actually predicted.
• The prediction contains a positive \((rT)^2\) piece from non-static gluons and quarks, and a negative \((rT)^2\) piece of order \(g^5\) from the Debye scale. Their partial cancellation is the reason the medium effect stays small until \(rT\approx0.3\) -- the figure is the evidence for that explanation.

Result. Agreement in shape at small \(rT\); failure once \(rT\) grows, which is the signal that the
hierarchy \(1/r\gg T\gg m_D\) has broken and one must resum, i.e. move to EQCD -- the next figure.

Plotted. \(-rF_S^{\rm sub}\) versus \(rT\) at several temperatures (curves offset vertically for
visibility), with hatched LO and solid NLO EQCD bands, scale varied \(\mu=\pi T,\,2\pi T,\,4\pi T\).

Why. This is the screening-regime counterpart of Fig. 7. The combination \(-rF_S^{\rm sub}\) is chosen
because a Yukawa-screened potential \(e^{-m_Dr}/r\) becomes a pure exponential once multiplied by \(r\), so
screening shows up as a straight line on a log axis rather than as a curved shape. Showing LO and NLO
side by side is the point of the plot: if NLO were not visibly closer to the data than LO, the EFT would
not be under control.

Result. NLO EQCD describes \(F_S^{\rm sub}\) up to \(rT\approx0.6\) for \(T\gtrsim600\) MeV, and the scale
band shrinks from LO to NLO. Beyond \(rT\approx0.6\) the description fails -- the magnetic scale.

Plotted. The octet piece of the Polyakov loop correlator (left) and the octet free energy (right)
versus \(r\), several temperatures, \(N_\tau=8\).

Why. Before one can predict the colour-averaged free energy from pNRQCD (Fig. 10), one needs the
octet channel, which is the part nobody measures directly. This figure isolates and displays it, so
that the reconstruction in the next figure is not a black box.

Result. Quantifies how large the repulsive octet contribution is and where it matters -- explaining,
in passing, why \(F_{Q\bar Q}\) departed from \(V_S\) earlier than \(F_S\) in Fig. 1.

Plotted. \(-r^2T F_{Q\bar Q}^{\rm sub}\) from \(N_\tau=12\) lattices (squares) against the pNRQCD
reconstruction (bands), at \(T=172\) MeV and \(T=666\) MeV. In the right panel three bands: full
reconstruction, one ignoring Casimir-scaling violation in the octet potential, one ignoring the octet
contribution altogether.

Why. This is the most stringent test in the paper, because the reconstruction formula uses only
lattice inputs -- the \(T=0\) static energy, the renormalized Polyakov loop, the measured octet shift --
plus one perturbative ingredient. So it predicts a measured quantity almost parameter-free. The three
bands are a deliberate ablation study: they show which physical ingredient each part of the agreement
comes from, rather than claiming agreement from a single tuned curve.

Result. The reconstruction works up to \(rT\approx0.3\), and both the octet contribution and its
Casimir-scaling violation are needed. Two very different temperatures are shown to prove that the
success is not accidental at one \(T\).

Plotted. Continuum data at a single high temperature against the weak-coupling expression to
next-to-next-to-next-to-leading order, three lines for \(\mu=\pi T,2\pi T,4\pi T\).

Why. One very high temperature is singled out because that is where the weak-coupling expansion has
the best chance, and where the deepest available order can be tested cleanly. The factor \(9=N_c^2\) in
the plotted combination normalizes away the leading colour factor so that the curve is directly a
statement about \(\alpha_s\).

Result. Agreement at short \(rT\), with the ordering of the three \(\mu\)-lines showing the residual scale
sensitivity -- the honest measure of how far the expansion can be trusted.

Plotted. Lattice data (including \(N_\tau=4\), aspect ratio 6, corrected for cutoff effects, as open
symbols) against LO and NLO EQCD bands with \(\mu\) varied.

Why. Fig. 8 tested EQCD on the singlet; this tests it on the gauge-invariant colour-averaged
quantity, where the LO prediction is a double-exponential \(e^{-2m_Dr}\) because the colour-averaged
correlator requires two screened gluon exchanges. Coarse \(N_\tau=4\) lattices with a large box are used
deliberately: the screening regime needs large physical volume, which is cheapest on coarse lattices,
so the price is cutoff effects, which is exactly why they are corrected and shown as open symbols.

Result. NLO EQCD reproduces \(F_{Q\bar Q}^{\rm sub}\) within errors for \(T\gtrsim1\) GeV. Together with
Fig. 8, this fixes the upper edge of the weak-coupling window at \(rT\approx0.6\).

Plotted. \(-rF_{Q\bar Q}^{\rm sub}\) (upper) and \(-rF_S^{\rm sub}\) (lower) versus \(r\) on a log scale at
\(N_\tau=4\), several temperatures.

Why. This is the figure that justifies the fit ansatz used in Figs. 14-15. Before quoting a screening
mass one must demonstrate that the data really do decay exponentially over a usable range -- otherwise
"the screening mass" is just a fit parameter of the wrong model. The straight portion of each curve is
the fit window.

Result. Clean exponential fall-off, with a slope that grows with \(T\); the singlet and colour-averaged
slopes differ by roughly a factor two, as the one- versus two-gluon-exchange picture requires.

Plotted. The asymptotic screening mass from exponential fits to \(-rF_S^{\rm sub}\), in units of \(T\),
versus \(T\); lines are the NLO Debye mass for three scales, multiplied by a constant. Open squares are
\(N_\tau=4\), aspect ratio 6.

Why. The purpose is to ask whether the measured screening is perturbative. The honest way to ask is
to check whether \(m_S/T\) has the same shape in \(T\) as \(m_D/T\) -- because at asymptotically high \(T\)
both must run like \(g(T)\). The multiplicative constant on the perturbative lines is the punchline, not a
fudge: a constant is needed, and its size measures the non-perturbative (magnetic) enhancement.

Result. \(m_S \approx (1.6\)--\(2.0)\,m_D\). Screening is qualitatively perturbative and quantitatively not.

Plotted. The same for \(-rF_{Q\bar Q}^{\rm sub}\), compared with \(2\times\) the NLO Debye mass (again
times a constant), plus a horizontal band from an independent 3D-EQCD lattice determination.

Why. Two independent cross-checks in one plot. First, internal consistency: the colour-averaged
channel must screen with roughly twice the singlet mass, because it needs two gluon exchanges -- so
\(m_{Q\bar Q}\approx2m_S\) is a prediction of the picture, not an input. Second, external consistency: the
3D-EQCD band comes from a completely different lattice calculation, so agreement with it validates the
whole chain from a direction the paper does not control.

Result. \(m_{Q\bar Q}\) sits about \(25\%\) above \(2m_S\) at high \(T\) and is consistent with the 3D-EQCD
band. Both checks pass.

Only Fig. 16 could be identified by number from the text available to me; the remaining appendix
figures are listed by content. If you have the PDF open, renumber them.

• Fig. 16: the ratio of bare single-quark free energies, and the screening functions for \(F_S^{\rm sub}\) and \(F_{Q\bar Q}^{\rm sub}\), for two light sea-quark masses (\(m_s/20\) vs \(m_s/5\)). Why: if the free energy depended on the light quark mass, no high-\(T\) conclusion would survive. Result: differences are consistent with statistics.
• Continuum-extrapolated static-quark entropy \(S_Q\) at high \(T\). Why: \(S_Q=-\partial F_Q/\partial T\) is a sharp, scheme-sensitive test of the Polyakov loop renormalization that fixed the additive constant in Fig. 1.

This is the largest block of plots, and the reason is arithmetic: the physics claim lives at
\(rT\approx0.3\), which at high \(T\) is only a few lattice spacings.

• Tree-level lattice static energy for Luescher-Weisz versus Wilson gauge action, normalized to the continuum. Why: to show that the improved action already removes most of the short-distance distortion, before any non-perturbative correction is applied.
• Six panels of the deviation of the static energy from the continuum estimate at \(\beta=7.280,\,7.373,\,7.596,\,7.825,\,8.000,\,8.200\). Why: the correction is measured, not computed, so it must be shown to be smooth and consistent across every \(\beta\) actually used.
• Interpolation of the correction factors \(\langle K\rangle\), and a comparison of \(\langle K\rangle\) obtained in two different renormalization schemes. Why: the correction must not depend on the scheme used to derive it; if it did, it would be putting physics in by hand.
• "Improvement of \(F_S\)" and "improvement of \(F_{Q\bar Q}\)", four panels each ($\beta = 9.67,\,8.2,\, 7.373,\,6.664\( at \)N_\tau=12$). Why: to display the size of the correction actually applied at each coupling -- the reader can see whether the corrected points move by a little or a lot.

• Scaling plots versus \(1/N_\tau^2\) at fixed \(rT\), separately for \(F_S^{\rm sub}\), \(F_{Q\bar Q}^{\rm sub}\), the effective coupling and \((V_S-F_S)/T\), at several \(rT\) values. Why: fixed \(rT\) is the only way to compare different lattices at the same physical point. Four different observables are shown because each enters a different conclusion, and an extrapolation that works for one need not work for another.
• Local versus global interpolation fits, log-log, for singlet and colour-averaged channels. Why: the interpolation in \((r,T)\) is a model. Showing two independent models agreeing in the overlap region is what makes the continuum result model-independent.
• Finite-volume correction at \(\beta=6.664\), \(N_\tau=4\), aspect ratio 4 versus 6. Why: the coarse large-volume lattices are the ones carrying the screening regime in Fig. 12, so the volume correction applied to them has to be exhibited.

• Volume dependence of \(m_S\) and of \(m_{Q\bar Q}\) at three temperatures (\(T=421,\,1687,\,5814\) MeV), two volumes, and for \(m_{Q\bar Q}\) also two source types (point and plane Polyakov loop sources). Why: an exponential fit at large \(r\) is exactly where a finite box bites, and where the fit range choice matters most. The horizontal reference lines for \(m_D/T\) and \(2m_D/T\) are drawn in the same figures so the reader can judge the fitted value against the perturbative expectation directly. Using two different sources is a further check that the extracted mass is a property of the medium and not of the operator.

Each figure closes one specific line of attack. Read the table as "if this figure did not exist, a
referee could say...".

figure	the objection it removes
1	--
2, 4	"your result depends on the additive renormalization constant" (a derivative kills it)
4, 5	"your singlet free energy is a Coulomb-gauge artifact"
3	"you never checked that the short-distance functional form is the perturbative one"
6	"the difference you claim is invisible on your own plot"
7	"you never compared with perturbation theory where it should work"
7, 8, 11, 12 (scale bands)	"you quote perturbation theory without an uncertainty"
9, 10	"you ignored the octet channel" / "the agreement is tuned" (ablation bands)
8 vs 12	"EQCD was only tested on one, gauge-dependent, observable"
13	"the exponential fit is a fit to a shape you never showed"
14, 15	"the screening mass may be pure lattice artifact" (compared to \(m_D\), to \(2m_S\), and to independent 3D-EQCD)
App. A	"your sea quarks are too heavy"
App. B	"your signal is at a few lattice spacings"
App. C	"your continuum limit is model-dependent"
App. D	"your fit range and box size drive the screening mass"

Structurally the paper has three layers, and the plot count is roughly the product of them:

(3 observables) x (2 EFT regimes) x (systematics that must be shown for each)

That is why fifteen main figures is not extravagance. The physics statement is a boundary --
"screening starts at \(rT\approx0.3\), weak coupling holds to \(rT\approx0.6\)" -- and a boundary can only
be established by showing where each description starts working and where it stops.

• Colour screening sets in at \(rT \approx 0.3\); below that, \(F_S\) is indistinguishable from the \(T=0\) static energy \(V_S(r)\).
• In the window \(0.3 \lesssim rT \lesssim 0.6\), weak-coupling EFTs work: pNRQCD for the short side, NLO EQCD for the screening side. Above \(rT\approx0.6\), non-perturbative magnetic physics takes over.
• The effective coupling peaks at \(r_{\max}\approx0.4/T\) with \(\alpha_{Q\bar Q}<0.5\) for \(T\gtrsim320\) MeV -- the quantitative basis for using weak-coupling pNRQCD for quarkonium above roughly \(300\) MeV.
• Asymptotic screening masses are perturbative in shape but not in size: \(m_S\approx(1.6\text{--}2.0)\,m_D\), and \(m_{Q\bar Q}\approx1.25\times 2m_D\), the latter consistent with independent 3D-EQCD lattice results.
• \(F_{Q\bar Q}\) feels the medium earlier than \(F_S\) purely because of the repulsive octet channel; the short-distance law is \(1/r^2\) (two-gluon exchange), not \(1/r\).
• Cyclic Wilson loops, sufficiently HYP smeared, agree with the Coulomb-gauge singlet free energy.