Towards a sensing model using random laser combined with diffuse reflectance
spectroscopy

Dongqin Ni1,2,* , Florian Klämpfl1,2, Michael Schmidt1,2, Martin Hohmann1,2

1 Institute of Photonic Technologies (LPT), Friedrich-Alexander-Universität
Erlangen-Nürnberg (FAU), Konrad-Zuse-Straße 3/5, 91052 Erlangen, Germany
2 Erlangen Graduate School in Advanced Optical Technologies (SAOT), Paul-
Gordan-Straße 6, 91052 Erlangen, Germany

*<EMAIL_ADDRESS>
## Abstract

The previous research proves that the random laser emission reflects not only
the scattering properties but also the absorption properties. The random laser
is therefore considered a potential tool for optical properties sensing.
Although the qualitative sensing using the random laser is extensively
investigated, a quantitative measurement is still rare. In this study, a
generalized mathematical quantitative model using random laser combined with
diffuse reflectance spectroscopy is proposed for optical sensing in turbid
media. This model describes the gain effect of the active medium and the
optical properties effect of the passive medium separately. Rhodamine 6G is
used as the active medium. Intralipid and ink are employed to demonstrate the
effect of the scattering and absorption, respectively. The peak wavelength
shift of the random laser is proved to be an ideal sensing parameter for this
sensing model. It is also revealed that the scaling parameters in the sensing
model are interrelated and can be simplified to one. With this combined model,
the direct sensing of optical properties in diverse turbid media is promising.

##  1 Introduction

The understanding of light-matter interactions and the characterization of
optical properties (OPs) are of paramount importance in various fields of
research [1, 2]. In biological and biomedical fields, a profound understanding
of OPs is essential for the interpretation of diagnostic and therapeutic
measurements, as well as for the development of medical technologies [2]. The
characterization of OPs in turbid media, such as biological tissues, typically
involves parameters such as absorption coefficient $\mu_{a}$, reduced
scattering coefficient $\mu_{s}^{\prime}$, and anisotropy factor $g$ [2].
While absorption and scattering strength are described by $\mu_{a}$ and
$\mu_{s}^{\prime}$, respectively, $g$ denotes the scattering angle dependence
[2]. The relation between the scattering parameters of $\mu_{s}^{\prime}$ and
$g$ is expressed by the equation:

$\mu_{s}^{\prime}=\mu_{s}(1-g),$ (1)

where $\mu_{s}$ represents the scattering coefficient.

The use of a random laser (RL) is emerging as a promising approach for the
direct characterization of OPs in turbid media [3]. Within a RL system, the
optical feedback for lasing stems from multiple scattering among random
scatterers [4]. Consequently, the emission properties of RLs reflect the
scattering properties of the random media, offering potential applicability in
optical sensing in turbid media, particularly in biological tissues [4]. For
example, RLs have demonstrated success in tissue differentiation [5],
cancerous tissue screening [6, 7], and the detection of optomechanical strain
in tissues [8]. These achievements are noteworthy, although the aforementioned
sensing applications are found solely on qualitative measurements, thereby
limiting their generalizability to diverse turbid systems. There is a need for
a quantitative model using RL emission properties to assess the universal OPs
in various turbid media.

Researchers have attempted to formulate a quantitative OPs sensing model based
on RL emission [9, 3]. The sensing parameters include the RL emission
intensity, peak wavelength, linewidth, lasing modes and lasing threshold [4].
In particular, the work of Tommasi et al. [9] demonstrated the measurement of
RL emission intensity to characterize the constant value of $\mu_{s}^{\prime}$
arising from microspheres with different diameters. They separated the active
gain medium from the passive scattering medium by using an isolated
transparent spherical cell. Both pump and emission light were guided through
one fiber. In the work from Hohmann et al. [3], the active gain medium was
mixed with the scattering medium. They reported the RL dependence on changes
in $\mu_{s}$, not only from the RL intensity, but also from the spectral peak
wavelength and linewidth. The responses to changes in $\mu_{s}$ are consistent
for all three parameters. Furthermore, it was empirically found that critical
alterations in the RL emission variations are induced by $\mu_{s}$ rather than
$\mu_{s}^{\prime}$. Specifically, the RL emission exhibited maximum intensity
when the laser cavity length was an integer multiple of the scattering mean
free path $l_{s}$ (equivalent to $1/\mu_{s}$), rather than the reduced
scattering mean free path $l_{s}^{\prime}$ (equivalent to
$1/\mu_{s}^{\prime}$) [3]. Consequently, direct measurement employing RLs may
provide a solution to the precise characterization of $\mu_{s}$. The question
remains in the optical characterization of $\mu_{a}$. It was explored in a
previous study [10], wherein an increase in $\mu_{a}$ resulted in a decrease
in RL intensity, a broadening of the linewidth, a blueshift of the peak
wavelength, and an increase in the lasing threshold [10]. It has been
demonstrated that the RL emission exhibits an inverse behavior on $\mu_{a}$ in
comparison to $\mu_{s}$ or $\mu_{s}^{\prime}$ [10]. The above three studies
showed the feasibility of the RL emission to sense the OPs of
$\mu_{s}^{\prime}$, $\mu_{s}$ and $\mu_{a}$. However, a mathematical
representation of the sensing model is still lacking.

In fact, the aforementioned RL behavior is analogous to the well-known diffuse
reflectance: diffuse reflectance intensity is directly proportional to
$\mu_{s}^{\prime}$ and inversely proportional to $\mu_{a}$ [11]. The
phenomenon that the RL responds similarly to the diffuse reflectance in the
context of OPs is not fortuitous. Pioneering RL research has evidenced that RL
generated from a turbid medium with scattering strength in the diffuse regime
can be modeled as light diffusion with gain [12]. In this study, when observed
in the backscattered direction, the RL can be hypothesized to be light diffuse
reflectance with gain. Since the modeling of the diffuse reflectance
spectroscopy (DRS) signal is well established in literature [11, 13], the RL
emission has high potential to be also modeled and applied for optical
sensing.

In addition, the utilization of RL emission for sensing compared to the
diffuse reflectance alone may facilitate deeper detection, since the typical
sensing depth using diffuse reflectance is limited, for example, 0.5 - 1.9 mm
in depth using a light source operating in the wavelength range of 350 - 1919
nm and applying a source-detector separation of 2500 µm [14]. This limitation
will be overcome in RLs by cascaded gain amplification along the light path,
leading to an extended active light path length in turbid media [4].
Furthermore, this effect may improve the sensitivity of RL-based OPs sensing
due to the enhanced laser emission compared to the lamp typically used in the
diffuse reflectance measurements. Therefore, the prospect of constructing a
mathematical RL-based OPs sensing model, adapted from the diffuse reflectance
model, with a higher sensitivity of a deeper sensing target, seems very
promising.

The aim of this study is therefore to use the RL spectral properties to
construct a mathematical DRS-RL model for quantitative OPs sensing in turbid
media. In particular, for the first time, this RL model is able to separately
describe the gain effect in addition to the OPs effect. This model can be
adapted to various turbid media, but in this study it is validated only in the
simplest turbid system consisting of Intralipid (IL) as scatterer and black
ink as absorber.

##  2 Theoretical consideration

###  2.1 DRS model

The diffuse reflectance is indeed a function of $\mu_{s}^{\prime}$, $\mu_{a}$,
phase function $p(\theta)$ (or $g$ factor in the multiple scattering regime in
this study) and geometry $G$ [15]:

$R_{d}=f(\mu_{s}^{\prime},\mu_{a},g,G).$ (2)

The phase function $p(\theta)$ describes the scattering angle dependence in a
single scattering event, analogous to the function of the $g$ factor in the
multiple scattering condition. One study demonstrated the non-negligible
influence of the phase function $p(\theta)$ on the diffuse reflectance in
certain probe geometry configurations, e.g. when the light collection is close
to the source [16]. The probe geometry $G$ includes the effect of the source-
detector separation and the effect of the numerical aperture of the fiber
probe [17, 18]. The symbol $f()$ denotes a nonspecific function for the
argument in bracket. The same function notation is applied in the following
equations in this study.

Among various applied mathematical models [17], the Zonios’s model of diffuse
reflectance in semi-infinite turbid media with fiber probes is the most
practical one [17]. The diffuse reflectance $R_{d}$ in the Zonios’s model
reads:

$R_{d}=\frac{\mu_{s}^{\prime}}{k_{1}+k_{2}\mu_{a}},$ (3)

where $k_{1}$ and $k_{2}$ are the scaling parameters depending on the optical
probe geometry $G$ [17]. The phase function was not formulated in this model.

Nevertheless, Equation 2 and 3 indicate that a direct measurement of the OPs
of $\mu_{s}^{\prime}$, $\mu_{s}$ and $\mu_{a}$ is promising ($\mu_{s}$ can be
calculated with known $\mu_{s}^{\prime}$ and $g$ according to the Equation 1).
The difficulty lies in the measurement of $g$. Although $g$ is a relevant
parameter of OPs, a mathematical representation that separates its effect from
the scaling parameters $k_{1}$ and $k_{2}$ in Equation 3 is still lacking in
the DRS literature.

###  2.2 DRS-RL sensing principle: sensing parameter and sensing curve

The emission from the RL sample is either non-lasing or lasing emission. Only
the RL signal in the lasing regime is utilized for sensing, because the RL
emission in the lasing regime is more stable and stronger due to the dominant
stimulated emission. To separate the non-lasing and lasing emissions, the peak
wavelength of the RL is employed. As sketched in Figure 1 (a), the peak
wavelength $\lambda_{p}$ always changes from blueshift to redshift at the
lasing threshold when the pump energy $E$ is increased. This changed trend of
the peak wavelength shift at the lasing threshold is universal and independent
of the OPs of the samples [10]. Moreover, this non-monotonic change of the
peak wavelength makes it a better lasing indicator than the other RL
parameters which only show monotonic changes at the lasing threshold [10].
Hence, the RL peak wavelength is chosen as the sensing parameter of the model.

Figure 1: Sensing principle: (a) raw data showing the response of the RL peak
wavelength $\lambda_{p}$ on the pump energy $E$. The lasing threshold is
determined upon the peak wavelength changes from blueshift to redshift.
Different colors represent the samples with different OPs of
$\mu_{s}^{\prime}$ and $\mu_{a}$. (b) pre-processed data to separate the non-
lasing and lasing regime. $\lambda_{p}$ and $E$ are replaced by the relative
peak wavelength $\Delta\lambda_{p}$ and the relative pump energy $\Delta E$,
both taking the values at the lasing threshold as a reference. Only the data
representing the lasing emission on the non-negative axes are further utilized
for sensing. (c) converged sensing curve. The $\Delta\lambda_{p}$ of all
samples after the correction of the $R_{d}$ effect, i.e.
$\Delta\lambda_{p}/R_{d}$, responds identically to $\Delta E$, due to the fact
that the remaining gain effect is identical in all samples.

In practice, to separate the non-lasing and lasing regimes, the peak
wavelength and pump energy are replaced by a relative peak wavelength
$\Delta\lambda_{p}=\lambda-\lambda_{th}$ and a relative pump energy $\Delta
E=E-E_{th}$, both of which take the values at the lasing threshold as a
reference. In other words, both $\Delta\lambda_{p}$ and $\Delta E$ are zero at
the lasing threshold, as shown in Figure 1 (b). Only the lasing regime where
the $\Delta\lambda_{p}$ and $\Delta E$ are non-negative is used for the
further sensing application. Meanwhile, the nonlinearity of the lasing effect
at the threshold is also excluded by using $\Delta\lambda_{p}$ instead of
$\lambda_{p}$ as the sensing parameter, which makes the formulation of the
gain effect much easier.

The hypothesis that the RL emission is the light diffuse reflectance with gain
is then mathematically expressed as:

$\Delta\lambda_{p}=f(R_{d},\Delta E).$ (4)

Since the response of the peak wavelength to the pump energy is OPs or diffuse
reflectance independent, the gain effect induced by $\Delta E$, i.e. $f(\Delta
E)$ can be formulated separately:

$\Delta\lambda_{p}=R_{d}\cdot f(\Delta E).$ (5)

When the effect of $R_{d}$ is excluded from $\Delta\lambda_{p}$ by division,
it is assumed that the $\Delta\lambda_{p}/R_{d}$ depends only on the $\Delta
E$, i.e. the gain effect. Mathematically, this relationship can be expressed
as:

$\Delta\lambda_{p}/R_{d}=f(\Delta E)$ (6)

and by applying the Zonios’s model from Equation 3:

$\Delta\lambda_{p}\cdot\frac{k_{1}+k_{2}\mu_{a}}{\mu_{s}^{\prime}}=f(\Delta
E).$ (7)

The right side of the Equation 6 and 7 can be interpreted as the ”OPs-
independent gain effect” that can be measured when the sample has no scatterer
and no absorber, but only a gain medium, i.e. a transparent active medium. In
this study, the measurement of the ”OPs-independent gain effect” was performed
on the active medium of the Rhodamine 6G (R6G) water solution, as shown in
Figure 2.

Figure 2: Response of peak wavelength to the pump energy applied to a R6G gain
medium.

Since the gain effect of the same gain medium is identical, the response of
$\Delta\lambda_{p}/R_{d}$ to $\Delta E$ is also identical for each sample.
This identical response is sketched in Figure 1 (c), where all RL peak
wavelength after correcting for the OPs effect converge to an identical curve.
This curve is referred to the ”OPs-independent sensing curve” in this study.

Since the ”OPs-independent gain effect” shows an exponential response in
Figure 2, the fit function of the ”OPs-independent sensing curve” is
predefined as an exponential function too:

$\Delta\lambda_{p}\cdot\frac{k_{1init}+k_{2init}\mu_{a}}{\mu_{s}^{\prime}}=-a\cdot
e^{-b\cdot\Delta E}+c,$ (8)

where $a$, $b$, $c$ are the fitting parameters, and $k_{1init}$, $k_{2init}$
are the initial assumptions of the scaling parameters of $k_{1}$, $k_{2}$.

In a short summary, the RL emission detected from the backscattered direction
is assumed to be the diffuse reflectance with gain. The RL peak wavelength
shift after correcting for the diffuse reflectance effect is assumed to be
identical due to the identical gain medium. Such a response can be expressed
mathematically as Equation 8 and illustrated graphically as the identical
sensing curve in Figure 1 (c). Given the sensing curve, the measured
$\Delta\lambda_{p}$ and $\Delta E$, and the optimized scaling parameters
$k_{1}$ and $k_{2}$, the OPs of $\mu_{s}^{\prime}$ and $\mu_{a}$ can be
derived from Equation 8.

##  3 Experimental validation

###  3.1 RL samples and experimental setup

The RL sample prepared in this study is a water-based liquid medium,
consisting of Rhodamine 6G (R6G; Sigma Aldrich, Germany) as the laser gain and
Intralipid (IL; Fresenius Kabi, Germany) as the scatterer. Indian black ink
(Royal Talens, The Netherland) was added as an external absorber. While the IL
and ink concentrations were varied to change the scattering and absorption
strength of the turbid media, the R6G concentration was kept at 2×10-4 g/ml to
achieve the optimal gain efficiency.

The method for calculating the OPs of scatters and absorbers was described in
the previous study [10, 19]. Two concentrations of IL scatterers of 5% and 9%
v/v were selected, leading to $\mu_{s}^{\prime}$ of 63.78 and 106.13 cm-1,
respectively. Since $\mu_{s}^{\prime}$ is wavelength dependent, the above
values took the average of the individual values calculated at different peak
wavelengths of the RL. The ink concentrations ranged from 0 to 0.09% v/v,
corresponding to $\mu_{a}$ from 0 to 8.64 cm-1. The OPs values of
$\mu_{s}^{\prime}$ and $\mu_{a}$ were selected according to the OPs of
biological tissues [2]. Meanwhile, the values satisfy the diffusion
approximation where $\mu_{s}^{\prime}$ is much larger than $\mu_{a}$.

Figure 3: RL experimental setup from our previous study [10]. The pump light
from a pulsed Nd:YAG laser is focused onto the sample. The back-scattered
emission light is collected from the excitation spot with an angle of 45∘, and
then guided by the optical fiber to a spectrograph. M: mirror; L: lens; NF:
notch filter; NDF: neutral density filter.

The experimental setup is shown in Figure 3 [10]. The liquid sample with a
volume of 75 mL was filled into a beaker, resulting in a cylindrical shape
with a diameter of 50 mm and a height of 38 mm. The pump light from a pulsed
Nd:YAG laser (Q-smart 450, Quantel) was focused onto the sample, and the
backscattered light was collected and directed to a spectrograph (Mechelle
Me5000 Echelle, Andor). For each measurement, a laser pulse with a wavelength
of 532 nm and a pulse duration of 5 ns was generated, and was focused on the
sample surface with a spot diameter of 0.25 mm. The emission light from the
same position of the excitation spot was detected with an angle of 45∘. The
spectrograph has a spectral range from 200 to 975 nm with a spectral resolving
power ($\lambda$/$\Delta\lambda$) of 6,000, i.e. a spectral resolution of 0.1
nm at the wavelength of 600 nm. For each sample, the laser pump energy $E$ was
varied from 1.82 mJ to 40.86 mJ in a trend of half-Gaussian distribution to
probe both the non-lasing and lasing regime of the RL emission. Each
measurement was repeated five times, and the averaged values were used.

###  3.2 Optimization of scaling parameters $k_{1}$ and $k_{2}$

The optimization of the scaling parameters $k_{1}$ and $k_{2}$ is estimated by
the convergence of the data points to the identical sensing curve. To assess
the convergence, the $R^{2}$ value, which estimates the goodness of the fit,
was applied as the optimization function. In practice, since the minimization
algorithms are more commonly available than maximization algorithms, the
minimum value of $1-R^{2}$ was evaluated instead of the maximum value of
$R^{2}$. Three scales were chosen as the initial assumption for $k_{1}$ and
$k_{2}$: 10, 1 and 0.1. These scales are comparable or smaller than the values
of $\mu_{s}^{\prime}$ and $\mu_{a}$ in Equation 8, so that the OPs rather than
the scaling parameters of $k_{1}$ and $k_{2}$ dominate the diffuse
reflectance. The tested values of $k_{1}$ and $k_{2}$ that lead to a minimum
value of the objective function of $1-R^{2}$ were returned as the optimal
values.

Two different optimization methods: direct search algorithm Nelder-Mead and
Markov chain Monte Carlo (MCMC) algorithm Metropolis-Hastings were applied to
avoid the method-induced optimization bias. The former one was implemented by
using the Python package of the ”scipy.optimize.minimize” with ”Nelder-Mead”
as the optimizer, $1-R^{2}$ as the objective function and $k_{1init}$,
$k_{2init}$ as the initial assumption. In the latter one, a Gaussian
distribution near the initial assumption $k$ was proposed ($k$ includes
$k_{1}$ and $k_{2}$). A random point $k^{\prime}$ from the proposed
distribution was selected, and accepted if $p(k)-p(k^{\prime})\geq\alpha$,
where $p()$ denotes the objective function of $1-R^{2}$ and $\alpha$ denotes a
uniform random number between 0 and 1.

##  4 Results

###  4.1 Data preprocessing

The response of the RL peak wavelength $\lambda_{p}$ on the pump energy $E$ is
shown in Figure 4 (a). For samples with the same IL concentration, increasing
ink concentration induces the higher lasing threshold $E_{th}$ as well as the
blueshift of the peak wavelength $\lambda_{p}$ in the lasing regime. For
samples with the same ink concentration, increasing the IL concentration leads
to an inverse response: a lower lasing threshold and a redshift of the peak
wavelength in the lasing regime.

Figure 4: Data preprocessing for the DRS-RL sensing. (a) response of the RL
peak wavelength to the pump energy that is applied to the samples. (b)
preprocessed data in the lasing regime. Different colors represent the samples
with different OPs. (c) data cloud after the correction of the OPs effect. The
data points are converged when dividing by $R_{d}$, the value of which is
calculated taking the arbitrary numbers of 8 and 2 for $k_{1}$ and $k_{2}$.

Figure 4 (b) shows the preprocessed data in the lasing regime. The peak
wavelength of each RL sample with different OPs increases with increasing pump
energy. The identical exponential response proves that the gain effect is
identical and independent of the OPs. In Figure 4 (c), after correcting for
the OPs effect by dividing the diffuse reflectance $R_{d}$, the data points
representing the peak wavelength shift are remarkably converged. It is
noteworthy that the $k_{1}$ and $k_{2}$ in $R_{d}$ take the arbitrary numbers
of 8 and 2 to show the convergence. The convergence can be further improved by
adjusting the scaling parameters of $k_{1}$ and $k_{2}$.

###  4.2 Model parameter optimization: $k_{1}$ and $k_{2}$

Table 4.2 summarizes the optimization results of the two optimizers applying
three different scales of initial assumptions of $k_{1}$ and $k_{2}$. Although
the returned optimal values of $k_{1}$ and $k_{2}$ are different, the values
of the fitting goodness $R^{2}$ are the same (the $R^{2}$ values differ only
in the 9th digit after the decimal). The bias caused by the selection of the
optimizers or initial assumptions is thus eliminated. Moreover, the optimal
convergence of the data points to the sensing curve is 94%, which is
independent of the optimizer and initial assumption.

Table 1: Optimization of the scaling parameters of $k_{1}$ and $k_{2}$.

—c—cc—cc—cc— gray!503,5,7 white!1001 $k_{1init}$, $k_{2init}$ 10, 10  1, 1
0.1, 0.1
Optimal $k_{1}$, $k_{2}$$R^{2}$ Optimal $k_{1}$, $k_{2}$$R^{2}$ Optimal
$k_{1}$, $k_{2}$$R^{2}$
Nelder-Mead 4.22, 1.520.94 0.43, 0.150.94 0.04, 0.020.94
Metropolis Hastings 8.17, 2.940.94 0.54, 0.190.94 0.14, 0.050.94

Figure 5: Linearity of the $k_{1}$ and $k_{2}$ optimal values derived by the
optimizer (a) Nelder-Mead and (b) Metropolis Hastings.

Furthermore, the relationship between the optimal $k_{1}$ and $k_{2}$ derived
from both optimizers is illustrated in Figure 5 for a comparison. The results
from both cases reveal the same linearity between $k_{1}$ and $k_{2}$:

$k_{2}=0.36\cdot k_{1}.$ (9)

###  4.3 DRS-RL sensing curve

Figure 6 shows the DRS-RL sensing curve after the normalization. The sensing
data points fit well to the exponential sensing curve with $R^{2}\approx
0.94$. The outlines located outside of the 95% prediction interval may result
from measurement errors.

Figure 6: The normalized DRS-RL sensing curve.

The sensing curve is formulated by applying the optimal values of $k_{1}$ and
$k_{2}$ in Table 4.2 to the fitting function in Equation 8. Although the
optimal values of $k_{1}$ and $k_{2}$ vary widely using different optimizers
and different initial assumptions, a unique fitting function is returned after
normalizing $\Delta\lambda_{p}/R_{d}$ to its maximum, as shown in Figure 6.
Taking into account the Equation 9, the sensing curve is expressed as

$\Delta\lambda_{p}\cdot\frac{k_{1}(1+0.36\mu_{a})}{\mu_{s}^{\prime}}|_{norm}=1.13\cdot(1-e^{-0.034\cdot\Delta
E}).$ (10)

A general scaling factor $k$ is further introduced to simplify the above
equation:

$k\cdot\frac{\Delta\lambda_{p}\cdot(1+0.36\mu_{a})}{\mu_{s}^{\prime}}=1-e^{-0.034\cdot\Delta
E},$ (11)

$k=9.59$ applying the measured ($\Delta\lambda_{p}$, $\Delta E$), the known
($\mu_{s}^{\prime}$, $\mu_{a}$) values and the fitting parameters, and taking
the average of the $k$ values calculated for each data point. Therefore, the
mathematical expression of the sensing curve is finally formulated as:

$9.59\cdot\frac{\Delta\lambda_{p}\cdot(1+0.36\mu_{a})}{\mu_{s}^{\prime}}=1-e^{-0.034\cdot\Delta
E}.$ (12)

From Equation 12, the RL sensing parameter of $\Delta\lambda_{p}$ is derived:

$\Delta\lambda_{p}=\frac{1}{9.59}\cdot\frac{\mu_{s}^{\prime}}{1+0.36\mu_{a}}\cdot(1-e^{-0.034\cdot\Delta
E}).$ (13)

This equation indicates that the RL peak wavelength shift in the lasing regime
can be represented by three combined components: (1) a general scaling factor,
(2) the OPs and (3) the gain effect induced by the pump energy. Furthermore,
the calculated values of $\Delta\lambda_{p}$ from Equation 13 are approximate
to the original values measured experimentally, with an average $R^{2}=0.9323$
as shown in Figure 7. The deviation is larger for samples with lower values of
$mu_{a}/\mu_{s}^{\prime}$. One possible cause is the additional gain
saturation as the RL medium has lower absorption loss and stronger feedback
caused by the stronger scattering strength. Despite this deviation, the
sensing curve and equation are representative for the experimental
measurements and can be applied for optical sensing.

Figure 7: The values of $\Delta\lambda_{p}$ calculated from the sensing curve
or from the Equation 13 (represented by the curves), compared with the
experimentally measured values of $\Delta\lambda_{p}$ (represented by the
dots).

###  4.4 Discussion and conclusion

Apart from the gain effect, RL in the backscattering direction responds
similarly as the diffuse reflectance changes on optical properties. This study
proposes a mathematical model of optical sensing that adapts the features of
RLs to the common DRS model, the Zonios’ model. A linear relationship between
the two scaling parameters in the Zonios’ model was revealed, and therefore
the two scaling parameters can be simplified to one general scaling factor. In
addition, it was shown that the pump energy induced gain effect in RL is
universal and independent of the OPs of the samples investigated. The RL gain
effect follows an exponential response of the transparent gain medium to the
pump energy. For the first time, the OPs effect and the gain effect of the RL
can be represented separately in a mathematical model for optical sensing.

The proposed optical sensing model was also experimentally validated. The
reconstruction of the RL emission quantity using the sensing model can achieve
an average $R^{2}$ of 93.23%. However, the validation only included the turbid
samples with two scattering coefficients. More samples with more scattering
coefficients are needed to test the sensing model. Besides, the Zonios’ model
was used because it is simplified and also representative of diffuse
reflectance. A limitation is that the phase function $p(\Theta)$ or the $g$
factor is not considered, although it has been shown that the $p(\Theta)$ or
$g$ factor has an influence on the diffuse reflectance, especially when the
source and detector separation is very close. Whether these factors also
affect the RL quantities is still questionable. Therefore, further
investigations such as the influence of $p(\Theta)$ or $g$ factor, the
influence of the setup geometry and their interrelationship are needed to
complete the proposed DRS-RL sensing model. Furthermore, the transformation of
the proposed sensing model from the liquid sample to a solid tissue phantom,
which is inharmonious, needs to be investigated for a more practical
application.

## Acknowledgment

This work was funded by the Deutsche Forschungsgemeinschaft (German Research
Foundation – DFG 414732368). In addition, the authors gratefully acknowledge
the funding of SAOT by the Bavarian State Ministry for Science and Art.

## References

  * 1. François Flory.  Optical properties of nanostructured materials: a review.  Journal of Nanophotonics, 5(1):052502, 2011.
  * 2. Steven L. Jacques.  Optical properties of biological tissues: A review.  Physics in Medicine and Biology, 58(14):5007–5008, 2013.
  * 3. Martin Hohmann, Moritz Späth, Dongqin Ni, Dominique Dörner, Benjamin Lengenfelder, Florian Klämpfl, and Michael Schmidt.  Random laser as a potential tool for the determination of the scattering coefficient.  Biomedical Optics Express, 12(9):5439, 2021.
  * 4. Dongqin Ni, Moritz Späth, Florian Klämpfl, and Martin Hohmann.  Properties and Applications of Random Lasers as Emerging Light Sources and Optical Sensors: A Review.  Sensors, 23(1):1–29, 2023.
  * 5. Martin Hohmann, Dominique Dörner, Fanuel Mehari, Chen Chen, Moritz Späth, Sebastian Müller, Heinz Albrecht, Florian Klämpfl, and Michael Schmidt.  Investigation of random lasing as a feedback mechanism for tissue differentiation during laser surgery.  Biomedical Optics Express, 10(2):807, 2019.
  * 6. R. C. Polson and Z. V. Vardeny.  Random lasing in human tissues.  Applied Physics Letters, 85:1289–1291, 2004.
  * 7. Neda Mogharari and Batool Sajad.  Random Laser Emission Spectra of the Normal and Cancerous Thyroid Tissues.  Iranian Journal of Science and Technology, Transaction A: Science, 43:2055–2060, 2019.
  * 8. Qinghai Song, Zhengbin Xu, Seung Ho Choi, Xuanhao Sun, Shumin Xiao, Ozan Akkus, and Young L Kim.  Detection of nanoscale structural changes in bone using random lasers.  Biomedical Optics Express, 1(5):1401, 2010.
  * 9. Federico Tommasi, Emilio Ignesti, Lorenzo Fini, Fabrizio Martelli, and Stefano Cavalieri.  Random laser based method for direct measurement of scattering properties.  Optics Express, 26(21):27615–27627, 2018.
  * 10. Dongqin Ni, Moritz Späth, Florian Klämpfl, Michael Schmidt, and Martin Hohmann.  Threshold behavior and tunability of a diffusive random laser.  Optics Express, 31(16):25747, 2023.
  * 11. A. J. Moy and J. W. Tunnell.  Diffuse Reflectance Spectroscopy and Imaging.  Elsevier Inc., 2016.
  * 12. Diederik S. Wiersma and Ad Lagendijk.  Light diffusion with gain and random lasers.  Physical Review E, 54(4):4256–4265, 1996.
  * 13. Sharmin Akter, Md Golzar Hossain, Izumi Nishidate, Hisanao Hazama, and Kunio Awazu.  Medical applications of reflectance spectroscopy in the diffusive and sub-diffusive regimes.  Journal of Near Infrared Spectroscopy, 26(6):337–350, 2018.
  * 14. Marcelo Saito Nogueira, Siddra Maryam, Michael Amissah, Huihui Lu, Noel Lynch, Shane Killeen, Micheal O’Riordain, and Stefan Andersson-Engels.  Evaluation of wavelength ranges and tissue depth probed by diffuse reflectance spectroscopy for colorectal cancer detection.  Scientific Reports, 11(1):1–17, 2021.
  * 15. Katherine W. Calabro and Irving J. Bigio.  Influence of the phase function in generalized diffuse reflectance models: review of current formalisms and novel observations.  Journal of Biomedical Optics, 19(7):075005, 2014.
  * 16. S. C. Kanick, V. Krishnaswamy, U. A. Gamm, H. J. C. M. Sterenborg, D. J. Robinson, A. Amelink, and B. W. Pogue.  Scattering phase function spectrum makes reflectance spectrum measured from Intralipid phantoms and tissue sensitive to the device detection geometry.  Biomedical Optics Express, 3(5):1086, 2012.
  * 17. George Zonios and Aikaterini Dimou.  Modeling diffuse reflectance from semi-infinite turbid media: application to the study of skin optical properties.  Optics Express, 14(19):8661, 2006.
  * 18. Ricky Hennessy, Will Goth, Manu Sharma, Mia K. Markey, and James W. Tunnell.  Effect of probe geometry and optical properties on the sampling depth for diffuse reflectance spectroscopy.  Journal of Biomedical Optics, 19(10):107002, 2014.
  * 19. Ben Aernouts, Robbe Van Beers, Rodrigo Watté, Jeroen Lammertyn, and Wouter Saeys.  Dependent scattering in Intralipid® phantoms in the 600-1850 nm range.  Optics Express, 22(5):6086, 2014.


# Next-to-leading order QCD corrections to diphoton-plus-jet production
through gluon fusion at the LHC

Simon Badger  Thomas Gehrmann  Matteo Marcoli  Ryan Moodie

###### Abstract

We compute the next-to-leading order (NLO) QCD corrections to the gluon-fusion
subprocess of diphoton-plus-jet production at the LHC. We compute fully
differential distributions by combining two-loop virtual corrections with one-
loop real radiation using antenna subtraction to cancel infrared divergences.
We observe significant corrections at NLO which demonstrate the importance of
combining these corrections with the quark-induced diphoton-plus-jet channel
at next-to-next-to-leading order (NNLO).

###### keywords:

NLO corrections

††journal: Physics Letters B

IPPP/21/31, ZU-TH 45/21

[tor]organization=Dipartimento di Fisica and Arnold-Regge Center, Università
di Torino, and INFN, Sezione di Torino,addressline=Via P. Giuria 1, city=10124
Torino, country=Italy [uzh]organization=Physik-Institut, Universität
Zürich,addressline=Wintherturerstrasse 190, city=CH-8057 Zürich,
country=Switzerland [dur]organization=Institute for Particle Physics
Phenomenology, Department of Physics, Durham University,addressline=South
Road, city=Durham, postcode=DH1 3LE, country=UK

##  1 Introduction

Recent breakthroughs in two-loop amplitude technology are opening up a new
range of precision two-to-three scattering problems. Diphoton-plus-jet
production has been one of the first predictions to appear at NNLO in QCD [1].
This progress is extremely timely given the continually improving experimental
measurements of diphoton signatures [2]. Predictions for pure diphoton
production have been known to NNLO accuracy for some time [3, 4, 5, 6, 7]. A
$q_{T}$-resummed calculation at order
$\text{N}^{3}\text{LL}^{\prime}+\text{NNLO}$ was presented recently [8]. Steps
towards N3LO are being taken with the completion of the three-loop amplitudes
[9]. Diphoton-plus-jet signatures are of particular importance at the LHC
since they form the largest background to Higgs production at high transverse
momenta. The extra jet is necessary to ensure a non-zero transverse momenta in
the diphoton system.

The recently computed NNLO corrections [1] of diphoton-plus-jet production
display a good perturbative convergence, except in regions where the loop-
mediated gluon-fusion process (which contributes to the cross section only
from NNLO onwards) is numerically sizable compared to other contributions. In
order to capture the full effects of the QCD corrections, it is important to
include loop-induced gluon-fusion channels from at least one order higher in
the perturbative series. These corrections are the subject of this article.

High precision two-to-three scattering problems have presented an enormous
challenge to the theoretical community. The development of new techniques and
methodology have been necessary to address several major bottlenecks that have
prevented predictions at NNLO in QCD from being completed.

One important ingredient is the two-loop amplitudes for which complete sets of
helicity amplitudes have recently been completed [10, 11, 12, 13]. These new
results have been achieved thanks to a complete understanding of the special
functions basis [14, 15, 16, 17] and a new range of reduction tools based in
finite field arithmetic [18, 19, 20]. The end products are fully analytic
formulae which can be evaluated efficiently over the phenomenologically
relevant phase-space [21, 22, 23, 10, 11, 12, 13].

Combining and integrating the amplitudes into differential cross sections
requires the subtraction of infrared divergences. To achieve this in a stable
and efficient way is an extremely hard problem and many solutions have been
proposed and applied in calculations up to NNLO. Such algorithms often scale
poorly with the number of external particles and only a handful of examples
for high multiplicity processes at NNLO currently exist [24, 25, 26, 1].

For the process considered in this article, the infrared divergences are only
at NLO. However, since the real radiation involves $2\to 4$ one-loop squared
amplitudes, the automated numerical algorithms are tested in extreme phase-
space regions. The leading order (LO) QCD contributions to the gluonic
subprocess were first considered in Ref. [27] based on the compact one-loop
five-gluon amplitudes [28].

Our paper is organised as follows. We first review the computational setup,
discussing the amplitude-level ingredients and antenna subtraction method used
to cancel infrared divergences. We then present results for the NLO
corrections to differential cross sections at the 13 TeV LHC. We study the
perturbative convergence in both transverse momentum and mass variables as
well angular distributions in rapidity and the Collins-Soper angle before
drawing our conclusions.

##  2 Computational setup

We consider the scattering process

$gg\to\gamma\gamma g+X$ (1)

at a hadron collider. As the process is loop-induced, the LO contribution is
at order $\alpha_{s}^{3}$ and involves the integration of a one-loop amplitude
squared. The NLO QCD corrections are computed by combining the two-loop
virtual corrections to the $2\to 3$ process with the $2\to 4$ processes with
an additional unresolved parton: $gg\to\gamma\gamma gg$ and $gg\to\gamma\gamma
q\bar{q}$. Pictorially, we can represent the parton level cross sections up to
NLO in QCD as,

$\displaystyle\sigma_{gg\to\gamma\gamma g+X}^{\rm NLO}=$ $\displaystyle\int
d\Phi_{3}\left|\raisebox{-17.07164pt}{\includegraphics[width=42.67912pt]{amppic-3g2a1L}}\right|^{2}+\int
d\Phi_{3}2{\rm
Re}\left(\raisebox{-17.07164pt}{\includegraphics[width=42.67912pt]{amppic-3g2a1L}}^{\dagger}\cdot\raisebox{-17.07164pt}{\includegraphics[width=51.21504pt]{amppic-3g2a2L}}\right)+$
$\displaystyle\int
d\Phi_{4}\left|\raisebox{-17.07164pt}{\includegraphics[width=42.67912pt]{amppic-4g2a1L}}\right|^{2}+\int
d\Phi_{4}\left|\raisebox{-17.07164pt}{\includegraphics[width=42.67912pt]{amppic-2q2g2a1L}}\right|^{2}+\mathcal{O}(\alpha_{s}^{5}),$
(2)

where $d\Phi_{n}$ represents the on-shell phase-space measure for $n$ massless
final state particles. The one-loop amplitude for $q\bar{q}gg\gamma\gamma$
indicates the loop contribution in which the photons couple to an internal
fermion loop. The observable process $pp\to\gamma\gamma j$ also includes
channels where the photons couple to an external quark pair. The expansion up
to the NNLO of $pp\to\gamma\gamma j$ includes terms up to
$\mathcal{O}(\alpha_{s}^{3})$ and so the contributions coming from Eq. (2) are
technically N3LO. However, due to the large gluon flux at high energy hadron
colliders, such contributions can be significant.

The one-loop amplitudes for the LO process and the real correction are finite,
since the corresponding tree-level processes vanish. The renormalised two-loop
five-particle amplitude contains explicit infrared divergences generated by
the integration over the loop momenta, while the one-loop six-particle
amplitudes exhibit a divergent behavior when a final-state parton becomes
unresolved. The divergences cancel in the final result, as established by the
KLN theorem, and a finite remainder of the virtual amplitudes can be defined
using QCD factorization [29]. In our calculation, this cancellation is
performed using the antenna subtraction method [30, 31, 32]. The method
extracts the infrared singular contributions from the real radiation
subprocess, and combines their integrated form with the virtual subprocess,
thus enabling their numerical integration using Monte Carlo methods, performed
here in the NNLOjet framework. The QCD structure of the process under
consideration is very similar to Higgs-plus-jet production in gluon fusion,
which has been computed previously [33, 34] using antenna subtraction, and
identical antenna subtraction terms are applied here.

The infrared-finite remainders of the two-loop amplitudes have recently been
computed [13] using a basis of pentagon functions [14, 15, 17], which permit
efficient and reliable numerical evaluation [17]. The full colour and helicity
summed expressions are obtained from the NJet amplitude library. Within NJet,
a dimension scaling test is performed for each phase-space point to assess the
accuracy of the evaluation. If the test is unsuccessful, the point is
recomputed in higher precision. We set a three digits accuracy threshold for
this test, which guarantees a stable result without significantly affecting
the performance.

The one-loop six-particle amplitudes are obtained using a combination of
implementations from the OpenLoops2 [35] generator and from the generalised
unitarity [36, 37, 38] approach within NJet [39]. We use an improved version
of OpenLoops2 in combination with the new extension Otter. Otter [40] is a
tensor integral library based on the on-the-fly reduction [41] of OpenLoops2
and on stability improvements described in Ref. [35]. This new version of
OpenLoops2 allows for a stable computation of the needed one-loop squared
amplitudes in deep infrared regions. Internally, Otter uses double-precision
scalar integrals that are provided by Collier [42, 43], as well as quad-
precision scalar integrals provided by OneLoop [44]. Minor modifications were
made in NJet to avoid de-symmetrisation over the two photons and allow for a
pointwise correspondence with the subtraction terms. To compute the one-loop
amplitude $gggg\gamma\gamma$, the OpenLoops implementation was generally more
efficient, but for exceptional phase-space points it was necessary to use the
high precision (32 digits) implementation within NJet. For the
$q\bar{q}gg\gamma\gamma$ channel, we used NJet, which allowed for a
straightforward selection of the required loop contribution. We note that this
amplitude is also available within OpenLoops2 and we checked that the two
implementations agree.

The amplitude-level ingredients have been validated in all relevant collinear
and soft limits by checking their convergence towards the respective antenna
subtraction terms.

##  3 Results

For the numerical evaluation of our NLO results on the gluon-induced diphoton-
plus-jet process, we apply the same kinematical cuts as were used for the NNLO
calculation [1] of the quark-induced processes. These represent a realistic
setup relevant for physics studies at the 13 TeV LHC. The cuts are as follows:

  * 1.

minimum photon transverse momenta and rapidities: $p_{T}(\gamma_{1})>30$ GeV,
$p_{T}(\gamma_{2})>18$ GeV and $|\eta(\gamma\gamma)|<2.4.$

  * 2.

smooth photon isolation criterion [45] with $\Delta R_{0}=0.4$, $E_{T}^{\rm
max}=10$ GeV and $\epsilon=1$.

  * 3.

minimal invariant mass of the photon pair: $m(\gamma\gamma)\geq 90$ GeV.

  * 4.

minimal separation of the photons: $\Delta R_{\gamma\gamma}>0.4$.

  * 5.

minimal transverse momentum of the photon pair: $p_{T}(\gamma\gamma)>20$ GeV.

We consider kinematical distributions in the following diphoton variables:
transverse momentum of the diphoton system $p_{T}(\gamma\gamma)$, pair
invariant mass $m_{\gamma\gamma}$, diphoton total rapidity $|y(\gamma\gamma)|$
and rapidity difference $\Delta y(\gamma\gamma)$, as well as Collins-Soper
angle $\left|\phi_{CS}(\gamma\gamma)\right|$ [46] and azimuthal decorrelation
$\Delta\phi(\gamma\gamma)$. For these distributions, no jet requirement is
applied, as done in Ref. [1], since the transverse momentum cut on the
diphoton system is already sufficient to avoid NNLO-like configurations where
all final-state QCD partons become unresolved.

Our numerical results use the NNLO set of the NNPDF3.1 parton distribution
functions [47] throughout, thus allowing a straightforward comparison with the
existing NNLO results [1] in the quark-initiated channels. The strong coupling
constant is evaluated using LHAPDF [48], with
$\alpha_{s}\left(m_{Z}\right)=0.118$. The electromagnetic coupling constant is
set to $\alpha=1/137.035999139$. Monte Carlo integration errors are below
$1\%$ on average and not displayed in the plots.

The uncertainty on our theory predictions is estimated by a seven-point
variation of the renormalisation and factorisation scales around a central
value, chosen in dynamical manner on event-by-event basis to be

$\mu_{F}=\mu_{R}=\frac{1}{2}m_{T}=\frac{1}{2}\left(m^{2}(\gamma\gamma)+p_{T}^{2}(\gamma\gamma)\right)^{1/2}$
(3)

Figure 1: Differential distributions in the transverse momentum
$p_{T}(\gamma\gamma)$ (left) and invariant mass $m(\gamma\gamma)$ (right) of
the diphoton system.

Figure 2: Differential distributions in the Collins-Soper angle
$\left|\cos\phi_{CS}(\gamma\gamma)\right|$ (left) the azimuthal decorrelation
$\Delta\phi(\gamma\gamma)$ (right) of the diphoton system.

Figures 2–3 display the theory predictions for the different single-
differential distributions in the diphoton variables. We observe the NLO
corrections to be very sizable, often being comparable in size to the LO
predictions. The corrections are largest at low $p_{T}(\gamma\gamma)$ or at
low invariant mass, Figure 2, where the NLO/LO ratio reaches 2 and NLO and LO
uncertainties fail to overlap, while the ratio is smoothly decreasing towards
values of 1.5 for large $p_{T}(\gamma\gamma)$ or $m(\gamma\gamma)$, with
overlapping scale uncertainty bands above $p_{T}(\gamma\gamma)=200$ GeV or
$m(\gamma\gamma)=175$ GeV.

Figure 3: Differential distributions in the diphoton rapidity difference
$\Delta y(\gamma\gamma)$ (left) and the diphoton total rapidity
$|y(\gamma\gamma)|$ (right).

The integrated cross section is dominated by the region of low
$p_{T}(\gamma\gamma)$ or low $m(\gamma\gamma)$, such that distributions that
are differential only in geometrical photon variables, Figures 2 and 3,
display typically near-uniform NLO/LO ratios of 2, and no overlap of the LO
and NLO scale uncertainty bands. Visually, the scale uncertainty bands at NLO
and LO appear to be of comparable width in all distributions. However, owing
to the large size of the NLO corrections, the relative scale uncertainty is
reduced from about 50% at LO to 30% at NLO.

Figure 4: Two-dimensional differential distributions in the diphoton invariant
mass $m(\gamma\gamma)$ and Collins-Soper angle
$\left|\phi_{CS}(\gamma\gamma)\right|$ (left) and in the diphoton rapidity
$|y(\gamma\gamma)|$ and transverse momentum $p_{T}(\gamma\gamma)$ (right).

Figure 5: Two-dimensional differential distributions in the diphoton
transverse momentum $p_{T}(\gamma\gamma)$ and invariant mass
$m(\gamma\gamma)$, in bins in $p_{T}(\gamma\gamma)$ (left) and for varying
lower $p_{T}(\gamma\gamma)$-cut (right)

By inspecting the two-dimensional differential distribution in
$m(\gamma\gamma)$ and $\left|\phi_{CS}(\gamma\gamma)\right|$, Figure 5 (left),
we observe that the relative magnitude of the NLO corrections decreases with
increasing $m(\gamma\gamma)$, while the corrections remain uniform in
$\left|\phi_{CS}(\gamma\gamma)\right|$ for all bins in $m(\gamma\gamma)$. The
two-dimensional differential distribution in $|y(\gamma\gamma)|$ and
$p_{T}(\gamma\gamma)$ also shows the decrease of the corrections towards
larger $p_{T}(\gamma\gamma)$. The decrease is more pronounced at forward
rapidity than at central rapidity.

Considering two-dimensional distributions in $p_{T}(\gamma\gamma)$ and
$m(\gamma\gamma)$, Figure 5, largely reproduces the features of the one-
dimensional distributions of Figure 2, both for distributions in bins of
$p_{T}(\gamma\gamma)$ or for varying lower cut in $p_{T}(\gamma\gamma)$. The
only novel feature is a non-uniform shape in $m(\gamma\gamma)$ for the highest
bin in $p_{T}(\gamma\gamma)$ (lowest curves in left Figure 5), which is
indicative of the onset of large logarithmic corrections in
$\log(m(\gamma\gamma)/p_{T}(\gamma\gamma))$ in this range.

The numerical size of the NLO corrections and the scale uncertainties at LO
and NLO are comparable to what was observed in inclusive Higgs boson
production in gluon fusion [49] or in the Higgs boson transverse momentum
distribution in gluon fusion [50, 51]. These processes are mediated through a
heavy top quark loop and are very similar to the diphoton-plus-jet production
considered here in terms of kinematics and initial-state parton momentum
range. The pathology of the NLO corrections observed here is thus not that
surprising after all; it does however indicate the potential numerical
importance of corrections beyond NLO.

The Born-level $gg\to\gamma\gamma g$ subprocess (corresponding to the LO in
our results) contributes to the full diphoton-plus-jet production as part of
the NNLO corrections. Corrections to this order were computed most recently
[1]: these were observed to be moderate and within the scale uncertainty of
the previously known NLO results for most of the kinematical range, where they
also led to a substantial reduction of the scale uncertainty at NNLO. Only at
low $p_{T}(\gamma\gamma)$ or low $m(\gamma\gamma)$, larger positive
corrections and an increased scale uncertainty were observed [1]. These
effects could be identified to be entirely due to the contribution of the
$gg\to\gamma\gamma g$, which only starts to contribute from NNLO onwards, and
it was anticipated in Ref. [1] that NLO corrections to the $gg\to\gamma\gamma
g$ (which form a subset of the N3LO corrections to the full diphoton-plus-jet
process) could help to stabilise the predictions in the relevant kinematical
ranges.

Our results demonstrate that this is not the case. The absolute scale
uncertainty on the gluon-induced process does not decrease from LO to NLO, and
the NLO correction is of about the size of the LO contribution. Consequently,
inclusion of the NLO corrections to the $gg\to\gamma\gamma g$ into the full
NNLO diphoton-plus-jet process will further enhance the predictions at low
$p_{T}(\gamma\gamma)$ or low $m(\gamma\gamma)$, thereby further elongating
them from the previously known order, and will leave the scale uncertainty
band largely unchanged.

##  4 Conclusions

In this article, we have presented the NLO QCD corrections to the diphoton-
plus-jet production in the gluon-fusion channel for the first time. The loop-
induced process requires the evaluation of six-point one-loop real emission
amplitudes and full-colour five-point two-loop virtual amplitudes. To the best
of our knowledge it is the first time that five-point two-loop amplitudes with
the full colour information have been integrated to provide fully differential
cross section predictions relevant for the LHC experiments.

Using a realistic set of kinematic cuts and simulation parameters, we find
significant corrections at NLO. This is particularly relevant at low values of
$p_{T}(\gamma\gamma)$ and $m(\gamma\gamma)$. Since angular observables such as
rapidity and the Collins-Soper angle are inclusive over the energy variables,
one observes significant NLO corrections across the full parameter range.
Double differential distributions further highlight this feature, which is
reminiscent of the perturbative convergence observed in other gluon-induced
processes such as inclusive Higgs production and the Higgs boson transverse
momentum distribution. The relative scale uncertainty is reduced by the higher
order corrections, although in absolute terms the scale uncertainty does not
decrease from LO to NLO in the low $p_{T}(\gamma\gamma)$ and $m(\gamma\gamma)$
regions.

This work demonstrates the importance of a combined prediction for quark-
induced and gluon-induced diphoton-plus-jet signatures for future precision
studies at the LHC.

## Acknowledgements

We would like to thanks Johannes Henn, Xuan Chen, Alexander Huss and Simone
Zoia for numerous discussions and interesting input in the course of this
project. Help with OpenLoops2 from Jean-Nicolas Lang and Federico Buccioni is
gratefully acknowledged. This research has been supported in part by the Swiss
National Science Foundation (SNF) under contract number 200020-175595 and by
the Swiss National Supercomputing Centre (CSCS) under project ID UZH10. This
project received funding from the European Union’s Horizon 2020 research and
innovation programme High precision multi-jet dynamics at the LHC (ERC
Condsolidator grant agreement No 772009). RM was supported by STFC grants
ST/S505365/1 and ST/P001246/1.

## References

  * [1] H. A. Chawdhry, M. Czakon, A. Mitov, R. Poncelet, NNLO QCD corrections to diphoton production with an additional jet at the LHC (2021).  arXiv:2105.06940.
  * [2] G. Aad, et al., Measurement of the production cross section of pairs of isolated photons in $pp$ collisions at 13 TeV with the ATLAS detector (2021).  arXiv:2107.09330.
  * [3] C. Anastasiou, E. W. N. Glover, M. E. Tejeda-Yeomans, Two loop QED and QCD corrections to massless fermion boson scattering, Nucl. Phys. B 629 (2002) 255–289.  arXiv:hep-ph/0201274.
  * [4] Z. Bern, A. De Freitas, L. J. Dixon, Two loop amplitudes for gluon fusion into two photons, JHEP 09 (2001) 037.  arXiv:hep-ph/0109078.
  * [5] S. Catani, L. Cieri, D. de Florian, G. Ferrera, M. Grazzini, Diphoton production at hadron colliders: a fully-differential QCD calculation at NNLO, Phys. Rev. Lett. 108 (2012) 072001, [Erratum: Phys.Rev.Lett. 117, 089901 (2016)].  arXiv:1110.2375.
  * [6] J. M. Campbell, R. K. Ellis, Y. Li, C. Williams, Predictions for diphoton production at the LHC through NNLO in QCD, JHEP 07 (2016) 148.  arXiv:1603.02663.
  * [7] T. Gehrmann, N. Glover, A. Huss, J. Whitehead, Scale and isolation sensitivity of diphoton distributions at the LHC, JHEP 01 (2021) 108.  arXiv:2009.11310.
  * [8] T. Neumann, The diphoton $q_{T}$ spectrum at N3LL′ \+ NNLO, Eur. Phys. J. C 81 (10) (2021) 905.  arXiv:2107.12478.
  * [9] F. Caola, A. Von Manteuffel, L. Tancredi, Diphoton Amplitudes in Three-Loop Quantum Chromodynamics, Phys. Rev. Lett. 126 (2021) 112004.  arXiv:2011.13946.
  * [10] B. Agarwal, F. Buccioni, A. von Manteuffel, L. Tancredi, Two-loop leading colour QCD corrections to $q\bar{q}\to\gamma\gamma g$ and $qg\to\gamma\gamma q$, JHEP 04 (2021) 201.  arXiv:2102.01820.
  * [11] H. A. Chawdhry, M. Czakon, A. Mitov, R. Poncelet, Two-loop leading-colour QCD helicity amplitudes for two-photon plus jet production at the LHC, JHEP 07 (2021) 164.  arXiv:2103.04319.
  * [12] B. Agarwal, F. Buccioni, A. von Manteuffel, L. Tancredi, Two-loop helicity amplitudes for diphoton plus jet production in full color (2021).  arXiv:2105.04585.
  * [13] S. Badger, C. Brønnum-Hansen, D. Chicherin, T. Gehrmann, H. B. Hartanto, J. Henn, M. Marcoli, R. Moodie, T. Peraro, S. Zoia, Virtual QCD corrections to gluon-initiated diphoton plus jet production at hadron colliders (2021).  arXiv:2106.08664.
  * [14] D. Chicherin, J. Henn, V. Mitev, Bootstrapping pentagon functions, JHEP 05 (2018) 164.  arXiv:1712.09610.
  * [15] T. Gehrmann, J. Henn, N. Lo Presti, Pentagon functions for massless planar scattering amplitudes, JHEP 10 (2018) 103.  arXiv:1807.09812.
  * [16] D. Chicherin, T. Gehrmann, J. Henn, P. Wasser, Y. Zhang, S. Zoia, All Master Integrals for Three-Jet Production at Next-to-Next-to-Leading Order, Phys. Rev. Lett. 123 (2019) 041603.  arXiv:1812.11160.
  * [17] D. Chicherin, V. Sotnikov, Pentagon Functions for Scattering of Five Massless Particles, JHEP 12 (2020) 167.  arXiv:2009.07803.
  * [18] A. von Manteuffel, R. M. Schabinger, A novel approach to integration by parts reduction, Phys. Lett. B744 (2015) 101–104.  arXiv:1406.4513.
  * [19] T. Peraro, Scattering amplitudes over finite fields and multivariate functional reconstruction, JHEP 12 (2016) 030.  arXiv:1608.01902.
  * [20] T. Peraro, FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs, JHEP 07 (2019) 031.  arXiv:1905.08019.
  * [21] S. Abreu, B. Page, E. Pascual, V. Sotnikov, Leading-Color Two-Loop QCD Corrections for Three-Photon Production at Hadron Colliders, JHEP 01 (2021) 078\.  arXiv:2010.15834.
  * [22] H. A. Chawdhry, M. Czakon, A. Mitov, R. Poncelet, Two-loop leading-color helicity amplitudes for three-photon production at the LHC, JHEP 06 (2021) 150\.  arXiv:2012.13553.
  * [23] S. Abreu, F. F. Cordero, H. Ita, B. Page, V. Sotnikov, Leading-color two-loop QCD corrections for three-jet production at hadron colliders, JHEP 07 (2021) 095\.  arXiv:2102.13609.
  * [24] M. Czakon, A. Mitov, R. Poncelet, Next-to-Next-to-Leading Order Study of Three-Jet Production at the LHC, Phys. Rev. Lett. 127 (15) (2021) 152001.  arXiv:2106.05331.
  * [25] H. A. Chawdhry, M. L. Czakon, A. Mitov, R. Poncelet, NNLO QCD corrections to three-photon production at the LHC, JHEP 02 (2020) 057.  arXiv:1911.00479.
  * [26] S. Kallweit, V. Sotnikov, M. Wiesemann, Triphoton production at hadron colliders in NNLO QCD, Phys. Lett. B 812 (2021) 136013.  arXiv:2010.04681.
  * [27] D. de Florian, Z. Kunszt, Two photons plus jet at LHC: The NNLO contribution from the g g initiated process, Phys. Lett. B 460 (1999) 184–188.  arXiv:hep-ph/9905283.
  * [28] Z. Bern, L. J. Dixon, D. A. Kosower, One loop corrections to five gluon amplitudes, Phys. Rev. Lett. 70 (1993) 2677–2680.  arXiv:hep-ph/9302280.
  * [29] S. Catani, The Singular behavior of QCD amplitudes at two loop order, Phys. Lett. B 427 (1998) 161–171.  arXiv:hep-ph/9802439.
  * [30] A. Gehrmann-De Ridder, T. Gehrmann, E. W. N. Glover, Antenna subtraction at NNLO, JHEP 09 (2005) 056.  arXiv:hep-ph/0505111.
  * [31] A. Daleo, T. Gehrmann, D. Maitre, Antenna subtraction with hadronic initial states, JHEP 04 (2007) 016.  arXiv:hep-ph/0612257.
  * [32] J. Currie, E. W. N. Glover, S. Wells, Infrared Structure at NNLO Using Antenna Subtraction, JHEP 04 (2013) 066.  arXiv:1301.4693.
  * [33] X. Chen, T. Gehrmann, E. W. N. Glover, M. Jaquier, Precise QCD predictions for the production of Higgs + jet final states, Phys. Lett. B 740 (2015) 147–150.  arXiv:1408.5325.
  * [34] X. Chen, J. Cruz-Martinez, T. Gehrmann, E. W. N. Glover, M. Jaquier, NNLO QCD corrections to Higgs boson production at large transverse momentum, JHEP 10 (2016) 066.  arXiv:1607.08817.
  * [35] F. Buccioni, J.-N. Lang, J. M. Lindert, P. Maierhöfer, S. Pozzorini, H. Zhang, M. F. Zoller, OpenLoops 2, Eur. Phys. J. C 79 (2019) 866.  arXiv:1907.13071.
  * [36] Z. Bern, L. J. Dixon, D. C. Dunbar, D. A. Kosower, One loop n point gauge theory amplitudes, unitarity and collinear limits, Nucl. Phys. B 425 (1994) 217–260.  arXiv:hep-ph/9403226.
  * [37] Z. Bern, L. J. Dixon, D. C. Dunbar, D. A. Kosower, Fusing gauge theory tree amplitudes into loop amplitudes, Nucl. Phys. B 435 (1995) 59–101.  arXiv:hep-ph/9409265.
  * [38] R. Britto, F. Cachazo, B. Feng, Generalized unitarity and one-loop amplitudes in N=4 super-Yang-Mills, Nucl. Phys. B 725 (2005) 275–305.  arXiv:hep-th/0412103.
  * [39] S. Badger, B. Biedermann, P. Uwer, V. Yundin, Numerical evaluation of virtual corrections to multi-jet production in massless QCD, Comput. Phys. Commun. 184 (2013) 1981–1998.  arXiv:1209.0100.
  * [40] J.-N. Lang et al., OTTER: On-The-fly TEnsor Reduction, To be published.
  * [41] F. Buccioni, S. Pozzorini, M. Zoller, On-the-fly reduction of open loops, Eur. Phys. J. C 78 (2018) 70.  arXiv:1710.11452.
  * [42] A. Denner, S. Dittmaier, L. Hofer, COLLIER - A fortran-library for one-loop integrals, PoS LL2014 (2014) 071.  arXiv:1407.0087.
  * [43] A. Denner, S. Dittmaier, L. Hofer, Collier: a fortran-based Complex One-Loop LIbrary in Extended Regularizations, Comput. Phys. Commun. 212 (2017) 220–238.  arXiv:1604.06792.
  * [44] A. van Hameren, OneLOop: For the evaluation of one-loop scalar functions, Comput. Phys. Commun. 182 (2011) 2427–2438.  arXiv:1007.4716.
  * [45] S. Frixione, Isolated photons in perturbative QCD, Phys. Lett. B 429 (1998) 369–374.  arXiv:hep-ph/9801442.
  * [46] J. C. Collins, D. E. Soper, Angular Distribution of Dileptons in High-Energy Hadron Collisions, Phys. Rev. D 16 (1977) 2219.
  * [47] R. D. Ball, et al., Parton distributions from high-precision collider data, Eur. Phys. J. C 77 (2017) 663.  arXiv:1706.00428.
  * [48] A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr, G. Watt, LHAPDF6: parton density access in the LHC precision era, Eur. Phys. J. C 75 (2015) 132.  arXiv:1412.7420.
  * [49] A. Djouadi, M. Spira, P. M. Zerwas, Production of Higgs bosons in proton colliders: QCD corrections, Phys. Lett. B 264 (1991) 440–446.
  * [50] D. de Florian, M. Grazzini, Z. Kunszt, Higgs production with large transverse momentum in hadronic collisions at next-to-leading order, Phys. Rev. Lett. 82 (1999) 5209–5212.  arXiv:hep-ph/9902483.
  * [51] V. Ravindran, J. Smith, W. L. Van Neerven, Next-to-leading order QCD corrections to differential distributions of Higgs boson production in hadron hadron collisions, Nucl. Phys. B 634 (2002) 247–290.  arXiv:hep-ph/0201114.