(135) $\displaystyle\frac{1}{\sigma_{1}}\frac{d\sigma_{g}}{ds}$
$\displaystyle=\Pi_{g}(\partial_{{\eta_{\ell}}},\partial_{{\eta_{h}}})s\left(\frac{sQ}{\mu_{s}}\right)^{{{\eta_{\ell}}}}\left(\frac{sQ}{\mu_{s}}\right)^{{\eta_{h}}}\frac{e^{-\gamma_{E}({{\eta_{\ell}}}+{{\eta_{h}}})}}{\Gamma(2+{{\eta_{\ell}}}+\eta_{h})}\frac{\sin(\pi{{\eta_{h}}})}{\sin(\pi({{\eta_{\ell}}}+{{\eta_{h}}}))}$
(136)

where

$\Pi_{g}(\partial_{{\eta_{\ell}}},\partial_{{\eta_{h}}})=\exp\Big{[}4C_{F}S(\mu_{h},\mu_{j})+4C_{F}S(\mu_{s},\mu_{j})+2C_{A}S(\mu_{h},\mu_{j})+2C_{A}S(\mu_{s},\mu_{j})\Big{]}\\\
\times\exp\Big{[}2A_{\gamma sg}(\mu_{s},\mu_{h})+2A_{\gamma
sqq}(\mu_{s},\mu_{h})+2A_{\gamma jg}(\mu_{j},\mu_{h})+4A_{\gamma
jq}(\mu_{j},\mu_{h})\Big{]}\\\ \times
H(Q,\mu_{h})\widetilde{j}_{q}\left(\partial_{\eta_{h}}+\ln\frac{Q\mu_{s}}{\mu_{j}^{2}}\right)\widetilde{j}_{\bar{q}}\left(\partial_{\eta_{h}}+\ln\frac{Q\mu_{sh}}{\mu_{j}^{2}}\right)\widetilde{j}_{g}\left(\partial_{\eta_{\ell}}+\ln\frac{Q\mu_{s\ell}}{\mu_{j}^{2}}\right)\widetilde{s}_{qq}(\partial_{\eta_{h}})\widetilde{s}_{g}(\partial_{\eta_{\ell}})$
(137)

and

$\displaystyle{{\eta_{\ell}}}$
$\displaystyle=\eta_{jg}+\eta_{sg}=2C_{A}A_{\Gamma}(\mu_{j},\mu_{s})$ (138)
$\displaystyle{{\eta_{h}}}$
$\displaystyle=2\eta_{j_{q}}+\eta_{sqq}=4C_{F}A_{\Gamma}(\mu_{j},\mu_{s})$
(139)

We have chosen the same jet scales for the light and heavy hemispheres
although one could also choose them to be different. Similarly, we have taken
the same soft scales for the left and right hemispheres.

One can read off from Eq. (135) that the large logs will be resummed for the
left-shoulder of heavy jet mass with the canonical scale choices

$\mu_{h}=Q,\quad\mu_{j}=\sqrt{r}\,Q,\quad\mu_{s}=r\,Q$ (140)

For thrust or the right shoulder of heavy jet mass, the canonical scale
choices are the same with $r$ replaced by $t$ or $s$ respectively. We have
verified that the expansion of the resummed distribution is independent of the
matching scales $\mu_{h},\mu_{j}$ and $\mu_{s}$ at order $\alpha_{s}$.

The quark channels have the same form as Eqs. (134) to (136) but with

$\Pi_{q}(\partial_{{\eta_{\ell}}},\partial_{{\eta_{h}}})=\exp\Big{[}(2C_{F}+C_{A})[S(\mu_{h},\mu_{j})+S(\mu_{s},\mu_{j})]+2C_{F}[S(\mu_{h},\mu_{j})+S(\mu_{s},\mu_{j})]\Big{]}\\\
\times\exp\Big{[}2A_{\gamma sq}(\mu_{s},\mu_{h})+2A_{\gamma
sqg}(\mu_{s},\mu_{h})+2A_{\gamma jq}(\mu_{j},\mu_{h})+2A_{\gamma
jq}(\mu_{j},\mu_{h})+2A_{\gamma jq}(\mu_{j},\mu_{h})\Big{]}\\\ \times
H(Q,\mu_{h})\widetilde{j}_{q}\left(\partial_{{{\eta_{h}}}}+\ln\frac{Q\mu_{s}}{\mu_{j}^{2}}\right)\widetilde{j}_{g}\left(\partial_{{{\eta_{h}}}}+\ln\frac{Q\mu_{s}}{\mu_{j}^{2}}\right)\widetilde{j}_{q}\left(\partial_{{{\eta_{\ell}}}}+\ln\frac{Q\mu_{s}}{\mu_{j}^{2}}\right)\widetilde{s}_{qg}(\partial_{{{\eta_{h}}}})\widetilde{s}_{q}(\partial_{{{\eta_{\ell}}}})$
(141)

and

$\displaystyle{{\eta_{\ell}}}$
$\displaystyle=\eta_{jq}+\eta_{sq}=2C_{F}A_{\Gamma}(\mu_{j},\mu_{s})$ (142)
$\displaystyle{{\eta_{h}}}$
$\displaystyle=\eta_{jq}+\eta_{jg}+\eta_{sqg}=(2C_{F}+2C_{A})A_{\Gamma}(\mu_{j},\mu_{s})$
(143)

The the final resummed distribution for thrust is

$\frac{d\sigma}{dt}=\frac{d\sigma_{g}}{dt}+2\frac{d\sigma_{q}}{dt}$ (144)

and similarly for heavy jet mass.

##  4 Analysis

In Section 3 derived a factorization formula for the left and right Sudakov
shoulders for heavy jet mass as well as the right Sudakov shoulder for thrust
(thrust has no left shoulder). We will now perform some cross checks on those
results. We first perform the fixed order expansion and compare to a numerical
computation of the exact NLO expression to verify the singular behavior. Then
we demonstrate that there are no non-global logarithms and discuss power
corrections.

###  4.1 Fixed-order expansions

First of all, we observe that the full resummed distributions are
renormalization-group invariance. This invariance has let us write the
evolution kernels in Eqs. (137) and (141) in a form that depends only on the
hard, jet and soft matching scales $\mu_{i}$, and not on $\mu$. The
cancellation of the $\mu$-dependence is non-trivial and requires the Casimirs
associated with the Sudakov double logs to cancel and the anomalous dimensions
to satisfy

$\gamma_{h}=\gamma_{jg}+2\gamma_{jq}+\gamma_{sqq}+\gamma_{sg}=\gamma_{jg}+2\gamma_{jq}+\gamma_{sqg}+\gamma_{sq}$
(145)

These relations can be checked explicitly using Eqs. (124), (113), (108) and
(110).

Expanding the resummed distributions to order $\alpha_{s}$ we find

$\frac{1}{\sigma_{1}}\frac{d\sigma}{dt}=3t+\frac{\alpha_{s}}{4\pi}\left\\{B_{1}t-\frac{3}{2}(2C_{F}+C_{A})\Gamma_{0}t\ln^{2}t\right.\\\
\left.\phantom{\frac{3}{2}}+\Big{[}3\gamma_{jg}+6\gamma_{jq}+2\gamma_{sg}+4\gamma_{sq}+2\gamma_{sqq}+4\gamma_{sqg}+3(C_{A}+2C_{F})\Gamma_{0}\Big{]}t\ln
t\right\\}$ (146)

for some $B_{1}$.The linear terms $3t$ and $B_{1}t$ are not predicted with NLL
resummation. So to be consistent we should remove all the terms linear in $t$.
This can be done to all orders by subtracting from the full resummed
distribution $\sigma(t)$ the boundary condition t $\sigma(1)$. That is, we
consider

$\frac{1}{\sigma_{1}}\frac{d\sigma^{{\text{sub}}}}{dt}\equiv\frac{1}{\sigma_{1}}\frac{d\sigma}{dt}-t\left[\frac{1}{\sigma_{1}}\frac{d\sigma}{dt}\right]_{t=1}$
(147)

which has only terms of the form $t\ln^{n}t$ to all orders in $\alpha_{s}$. We
use an analogous definition with $t$ replaced by $r$ or $s$ for the subtracted
form of the heavy jet mass distribution. Plugging in the anomalous dimensions

$\frac{1}{\sigma_{1}}\frac{d\sigma^{{\text{sub}}}}{dt}=\frac{\alpha_{s}}{4\pi}\left\\{-6(2C_{F}+C_{A})t\ln^{2}t\right.\\\
\left.+\Big{[}6C_{F}(1-4\ln 3)+C_{A}(1-12\ln 3)+4n_{f}T_{F}\Big{]}t\ln
t\right\\}+{\mathcal{O}}(\alpha_{s}^{2})$ (148)

This is shown in comparison to the NLO calculation in Fig. 6.

For the left shoulder of heavy jet mass, the expansion gives

$\frac{1}{\sigma_{1}}\frac{d\sigma^{{\text{sub}}}}{dr}=\frac{\alpha_{s}}{4\pi}\left\\{-\frac{1}{2}(2C_{F}+C_{A})\Gamma_{0}r\ln^{2}r+\Big{[}(C_{A}+2C_{F})\Gamma_{0}+\gamma_{jg}+2\gamma_{jq}+2\gamma_{sg}+4\gamma_{sq}\Big{]}r\ln
r\right\\}\\\
=\frac{\alpha_{s}}{4\pi}\left\\{-2(2C_{F}+C_{A})r\ln^{2}r+\Big{[}2C_{F}\left(1+4\ln\frac{4}{3}\right)+C_{A}\left(\frac{1}{3}+4\ln\frac{4}{3}\right)+\frac{4}{3}n_{f}T_{F}\Big{]}r\ln
r\right\\}$ (149)

This agrees with our fixed-order computation in Section 2 and with the leading
shoulder logarithms at NLO as can be seen in Fig. 6.

Breaking down the expression in Eq. (149) the anomalous dimensions which
appear are $\gamma_{jg}+2\gamma_{sg}$ from the gluon channel and
$\gamma_{jq}+2\gamma_{sq}$ from the quark and antiquark channels. So in each
channel only anomalous dimensions associated with light-hemisphere side are
contributing logarithms as order $\alpha_{s}$. This is a somewhat remarkable
feature of the factorization formula: although both sides contribute 1-loop
anomalous dimensions, as is required for renormalization-group invariance, Eq.
(145) only one side contributes logarithms. Mechanically, what happens that

$(\gamma_{h}\partial_{{\eta_{h}}}+\gamma_{\ell}\partial_{{\eta_{\ell}}})\left[r^{{{\eta_{\ell}}}+{{\eta_{h}}}}-1\right]\frac{\sin(\pi{{\eta_{\ell}}})}{\sin(\pi({{\eta_{\ell}}}+{{\eta_{h}}}))}=\gamma_{\ell}\,r\ln
r$ (150)

So the
$\frac{\sin(\pi{{\eta_{\ell}}})}{\sin(\pi({{\eta_{\ell}}}+{{\eta_{h}}}))}$
factor replaces the full anomalous dimension $\gamma_{\ell}+\gamma_{h}$ with
just $\gamma_{\ell}$.

For the right shoulder of heavy jet mass

$\frac{1}{\sigma_{1}}\frac{d\sigma^{{\text{sub}}}}{ds}=\frac{\alpha_{s}}{4\pi}\left\\{-(2C_{F}+C_{A})\Gamma_{0}s\ln^{2}s+\Big{[}2(C_{A}+2C_{F})\Gamma_{0}+2\gamma_{jg}+4\gamma_{jq}+2\gamma_{sqq}+4\gamma_{sqg}\Big{]}s\ln
s\right\\}\\\
=\frac{\alpha_{s}}{4\pi}\left\\{-4(2C_{F}+C_{A})s\ln^{2}s+\Big{[}4C_{F}\left(1-4\ln{6}\right)+\frac{2C_{A}}{3}\left(1-12\ln{6}\right)+\frac{8}{3}n_{f}T_{F}\Big{]}s\ln
s\right\\}$ (151)

In this case, only the anomalous dimensions in the heavy hemisphere contribute
at NLO. This distribution is also shown in Fig. 6 and compared to the exact
NLO calculation.

Figure 6: Comparison of the resummed distribution expanded to NLO (colored
curves) to the exact NLO distribution in the Sudakov shoulder region (blue
histograms) for heavy jet mass (left) and thrust (right). We include in the
prediction an offset and a linear term which are fit separately on either side
of the peak. The $C_{F}n_{f}T_{F}$ color structure has been scaled up by a
factor of 10 for clarity. The NLO histograms is this figure were computed
using event2 Catani:1996jh ; Catani:1996vz  with a cutoff of $10^{-12}$ and 12
trillion events and normalized to born cross section $\sigma_{0}$.

###  4.2 Non-global logarithms

When observables are sensitive to emissions only in a restricted region of
phase space, there can be an incomplete cancellation of virtual and real
emissions leading to non-global logarithms Dasgupta:2001sh . The classic
example is the light-hemisphere mass in $e^{+}e^{-}$ collisions. For the light
hemsiphere mass, emissions into the heavy hemisphere do not affect the value
of the light hemisphere mass, making it non-global. Generally, to be able to
resum logarithms using a factorization formula, one would like the condition
that the observable be small to force soft and collinear kinematics. This does
not happen with the light hemisphere mass, for example, since demanding it be
small does not prevent additional hard emissions into the heavy hemisphere.
For light-jet mass, the leading non-global logarithm contributes to the cross
section at order $d\sigma/d\rho_{\ell}\sim\alpha_{s}^{2}\ln^{2}\rho_{\ell}$ so
it is the same order as terms in NLL resummation. The leading logarithmic
series of non-global logs for the light jet mass and related observables is
understood and can be resummed Banfi:2002hw ; Schwartz:2014wha . Progress has
also been made on systematic higher-order resummation of non-global logarithms
Becher:2015hka ; Larkoski:2015zka . Thus, if there were non-global logs in the
Sudakov shoulders it would not pose an insurmountable obstacle. Nevertheless,
we will show that for the Sudakov shoulders of thrust and heavy jet mass, non-
global logs are absent.

For the right shoulder of thrust, the constraint in Eq. (82) is of the form
$t>x+y$ where $x$ and $y$ represent contributions to the mass of the heavy or
light hemispheres respectively from soft and collinear radiation near the
trijet region (i.e.
$x=m_{2}^{2}+m_{3}^{2}+2p_{2}k_{2}+2p_{3}k_{3}+2v_{\bar{1}}k_{\bar{1}}$ and
$y=m_{1}^{2}+2p_{1}k_{1}+v_{\bar{2}}^{\prime}k_{\bar{2}}+v_{\bar{3}}^{\prime}k_{\bar{3}}$
when the light jet is in the 1 direction). Since the constraint imposes a
lower bound on $t$, demanding $t\ll 1$ does not force $x$ and $y$ to be small,
suggesting that the Sudakov shoulder for thrust might be non-global. However,
we can rewrite the core convolution integral in Eq. (131) as

$\displaystyle\int_{0}^{\infty}dx\int_{0}^{\infty}dy$
$\displaystyle~{}x^{a-1}y^{b-1}(x+y-t)\Theta(x+y-t)$ (152)
$\displaystyle=\int_{0}^{\infty}dz\left[\int_{0}^{\infty}dx\int_{0}^{\infty}dy\
x^{a-1}y^{b-1}\delta(x+y-z)\right](z-t)\Theta(z-t)$ (153)
$\displaystyle=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}\left[\underbrace{\int_{0}^{\infty}dz\
z^{a+b-1}(z-t)}_{\text{scaleless}}-\underbrace{\int_{0}^{t}dz\
z^{a+b-1}(z-t)}_{\text{global}}\right]$ (154)
$\displaystyle=t^{1+a+b}\frac{\Gamma(a)\Gamma(b)}{\Gamma(2+a+b)}$ (155)

The first integral in brackets in Eq. (154) is divergent but either
independent of $t$ or linear in $t$, so it is smooth across $t=0$ and does not
generate Sudakov shoulder logarithms. The remaining integral has $z=x+y<t$, so
taking $t\ll 1$ does force $x,y\ll 1$. We conclude that the right shoulder of
thrust should be free of non-global logarithms. The actual divergence is an
artifact of expanding the phase space limits to leading power. In the full
theory, the divergences would cut off by the hard scale $Q$, but still would
not generate logarithms of $t$.

It is also worth noting that the scaleless integral in Eq. (154) does generate
a divergent term proportional to $t$. This would be the same order as terms in
the NNLL resummation of the Sudakov shoulder. The presence of such a term does
not imply that the factorization formula is valid only to NLL. Indeed, this
divergent contribution is smooth across $t=0$, suggesting that it contributes
similarly to the left and right sides of $\tau=\frac{1}{3}$ and therefore does
not give a discontinuity or a kink at $\tau=\frac{1}{3}$. In any case, since
we are only working to NLL in this paper, we can safely ignore it.

For heavy jet mass, the analogous constraint is in Eq. (78) which corresponds
to $x<r+y$ for the left shoulder or $x+s<y$ for the right shoulder, as in Eqs.
(132) and (133). We can rewrite Eq. (132) as

$f(r)=\int_{0}^{\infty}dx\int_{0}^{\infty}dy\
x^{a-1}y^{b-1}(r+y-x)\theta(r+y-x)=\frac{1}{a(a+1)}\int_{0}^{\infty}dy\
(r+y)^{a+1}y^{b-1}$ (156)

This integral is both UV and IR divergent (for $a,b>0$) and gets contributions
from all scales, suggesting, again, that it may generate non-global
logarithms. To separate out the UV and IR divergences, we can take two
derivatives with respect to $r$, leaving an integral which is UV finite for
$a,b>0$. We also introduce a new scale $R$ to separate small $r$ from large
$r$. Then we have

$f^{\prime\prime}(r)=\underbrace{\int_{0}^{R}dy\
(r+y)^{a-1}y^{b-1}}_{\text{global}}+\underbrace{\int_{R}^{\infty}dy\
(r+y)^{a-1}y^{b-1}}_{\text{regular in $r$}}$ (157)

Since we are interested in the region with $r\ll 1$, we can take $0<r<R\ll 1$.
Then the first integral in Eq. (157)is global, since it gets contributions
only from the region where $y\ll 1$ and $x\ll 1$ (we had integrated $x$ from 0
to $r+y\ll 1$ in Eq. (156). The second integral in Eq. (157) is regular as
$r\to 0$. Thus it does not contribute to any discontinuities or kinks near the
shoulder, at $r=0$. As with thrust, it may contribute terms linear in $r$, but
will not give any Sudakov shoulder logs. So only the soft and collinear
regions should contribute to the Sudakov shoulder logs for heavy jet mass, as
with thrust, and there are no non-global logarithms.

To complete the computation, as far as the Sudakov shoulder logs are
concerned, we have

$\displaystyle f^{\prime\prime}(r)$ $\displaystyle\cong\int_{0}^{R}dy\
(r+y)^{a-1}y^{b-1}$ (158) $\displaystyle=\int_{0}^{\infty}dy\
(r+y)^{a-1}y^{b-1}-\underbrace{\int_{R}^{\infty}dy\
(r+y)^{a-1}y^{b-1}}_{\text{regular in $r$}}$ (159)
$\displaystyle\cong\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}\frac{\sin(\pi
a)}{\sin(\pi(a+b))}r^{a+b-1}-R^{1-a-b}\frac{1}{{a+b-1}}$ (160)

where we have taken $R\gg r$ to simplify the second integral. Integrating
twice with respect to $r$ then gives

$f(r)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(2+a+b)}\frac{\sin(\pi
a)}{\sin(\pi(a+b))}r^{1+a+b}-\frac{r^{2}}{2}R^{a+b-1}\frac{1}{{1-a-b}}+c_{1}r+c_{0}$
(161)

At small $a$ and $b$ (these are proportional to $\alpha_{s}$), the second term
on the right-hand side is suppressed by a factor of $\frac{r}{R}$ compared to
the first term, so it only gives power corrections and no Sudakov shoulder
logs, as anticipated. We should fix the integration constants $c_{1}$ and
$c_{2}$ so that the expansion of $f(r)$ at small $a$ and $b$ only has terms of
the form $r\ln^{n}r$ with $n>0$. The constant term we can simply discard,
$c_{0}=0$. To fix $c_{1}$ we should set $c_{1}=-f(1)$. This corresponds to
integrating $f^{\prime}(r)$ from $1$ to $r$. These integration constants were
used in Eq. (147).

In summary, the heavy jet mass distribution at a value of
$\rho\approx\frac{1}{3}$ does get contributions from phase space regions with
jets whose masses are not small. In this sense it is similar to light jet mass
near $\rho_{\ell}=0$ which gets contributions from phase space regions where
$\rho$ is not small. However, the contributions corresponding to heavy jets
for the Sudakov shoulder do not generate large logarithms. This is because the
phase space regions with heavy jets can contribute to both
$\rho\lesssim\frac{1}{3}$ and $\rho\gtrsim\frac{1}{3}$ and are smooth across
$\rho=\frac{1}{3}$. All the contributions to the distribution that are not
smooth across $\rho=\frac{1}{3}$ come from the regions with one nearly-
massless jet in the light hemisphere and two nearly massless jets in the heavy
hemisphere. There is no analog of this continuity argument for light jet mass,
which cannot have $\rho_{\ell}<0$. Thus, the Sudakov shoulders of heavy jet
mass (and thrust) are free of non-global logarithms.

###  4.3 Power corrections

In resummed distributions, there are typically different types of power
corrections. For threshold resummation, near $\rho=0$ for example, there can
be power corrections of order $\frac{{\Lambda_{\text{QCD}}}}{Q}$ associated
with the strong dynamics of QCD. There can also be hard power corrections,
suppressed by additional powers of $\rho=\frac{m_{H}^{2}}{Q^{2}}$ where
$m_{H}$ is the mass of the heavy jet. The $\Lambda_{\text{QCD}}$ power
corrections are often modeled with parameters fit to data. This allows for
predictivity closer to threshold than with just the resummed distribution
alone, although one cannot get too close to threshold since more and more non-
perturbative parameters than become relevant. The hard power corrections are
typically accounted for in matching to an exact fixed order expression at
large $\rho$.

For Sudakov shoulder resummation, it is not clear whether
$\frac{{\Lambda_{\text{QCD}}}}{Q}$ power corrections are important near the
trijet threshold. On the one hand, the resummed distribution involves
evaluating $\alpha_{s}$ at scales such as $\mu_{s}=Qr$ which can reach
${\Lambda_{\text{QCD}}}$ for small enough $r$. On the other hand, the shoulder
is intrinsically perturbative, associated with fixed-order phase space
boundaries, so one might expect that it might be invisible to non-perturbative
physics.

The hard power corrections for the Sudakov shoulder are more interesting. In
Section 4.2 we argued that at leading power all the non-analytic behavior near
the shoulder is determined by soft and collinear physics. That is, there are
no non-global logarithms. One can see this from Eq. (161). The quantities $a$
and $b$ are to be replaced by ${{\eta_{\ell}}}$ and ${{\eta_{h}}}$ in the
resummed distribution, which are parameterically of the form
$\eta\sim\alpha_{s}\Gamma_{0}\ln r$. Thus at small $\alpha_{s}$, all the terms
of the form $r\ln^{n}r$ will come from the expansion of the first term on the
right-hand side in Eq. (161). On the other hand if $r$ is sufficiently small
then $a+b$ can be of order $1$. As $a+b$ nears $1$ a pole from the
$\sin^{-1}(\pi(a+b))$ factor in Eq. (161) is approached. However, when
$a+b\approx 1$, power-suppressed term is no-longer power suppressed. Indeed,
it has precisely the behavior needed to remove the singular behavior from the
leading power term. We show this in Fig. 7.

Figure 7: The solid red curve shows
$r^{1+a+b}\frac{\Gamma(a)\Gamma(b)}{\Gamma(2+a+b)}\frac{\sin(\pi
a)}{\sin(\pi(a+b))}$ with $a=b=-6\frac{\alpha_{s}}{\pi}\ln r$ and
$\alpha_{s}=0.1$. The pole is at $a+b=1$. The dashed curve shows this same
function once the power-suppressed term
$-\frac{r^{2}}{2}R^{a+b-1}\frac{1}{1-a-b}$ is added with $R=1$. The power
suppression is apparent at small $r$, where the difference between the two
curves is negligible. The pole at $a+b=2$ is not canceled and appears as the
spike at $r\approx 0.0003$.

To see what is happening analytically, notating that the leading power
expression scales like $r^{1+a+b}$ we can write the subleading power
expression as

$r^{2}R^{a+b-1}=r^{1+a+b}\left(\frac{r}{R}\right)^{1-a-b}$ (162)

For $a+b\ll 1$ there is a $\frac{r}{R}$ linear power suppression. However for
$a+b\sim 1$ there is no power suppression at all; for $a+b=1$ this expressions
reduces to $r^{2}$ as does $r^{1+a+b}$. In effect, the scaling dimensions of
the leading power and subleading power pick up such large anomalous
contributions that their relative scaling changes.

Taking $R\to\infty$ gives the leading contribution. The first subleading power
contribution in this limit cancels the pole at $a+b=1$. To cancel subsequent
poles, one can use the exact integrated form of Eq. (157). This effectively
replaces $f(r)$ in Eq. (161) by

$f(r)=r^{1+a+b}e^{-ib\pi}\frac{1}{a(a+1)}B_{-\frac{R}{r}}(b,2+a)+c_{1}r+c_{0}$
(163)

where $B_{z}(x,y)$ is the incomplete Euler $\beta$ function and $c_{0}$ and
$c_{1}$ are again integration constants to be fixed with physical boundary
conditions.

##  5 Discussion

Next, we want to evaluate the resummed distribution numerically and compare to
fixed order, to see the effect of the higher order logarithms. There are a
number of issues which complicate the analysis, compared with typical
threshold resummation of large logarithms.

First, the relevant domain of the observable is rather small for Sudakov
shoulders. For example, for thrust in the threshold limit, although the
logarithms are largest at small $\tau$, power corrections and subleading logs
are also large there. Typical fits restrict $\tau\gtrsim 0.1$ where
perturbative control is best. For example, with $Q=92$ GeV, Ref. Becher:2008cf
used $0.1<\tau<0.24$ for their $\alpha_{s}$ fits to thrust while Ref.
Abbate:2010xh  took $\tau\gtrsim\frac{Q}{6\,\text{GeV}}=0.066$. For the right
shoulder of thrust which begins at the 3-parton maximum $\tau=\frac{1}{3}$ if
one excludes the region up to $\frac{1}{3}+0.1=0.43$ there is no cross section
or phase space left! Moreover, the 4-particle phase space forces $\tau\lesssim
0.42$, so there is another Sudakov shoulder at this thrust value whose logs
must be resummed separately. So it is not clear if there is a region on the
right shoulder where the resummed formula might even be valid. For the left
shoulder, in contrast, one can exclude the region with
$r=\frac{1}{3}-\rho<0.1$ which still leaves a region of
$0.1\lesssim\rho\lesssim 0.23$ in which Sudakov shoulder logarithms might be
important and renormalization-group improved perturbation theory could be
valid.

Second, in the threshold region, the logarithms of thrust are of the form
$\frac{d\sigma}{d\tau}\sim\alpha_{s}^{n}\frac{\ln^{m}\tau}{\tau}$. In
contrast, the logarithms near the shoulder region are of the form
$\frac{d\sigma}{dt}\sim\alpha_{s}^{n}t\ln^{m}t$. So they are suppressed
effectively by $t^{2}$ compared to the threshold region. The thrust and heavy
mass distributions are indeed finite at the trijet threshold to all orders
while they are divergent at the dijet threshold. Despite this additional
suppression, the logarithms are noticeable, as can be been in Fig. 2.

Third, in the important left-shoulder region for heavy jet mass, the resummed
distribution has usually singular behavior. Let us recall the form of the
resummed heavy jet mass distribution in the region $r=\frac{1}{3}-\rho\ll 1$
when the light hemisphere has a gluon jet from Eq. (135)

$\frac{1}{\sigma_{1}}\frac{d\sigma_{g}}{dr}=\Pi_{g}(\partial_{{\eta_{\ell}}},\partial_{{\eta_{h}}})r\left(\frac{rQ}{\mu_{s\ell}}\right)^{{{\eta_{\ell}}}}\left(\frac{rQ}{\mu_{sh}}\right)^{{\eta_{h}}}\frac{e^{-\gamma_{E}({{\eta_{\ell}}}+{{\eta_{h}}})}}{\Gamma(2+{{\eta_{\ell}}}+\eta_{h})}\frac{\sin(\pi{{\eta_{\ell}}})}{\sin(\pi({{\eta_{\ell}}}+{{\eta_{h}}}))}$
(164)

with ${{\eta_{h}}}$ and ${{\eta_{\ell}}}$ in Eqs.(143) and (142). This
expression has singularities whenever
${{\eta_{\ell}}}+{{\eta_{h}}}\in\mathbb{Z}$.

Choosing canonical scales as in Eq. (140) at leading logarithmic level gives

${{\eta_{\ell}}}+{{\eta_{h}}}=-\frac{\alpha_{s}}{2\pi}(C_{A}+2C_{F})\Gamma_{0}\ln\frac{\mu_{j}}{\mu_{s}}=-\frac{\alpha_{s}}{4\pi}(C_{A}+2C_{F})\Gamma_{0}\ln
r$ (165)

The singularity ${{\eta_{\ell}}}+{{\eta_{h}}}=0$ occurs when
$\mu_{j}=\mu_{s}$, which happens at $r=1$. At $r=1$ there are no logarithms,
so this singularity is entirely removed by the subtraction in Eq. (147). That
is,

$\frac{d\sigma^{\text{sub}}}{dr}=\frac{d\sigma}{dr}-r\left[\frac{d\sigma}{dr}\right]_{r=1}$
(166)

is regular at $r=1$. Note, however, if the soft and jet scales meet at some
lower scale, this singularity may be reintroduced.

The singularity at ${{\eta_{\ell}}}+{{\eta_{h}}}=1$ is more troublesome.
Similar singularities have been seen in other processes, such as Drell-Yan or
Higgs production at small $p_{T}$ Catani:1996yz ; Frixione:1998dw ;
Becher:2010tm ; Monni:2016ktx ; Banfi:2018mcq  or the jet shape Seymour:1997kj
; Cal:2019hjc . Writing $L=\ln\frac{1}{r}$, resummation at order NLL is meant
to get right all terms of order $\alpha_{s}^{n}L^{j}$ with $j\geq 2n-1$ in
$R(r)$ or equivalently all terms of order $\alpha_{s}^{n}L^{j}$ with $j\geq n$
in the exponent, i.e. in $\ln R(r)$. In the notation of Catani:1992ua , we can
write

$\ln\frac{d\sigma^{\text{sub}}}{dr}=Lg_{1}(\alpha_{s}L)+g_{2}(\alpha_{s}L)+\cdots=\sim\cdots-\ln\sin(\pi({{\eta_{\ell}}}+{{\eta_{h}}}))+\cdots$
(167)

with $g_{1}(\alpha_{s}L)$ and $g_{2}(\alpha_{s}L)$ completely fixed by the
expansion and reorganization of our resummed expression. Normally, when
$\alpha_{s}L\sim 1$ then we must go to higher order in RG-improved
perturbation theory; at NNLL level, we would have additionally
$Lg_{3}(\alpha_{s}L)$ which would extend the validity of the theoretical
prediction. Here, instead we find a singularity in the exponent: $\ln
R^{\text{sub}}$ is infinite at $\alpha_{s}L\sim 1$ due to the
${{\eta_{\ell}}}+{{\eta_{h}}}=1$ singularity. Therefore, going beyond NLL
would not allow us to make perturbative predictions beyond where the
singularity occurs. Instead, the singularity is canceled by including
subleading power effects, as discussed in Section 4.3 and shown in Fig. 7.

The singularity at $\alpha\ln r\sim 1$ is reminiscent of the Landau pole in
QCD. There, already at 1-loop one can see a pole in the running coupling at
$\mu=\Lambda_{\text{QCD}}$. With 2-loop or higher order running, the precise
location of $\Lambda_{\text{QCD}}$ moves around, but cannot be surpassed. Thus
what we see here is a kind of Sudakov Landau pole. Using the LL form in Eq.
(165) it occurs at

$r=\exp\left[-\frac{4\pi}{(C_{A}+2C_{F})\alpha_{s}\Gamma_{0}}\right]=\exp\left[-\frac{3\pi}{17\alpha_{s}}\right]$
(168)

For $\alpha_{s}=0.119$ this gives $r\approx 0.01$. Using the NLL expressions
for ${{\eta_{\ell}}}$ and ${{\eta_{h}}}$ with canonical scale choices, Eq.
(140) the pole ascends to $r\approx 0.06$. Thus we cannot expect the leading-
power NLL resummed distribution to be predictive between $0.27<\rho<0.39$.
This essentially excludes the entire region on the right shoulder, but leaves
the region with $\rho\lesssim 0.27$ as potentially viable for a precise
prediction. To stay well away from the singular region however, one must take
$\rho$ smaller, $\rho\lesssim 0.2$ where the logarithms are no longer
particularly large.

We emphasize that excluded range is larger than that associated with strong
coupling. With 2-loop running and $\alpha_{s}(m_{Z})=0.119$, we find
$\alpha_{s}(m_{s})=1$ at $r=0.005$. Thus the singularity comes in at a factor
of 10 larger values of $r$ than where the soft scale probes strong dynamics.
This is because the singularity is associated with the cusp anomalous
dimension, not the QCD $\beta$ function: the two Landau poles are unrelated.

Because of the Sudakov Landau pole in the resummed distribution it is
difficult to make quantitative predictions, particularly at the NLL level,
without a better understanding of the power corrections. There are a number of
approaches that could be applied to ameliorate the problem. In Frixione:1998dw
, a similar pole in the Drell-Yan spectrum at small $p_{T}$ (at
$q^{\star}=m_{Z}\exp(-\frac{3\pi}{8\alpha_{s}})$ Becher:2011xn ) was shown to
be associated with power-suppressed region of small impact parameter, but
could be softened with higher-order resummation. In Monni:2016ktx  it is
argued that one could also do resummation in momentum space directly with a
modified expansion of the Sudakov radiator. Related ideas can be found in
Ebert:2016gcn ; Cal:2019hjc . It will be important to understand which of
these approaches might apply for Sudakov shoulder resummation, but we do not
attempt a complete analysis here.

At the LL level, however, because the Sudakov Landau pole is very close to the
shoulder, we can at least begin to get a quantitative feel of how important
resummation is. Consider the LL distribution using canonical scales in Eq.
(140). When the jet in the light hemisphere is a gluon, it has the form as in
Eq. (164) with $\Pi_{g}$ from Eq. (137) becoming

$\Pi_{g}=e^{-\frac{\alpha_{s}}{8\pi}\Gamma_{0}(C_{A}+2C_{F})\log^{2}r}\left[1-\frac{\alpha_{s}}{8\pi}\Gamma_{0}\left(C_{A}\partial^{2}_{{\eta_{\ell}}}+2C_{F}\partial_{{\eta_{h}}}^{2}\right)\right]$
(169)

Note that we include every term with $\Gamma_{0}$ in it for leading-log
resummation, not just the exponential prefactor. Including only the prefactor
would give the double-logarithmic approximation, as used in previous work on
resummation of the $C$ parameter Sudakov shoulder Catani:1997xc . We subtract
off from the resummed distribution $r$ times its $r\to 1$ limit as done in Eq.
(147). Note that this subtraction must be done before setting canonical
scales. We then match to the fixed order LO+NLO calculation by subtracting
from the resummed distribution its expansion to order $\alpha_{s}$. In this
case, the matching subtraction is

$\frac{1}{\sigma_{1}}\frac{d\sigma^{\text{match}}}{dr}=\frac{\alpha_{s}}{4\pi}(C_{A}+2C_{F})\Gamma_{0}(r\ln
r-\frac{1}{2}r\ln^{2}r)$ (170)

Finally we include the subtraction of the first subleading power contribution.
For the gluon channel this amounts to subtracting:

$\frac{1}{\sigma_{1}}\frac{d\sigma_{g}^{R}}{dr}=\frac{r^{2}}{2R}\Pi_{g}(\partial_{{\eta_{\ell}}},\partial_{{\eta_{h}}})\left(\frac{RQ}{\mu_{s\ell}}\right)^{{{\eta_{\ell}}}}\left(\frac{RQ}{\mu_{sh}}\right)^{{\eta_{h}}}\frac{e^{-\gamma_{E}({{\eta_{\ell}}}+{{\eta_{h}}})}}{\Gamma({{\eta_{\ell}}})\Gamma({{\eta_{h}}})}\frac{1}{(1-{{\eta_{\ell}}}-\eta_{h})}$
(171)

from the resummed distribution. The resulting LL resummed and matched result
at $Q=m_{Z}$ and $\alpha_{s}(Q)=0.119$ with $R=\frac{1}{3}$ for the left
shoulder is shown in Fig. 8. To be clear, this is not the complete power
correction, but amounts to integrating the leading power soft and collinear
matrix elements outside of their formal region of validity upto the kinematic
limit of $r=\frac{1}{3}$.

Figure 8: Resummation of the left Sudakov shoulder for heavy jet mass at
leading-logarithmic level (upper curve) compared to NLO (lower curve). The
strong coupling constant is fixed to $\alpha_{s}=0.119$

##  6 Conclusion

Thrust $\tau$ and heavy jet mass $\rho$ are two of the most important
observables at $e^{+}e^{-}$ colliders. They have been used for decades for
tests of precision QCD and measurements of $\alpha_{s}$. At leading order in
perturbation theory, both $\rho$ and $\tau$ are phase-space limited to be less
than $\frac{1}{3}$ and have a non-vanishing slope as $\frac{1}{3}$ is
approached. At next-to-leading order, thrust behaves like
$\alpha_{s}^{2}(\tau-\frac{1}{3})\ln^{2}(\tau-\frac{1}{3})$ for
$\tau>\frac{1}{3}$ so that the slope diverges as $\frac{1}{3}$ is approached
from the right. This behavior is called a right Sudakov shoulder. Heavy jet
mass has a slope which diverges as $\rho$ nears $\frac{1}{3}$ both from the
left and the right: it has two Sudakov shoulders. The left shoulder of heavy
jet mass is particularly important as the large logarithms can extend well
into the region where $\alpha_{s}$ fits are typically done
($0.1\lesssim\rho\lesssim 0.24$). Thus understanding and resumming its Sudakov
shoulders could be very important for improving agreement of theoretical
predictions with data and subsequent extractions of $\alpha_{s}$. We also
point out that it has been noted recently in the literature that in the
context of other event shape observables such as fractional moments of energy-
energy correlation Banfi:2018mcq  or projected energy correlators Chen:2020vvp
, one must resort to a joint resummation of the Sudakov shoulders and endpoint
peaks.

We derived a factorization formula for both thrust and heavy jet mass in the
Sudakov shoulder region. The basic mechanism for generating Sudakov shoulder
logs is when a soft or collinear emission goes into one hemisphere a global
constraint such as $m^{2}<\rho-\frac{1}{3}$ transfers large logs from the
emissions to the shoulder. Although the constraint seems non-local, involving
both hemispheres, and therefore might violate factorization, we show that it
does not. Moreover regions of large jet mass do not contribute Sudakov
shoulder logs, showing that there is no non-global log contribution in the
shoulder region. We checked our factorization formula by expanding to NLO and
comparing to the exact numerical NLO calculation very close to the shoulder
region. As can be seen in Fig. 6 the agreement is excellent.

The calculation involves some unusual ingredients. Since the emissions come
off a trijet configuration with 2 quarks and 1 gluons, there is no azimuthal
symmetry (unlike the threshold case), and the polarization of the gluon
affects the spectrum. At leading order, only a uniform azimuthal angle
integral was needed for the resummed expression, but in general polarized
splitting function may be necessary. We also saw the appearance of Gieseking’s
constant, a transcendentally-two number. Although it also drops out of the NLL
expression, at higher orders it or related constants may be involved.

The resummed distribution for heavy jet mass has a term of the form
$\sin^{-1}(\pi\eta)$ with $\eta\sim\alpha_{s}\Gamma_{0}\ln r$, where
$r=\frac{1}{3}-\rho$. The expansion near $\alpha_{s}=0$ (or $\ln r=0$)
produces the leading and next-to-leading logarithmic series: terms like
$\alpha^{n}\ln^{2n}r$. However, there is also a pole at $\eta=1$. This pole in
the resummed distribution is not due to the running coupling – it is present
even with $\beta(\alpha_{s})=0$ – but due to the cusp anomalous dimension.
Thus it is a kind of Sudakov Landau pole. Similar behavior has been seen
before, in the Drell-Yan process at small $p_{T}$, for example Frixione:1998dw
; Catani:1996yz ; Monni:2016ktx . In both cases there is a connection between
the pole and subleading power effects (subleading in $r$ for the shoulder, or
in impact parameter $b$ for Drell-Yan). We show that subleading power terms
can in fact cancel the $\eta=1$ pole but do not affect the NLL series. This
implies that a better understanding of power corrections will be necessary to
establish proper theoretical uncertainty on the resummed distribution. There
are many approaches that may help improve the convergences of the resummed
distribution Becher:2007ty ; Ebert:2016gcn ; Monni:2016ktx .

Although our results are only valid to NLL level, the factorization formula
applies to all orders. In fact, since the anomalous dimensions of the jet and
hard functions are known to two loops, and therefore the soft function
anomalous dimension as well by renormalization-group invariance, NNLL
resummation should be possible. At NNLL level, terms linear in $\rho$ or
$\tau$ are determined. The slope can be discontinuous from the left to right
side of the shoulder, as it is already at LO. This discontinuity should be
computable. However, because there is also a linear term in the distribution
not associated with the shoulder, confirming the predictions at NNLL will be
challenging. Nevertheless, pushing the limits of Sudakov shoulder resummation,
not just for $e^{+}e^{-}$ event shapes but for collider observables more
broadly, provides opportunities to improve our understanding of precision QCD.

## Acknowledgements

The authors would like to thank Vicent Mateu for valuable conversations that
helped inspire this work as well as Thomas Becher and Pier Monni for valuable
feedback. This study was supported by the U.S. Department of Energy under
contract DE-SC0013607.

##  Appendix A Calculation of the one-loop trijet soft function

Before using rotational invariance all the integrals needed for the one-loop
trijet soft functions are of the form

$I_{n_{a},n_{b},n_{c},n_{d},N}(q)=\int
d^{d}k\frac{{n_{a}}\cdot{n_{b}}}{({n_{a}}\cdot k)({n_{b}}\cdot
k)}\delta(k^{2})\theta(k^{0})\delta\left(q-\frac{2}{3}N\cdot
k\right)\theta({n_{c}}\cdot k-{\bar{n}_{c}}\cdot k)\theta({n_{d}}\cdot
k-{\bar{n}_{d}}\cdot k)$ (172)

where the $N$ is selected from the direction vectors

$\displaystyle n_{1}$ $\displaystyle=(1,0,0,1),\quad
n_{2}=\left(1,0,\frac{\sqrt{3}}{2},-\frac{1}{2}\right),\quad
n_{3}=\left(1,0,-\frac{\sqrt{3}}{2},-\frac{1}{2}\right)$ (173) $\displaystyle
N_{1}$ $\displaystyle=(2,0,0,-2),\quad N_{2}=(2,0,\sqrt{3},3),\quad
N_{3}=(2,0,-\sqrt{3},3)$ (174)

and $n_{a}\cdots n_{d}$ only from the $n_{j}$ in Eq. (173). The
$\theta$-functions restrict the phase space to one of the sextants in Fig. 4.

The Wilson lines in the trijet configuration have an $S_{3}$ symmetry which
includes a $Z_{3}$ rotational invariance and a reflection symmetry. We can use
the rotational invariance to rotate the Wilson lines so that they always point
in the $n_{1}$ and $n_{2}$ directions. Thus we only need to consider integrals
as in Eq. (100):

$I_{n_{a},n_{b},N}(q)=\int d^{d}k\frac{{n_{1}}\cdot{n_{2}}}{({n_{1}}\cdot
k)({n_{2}}\cdot k)}\delta(k^{2})\theta(k^{0})\delta\left(q-\frac{2}{3}N\cdot
k\right)\theta({n_{a}}\cdot k-{\bar{n}_{a}}\cdot k)\theta({n_{b}}\cdot
k-{\bar{n}_{b}}\cdot k)$ (175)

When $N=n_{i}$ is one of the Wilson line directions, then only two sextants
are relevant, $I_{1}(q)$ and $I_{2}(q)$ from Fig. 5:

$\displaystyle I_{1}(q)=I_{n_{2},n_{3},n_{1}}(q)\quad
I_{2}(q)=I_{n_{1},n_{2},n_{3}}(q)$ (176)

When $N=N_{1}=2\bar{n}_{1}$, there are two configurations relevant, with both
Wilson lines adjacent to the $n_{\bar{1}}$ measurement region or just one of
them adjacent two it. The two integrals are, as in Fig. 5,

$\displaystyle I_{3}(q)=I_{\bar{n}_{1},\bar{n}_{2},2\bar{n}_{3}}(q),\quad
I_{4}(q)=I_{\bar{n}_{1},\bar{n}_{3},2\bar{n}_{2}}(q)$ (177)

where $2\bar{n}_{3}$ comes from rotating $N_{1}$ as the Wilson lines are
rotated to the $n_{1}$ and $n_{2}$ directions.

The remaining integrals involve $N_{2}$ and $N_{3}$. These vectors are not
lightlike, but they are related by a $Z_{2}$ symmetry. (Recall that the origin
of the asymmetry between $N_{1}$ and $N_{2}/N_{3}$ is that $n_{1}$ points to
the light-hemisphere which affects the soft projections in the factorization
formula.) Since $N_{2}$ and $N_{3}$ are related by a reflection in the $y$
direction, which is a symmetry of the Wilson lines, if we know the integral
for all Wilson line configurations for $N_{2}$ we know it for $N_{3}$ as well.
So there are 3 possibilities, corresponding to the location of the 3
measurement regions with respect to the Wilson line. Rotating the $N_{2}$
Wilson line by $\frac{2\pi}{3}$ and $\frac{4\pi}{3}$ gives

$N_{2}^{\prime}=(2,0,2\sqrt{3},0),\quad
N_{2}^{\prime\prime}=(2,0,-\sqrt{3},-3)$ (178)

Thus the last 3 integrals we need are

$I_{5}(q)=I_{\bar{n}_{1},\bar{n}_{2},N_{3}}(q),\quad
I_{6}(q)=I_{\bar{n}_{1},\bar{n}_{3},N_{3}^{\prime}}(q),\quad
I_{7}(q)=I_{\bar{n}_{2},\bar{n}_{3},N_{3}^{\prime\prime}}(q)$ (179)

To perform the integrals, we parameterize the phase space with lightcone
components in some direction $n_{1}$:

$k^{\mu}=k_{+}\frac{n_{1}^{\mu}}{2}+k_{-}\frac{\bar{n}_{1}^{\mu}}{2}+k_{\perp}^{\mu}$
(180)

so that

$d^{d}k=\frac{1}{2}d\Omega_{d-2}k_{\perp}^{d-3}dk_{\perp}dk_{+}dk_{-}=\frac{1}{2}d\Omega_{d-3}\sin^{d-4}\theta
d\theta k_{\perp}^{d-3}dk_{\perp}dk_{+}dk_{-}$ (181)

and

$\delta(k^{2})=\delta(\vec{k}_{\perp}^{2}-k_{+}k_{-})$ (182)

The integral over $k_{\perp}$ can be calculated using the $\delta$ function
and then $k_{+}$ and $k_{-}$ recsaled by $q$ to obtain the $q$ dependence
$q^{d-5}$ as expected by dimensional analysis. We also introduce
$t=\sqrt{\frac{1-\cos\theta}{1+\cos\theta}}$ to rationalize
$\sin\theta=\sqrt{1-\cos^{2}\theta}$.

With these preliminariles, the integral $I_{1}(q)=I_{n_{2},n_{3},n_{1}}(q)$
takes the form in $d=4-2\epsilon$ dimensions

$\displaystyle I_{1}$
$\displaystyle=\frac{3^{1-2\epsilon}\Omega_{1-2\epsilon}}{q^{1+2\epsilon}}\int_{0}^{\infty}\frac{dt}{t^{2\epsilon}(1+t^{2})^{1-2\epsilon}}\int_{0}^{\infty}\frac{dk_{+}}{k_{+}^{\epsilon}}\frac{1}{1+3k_{+}-\frac{2\sqrt{3k_{+}}(1-t^{2})}{1+t^{2}}}$
(183)
$\displaystyle\times\theta\left(1+k_{+}\right)\theta\left(-1+k_{+}-\frac{2\sqrt{3k_{+}}(1-t^{2})}{1+t^{2}}\right)\theta\left(-1+k_{+}+\frac{2\sqrt{3k_{+}}(1-t^{2})}{1+t^{2}}\right)$
(184)

The $\theta$-functions impose that

$0<t<\infty,\quad
k_{+}>\frac{7-10t^{2}+7t^{4}}{(1+t^{2})^{2}}+\frac{4\sqrt{3(1-t^{2})^{2}(1-t^{2}+t^{4})}}{(1+t^{2})^{2}}\equiv
k_{c}$ (185)

To handle the UV divergence as $k_{+}\to\infty$, we can add and subtract the
integral $I_{2}^{\text{div}}$ over the integrand expanded at large $k^{+}$.
This subtraction term requires the integral

$\int_{0}^{\infty}\frac{dt}{t^{2\epsilon}(1+t^{2})^{1-2\epsilon}}\int_{k_{c}}^{\infty}\frac{dk_{+}}{k_{+}^{\epsilon}}\frac{1}{3k_{+}}=\frac{\pi}{6}\frac{1}{\epsilon}+\frac{1}{9}\left(-5\kappa+\pi\ln\frac{64}{3}\right)+\mathcal{O}(\epsilon)$
(186)

where $\kappa=\text{Im}\,\text{Li}_{2}\,e^{\frac{\pi i}{3}}$ is Gieseking’s
constant. Then $I_{1}^{\text{fin}}=I_{1}-I^{\text{div}}_{1}$ is finite and can
be expanded in $\epsilon$ and integrated order by order. Eventually, we arrive
at

$I_{1}(q)=\pi^{-\epsilon}e^{-\gamma_{E}\epsilon}\left(\frac{2}{3}\right)^{2\epsilon}\frac{1}{q^{1+2\epsilon}}\left(\frac{1}{\epsilon}+\ln
3-\frac{7}{2}\ln 2-\frac{3}{2\pi}\kappa+\mathcal{O}(\epsilon)\right)$ (187)

The calculation for other soft integrals are similar, and the results up to
order $\mathcal{O}(\epsilon)$ are summarized as followed:

$\displaystyle I_{1}(q)$
$\displaystyle=\mathcal{N}\bigg{[}\frac{1}{\epsilon}+\ln 3-\frac{7}{2}\ln
2-\frac{3}{2\pi}\kappa+\epsilon\bigg{(}\frac{18}{5\pi}c_{1}+\frac{8}{\pi}c_{2}-\frac{103}{180}\pi^{2}+\frac{3}{\pi}\kappa\ln
2+\frac{10}{3}\ln^{2}2$ $\displaystyle-\frac{7}{2}\ln 2\ln
3+\frac{17}{40}\ln^{2}3+\frac{5}{12}\text{Li}_{2}\left(\frac{1}{4}\right)\bigg{)}+\mathcal{O}(\epsilon^{2})\bigg{]}$
(188) $\displaystyle I_{2}(q)$
$\displaystyle=\mathcal{N}\bigg{[}{\frac{3}{\pi}\kappa-\ln
2}+\epsilon\left(2\ln^{2}2+\frac{3}{2\pi}c_{3}-\frac{6}{\pi}\kappa\ln
2\right)+\mathcal{O}(\epsilon^{2})\bigg{]}$ (189) $\displaystyle I_{3}(q)$
$\displaystyle=\mathcal{N}\left[{\frac{3}{\pi}\kappa+\ln
2}+\epsilon\left(-2\ln^{2}2-\frac{6}{\pi}\kappa\ln
2+\frac{3}{2\pi}c_{5}\right)+\mathcal{O}(\epsilon^{2})\right]$ (190)
$\displaystyle I_{4}(q)$
$\displaystyle=\mathcal{N}\bigg{[}{-\frac{3}{2\pi}\kappa+\frac{3}{2}\ln
2}+\epsilon\left(-3\ln^{2}2+\frac{3}{2\pi}c_{4}+\frac{3}{\pi}\kappa\ln
2\right)+\mathcal{O}(\epsilon^{2})\bigg{]}$ (191) $\displaystyle I_{5}(q)$
$\displaystyle=\mathcal{N}\bigg{[}{\frac{3}{\pi}\kappa+\ln
2}+\epsilon\left(-2\ln^{2}2+\frac{3}{2\pi}c_{7}-\frac{6}{\pi}\kappa\ln
2\right)+\mathcal{O}(\epsilon^{2})\bigg{]}$ (192) $\displaystyle I_{6}(q)$
$\displaystyle=\mathcal{N}\bigg{[}{-\frac{3}{2\pi}\kappa+\frac{3}{2}\ln
2}+\epsilon\left(-3\ln^{2}2+\frac{3}{2\pi}c_{6}+\frac{3}{\pi}\kappa\ln
2\right)+\mathcal{O}(\epsilon^{2})\bigg{]}$ (193) $\displaystyle I_{7}(q)$
$\displaystyle=\mathcal{N}\bigg{[}{-\frac{3}{2\pi}\kappa+\frac{3}{2}\ln
2}+\epsilon\left(-3\ln^{2}2+\frac{3}{2\pi}c_{8}+\frac{3}{\pi}\kappa\ln
2\right)+\mathcal{O}(\epsilon^{2})\bigg{]}$ (194)

Here the normalization factor is

$\mathcal{N}=\pi^{-\epsilon}e^{-\gamma_{E}\epsilon}\left(\frac{2}{3}\right)^{2\epsilon}\frac{1}{q^{1+2\epsilon}}$
(195)

with

$c_{1}=\text{Im}\left[\text{Li}_{3}\left(\frac{i}{\sqrt{3}}\right)\right],\quad
c_{2}=\text{Im}\left[\text{Li}_{3}\left(1+i\sqrt{3}\right)\right],$ (196)

and

$\displaystyle c_{3}$ $\displaystyle=-0.949789416853385,\quad
c_{4}=-0.305492372030520,\quad c_{5}=-1.34784474998184,$ $\displaystyle c_{6}$
$\displaystyle=4.39528715012114,\quad c_{7}=12.8006547420728,\quad
c_{8}=3.37308464401608$ (197)

As a cross check, we can add the three sextants in one hemisphere with the
same axis projection and compare to the hemisphere soft function
Fleming:2007xt ; Becher:2009th ; Chien:2010kc ; Kelley:2011ng ; Hornig:2011iu
.

$I_{\bar{n}_{1},\bar{n}_{1},n_{1}}(q)=I_{1}(q)+I_{\bar{n}_{1},\bar{n}_{3},n_{1}}(q)+I_{\bar{n}_{1},\bar{n}_{2},n_{1}}(q)$
(198)

The extra two soft integrals are

$\displaystyle I_{\bar{n}_{1},\bar{n}_{3},n_{1}}(q)$
$\displaystyle=\mathcal{N}\bigg{[}{\frac{3}{2}\ln
2-\frac{3}{2\pi}\kappa}+\epsilon\bigg{(}\frac{71}{540}\pi^{2}+\frac{18}{5\pi}c_{1}-\frac{4}{\pi}c_{2}+\frac{3}{\pi}\kappa\ln
2-\frac{8}{3}\ln^{2}2$ $\displaystyle+\frac{3}{2}\ln 2\ln
3-\frac{3}{40}\ln^{2}3-\frac{13}{12}\text{Li}_{2}\left(\frac{1}{4}\right)\bigg{)}+\mathcal{O}(\epsilon^{2})\bigg{]}$
$\displaystyle I_{\bar{n}_{1},\bar{n}_{2},n_{1}}(q)$
$\displaystyle=\mathcal{N}\bigg{[}{\ln
2+\frac{3}{\pi}\kappa}+\epsilon\bigg{(}\frac{119}{270}\pi^{2}-\frac{36}{5\pi}c_{1}-\frac{4}{\pi}c_{2}-\frac{6}{\pi}\kappa\ln
2-\frac{2}{3}\ln^{2}2$ $\displaystyle+\ln 2\ln
3+\frac{3}{20}\ln^{2}3+\frac{1}{6}\text{Li}_{2}\left(\frac{1}{4}\right)\bigg{)}+\mathcal{O}(\epsilon^{2})\bigg{]}$
(199)

which leads to the same hemisphere soft function as in Ref. Becher:2009th .

Another interesting fact is that the divergent part of our trijet soft
function does not depend on the projection vector $N$. The soft function
integral is

$I_{n_{a},n_{b},n_{c},n_{d},N}(q)\sim\int
d^{d}k\frac{{n_{a}}\cdot{n_{b}}}{({n_{a}}\cdot k)({n_{b}}\cdot
k)}\delta(k^{2})\delta\left(q-\frac{2}{3}N\cdot
k\right)\times\left[\cdots\right]\\\ $ (200)

If we rotate the projection vector $N$ to another direction $N^{\prime}$, then
the $\delta$-function transforms as

$\delta\left(q-\frac{2}{3}N^{\prime}\cdot
k\right)=\delta\left(q-\frac{2}{3}(N\cdot k)\frac{N^{\prime}\cdot k}{N\cdot
k}\right)=\frac{N\cdot k}{N^{\prime}\cdot k}\delta\left(\frac{N\cdot
k}{N^{\prime}\cdot k}q-\frac{2}{3}N\cdot k\right)$ (201)

Then after rescaling

$k\to\frac{N\cdot q}{N^{\prime}\cdot q}k$ (202)

the integral becomes

$I_{n_{a},n_{b},n_{c},n_{d},N}(q)\sim\left(\frac{N\cdot q}{N^{\prime}\cdot
q}\right)^{2\epsilon}q^{d-5}\int
d^{d}k\frac{{n_{a}}\cdot{n_{b}}}{({n_{a}}\cdot k)({n_{b}}\cdot
k)}\delta(k^{2})\delta\left(q-\frac{2}{3}N\cdot
k\right)\times\left[\cdots\right]\\\ $ (203)

So the effect of using different $N$’s only shows up at order $\epsilon$ in
the expansion. This only affects the anomalous dimension for the integrals
which have soft-collinear divergences. This is only $I_{1}(q)$, since that is
the only integral where a Wilson line is in within the integration region.
However, for $I_{1}(q)$ the projection is on $n_{1}$ (the Wilson line
direction) in both thrust and heavy jet mass. For the others, using the
definitions in Eqs.(176), (177) and (179), we see that at NLL level after
rescaling the projection vector $N$

$I_{5}(q)=I_{3}(q),\quad I_{6}(q)=I_{4}(q),\quad\text{and}\quad
I_{7}(q)=I_{6}(q)$ (204)

where we need a reflection with respect to $z$-axis to see the third equation.
This agrees with our explicit calculations in Eqs.(A)-(194) to order
$\epsilon^{0}$. Note however, that the $\epsilon^{1}$ terms differ, as
expected. In summary, at the NLL level, the thrust and heavy jet mass trijet
soft function can be taken to be the same. For NNLL resummation and beyond,
they will generically be different.

## References

  * (1) S. Catani and B. R. Webber, “Infrared safe but infinite: Soft gluon divergences inside the physical region,” JHEP 10 (1997) 005, arXiv:hep-ph/9710333.
  * (2) M. H. Seymour, “Jet shapes in hadron collisions: Higher orders, resummation and hadronization,” Nucl. Phys. B 513 (1998) 269–300, arXiv:hep-ph/9707338.
  * (3) G. Luisoni, P. F. Monni, and G. P. Salam, “$C$-parameter hadronisation in the symmetric 3-jet limit and impact on $\alpha_{s}$ fits,” Eur. Phys. J. C 81 no. 2, (2021) 158, arXiv:2012.00622 [hep-ph].
  * (4) E. Farhi, “A QCD Test for Jets,” Phys. Rev. Lett. 39 (1977) 1587–1588.
  * (5) G. P. Salam and D. Wicke, “Hadron masses and power corrections to event shapes,” JHEP 05 (2001) 061, arXiv:hep-ph/0102343.
  * (6) Y.-T. Chien and M. D. Schwartz, “Resummation of heavy jet mass and comparison to LEP data,” JHEP 08 (2010) 058, arXiv:1005.1644 [hep-ph].
  * (7) A. Gehrmann-De Ridder, T. Gehrmann, E. W. N. Glover, and G. Heinrich, “NNLO corrections to event shapes in e+ e- annihilation,” JHEP 12 (2007) 094, arXiv:0711.4711 [hep-ph].
  * (8) S. Catani and M. H. Seymour, “The Dipole formalism for the calculation of QCD jet cross-sections at next-to-leading order,” Phys. Lett. B 378 (1996) 287–301, arXiv:hep-ph/9602277.
  * (9) S. Catani and M. H. Seymour, “A General algorithm for calculating jet cross-sections in NLO QCD,” Nucl. Phys. B 485 (1997) 291–419, arXiv:hep-ph/9605323. [Erratum: Nucl.Phys.B 510, 503–504 (1998)].
  * (10) S. Frixione, P. Nason, and G. Ridolfi, “Problems in the resummation of soft gluon effects in the transverse momentum distributions of massive vector bosons in hadronic collisions,” Nucl. Phys. B 542 (1999) 311–328, arXiv:hep-ph/9809367.
  * (11) T. Becher and M. Neubert, “Drell-Yan Production at Small $q_{T}$, Transverse Parton Distributions and the Collinear Anomaly,” Eur. Phys. J. C 71 (2011) 1665, arXiv:1007.4005 [hep-ph].
  * (12) P. F. Monni, E. Re, and P. Torrielli, “Higgs Transverse-Momentum Resummation in Direct Space,” Phys. Rev. Lett. 116 no. 24, (2016) 242001, arXiv:1604.02191 [hep-ph].
  * (13) P. Nogueira, “Automatic feynman graph generation,” Journal of Computational Physics 105 no. 2, (1993) 279–289. https://www.sciencedirect.com/science/article/pii/S0021999183710740.
  * (14) J. A. M. Vermaseren, “New features of FORM,” arXiv:math-ph/0010025.
  * (15) T. Hahn, “Generating feynman diagrams and amplitudes with feynarts 3,” Computer Physics Communications 140 no. 3, (2001) 418–431. https://www.sciencedirect.com/science/article/pii/S0010465501002909.
  * (16) R. Mertig, M. Böhm, and A. Denner, “Feyn calc - computer-algebraic calculation of feynman amplitudes,” Computer Physics Communications 64 no. 3, (1991) 345–359. https://www.sciencedirect.com/science/article/pii/001046559190130D.
  * (17) V. Shtabovenko, R. Mertig, and F. Orellana, “New developments in feyncalc 9.0,” Computer Physics Communications 207 (2016) 432–444. https://www.sciencedirect.com/science/article/pii/S0010465516301709.
  * (18) R. K. Ellis, W. J. Stirling, and B. R. Webber, QCD and collider physics, vol. 8.  Cambridge University Press, 2, 2011.
  * (19) T. Becher and M. D. Schwartz, “Direct photon production with effective field theory,” JHEP 02 (2010) 040, arXiv:0911.0681 [hep-ph].
  * (20) M. D. Schwartz, “Precision direct photon spectra at high energy and comparison to the 8 TeV ATLAS data,” JHEP 09 (2016) 005, arXiv:1606.02313 [hep-ph].
  * (21) T. Becher, C. Lorentzen, and M. D. Schwartz, “Resummation for W and Z production at large pT,” Phys. Rev. Lett. 108 (2012) 012001, arXiv:1106.4310 [hep-ph].
  * (22) T. Becher, C. Lorentzen, and M. D. Schwartz, “Precision Direct Photon and $W$-Boson Spectra at High $p_{T}$ and Comparison to LHC Data,” Phys. Rev. D 86 (2012) 054026, arXiv:1206.6115 [hep-ph].
  * (23) A. Banfi, G. Marchesini, and G. Smye, “Away from jet energy flow,” JHEP 08 (2002) 006, arXiv:hep-ph/0206076.
  * (24) M. D. Schwartz and H. X. Zhu, “Nonglobal logarithms at three loops, four loops, five loops, and beyond,” Phys. Rev. D 90 no. 6, (2014) 065004, arXiv:1403.4949 [hep-ph].
  * (25) T. Becher, M. Neubert, and B. D. Pecjak, “Factorization and Momentum-Space Resummation in Deep-Inelastic Scattering,” JHEP 01 (2007) 076, arXiv:hep-ph/0607228.
  * (26) M. D. Schwartz, “Resummation and NLO matching of event shapes with effective field theory,” Phys. Rev. D 77 (2008) 014026, arXiv:0709.2709 [hep-ph].
  * (27) T. Becher and M. D. Schwartz, “A precise determination of $\alpha_{s}$ from LEP thrust data using effective field theory,” JHEP 07 (2008) 034, arXiv:0803.0342 [hep-ph].
  * (28) S. Fleming, A. H. Hoang, S. Mantry, and I. W. Stewart, “Top Jets in the Peak Region: Factorization Analysis with NLL Resummation,” Phys. Rev. D 77 (2008) 114003, arXiv:0711.2079 [hep-ph].
  * (29) T. Becher, M. Neubert, and G. Xu, “Dynamical Threshold Enhancement and Resummation in Drell-Yan Production,” JHEP 07 (2008) 030, arXiv:0710.0680 [hep-ph].
  * (30) I. W. Stewart, F. J. Tackmann, and W. J. Waalewijn, “N-Jettiness: An Inclusive Event Shape to Veto Jets,” Phys. Rev. Lett. 105 (2010) 092002, arXiv:1004.2489 [hep-ph].
  * (31) T. T. Jouttenus, I. W. Stewart, F. J. Tackmann, and W. J. Waalewijn, “The Soft Function for Exclusive N-Jet Production at Hadron Colliders,” Phys. Rev. D 83 (2011) 114030, arXiv:1102.4344 [hep-ph].
  * (32) I. Feige, M. D. Schwartz, I. W. Stewart, and J. Thaler, “Precision Jet Substructure from Boosted Event Shapes,” Phys. Rev. Lett. 109 (2012) 092001, arXiv:1204.3898 [hep-ph].
  * (33) R. K. Ellis, D. A. Ross, and A. E. Terrano, “The Perturbative Calculation of Jet Structure in e+ e- Annihilation,” Nucl. Phys. B 178 (1981) 421–456.
  * (34) T. Becher and M. Neubert, “On the Structure of Infrared Singularities of Gauge-Theory Amplitudes,” JHEP 06 (2009) 081, arXiv:0903.1126 [hep-ph]. [Erratum: JHEP 11, 024 (2013)].
  * (35) A. H. Hoang and S. Kluth, “Hemisphere Soft Function at O(alpha(s)**2) for Dijet Production in e+ e- Annihilation,” arXiv:0806.3852 [hep-ph].
  * (36) M. Dasgupta and G. P. Salam, “Resummation of nonglobal QCD observables,” Phys. Lett. B 512 (2001) 323–330, arXiv:hep-ph/0104277.
  * (37) T. Becher, M. Neubert, L. Rothen, and D. Y. Shao, “Effective Field Theory for Jet Processes,” Phys. Rev. Lett. 116 no. 19, (2016) 192001, arXiv:1508.06645 [hep-ph].
  * (38) A. J. Larkoski, I. Moult, and D. Neill, “Non-Global Logarithms, Factorization, and the Soft Substructure of Jets,” JHEP 09 (2015) 143, arXiv:1501.04596 [hep-ph].
  * (39) R. Abbate, M. Fickinger, A. H. Hoang, V. Mateu, and I. W. Stewart, “Thrust at $N^{3}LL$ with Power Corrections and a Precision Global Fit for $\alpha_{s}(mZ)$,” Phys. Rev. D 83 (2011) 074021, arXiv:1006.3080 [hep-ph].
  * (40) S. Catani, M. L. Mangano, P. Nason, and L. Trentadue, “The Resummation of soft gluons in hadronic collisions,” Nucl. Phys. B 478 (1996) 273–310, arXiv:hep-ph/9604351.
  * (41) A. Banfi, B. K. El-Menoufi, and P. F. Monni, “The Sudakov radiator for jet observables and the soft physical coupling,” JHEP 01 (2019) 083, arXiv:1807.11487 [hep-ph].
  * (42) P. Cal, F. Ringer, and W. J. Waalewijn, “The jet shape at NLL’,” JHEP 05 (2019) 143, arXiv:1901.06389 [hep-ph].
  * (43) S. Catani, L. Trentadue, G. Turnock, and B. R. Webber, “Resummation of large logarithms in e+ e- event shape distributions,” Nucl. Phys. B 407 (1993) 3–42.
  * (44) T. Becher, M. Neubert, and D. Wilhelm, “Electroweak Gauge-Boson Production at Small $q_{T}$: Infrared Safety from the Collinear Anomaly,” JHEP 02 (2012) 124, arXiv:1109.6027 [hep-ph].
  * (45) M. A. Ebert and F. J. Tackmann, “Resummation of Transverse Momentum Distributions in Distribution Space,” JHEP 02 (2017) 110, arXiv:1611.08610 [hep-ph].
  * (46) H. Chen, I. Moult, X. Zhang, and H. X. Zhu, “Rethinking jets with energy correlators: Tracks, resummation, and analytic continuation,” Phys. Rev. D 102 no. 5, (2020) 054012, arXiv:2004.11381 [hep-ph].
  * (47) R. Kelley, M. D. Schwartz, R. M. Schabinger, and H. X. Zhu, “The two-loop hemisphere soft function,” Phys. Rev. D 84 (2011) 045022, arXiv:1105.3676 [hep-ph].
  * (48) A. Hornig, C. Lee, I. W. Stewart, J. R. Walsh, and S. Zuberi, “Non-global Structure of the $\mathcal{O}({\alpha}^{2}_{s})$ Dijet Soft Function,” JHEP 08 (2011) 054, arXiv:1105.4628 [hep-ph]. [Erratum: JHEP 10, 101 (2017)].


# Online control of the gain drift with temperature of SiPM arrays used for
the readout of LaBr3:Ce crystals

M.Bonesini,11footnotetext: Corresponding author. R.Bertoni  M.Prata  and
M.Rossella

###### Abstract

LaBr3:Ce crystals have been introduced for radiation imaging in medical
physics, with photomultiplier or single SiPM readout. An R$\&$D was pursued
with 1/2" and 1" LaBr3:Ce crystals, from different producers, to realize
compact large area detectors (up to some cm2 area) with SiPM array readout,
aiming at high light yields, good energy resolution, good detector linearity
and fast time response for low-energy X-rays. A natural application was found
inside the FAMU project at RIKEN-RAL muon facility, that aims at a precise
measure of the proton Zemach radius to solve the so-called "proton radius
puzzle", triggered by the recent measure of the proton charge radius at PSI.
The goal is the detection of characteristic X-rays around 130 keV. Other
applications may be foreseen in medical physics, such as PET, and gamma-ray
astronomy. A limiting factor is the gain drift of SiPM arrays with
temperature, that give a deterioration of the detector’s FWHM energy
resolution. To solve this problem, a custom NIM module, based on CAEN A7585
digital power supply, was developed. Test results of the correction of gain
drift with temperature for SiPM arrays from Advansid, Sensl, Hamamatsu will be
presented. At the137 Cs peak, an energy resolution better than $3\%$ was
obtained for a typical LaBr3:Ce crystal, using Hamamatsu S13461 arrays. This
compares well with best available results obtained with a PMT.

##  1 Introduction

Silicon photomultipliers (SiPM) are a valuable alternative to conventional
photomultipliers (PMTs) for the readout of scintillation detectors. They may
be arranged as an array. As readout devices they have a high reliability, a
low sensitivity to external magnetic fields and can operate at voltages
significantly lower than the ones used for PMTs. Using a SiPM array it is
possible to obtain energy resolutions comparable to what obtained with PMTs.
However, SIPMs have a relevant problem: in addition to an increased noise
level, their gain drifts significantly as a function of temperature. This
feature prevents their use in conditions with a changing temperature, as
homeland security and military applications. To overcome this problem, SiPMs
must be operated at a fixed gain if temperature changes. Thus the operating
voltage $V_{op}=V_{bd}+V_{ov}$, with $V_{bd}$ breakdown voltage and $V_{ov}$
overvoltage, must be modified as a function of temperature according to:

$V_{bd}(T)=V_{bd}(T_{0})\times(1+\beta(T-T_{0})$

with T working temperature, $T_{0}$ reference temperature (typically 25∘ C)
and $\beta=\Delta V_{bd}/\Delta T$ temperature coefficient of the used SiPM,
as explained in reference [1]. For an online hardware correction, a custom NIM
module was developed, based on CAEN A7585D digital power supplies, with
temperature feedback. Up to eight channels may be powered by a single 2-slots
NIM module, as shown in figure 1.

Figure 1: Left panel - top, A7585D functional layout (courtesy of CAEN srl);
bottom, exploded view of one custom NIM module. Right panel: front panel of
one custom NIM modules: a) is the BNC connector for a channel HV, b) the 3.5
mm jack stereo for the connection to a temperature TMP37 sensor, c) the USB
interface connection.

The control of this module was implemented via either a FDTI USB-I2C converter
or an Arduino Nano chip. Three modules (with which up to 24 channels may be
powered) were realized and are linkable in daisy chain, via the I2C protocol.
Our application has developed inside the FAMU experiment at RAL, to measure
with high precision the proton Zemach radius [2]: its main requirement was the
measure of X-rays around 100 keV with compact crystal detectors.

##  2 Laboratory tests

The crystal detectors are mounted inside an ABS holder, realized with a 3D
printer and are fully described in reference [3]. The SiPM array temperature
is measured by a TMP37 temperature sensor, mounted on the PCB where the array
socket is soldered. A 3.5 mm stereo jack cable connects the sensor to the
custom NIM module for online temperature correction. Laboratory tests were
done putting the detectors under test inside a Memmert IPV-30 climatic
chamber, where the temperature could be stabilized with a precision of $\sim$
0.1 ∘C. Detectors were powered at their nominal operating voltage $V_{op}$.
The summed analogue signal from the cells of a SiPM array is directly fed into
a CAEN V1730 fast digitizer (500 MHz bandwidth, 14 bit resolution) and is
acquired by a custom DAQ system. In the reported laboratory tests, a 137Cs
exempt source was used and the dependence of the photo peak at 662 keV, as a
function of temperature is reported in figure 2. The same crystal is used in
all the tests. While Advansid and Hamamatsu SiPM arrays show a larger
variation of response (up to 70 $\%$), the effect is smaller for SENSL devices
($\sim 30\%$). After correction, in all cases the variation is reduced to
$\sim 5\%$.

Figure 2: Dependence of the photo peak position at 662 keV for a typical 1/2"
detector, with and without temperature correction. Different SiPM arrays have
been used for the readout.

The response of a typical 1/2" crystal, inside the climatic chamber, with a
temperature excursion from 20 ∘C to 30 ∘C, is shown in the left panel of
figure 3. Without the online correction, the resolution of the photo peak at
662 keV is sensibly degraded. The variation of FWHM energy resolution of the
same crystal, with different readout, is shown instead in the right panel.

Figure 3: Left panel: Cs137 spectra recorded by a LaBr3:Ce 1/2" detector read
by an Advansid ASD-NUV3S SiPM array during a temperature scan between 20 ∘C
and 30 ∘C. Top (bottom) is without (with) temperature correction. Fit results
to the 662 keV photo peak, with a simple gaussian, are reported. Right panel:
dependence of the energy resolution at 662 keV on temperature, for a typical
1/2" LaBr${}_{3}:$Ce detector read by different SiPM arrays, with/without
temperature correction.

##  3 Conclusions

A custom NIM module, based on CAEN A7585D digital power supply, has been
developed to allow an effective online correction of the effect of temperature
on the gain of SiPM arrays. The effect, depending on the type of SiPM array,
is reduced from $30-70\%$ to $\sim 5\%$.

## References

  * [1] N. Dinu et al., _Temperature and bias voltage dependence of the MPPC detectors_ ,IEEE NSS-MIC Symposium (2010), Piscateway, NJ.
  * [2] A. Vacchi et al., _Measuring the size of the proton_ , SPIE Newsroom (2012), DOI:10.1117/2/1201207.004274;
C. Pizzolotto et al., _The FAMU experiment: muonic hydrogen high precision
spectroscopy studies_ , EPJ A56 (2020) 7,185 ;
M. Bonesini, _The FAMU experiment at RIKEN RAL for a precise measurement of
the proton radius_ , PoS EPS-HEP2019 132.

  * [3] M. Bonesini et al., _Ce:LaBr 3 crystals with SiPM array readout and temperature control for the FAMU experiment at RAL_, JINST 15 (2020) 05, C05065.