\section{Introduction}
\input{figures/title_figure}

Radiologists interpret a vast amount of imaging data and summarize their insights as medical reports. This documentation accumulates large databases of radiological imaging and accompanying findings, i.e.,  millions of collected chest radiographs annually~\cite{nhs}. Computer-aided-diagnosis (CAD) systems utilize these databases to streamline the clinical workflow and save time~\cite{qin2018computer,kim2019deep,jaiswal2019identifying}.
Modern CAD tools often rely on deep learning models~\cite{qin2018computer} using large-scale data sets such as MIMIC-CXR~\cite{wang2017chestx,johnson2019mimic,irvin2019chexpert,bustos2020padchest} for training.
Training for such tasks requires hand-designed supervision, typically by extracting a fixed set of predefined labels from the reports using rule- or deep learning-based models~\cite{wang2017chestx,irvin2019chexpert,smit2020chexbert}.
Such training typically requires hand-designed supervision in the form of extracting a set of labels from the reports using rule- or deep learning-based models~\cite{wang2017chestx,irvin2019chexpert,smit2020chexbert}.
While these tools can deliver acceptable performance on a subset of diseases~\cite{wu2020comparison}, they lack generalization capabilities for diseases that were not part of the fixed label-set used for training.
To add disease classes requires substantial effort annotating the data with extra labels and retraining the system. To circumvent this, one can approach training in a class-agnostic manner, however, it becomes unclear how models can still be applied to classify diseases.

Recent methods based on contrastive language-image pre-training~\cite{radford2021learning,jia2021scaling,pham2021combined,zhang2020contrastive,wang2021self} indicate that by large-scale multi-modal representation learning, object recognition can be detached from prior fixed-set, hand-designed class definitions. These models learn joint feature spaces between images and textual descriptions and utilize text prompts to transform recognition from learned fix-set classification to a matching task between text and image embeddings.
Radiological reports, in contrast to natural-image captions, have an inherently different structure, as they encompass multiple distinct sentences such that their entirety describes all relevant information. This shift makes a direct application of existing methods non-trivial.

In this work, we see our contributions as the following:
\begin{enumerate}
    \item We address training through report supervision by considering radiological reports in one of two ways: The local level, assuming each sentence conveys a distinct concept relevant for the patient, and secondly, the global report view, which encodes the entirety of the findings.
    \item We propose a novel inference setting that allows us to query any desired finding, and the CAD system generates a binary decision regarding its presence in the given radiological imagery.
    \item We provide an extensive study on various factors impacting the performance of multi-modal training and inference.
\end{enumerate}

\subsection{Related Work}

\iffalse

https://arxiv.org/pdf/2010.08670.pdf
https://arxiv.org/pdf/1901.07441.pdf
https://arxiv.org/abs/2011.10566
https://papers.nips.cc/paper/2020/file/d89a66c7c80a29b1bdbab0f2a1a94af8-Paper.pdf
https://arxiv.org/pdf/2201.11736.pdf
https://openaccess.thecvf.com/content_CVPR_2020/papers/He_Momentum_Contrast_for_Unsupervised_Visual_Representation_Learning_CVPR_2020_paper.pdf
https://openaccess.thecvf.com/content/CVPR2021W/ECV/papers/Cheng_Data-Efficient_Language-Supervised_Zero-Shot_Learning_With_Self-Distillation_CVPRW_2021_paper.pdf
https://proceedings.neurips.cc/paper/2021/file/505259756244493872b7709a8a01b536-Paper.pdf
https://proceedings.neurips.cc/paper/2020/file/d27b95cac4c27feb850aaa4070cc4675-Paper.pdf
https://arxiv.org/pdf/2006.13979.pdf
https://arxiv.org/abs/2011.02803
https://arxiv.org/pdf/2103.00020.pdf
https://arxiv.org/pdf/2112.12750.pdf
https://arxiv.org/pdf/2011.01403.pdf
https://arxiv.org/pdf/2201.11736.pdf
https://proceedings.mlr.press/v158/endo21a/endo21a.pdf
https://arxiv.org/pdf/2010.00747.pdf
https://proceedings.neurips.cc/paper/2020/file/d89a66c7c80a29b1bdbab0f2a1a94af8-Paper.pdf
https://arxiv.org/pdf/2112.02889.pdf
https://arxiv.org/pdf/1912.06430.pdf
https://arxiv.org/pdf/2110.05208.pdf
https://arxiv.org/pdf/2111.03930.pdf
https://arxiv.org/pdf/2111.14745.pdf
\fi
\section{Global-Local Contrastive Learning}
We illustrate our method for report-based training of a vision model and inference protocol specifically designed for disease recognition in Figure \ref{fig:method}.
To tackle the complexity of medical reports, we split representations into a sentence- and report-level from a shared visual and language encoder. We consider embeddings for both the presence and absence of a pathology for its prediction.
\subsection{Model Overview}
Contrastive language and image pretraining (CLIP) has shown immense potential for object recognition in natural images~\cite{radford2021learning} through learning from image-caption pairs. In contrast to textual descriptions in CLIP-based models, medical reports typically consists of multiple sentences focusing on different parts in the image. As each sentence contains specific subset of information, we aim to capture sentence- (local) and report-level (global)  context in our representations.
Thus, our model builds on separate image- and text encoders $\phi$ and $\theta$, which embed an image $I$ via $z_I = \phi(I)$ and a sentence by $z_s = \theta(s)$, respectively.

In training, for a given report $R$, we capture the local context by splitting the full report into its sentences $R = \{s_1,\dots,s_n\}$ and subsequently extract sentence-level embeddings $z_{s_i} = \theta(s_i), s_i \in R$. To generate global embeddings that contain the full information of the whole report, the sentence-level embeddings are aggregated through attention pooling: $z_{R} = \text{Attn}([z_{s_1},\dots,z_{s_{|R|}}])$~\cite{NIPS2017_3f5ee243}.
To embed $z_I$, $z_s$ and $z_R$ into shared multi-modal representations, we project sentences and reports via linear transformations $p^S$ and $p^R$ into two feature spaces.
As the image encoder has access to global image information for report-level prediction as well as to local image patterns for selective sentence-dependent prediction, we project $z_I$ twice: into a global representation $p^{G}(z_I)$ which shall align with $p^R(z_R)$ and a representation for local patterns $p^{L}(z_I)$ for alignment with $p^S(z_s)$.

During training, we are provided with a dataset of image-report pairs $(I_i, R_i) \in {(I_1, R_1), \dots, (I_N, R_N)}$.
For brevity and clarity in subsequent formulas, we will write projections, \textit{e.g.} the global projection of an image $I_i$ as $p^G_i$ instead of $p^G_i(\theta(I_i))$, for the projection of the $k^{th}$ sentence from report $R_i$, we write $p^S_{ik}$.

\subsection{Training Objectives}

\noindent\textbf{Local Contrast:} While radiological reports describe the assessment of a patient's health, not every sentence is directly linked to specific findings, some sentences mention clinical procedures or required follow-up examinations.
However, we can assume that all clinically relevant information is present in a subset of sentences in the report due to the doctors' obligation to document the findings.
This property is the core of the multiple-instance learning (MIL) assumption.
Therefore, it might seem natural to choose MILNCE~\cite{miech2020end} as the MIL-based objective for integrate sentences in training, yet, this assumption only holds when normalizing over sentences, as not every sentence has to match the image. However, if a sentence fits an image, it should match strictly that image, thus, we hold the regular formulation when normalizing over images.
As such, MILNCE does not quite fit this use-case and we redesign it by splitting its symmetry:
\begin{equation}
    \begin{aligned}
        \mathcal{L}_{L}(I_i, R_i) =  & - \log \frac{
    \sum^{n}_{k=1} \exp(\sigma(p^L_i, p_{ik}^{S})/\tau_L) }
    {\sum^N_{j=1} \sum^{n}_{m=1}  \exp(\sigma(p_i^{L}, p_{jm}^{S}) /\tau_L)  } \\&- \sum^{n}_{k=1} log \frac{
    \exp(\sigma(p_i^{L}, p_{ik}^{S})/\tau_L) }
    {\sum^N_{j=1} \exp(\sigma(p_j^{L}, p_{ik}^{S}) /\tau_L)  },
    \end{aligned}
\end{equation}
with $\tau_L$ being a learned parameter and $\sigma(\cdot,\cdot)$ denoting the cosine similarity.

\noindent\textbf{Global Contrast:} For our batch we assume that an image-report pair is unique and formulate the following objective leveraging the attention-fused  reports via:
\begin{equation}
    \mathcal{L}_{G}(I_i, R_i) =  - \log \frac{
    \exp(\sigma(p_i^{G}, p_{i}^{R})/\tau_G) }
    {\sum^N_{j=1} \exp(\sigma(p_i^{G}, p_{j}^{R}) /\tau_G)  } \\- log \frac{
    \exp(\sigma(p_i^{G}, p_{i}^{R})/\tau_G) }
    {\sum^N_{j=1} \exp(\sigma(p_j^{G}, p_{i}^{R}) /\tau_G)  }
\end{equation}

\noindent\textbf{Self-Supervision:} CLIP has been established as a data-hungry algorithm~\cite{radford2021learning,pham2021combined}. Several recent methods combine intrinsic supervision signals with the CLIP objective to make full use of the available data~\cite{li2021supervision,mu2021slip}.
As we have access to severely smaller datasets in the medical domain as compared to the natural image domain, we follow Li~et~al.~\cite{li2021supervision} and integrate SimSiam~\cite{chen2021exploring}. For this, we generate two augmented versions of the input image $A_1(I)$ and $A_2(I)$ and add a three-layer encoder-head $p^{E}$ and a two-layer prediction-head $p^{P}$ on top of the visual backbone $\phi$ to enforce similarity between the two views:
\begin{equation}
    \mathcal{L}_{S}(A_1(I),A_2(I)) = - \sigma(p^{P}_{A_1(I)}, \text{detach}(p^{E}_{A_2(I)}))  - \sigma(p^{P}_{A_2(I)}, \text{detach}(p^{E}_{A_1(I)}))
\end{equation}
Furthermore, we utilize the augmented images used for the self-supervised objective to mirror our text-image objectives to the augmented samples.
\begin{equation}
    \mathcal{L}_{M}(I_i,R_i) = \mathcal{L}_G(A_1(I_i), R_i) + \mathcal{L}_L(A_1(I_i), R_i) + \mathcal{L}_G(A_2(I_i), R_i) + \mathcal{L}_L(A_2(I_i), R_i)
\end{equation}

Our final objective for report-based contrastive learning amounts to:
\begin{equation}
    \mathcal{L}(I_i,R_i) =
    \lambda_1 * \mathcal{L}_{L}(I_i,R_i) +
    \lambda_2 * \mathcal{L}_{G}(I_i,R_i) +
    \lambda_3 * \mathcal{L}_{S}(I_i,R_i) +
    \lambda_4 * \mathcal{L}_{M}(I_i,R_i)
\end{equation}
with $\lambda_1 = \lambda_2 = \lambda_3 = 0.5$ and $\lambda_4 = 0.25$.

\subsection{Model Inference}
For fixed set classification models, the inference process is straightforward: A given image $I$ passed through a network with an activation in the final layer, returning class-wise pseudo probabilities. When model architecture and training procedure do not permit a classification layer, methods often resort to zero-shot-like inference~\cite{huang2021gloria,radford2021learning,wang2021self} where a nearest neighbor search in semantic space is conducted~\cite{frome2013devise}.
In our considered design a text-based query is used to infer the presence or absence of a given disease. In similar CLIP-like models, the text embeddings (e.g., of the disease names) can be matched to an image embedding. The query with the maximum similarity can then be retained as matched.

Such a matching based disease discovery is feasible for detecting single disease class. The underlying assumption of having exclusively one dominant class to predict does not hold for chest radiographs as pathologies are not mutually exclusive. Similarly, modeling co-occurrence as individual classes is also infeasible due to the exponentially rising number of possible class combinations. As such, inference for a multi-label classification needs to be formulated for such contrastively trained methods.

We perform this, by querying an image with class-related textual prompts and interpreting their similarity scores as prediction probabilities for the respective disease class. In practice, we notice that a single query for class presence is ambiguous since the text embeddings of words and their negations may fall close to one another in the feature space. Due to this proximity of opposing semantics a query could be mistaken with the negation of its class.

To overcome this issue, we propose to perform inference over two sets of queries ($q_c^p,q_c^n$) for each class. While the query embedding $q_c^p$  indicates the occurrence of a class $c$, $q_c^n$ indicates its absence, \textit{e.g.} \textsc{opacities consistent with pneumonia} for presence as opposed to \textsc{the lungs are clear} for its absence.
Then, the cosine similarity between the image and both queries $\sigma(p^I,q_c^p)$ and $\sigma(p^I,q_c^n)$ is computed with the final prediction wrt. class $c$ being defined as:
\begin{equation}
    P(c,I) = \frac{\exp(\sigma(p^I,q_c^p)/\tau)}{\exp(\sigma(p^I,q_c^p)/\tau) + \exp(\sigma(p^I,q_c^n)/\tau)},
\end{equation}
where $\tau$ is the respective learned scaling factor depending on the used projection.

\section{Prompt Engineering}
Several works on zero-shot classification perform their inference by extracting the features of the class name through a word embedding model~\cite{wang2019survey}. While sufficient for most zero-shot applications a lot of context regarding the class is lost.  In order to effectively utilize language-vision models it is necessary to align the downstream task to the training~\cite{radford2021learning}.
As such we model a set of positive and negative prompts applicable for pathologies to enrich our matching process between visual and textual projections.
While for our basic approach we consider (`\{class\}', `No \{class\}') prompts, we found that a more detailed prompt design can overall deliver improved performance.  As such we consider a set of prompts following the templates (`\{adverb\}\{indication\_verb\} \{effect\}$^*$ \{location\}$^*$ \{class\_synonym\}', `\{adverb\} \{indication\_verb\} \{absence\} \{class\_synonym\}'). Hereby, we utilize all combinations over a small set of categories to gather a variety of different settings. During inference, features of all queries of the same set are averaged.

\noindent\textbf{Prompt-based Dataset Extension: } Despite medical reports being the more common resource in the practical field, currently the majority of large-scale datasets are only publically available with fixed sets of labels. In order to investigate the effect of additional data in training of our method, we reverse our proposed prompt engineering to generate synthetic reports for the datasets PadChest and ChestX-Ray14 based on their class-labels.
Through this procedure, we are able to sample sentences indicating presence or absence of a class and generate more than 200k added image-report pairs.

\section{Experiments}
\subsection{Experimental Setup}
\noindent\textbf{Datasets:}
\begin{itemize}
    \item \textbf{MIMIC-CXR: } It contains 377,110
        chest X-rays taken from 65179 patients with 14 disease labels and 227,835 reports.
        We use the splits provided by~\cite{johnson2019mimic}. Unless further specified all models were trained on this dataset.
\item \textbf{CheXpert: } It contains 224,316
    chest X-rays taken from 65,240 patients with 14 disease labels. The labels are shared with MIMIC-CXR.
    We only consider the validation split provided by~\cite{irvin2019chexpert}.
\item \textbf{ChestX-ray14: } It contains 112,120
frontal-view chest X-rays taken from 30,805 patients with 14 disease labels.
    We use the splits provided by~\cite{seibold2020self}.
    \item \textbf{PadChest: } It consists  160k
    chest X-rays of 67k patients with 174 findings.
\end{itemize}
\noindent\textbf{Evaluation Setup: } We evaluate the multi-label classification ability of all networks via the Area Under the ROC-curve (AUROC) and show the performance over MIMIC-CXR, CheXpert and ChestX-Ray14. For all experiments expect Table~\ref{tab:final} we consider validation performance. Labels with value -2 and -1 are ignored for the calculation of the metric as their state is not certain. For all ablations, we use the "basic"-prompting scheme, while for further experiments the "detailed"-scheme is used.

\noindent\textbf{Implementation Details: } For all experiments we use the same ResNet50 and Transformer as Redford~et~al.~\cite{radford2021learning} as backbones. We optimize with AdamW~\cite{loshchilov2017decoupled}, a learning rate of 0.0001 and a cosine schedule. We trained classification models with a learning rate of 0.0005 as this has shown slightly better performance. During training, we resize the images to the inference size of 320 × 320 and randomly crop by 288x288. For specific further augmentations, we follow SimSiam~\cite{chen2021exploring}.

\subsection{Results}
\noindent\textbf{Ablation - Effect of Heads:}
We investigate the impact of both prediction heads during inference. We start by showing the individual head performance and then go over to different fusion approaches on the left of Table\ref{tab:head}. For feature fusion we consider the concatenation of local and global features of the same modality. For score-fusions, we calculate scores as described above and aggregate the class predictions based on their maximum or average.

We see that for our method both the global and local head show nearly the same performance. While performing max-score fusion the across-dataset-performance drops by 0.4\%, where mean-score fusion improves by 0.6\%.
\input{figures/data_performance_plot}
\input{tables/mimic}

\noindent\textbf{Ablation - Effect of Losses:}
On the right of Table\ref{tab:head}, we show the impact of the objective functions. We see that adding local contrast improves the model by 2\%. Adding the self-supervised and mirrored objective worsen performance by 0.45\%. It can be noted that the self-supervised loss achieves the best performance on ChestX-Ray14 by more than 2\%.
 whereas adding both simultaneously improves across dataset performance by 0.6\%.

\noindent\textbf{Multi-label Inference and Prompt Engineering:} We show the impact of our proposed inference scheme and prompts in Figure~\ref{fig:prompts}. We see that performance overall improves with significant improvements for some classes such as fractures which were unable to be categorized just using cosine-similarity. When the detailed prompt the mean performance further improves.

\noindent\textbf{Data size Impact: }
 We show the impact of using additional prompt-based reports during training  in Figure~\ref{fig:data}. We see that including artificial training data for ChestX-Ray14 significantly improves its validation performance.  In general it seems that while for some classes performance seems to worsen, the overall performance improves when adding additional data.

\noindent\textbf{Comparison with Other Approaches: }
We compare against the same vision network trained with label supervision on its respective dataset. All other methods were trained using the MIMIC-CXR dataset. SLIP~\cite{mu2021slip} refers to a version of CLIP, which incorporates self-supervision in form of a SIMCLR-like objective~\cite{chen2020simple}. $MILNCE_{local}$ refers to our local branch trained with the MILNCE~\cite{miech2020end} objective alone. LoCo and GloCo refer to our method trained with either just the local or global objective respectively. We evaluate using the "detailed"-prompt scheme.  We show the results in Table~\ref{tab:final}.

We see that our formulation of the local contrastive loss outperforms the MILNCE version across all datasets. Our proposed method outperforms the considered contrastive language-image pretraining baselines in the form of CLIP and SLIP and manages to achieve similar performance as the supervised ResNet for domains similar to MIMIC, however, underperforms for the ChestX-Ray14 dataset. When adding the additional report datasets of PadChest and ChestX-Ray14 we manage to beat label-supervised performance across all datasets.

\section{Conclusion}
In this paper, we proposed an approach to make networks less reliant to label supervision through contrastive language-image pre-training on report level. In order to still maintain competitive levels of performance we introduced a novel way of constructing inference. Doing so we are able to offset issues stemming from explicit class similarities.
We show that despite using unstructured medical report supervision, we perform on par with explicit label supervision through a sophisticated inference setting across different datasets.

\section{Acknowledgements}
The present contribution is supported by the Helmholtz Association under the joint research school “HIDSS4Health – Helmholtz Information and Data Science School for Health” and by the Helmholtz Association Initiative and
Networking Fund on the HAICORE@KIT partition.


\section{Introduction.}

\setcounter{equation}{0} \renewcommand{\theequation}{\arabic{section}
\arabic{equation}}

The first mathematically rigorous studies of infinite statistical systems were the works of Bogolyubov \cite {Bog46}
in the canonical ensemble and Ruelle \cite{Ru63, Ru69} in the grand canonical ensemble.
In the monograph \cite {Bog46} Bogolyubov developed a theory
of expansions in degrees of density, which, although not yet rigorously mathematically justified,
but was based on the fundamental principles of statistical mechanics and included all
known at that time expressions of Ursel-Mayer theory on the basis of his derived integro-differential
equations for $ m $ -particle distribution functions. Really rigorous result was formulated in \cite{BogKh49}
(see details in \cite{Kh56-1, Kh56-2}), but his description covered only non-negative
interaction potentials of the particles.
It was only in 1969 that Bogolyubov, Petrina, and Khatset \cite{BPH69} succeeded in generalizing these results
to the case of stable regular potentials.
Following Ruelle's operator method, they derived a recurrent relation, which, in the thermodynamic limit,
 turned into the Kirkwood-Salzburg equation.
But the role of activity $z$  is played  some function of the density $ a(\varrho) $, which is determined
by series in the set of all distribution functions.
This essentially means that the equation is an infinite system of nonlinear equations,
and its iterations are not a series expansion in the density.
A little later, these equations were derived by Pogorelov \cite {Po75} directly
from Bogolyubov's integro-differential
equations (see \cite {Bog46}) and a theorem  on the existence of a unique solution
(as of nonlinear operator equation) for small values of density  was proved.

The virial expansion is most often understood as the expansion of pressure in terms of the density \cite{MaMa40}
and obtained, by formally inverting the density as a power series of the fugacity and
composing the latter with the expansion for the pressure (Mayer's   power series in the fugacity).
In addition, we work with expansions of correlation functions in the thermodynamic limit.
Existence of a thermodynamic limit for  correlation functions in a canonical ensemble
"is more subtle than for thermodynamic quantities like pressure and free energy which are on a logarithmic
scale" \cite{KuTs18}. One more thing, the terms of the Mayer's  cluster expansions can be expressed in terms of forest-graphs
(see, i.e. \cite{ Ru69, MP77, Re2021}, while those of the virial expansion are usually
written in terms of 2-connected diagrams (see, for instance, the classical reference \cite{UF63}).

In this article we construct first an expansion for correlation functions by some special type of forest-graphs,
 but with huge number of terms in every order of density. Due to the fact that this series consists of two parts
 that have opposite signs, it is possible to reduce the contributions and, as a result, obtain its representation by
 some set of graphs, connected components of which are 2-connected graphs. Such set of graphs was also defined
in the article of Jansen \cite{Ja21} (see sect. 2).

The well-known proofs of convergence for the virial expansion are based on the above mentioned
inversion. Obviously the first result goes back to the work of Lebowitz and Penrose \cite{LePe64}
and some alternative approaches were given by Groeneveld \cite{Gro67}  (see also \cite{PT12, RT15, NF20},
which was continued in a more general and modern presentation in the work \cite{Ja21}.
The present work apparently has quite a few points of contact with the work \cite{Ja21}.

The article is organized as follows. In Section 2 we give a short mathematical introduction of
concepts and formulas in which the presentation of the article is carried out.
In Section 3, we briefly show the derivation of equations for correlation functions in the canonical ensemble.
The fourth section is devoted to the representation of the solution in the form of the Lebesgue-Poisson integral
of some kernel, for which we write the recurrent nonlinear equation. And finally, in the last section,
the question of its solution in the form of an expansion in terms of powers of density and
the convergence of the integral representation for correlation functions are discussed.

\section{On some mathematical concepts}

\subsection{Configuration spaces and Lebesgue-Poisson measure.}
Let $ \sigma $ be a Lebesgue measure in $ \R^d $. {\it Configuration space} $ \Ga \;: = \; \Ga_{\R^d} $
consists of all locally finite subsets of the space $\R^d $.
This definition is quite natural from the point of view of physics, because
an infinite number of particles cannot be in a bounded volume.

Denote the set of all finite configurations of the space  $\Ga $ by $ \Ga_0 $, and by $ \Ga^{(n)} $
configuration space with a fixed number of points.

If all such configurations are in some bounded set $ \La\subset\R^d $,
then the corresponding space will be $ \Ga^{(n)}_\La $.

State of an ideal gas in equilibrium statistical mechanics
is described by the {\it Poisson} measure  $ \pi_{z \sigma} $ on
configuration space $ \Ga $, where $ z> 0 $ is the activity (physical parameter that
associated with the density of particles in the system).

But we need  so-called Lebesgue-Poisson measure, i.e. the Poisson measure
on the space of finite configurations, which in some sense is absolutely continuous
with respect  to the Lebesgue measure in contrast to the $ \pi_{\sigma} $ measure for $ \Ga $. Integral for
Lebesgue-Poisson measure in space
$ \Ga_0 $ (or $ \Ga_\Lambda $) is determined by the formula:
\begin{align}\label{IL-P
\int_{\Ga_X}F(\ga)\la_{z\si}(d\ga)\;&:=\;\sum_{n=0}^\infty\frac{z^n}{n!}\int_X\cdots\int_X
F(\{x_1,...,x_n\})\si(dx_1)\cdots\si(dx_n)=\notag\\
&=\sum_{n=0}^\infty\frac{z^n}{n!}\int_X\cdots\int_X
F_n(x_1,...,x_n)dx_1\cdots dx_n,
\end{align}
for all measurable functions $ F = \{F_n \}_{n \geq 0} $, $ F_n \in L^{1}(X^n) $, and
  $ \Gamma_X \in \{\Gamma_0, \Gamma_\Lambda $).
  Such  measure not only simplifies the recording of cumbersome expressions and combinatorial proofs,
  but also allows  to get some new results based on properties
  infinite divisibility (see, for example, \cite{Re93, ReSh97, Re98, PR07}.

In the language of integrals with the measure $\la_{z\sigma}$ we give one well-known and important identity,
which will be used very often below (see, for example, \cite{Ku99}).

\begin{lemma}\label{lem2}
 For all measurable functions $G:\Ga_0\mapsto \R$ and $H: \Ga_0\times\Ga_0\mapsto \R$,  for which
 $G(\xi\cup\ga)H(\xi,\ga)\in L^1(\Ga_0\times\Ga_0, \la_{\si}\otimes\la_\si) $, the following equality is true:
\begin{equation}\label{EgLem2
\int_{\Ga_0}G(\ga)\sum_{\xi\subseteq\ga}H(\xi,\ga\sm\xi)\la_\si(d\ga)=\int_{\Ga_0}\int_{\Ga_0}G(\xi\cup\ga)H(\xi,\ga)\la_{\si}(d\ga)\la_\si(d\xi).
 \end{equation}
\end{lemma}
{\it Proof}.
 Let
 $\xi\restriction\Ga_0^{(k)}=\{x_1,...,x_k\}:=\{x\}_1^k$,
 $\ga\restriction\Ga_0^{(m)}=\{x_{k+1},...,x_{k+m}\}:=\{x\}_{k+1}^{k+m}$. Then
\begin{align*}
&\int_{\Ga_0}\int_{\Ga_0}G(\xi\cup\ga)H(\xi,\ga)\la_{\si}(d\ga)\la_\si(d\xi)=\\
&=\sum_{k,m=0}^\infty\frac{1}{k!m!}\int_{\R^{dk}}\int_{\R^{dm}}
G(\{x\}_1^{k+m})H(\{x\}_1^k,\{x\}_{k+1}^{k+m})\si(dx)^{k+m}=\\
&=\sum_{n=0}^\infty\frac{1}{n!}\sum_{k=0}^n {n \choose k}\int_{\R^{dn}}
G(\{x\}_1^{n})H(\{x\}_1^k,\{x\}_{k+1}^{n})\si(dx)^{n}=\\
&=\int_{\Ga_0}G(\ga)\sum_{\xi\subseteq\ga}H(\xi,\ga\sm\eta)\la_\si(d\ga).
\end{align*}
\hfill $\blacksquare$

Note that the introduction of the measure \eqref{IL-P} and the application of the lemma \ref{lem2} greatly simplifies the calculations of many combinatorial statements.

\subsection{Algebra of generating functionals.}

For any function
$\psi\in C_0(\Gamma_0)$ define  {\it generating functional }:
\begin{align}\label{GFf
\w{F}_\psi(j)&\;=\;\sum_{N=0}^{\infty}\frac{1}{N!}\int_{\R^{dN}}(dx)^{N}j(x_1)\cdots j(x_N)\psi(x_1,\ldots, x_N)=\notag \\
&\;=\;\int_{\Gamma_0}e(j;\gamma)\psi(\gamma)\lambda_\sigma(d\gamma),
\end{align}
where functional $e(j;\gamma)$ is defined by the formula
\begin{align}\label{I-lam4
e(j;\gamma):=\;&\begin{cases}
 1,\;&\gamma=\emptyset,\\
  \prod_{x\in\gamma}j(x),
\;&\gamma\in\Ga_0\sm\{\emptyset\},
  \end{cases}
\end{align}

 and $j\in C_0(\R^d)$ is bounded (take for simplicity  $j\leq 1$), continuum, nonnegative  function.

It is easy to calculate using the lemma \ref{lem2} that
\begin{equation}\label{odprod
  \w{F}_{\psi_1}(j)\w{F}_{\psi_2}(j)\;=\;\w{F}_{\psi_1*\psi_2}(j)=\int_{\Gamma_0}e(j;\gamma)(\psi_1*\psi_2)(\gamma)\lambda_\sigma(d\gamma),
\end{equation}
where $\psi_1*\psi_2$ is product $\psi_1$ and $\psi_2$ in commutative algebra $\mathcal{A}$,
which was introduced by Ruelle (see \cite{Ru69}, Ch. 4). Formula \eqref{odprod}
establishes an unambiguous correspondence between the algebra of generating functionals with operation of ordinary multiplication
and the algebra $ \mathcal {A} $. Recall that for
$\gamma=\{x_1, x_2, ... x_n\}$
\begin{equation}\label{odprodRu
(\psi_1*\psi_2)(\gamma)= \sum_{\xi\subseteq\gamma}\psi_{1}(\xi)\psi_{2}(\gamma\setminus\xi).
\end{equation}

\section{Description of the system of interacting particles in the canonical ensemble.}

\subsection{Interaction between particles}
We  consider the 2-particle interaction, which is described by the potential
$V_2(x,y)=\phi(|x-y|),\;\phi(0)=+\infty$. For any configuration $\ga\in\Ga_0$
an interaction energy  is
\begin{equation}\label{U
U(\gamma)=U_{\phi}(\gamma):=\sum_{\{x,y\}\subset\gamma}\phi(|x-y|),
\end{equation}
and interaction $W(\eta; \gamma)$ between particles of configuration $\eta\in\Ga_0$ with particles of
configuration $\ga\in\Ga_0$ is
\begin{equation}\label{W
W(\eta; \gamma):=\sum_{\substack{x\in\eta \\
y\in\gamma}}\phi(|x-y|).
\end{equation}

We impose the following conditions on the potential of interaction:

{\bf (A)}: 1. {\it Stability:}
\begin{equation}\label{stU
U(\gamma)\;\geq\;-B|\gamma|,\;B\,\geq 0,\;\gamma\in\Gamma_0,
\end{equation}
\qquad  2. {\it  Regularity:}
\begin{equation}\label{reg
C(\beta)\;=\;\int_{\R^d}dx|e^{-\beta\phi(|x|)}\;-\;1|\;<\;+\infty, \beta=\frac{1}{kT}.
\end{equation}

\subsection{Derivation of equations for correlation functions.}

The correlation functions of the canonical ensemble of the system $N$ of particles in some limited volume $\Lambda$ in configuration
$\eta=\{x_1,\ldots, x_m\}$ are determined by the sequence

\begin{align}\label{corfCA-2
\rho^{(N)}_{\La}(\eta)\;:=\;\begin{cases}
 & 1,\;\;\text{at}\; \eta = \emptyset,\\
 & \varrho=\frac{N}{V}\;\;\text{at}\;\; |\eta|=1, V=\sigma(\Lambda),\\
 & 0,\;\text{at};|\eta|=m\;>\; N,
  \end{cases}
\end{align}
and for $2\leq m\leq N$,
\begin{equation}\label{corfCA-3
\rho^{(N)}_{\La}(\eta)\;:=\;\frac{1}{Z_\Lambda(\beta, N)}\int_{\Ga^{(N-m)}_\La}\frac{d\gamma}{(N-m)!}
 e^{-\beta U(\eta\cup\ga)},
\end{equation}

where
\begin{equation}\label{SSCA
Z_{\La}(\beta, N)\;=\;\frac{1}{N!}\int_{\Ga^{(N)}_\La}
e^{-\beta U(\ga)}d\gamma,
\end{equation}
and $\varrho$ is the density of the particles in the system.

To derive equation for correlation functions using lemma \ref{lem2}
we write down the right hand part of \eqref{corfCA-3} by the integral \eqref{IL-P}:
\begin{equation}\label{corfCA1
\rho^{(N)}_{\La}(\eta)\;=\;\frac{1}{Z_{\La}(\beta, N)}\int_{\Ga_\La}
e^{-\beta U(\eta\cup\ga)}\1_{N-m}(\gamma)\lambda_\sigma(d\gamma), \;|\eta|\leq N,
\end{equation}
where
\begin{align}\label{unit
\1_{N-m}(\gamma)\;:=\;\begin{cases}
 &1,\; |\gamma|=N-m,\\
 & 0,\;|\gamma|\neq N-m.
  \end{cases}
\end{align}

Introduce a notation $\eta^{(\hat{x})}:=\eta\sm\{x\}$. Let in configuration $\eta=\{x_1,..., x_m\}$ the coordinate
$x_1$ such that, for which the following inequality is true:
\begin{equation}\label{St
  W(x_1; \eta^{(\hat{x}_1)})=\sum_{k=2}^{m}\phi(|x_1-x_k|)\;\geq\;-B,\;B\geq 0 .
\end{equation}
Then
\begin{equation}\label{SSCA
  U(\eta\cup\ga)=W(x_1; \eta^{(\hat{x}_1)})+W(x_1; \gamma)+U(\eta^{(\hat{x}_1)}\cup\ga)
\end{equation}
and
\begin{equation}\label{corf6
 e^{-\be W(x_1;\ga)}=\prod_{y\in\ga}\left[(e^{-\be \phi(|x_1-y|)}-1)+1\right]=\sum_{\xi\subseteq\ga}K(x_1;\xi),
\end{equation}
where
\begin{align}\label{corf7
K(x_1;\xi)\;:=\;\begin{cases}
 &\prod_{y\in\xi}(e^{-\be \phi(|x_1-y|)}-1),|\xi|\geq 1,\\
 & 1,\;\xi=\emptyset.
  \end{cases}
\end{align}
Then
\begin{align}\label{corfCA12-1
\rho^{(N)}_{\La}(\eta)\;=\;&\frac{e^{-\beta W(x_1; \eta^{(\hat{x}_1)})}}{Z_{\La}(\beta, N)}\int_{\Ga_\La}
e^{-\beta U(\eta^{(\hat{x}_1)}\cup\ga)} \times\notag\\
&\times\1_{N-m}(\gamma)\sum_{\xi\subseteq\ga}K(x_1;\xi)\1_{\Lambda}(\gamma\setminus\xi)\lambda_\sigma(d\gamma).
\end{align}
Apply to the integral in  \eqref{corfCA12-1} the lemma \ref{lem2} with
\begin{equation}\label{GH1
   G(\gamma)=e^{-\beta U(\eta^{(\hat{x}_1)}\cup\ga)}\1_{N-m}(\gamma)\;\;\text{and}\;\;H(\xi;\gamma\setminus\xi)= K(x_1;\xi)\1_{\Lambda}(\gamma\setminus\xi).
\end{equation}
Then
\begin{align}\label{corfCA14
\rho^{(N)}_{\La}(\eta)\;&=\;\frac{e^{-\beta W(x_1; \eta^{(\hat{x}_1)})}}{Z_{\La}(\beta, N)}\int_{\Ga_\La}\int_{\Ga_\La}
e^{-\beta U(\eta^{(\hat{x}_1)}\cup\xi\cup\ga)} \times\notag\\
&\times\1_{N-m}(\xi\cup\gamma)K(x_1;\xi)\1_{\Lambda}(\gamma)\lambda_\sigma(d\gamma)\lambda_\sigma(d\xi).
\end{align}
Rewrite the right hand side in the following way
\begin{align}\label{corfCA13
&\rho^{(N)}_{\La}(\eta)=e^{-\beta W(x_1; \eta^{(\hat{x}_1)})}\frac{Z_{\La}(\beta, N-1)}{Z_{\La}(\beta, N)}\int_{\Ga_\La}K(x_1;\xi)\times\\
&\times\left[\frac{1}{Z_{\La}(\beta, N-1)}\int_{\Ga_\La}e^{-\beta U(\eta^{(\hat{x}_1)}\cup\xi\cup\ga)} \1_{N-m-|\xi|}(\gamma) \lambda_\sigma(d\gamma)\right]\lambda_\sigma(d\xi).\notag
\end{align}
Taking into account \eqref{corfCA1}($N-m-|\xi|= N-1-(m-1+|\xi|)$) we obtain:
\begin{equation}\label{corfCA123
 \rho^{(N)}_{\La}(\eta)= e^{-\beta W(x_1; \eta^{(\hat{x}_1)})} a_N^\Lambda
 \int_{\Ga_\La}K(x_1;\xi)\rho^{(N-1)}_{\La}(\eta^{(\hat{x}_1)}\cup\xi)\lambda_\sigma(d\xi),
\end{equation}
where $|\eta|=m$\;\;$\eta=\{x_1,..., x_m\}$, $1<m<N$, and
\begin{equation}\label{coef
 a_N^\Lambda\;=\; \frac{Z_{\La}(\beta, N-1)}{Z_{\La}(\beta, N)}.
\end{equation}

For  $m=1$
\begin{equation}\label{I-1
 \rho^{(N)}_{\La}(\{x_1\})=\frac{N}{V}=a_N^\Lambda
 \int_{\Ga_\La}K(x_1;\xi)\rho^{(N-1)}_{\La}(\xi)\lambda_{\sigma}(d\xi),\;\;V=\sigma(\Lambda).
\end{equation}

Note, also, that
\begin{equation}\label{I-2-3
\rho^{(N)}_\La(\eta)= 0, \;\;\text{if}\;\;  |\eta|\;\geq\;N.
\end{equation}
In the thermodynamic limit
\begin{equation}\label{lim
\Lambda\uparrow\R^d,\;\;\lim_{\substack{V\rightarrow\infty \\
N\rightarrow\infty}}\frac{V}{N}=v,\;\;\frac{1}{v}=\varrho,
\end{equation}
 existing of which was proved in \cite{BPH69} we have
the following equation
\begin{equation}\label{limeq
 \rho(\eta)\;=\;\varrho a(\varrho) e^{-\beta W(x_1; \eta^{(\hat{x}_1)})}
 \int_{\Ga_0}K(x_1;\xi)\rho(\eta^{(\hat{x}_1)}\cup\xi)\lambda_{\sigma}(d\xi),
\end{equation}
where
\begin{equation}\label{ratioPF
\lim_{\substack{V=N/\varrho \\
N\rightarrow\infty}}\frac{Z_\Lambda(\beta, N-1)}{Z_\Lambda(\beta, N)}\;=\;\varrho a(\varrho),\;\;V=\sigma(\Lambda)
\end{equation}
and
\begin{equation}\label{ratioPF1
 a(\varrho)=\frac{1}{\w{Q}(\rho)}, \;\w{Q}(\rho)\;=\;\int_{\Ga_0}K(x_1;\xi)\rho(\xi)\lambda_{\sigma}(d\xi).
\end{equation}

\section{ Construction of expansion for correlation functions in the density parameter $ \varrho $.}

Solution of the equation \eqref {limeq} by the iteration method (Ruelle operator method \cite {Ru69}) or
by the Minlos-Poghosyan method \cite {MP77} (about the connection between these methods see \cite {Re2021})
leads to the series in the parameter $ \varrho a(\varrho) $. So, in fact, these equations are nonlinear,
because $ a(\varrho) $ depends on all correlation functions (see \eqref{limeq} and \eqref{ratioPF1}).
In this section, we construct expansion in the  parameter density $ \varrho $, considering the equation \eqref{limeq}
as nonlinear equation.

\subsection{Solution of equations for correlation functions.}

Rewrite the equation \eqref{limeq}, \eqref{ratioPF1} as follows:

\begin{equation}\label{limeq-1
 \w{Q}(\rho_j)\rho_j(\eta)\;=\;\varrho  e^{-\beta W(x_1; \eta^{(\hat{x}_1)})}
 j(x_1)\int_{\Ga_0}K(x_1;\xi)\rho_j(\eta^{(\hat{x}_1)}\cup\xi)\lambda_{\sigma}(d\xi),
\end{equation}
\begin{equation}\label{ratioPF1-1
 \w{Q}(\rho_j)\;=\;\int_{\Ga_0}K(x_1;\xi)\rho_j(\xi)\lambda_{\sigma}(d\xi),
\end{equation}
where, following \cite{MP77}, a continuous  positive function
$ j: \R^d \rightarrow \R_+ $ is introduced and $ \rho (\eta) $ becomes functional
such as $\rho(\eta)\;=\;\rho_j(\eta)|_{j=1}.$

We will look for the solution of these equations in the form:
\begin{equation}\label{SSCA-2
\rho_j(\eta)\;=\;e(j;\eta)\int_{\Gamma^0}T(\eta|\gamma)e(j;\ga)\lambda_\sigma(d\gamma),
\end{equation}
where $e(j;\eta)$ is given by \eqref{I-lam4}.

The idea of finding a solution is to consider the left side of the equation \eqref {limeq-1} as the product of two generating
functionals (see \eqref{odprod}) and compare it with the right part, which is also written in the form of a generating functional
of the form \eqref {GFf}. So, substituting in \eqref {ratioPF1-1} the corresponding expression for $ \rho_j (\xi) $ we have:
\begin{equation}\label{ratioPF1-1a
 \w{Q}(\rho_j)\;=\;\int_{\Ga_0}\int_{\Ga_0}e(j;\xi\cup\ga)K(x_1;\xi)T(\xi|\gamma)\lambda_{\sigma}(d\gamma)\lambda_{\sigma}(d\xi).
\end{equation}
Apply to the integral on the right side of the \eqref {ratioPF1-1a} lemma \ref {lem2} with
\begin{equation}\label{GH1
   G(\xi\cup\gamma)=e(j;\xi\cup\ga)\;\;\text{and}\;\;H(\xi;\gamma)= K(x_1;\xi)T(\xi|\gamma).
\end{equation}
Then
\begin{equation}\label{ratioPF1-2
 \w{Q}(\rho_j)\;=\;\w{F}_{\psi_1}(j)\;=\;\int_{\Ga_0} e(j;\ga)\psi_1(\gamma)\lambda_{\sigma}(d\gamma),
\end{equation}
with
\begin{equation}\label{ratioPF1-2a
 \psi_1(\gamma)\;=\;\sum_{\xi\subseteq\gamma}K(x_1;\xi)T(\xi|\gamma\setminus\xi).
\end{equation}
Write down  \eqref{SSCA-2} in the form of generating functional:
\begin{equation}\label{SSCA-3
\rho_j(\eta)\;=\; \w{F}_{\psi_2}(j)\;=\;\int_{\Ga_0} e(j;\ga)\psi_2(\gamma)\lambda_{\sigma}(d\gamma),
\end{equation}
with
\begin{equation}\label{ratioPF1-2b
 \psi_2(\gamma)\;=\;e(j;\eta) T(\eta|\gamma).
\end{equation}
Left hand side \eqref{limeq-1} is:
\begin{equation}\label{limeq-1-left
 \w{Q}(\rho_j)\rho_j(\eta)\;=\; \w{F}_{\psi_1}(j)\w{F}_{\psi_2}(j)\;=\;\w{F}_{\psi_1*\psi_2}(j)\;
 =\;\int_{\Ga_0} e(j;\ga)(\psi_1*\psi_2)(\gamma)\lambda_{\sigma}(d\gamma),
\end{equation}
where
\begin{align}\label{help}
(\psi_1*\psi_2)(\gamma)&=\sum_{\xi\subseteq\gamma}\psi_1(\xi)\psi_2(\gamma\setminus\xi)=\sum_{\xi\subseteq\gamma}\psi_2(\gamma\setminus\xi)\psi_1(\xi)=\\
&=e(j;\eta)\sum_{\xi\subseteq\gamma} T(\eta|\gamma\setminus\xi)\sum_{\zeta\subseteq\xi}K(x_1;\zeta)T(\zeta|\xi\setminus\zeta)\notag
\end{align}

Right hand side \eqref{limeq-1} after substitution of $\rho_j(\eta^{(\hat{x}_1)}\cup\xi)$
in the form \eqref{SSCA-2} will look as:
\begin{equation}\label{limeq-1-right0
\varrho  e^{-\beta W(x_1; \eta^{(\hat{x}_1)})}
 e(j;\eta)\int_{\Ga_0}\int_{\Ga_0}e(j;\xi\cup\ga)K(x_1;\xi)T(\eta^{(\hat{x}_1)}\cup\xi|\gamma)\lambda_{\sigma}(d\gamma)
\lambda_{\sigma}(d\xi).
\end{equation}
Let us use lemma \ref{lem2} again with
\begin{equation}\label{GH1
   G(\xi\cup\gamma)=e(j;\xi\cup\ga)\;\;\text{and}\;\;H(\xi;\gamma)= K(x_1;\xi)T(\eta^{(\hat{x}_1)}\cup\xi|\gamma).
\end{equation}
Then right hand side of \eqref{limeq-1} will be as
\begin{equation}\label{limeq-1-right
\varrho  e^{-\beta W(x_1; \eta^{(\hat{x}_1)})}
 e(j;\eta)\int_{\Ga_0} e(j;\ga)\sum_{\xi\subseteq\gamma} K(x_1;\xi)T(\eta^{(\hat{x}_1)}\cup\xi|\gamma\setminus\xi)\lambda_{\sigma}(d\gamma).
\end{equation}
We equate the functions near $ e(j;\ga) $ under the integrals of the expressions \eqref {limeq-1-left} and \eqref {limeq-1-right} and, reducing the multiplier
   $ e(j; \eta) $ we get the relation:
\begin{align}\label{help}
 &\sum_{\xi\subseteq\gamma} T(\eta|\gamma\setminus\xi)\sum_{\zeta\subseteq\xi}K(x_1;\zeta)T(\zeta|\xi\setminus\zeta)\;=\\
 &=\;\varrho  e^{-\beta W(x_1; \eta^{(\hat{x}_1)})}\sum_{\xi\subseteq\gamma} K(x_1;\xi)T(\eta^{(\hat{x}_1)}\cup\xi|\gamma\setminus\xi).\notag
\end{align}
Select the term with $ \xi = \emptyset $ in the left side of this relation and take into account that $ T(\emptyset | \emptyset) = 1 $ and $ T (\emptyset | \gamma) = 0 $,
if $ \gamma \neq \emptyset) $, we obtain the following nonlinear recurrent relation:
\begin{align}\label{final}
T(\eta|\gamma)&= \varrho  e^{-\beta W(x_1; \eta^{(\hat{x}_1)})}\sum_{\xi\subseteq\gamma} K(x_1;\xi)T(\eta^{(\hat{x}_1)}\cup\xi|\gamma\setminus\xi)-\notag\\
&-\sum_{\emptyset\neq\xi\subseteq\gamma} T(\eta|\gamma\setminus\xi)\sum_{\emptyset\neq\zeta\subseteq\xi}K(x_1;\zeta)T(\zeta|\xi\setminus\zeta).
\end{align}

{\bf  Initial conditions:}
\begin{equation}\label{incon0
    T(\emptyset|\emptyset)\;=\;1,\;\;\;T(\emptyset|\gamma)\;=\;0, \text{if}\;\;\gamma\neq\emptyset,
\end{equation}
\begin{equation}\label{incon1
    T(\{\eta|\gamma)\;=\;0,\text{if}\;\;\eta\cap\gamma\neq\;\emptyset.
\end{equation}

\subsection{Graphic interpretation of solutions.}

 As in the case of the Kirkwood-Salzburg equations for a large canonical ensemble (see
   details in \cite {Re2021}) the analytical expression of the kernel $ T (\eta|\gamma) $ can be given as
  contributions of forest-graphs.   The number $ N (m | n), m = |\eta|, n = |\gamma| $
all forest graphs can be calculated by solving the following recurrent equation:
\begin{equation}\label{Numberforests-vir
N(m|n)=\sum_{k=0}^n {n \choose k}N(m+k-1|n-k)+\sum_{k=1}^n {n \choose k}N(m|n-k)\sum_{l=1}^k {k \choose l}N(l|k-l).
\end{equation}
If we reject the nonlinear term, the recurrent equation can be solved exactly: $ N (m | n) = m (m + n)^{n-1} $ (see \cite {Re2021}).
This formula makes it possible to prove the convergence of the integral \eqref {SSCA-2} at $ j = 1 $ and small density values.
But the nonlinear term apparently significantly increases the number of forest graphs. But according to the formula \eqref {final},
each contribution of the forest graph, which contains a tree with the contribution of
the factor $ K (x_1; \xi), \xi \in \gamma $, which appears from the second sum in \eqref {final} is included with a minus sign. So, for example calculating
\begin{equation}\label{T11
    T(x_1|y)\;=\; \varrho^2\left[ K(x_1;y)\;-\;K(x_1;y)\right]=0.
\end{equation}
Therefore, it is necessary to analyze this recurrent relationship in more detail.
It is easy to get the following proposition.

\begin{prop}\label{dens}
From initial conditions \eqref{incon0}-\eqref{incon1} it follows that
\begin{equation}\label{den1
   T(x_1|\emptyset)\;=\;\varrho,\;\;T(x_1|\gamma)\;=\;0\;\text{for}\;\gamma\neq\emptyset.
\end{equation}
It gives the needed equality for density $\rho(x_1)\;=\;\varrho$. \\
In addition   $T(\eta|\emptyset)\;=\;\varrho^{|\eta|}\exp[-\beta U(\eta)]$.
\end{prop}
The {\it proof} follows from \eqref{final} by induction procedure.
To obtain the solution of \eqref{final} in the language of graph theory \cite {Ste64, Gro67}
define the set graphs $\mathcal{D}(\eta; \gamma)$.

\begin{remark}\label{graph}
The set  $\mathcal{D}(\eta; \gamma)$ is the set of rooted graph,
connected components of which are 2-connected graphs with respect to configuration $\gamma$.
Any graph $G\in\mathcal{D}(\eta; \gamma) $ has $m=|\eta|$ vertices of configuration $\eta$ and $n=|\gamma|$ vertices of
configuration $\gamma$. Every vertex of $\eta$ can be free of lines or can connect only with vertices $y$ of configuration
$\gamma$. Every vertex $y\in\gamma$ connects with vertices $x$ or $y$ by at least two lines.
This definitions is very close to definitions of the set of white and black vertices $\mathcal{D}(W; B)$ in \cite {Ja21}, Section 2.
The role of the roots of the graph is performed by the configuration $\eta$.
A singleton $G\in\mathcal{D}(x_1; \emptyset)$.  The contribution of any vertex is $\varrho$. The contribution of line $w(l_{xy})$,
which connected two vertices $x$ and $y$ or $y_i$ and $y_j$ is $K(x;y)$.
\end{remark}

\begin{remark}\label{graph1}
According to \eqref{SSCA-2}, integration is performed at the points of configuration $\gamma$, and according to equation \eqref{final},
 in the kernels $T(\theta; \xi)$, the configuration $\xi\subseteq\gamma$, and the configuration $\theta$ can contain both variables
 from $\eta$ and $\gamma$. Therefore, the corresponding graphs of sets $\mathcal{D}(\theta; \xi)$ can contain only black vertices.
\end{remark}

The main statement  of the section is

\begin{prop}\label{dens}
The solution of  the recurrence relation \eqref{final}  can be rewritten by contributions from the graphs of the set $\mathcal{D}(\eta; \gamma)$:
\begin{equation}\label{den1
   T(\eta|\gamma)\;=\;\varrho^{|\eta|+|\gamma|}\exp[-\beta U(\eta)] \sum_{G\in\mathcal{D}(\eta; \gamma) }w_G(L; \eta, \gamma),
\end{equation}
where $L$ is a configuration of lines(edges).
\end{prop}
The main points of {\it proof}. The proof is based on the induction procedure.
It is very easy to see the cancellation of the contributions of the corresponding graphs, for example, for kernels
$T(x_1, x_2|y_1)$ and $T(x_1, x_2|y_1, y_2)$. To facilitate these calculations, it is better to consider instead of the equation
\eqref{final} the equation for kernels that are integrated in variables $y_1,..., y_n$:

\begin{equation}\label{kern2
   \widehat{T}(\eta|n):=\int_{\Gamma^{(n)}} T(\eta|\gamma)d\gamma =\int_{\R^{dn}} (dy)^nT(\eta|\{y\}_1^n),
\end{equation}
especially since this is the value that must be substituted into the formula for the solution \eqref{SSCA-2}.
\begin{align}\label{final2}
\widehat{T}(\eta|n)&=\varrho  e^{-\beta W(x_1; \eta^{(\hat{x}_1)})}\sum_{k=0}^n {n \choose k}
\int (dy)^k K(x_1;\{y\}_1^k)\widehat{T}(\eta^{(\hat{x}_1)}\cup\{y\}_1^k)|n-k)-\notag \\
&-\sum_{k=1}^n {n \choose k}\sum_{l=1}^k {k \choose l}\int (dy)^l K(x_1;\{y\}_1^l)\widehat{T}(\eta|n-k)\widehat{T}(\{y\}_1^l|k-l).
\end{align}
For $k=0$ in the first line of formula \eqref{final2}, the term $\widehat{T}(\eta^{(\hat{x}_1)}|n)$ corresponds to the analytical
contribution of all graphs $\mathcal{D}(\eta; \gamma)$, in which the vertex $x_1$ has no edges.
For $k\geq 1$ the set of graphs for $\widehat{T}(\eta^{(\hat{x}_1)}\cup\{y\}_1^k)|n-k)$ is arranged in such a way
that $k$ vertices from $n$ vertices of the configuration $\gamma$ are transformed into vertices, which are connected
to the graphs of $\widehat{T}(\eta^{(\hat{x}_1)}|n-k)$ either by means of lines or as free vertices.
If a vertex $y_i$ is already connected by at least one line, then joining a line with vertex $x_1$
leaves this graph in the set $\mathcal{D}(\eta; \gamma)$. If the vertex did not have edges,
then operation  $\int (dy)^l K(x_1;\{y\}_1^l)$ does not ensure that it belongs to the set
$\mathcal{D}(\eta; \gamma)$, and such a diagram should be reduced.This is ensured by the presence of
a similar contribution in the second line of the formula \eqref{final2}.
This variant occurs in $l$ of $k$ cases for  each of the graphs in  $\widehat{T}(\eta|n-k)$.
\hfill $\blacksquare$

\subsection{ Conclusion on the convergence of the expansion.}

Existence of a unique solution of the nonlinear equation \eqref{limeq} as an operator equation in the Banach space $ E_\xi $
was proved in \cite{Po75}. In a recent work of Jansen \cite{Ja21} a generalized version of Groeneveld’s convergence criterion for the
virial expansion and generating functionals for weighted 2-connected graphs was proven. In fact, the present paper turned out
to be very close to article \cite{Ja21}, and here one can apply the criteria for the convergence of expansions.
Therefore, we do not yet discuss this issue in this preliminary publication. Although perhaps a more thorough study
of the recurrence relation \eqref{final2} may introduce some adjustments.

\vspace{0.2cm}

{\bf AVAILABILITY OF DATA}

\vspace{0.2cm}

Data available on request from the
author

\vspace{0.2cm}

\small{


\section{Introduction} \label{sec:intro}

The study of rare isotopes has unveiled an array of surprising exotic phenomena, arising from the subtle interplay of the different forces acting on loosely bound nucleons. In some instances, these nucleons are allowed to reside at unusually large distances from the center of mass of the system, creating a halo. In these so-called halo nuclei \cite{hove2018}, there is a strong decoupling of the core degrees of freedom and the degrees of freedom associated with the halo nucleons.
Halo nuclei can be found for many different isotopic chains, both at the neutron-rich limit as well as the proton-rich extreme. For example, in the Carbon isotopic chain, one finds the exotic $^{9}$C \cite{fukui2012} in the proton-rich side and the two-neutron halo $^{22}$C \cite{tanaka2010,pinilla2016,nagahisa2018} on the neutron dripline; in the Boron isotopic chain, while $^8$B exhibits a well studied proton halo (e.g.,  \cite{rangel2016,mazzocco2019,wang2021}), $^{19}$B has a neutron halo \cite{cook2020,casal2020}.

Halo nuclei are primarily studied through reactions and, over recent decades, a variety of experimental programs at many rare isotope facilities worldwide have unveiled their properties (e.g., \cite{davids2001,nakamura2009,chen2016,heine2017,cook2020}). Breakup reactions offer a unique probe into these systems: since the halo nucleus breaks up so easily, the cross sections for the process are large. From the analysis of breakup experiments we learn about the halo properties of the ground state (particularly the asymptotic normalization coefficient (ANC)) as well as the low-lying continuum including resonances (e.g., \cite{capel2006}). From breakup measurements, we can also infer capture reactions for astrophysics \cite{nunes2020}.

Reaction theory is needed in order to connect the breakup measurement to the quantities of interest, be it specific properties of the projectile or a capture reaction of astrophysical interest. One challenge faced when modeling breakup is that, for most cases, it  becomes computationally intensive to evaluate observables given theory parameters.
It has been shown that the simple perturbative approaches are unreliable \cite{capel2012} and instead there are large coupled equations that need to be solved, requiring a significant amount of run time (e.g., \cite{deltuva2009,upadhyay2012,hlophe2019}).
As a consequence, the current predictions include no uncertainty quantification.

Progress on uncertainty quantification in reaction theory has offered important insights, particularly concerning the effective potential, a strong contributor to the theoretical error \cite{lovell2017,lovell2018}. Many of the recent uncertainty quantification studies involve a Bayesian analysis of elastic scattering within the optical model, and at times, the propagation of uncertainties on the optical potential to specific reaction channels that can be modeled in a simple perturbative description \cite{king2019,catacora2019,lovell2020,catacora2021,whitehead2021}. It is important to go beyond the first-order reaction models covered in these pioneering works, and enable a Bayesian analysis of the state-of-the-art reaction theory models \cite{deltuva2009,upadhyay2012,hlophe2019}. Obviously, from the computational cost, it is not feasible to do statistical calculations (such as Markov chain Monte Carlo) directly on the model.
Emulators that can efficiently leverage model outputs are called for to solve this problem.

Recently, emulators for two-body elastic scattering have been developed based on the eigenvector continuation method \cite{furnstahl2020,drischler2021,melendez2021}. To describe the breakup dynamics one needs at least a three-body theory, as is the case for the Continuum Discretized Coupled Channel method (CDCC) \cite{cdcc}, the working horse in most analyses of breakup reactions.
Despite some preliminary work \cite{zhang2021}, it is not straightforward to extend such an eigenvector continuation method to the general three-body problem. Instead, a standard emulation by Gaussian Processes (GPs) (see \cite{santner2018design, gramacy2020surrogates}) should be considered. GP emulators often outperform other statistical learning tools for these types of tasks \citep{myren2021comparison}.

In this work, we implement an emulator based on GPs for breakup cross sections and couple this to a Bayesian framework. We focus on the breakup reaction of $^8$B$ + ^{208}$Pb$ \rightarrow ^7$Be$ + p + ^{208}$Pb at $80$ MeVA, which has been studied extensively  \cite{davids2001,capel2007}. The emulator is trained with CDCC calculations and the sharpness and accuracy of the uncertainty quantification of the emulator is studied. We then use the emulator for calibration, and extract posterior distributions for the parameters of the effective interaction between the halo nucleon and the core. Lastly, we obtain credible intervals for relevant reaction observables.

\section{Model description} \label{sec:model}

We now detail our reaction model and Bayesian calibration.

\subsection{Modeling breakup reactions}

Breakup reactions of halo nuclei of the type $a+t \rightarrow c+v+t$ are typically formulated in terms of a three-body problem $c+v+t$ with the projectile assumed to have two-body structure $a=c+v$. The three-body Hamiltonian
\beq
H_{3B}=T_r + T_R + V_{cv}(\ve r)+U_{ct}(\ve{R}_{ct})+U_{vt}(\ve{R}_{vt})
\eeqn{eq:h3b}
is written in terms of the respective pairwise interactions, the interaction describing the internal states of the projectile $V_{cv}$ and the optical potentials describing the scattering of the core (valence nucleons) and the target $U_{ct}$ ($U_{vt}$). Here, $T_r$ and $T_R$ are the kinetic energy operators associated with the Jacobi coordinates: $\ve r$ is the relative coordinate between the core and the valence nucleon and  $\ve R$ is the relative coordinate connecting the center of mass of the projectile and the target.

The CDCC method \cite{cdcc} expands the full three-body wavefunction $\Psi(\ve r, \ve R)$ in terms of the eigenstates of the projectile $\varphi_i(\ve r)$ defined by
\beq
\left[ T_r +  V_{cv}(\ve r) \right] \varphi_{i}(\varepsilon,\ve{r})=\varepsilon\varphi_{i}(\varepsilon,\ve{r}).
\eeqn{eq:hproj}
For halo nuclei, typically there is one bound state $i=0$ and the remaining states lie in the continuum. Because the continuum solutions of Eq.~\ref{eq:hproj} are not square integrable, the CDCC method discretizes the continuum into energy bins and averages the scattering states over the energy or momentum as described in \cite{book}. When this expansion is introduced in the three-body Sch\"odinger equation $H_{3B} \Psi = E \Psi$, one obtains the so-called CDCC equations, which consist of tightly coupled second-order differential equations that must be solved with the appropriate boundary conditions.
For more details, please consult \cite{book,nunes99}.

There are three pairwise interactions in this model, introducing close to $30$ parameters. Clearly breakup data alone would not be sufficient to inform all these parameters. Since this is the first Bayesian study for this type of reaction, here we choose to focus only on the uncertainty from the core-valence effective interaction, and thus fix the optical potentials to those used in a previous study \cite{mortimer2002,capel2007}. The effective interaction $V_{cv}$ has been shown to have a very important effect in the breakup cross section distributions, renormalizing the cross section according to the ANC of the ground state \cite{capel2006} and thus we expect the breakup data will be highly informative for the $V_{cv}$ parameters.
This $V_{cv}$ interaction is typically parameterized by a Woods-Saxon form \cite{capel2007}.
We take the radius  and diffuseness for the spin-orbit term to be the same as that for the central term. We then vary the Coulomb radius, the Woods-Saxon radius, the diffuseness, and the depth of the spin-orbit term ($R_C, R_{ws}, a_{ws}, V_{so}$)  by adjusting the depth of the central interaction $V_{ws}$ to reproduce the binding energy of the $^8$B system with $\varepsilon = 0.137$ MeV.  For simplicity, the same interaction is used to produce both the $^{8}$B ground state and all the continuum states.

For a combination of parameters ($R_C, R_{ws}, a_{ws}, V_{so}$) in the physical intervals defined in Table~\ref{tab:parameters}, the corresponding CDCC cross sections are obtained using {\sc frescox} \cite{fresco}. Given that the calculations of \cite{capel2007} reproduce well the data, the parameters used therein, and in particular the core-valence interaction in the projectile ($R_C=R_{ws}=2.391$ fm, $a_{ws}=0.52$ fm, and $V_{so}=4.898$ MeV), are used to produce our mock observable data, namely the breakup angular distributions and the breakup energy distributions. In the next subsection, we discuss the implementation of the emulators for these breakup calculations.
\begin{table}[tbh!]
    \centering
        \caption{Model parameters and their ranges. \label{tab:parameters}}
    \begin{tabular}{c|| c l}

    Parameter & Label & Range $[\underline{{\rho}_i}, \overline{{\rho}_i}]$ \\[0.5ex]
    \hline\hline
     $R_{C}$ & ${\rho}_1$ & [2, 3] (fm) \\

     $R_{ws}$ & ${\rho}_2$ & [2, 3] (fm) \\
     $a_{ws}$ & ${\rho}_3$ & [0.4, 0.9] (fm) \\
     $V_{so}$ & ${\rho}_4$ & [2, 8] (MeV)
 \\%[1ex]

    \end{tabular}
\end{table}

\subsection{Bayesian calibration of breakup reactions}
\label{sec:bayesian_calibration}

After reviewing Bayesian calibration principles and notation, we present our methodology for emulation of CDCC calculations.

\subsubsection{Bayesian calibration}

We represent a CDCC simulation with a mathematical function ${\pmb{\sigma}}(\cdot)$ that takes values of parameters ${\pmb{\rho}}$ and returns output ${\pmb{\sigma}}({\pmb{\rho}}) \in {\mathbb{R}}^d$.
We note that ${\pmb{\sigma}}({\pmb{\rho}})$ represents multidimensional output observables, which include both cross section angular distributions and energy distributions. ${\pmb{\sigma}}({\pmb{\rho}})$ should not be confused with total or reaction cross sections. We collect the four parameters in Table~\ref{tab:parameters} in the vector ${\pmb{\rho}} = ({\rho}_1, \ldots, {\rho}_4)$. To emulate and then calibrate the breakup reactions, we use breakup angular $\angular$ and energy $\frac{{\rm d} \sigma}{\rm d E}$ distributions. We denote the dimension of the $\angular$ and $\frac{{\rm d} \sigma}{\rm d E}$ outputs by $d_\Omega$ and $d_E$, respectively, with $d = d_\Omega + d_E$.

Next we describe the model space needed for convergence of the CDCC calculations.
The model space is optimized per observable and we find that the bin discretization needed for producing converged $\angular$ is different than that needed for converged $\frac{{\rm d} \sigma}{\rm d E}$.
For the angular distribution, the continuum is discretized up to a maximum relative energy of $\varepsilon_{\rm max} = 10$ MeV, and all spin/parity states are included up to a relative angular momentum between core and valence of $\ell=3$.
For each partial wave, the continuum is discretized into energy bins evenly spaced from 0 MeV to 3 MeV with a step size of 0.25 MeV, and from 3 MeV to 10 MeV with a step size of 1 MeV.
The energy distributions also contain all partial waves up to $\ell=3$, but we take the maximum relative energy of $\varepsilon_{\rm max} = 3$ MeV with a uniform grid with a step size of $0.20$ MeV. This finer discretization is needed to capture the details of the peak of the distribution.

A calibration process uses observations from the real system (here, the mock data), denoted by vector ${\mathbf{y}} = (y_1, \ldots, y_d) = ( \{ \angular_{\rm mock} \}, \{ \frac{{\rm d} \sigma}{\rm d E}_{\rm mock} \} )$ with $d=d_{\Omega} + d_E$, to constrain the uncertainty on the input parameters ${\pmb{\rho}}$.
This is done using a statistical model of the form
\begin{equation}
    {\mathbf{y}} = {\pmb{\sigma}}({\pmb{\rho}}) + {\pmb{\epsilon}},
\end{equation}
where ${\pmb{\epsilon}}\sim {\rm MVN}(\mathbf{0}, {\pmb{\Sigma}})$ denotes the residual error following a multivariate normal distribution (MVN) with mean $\mathbf{0}$ and covariance matrix ${\pmb{\Sigma}}$.

In the Bayesian framework, the model parameters are viewed as random variables, and the posterior probability density ${p}({\pmb{\rho}}|{\mathbf{y}})$ indicates the posterior probability of the parameters ${\pmb{\rho}}$ given the observations ${\mathbf{y}}$. Based on Bayes' rule, the posterior density has the form
\begin{equation}
    {p}({\pmb{\rho}}|{\mathbf{y}}) = \frac{{p}({\mathbf{y}}|{\pmb{\rho}}) {p}({\pmb{\rho}})}{{p}({\mathbf{y}})} \propto {p}({\mathbf{y}}|{\pmb{\rho}}) {p}({\pmb{\rho}}), \label{eq:posterior}
\end{equation}
where ${p}({\pmb{\rho}})$ is the prior probability for the parameter ${\pmb{\rho}}$ and ${p}({\mathbf{y}}|{\pmb{\rho}})$ is the likelihood function indicating how the model output of a set of parameters ${\pmb{\sigma}}({\pmb{\rho}})$ agrees with the observations ${\mathbf{y}}$. Assuming ${\pmb{\epsilon}} \sim {\rm MVN}(\mathbf{0}, {\pmb{\Sigma}})$, the likelihood satisfies
\begin{equation}
   {p}({\mathbf{y}}|{\pmb{\rho}}) \propto |{\pmb{\Sigma}}|^{-1/2} \exp\bigg(-\frac{1}{2} ({\mathbf{y}} - {\pmb{\sigma}}({\pmb{\rho}}))^\top {\pmb{\Sigma}}^{-1} ({\mathbf{y}} - {\pmb{\sigma}}({\pmb{\rho}} )) \bigg). \label{eq:truelike}
\end{equation}

Markov chain Monte Carlo (MCMC) techniques \cite{Gelman2004} (e.g., Metropolis-Hastings algorithm) are generally used to sample from the posterior distribution in Eq.~\ref{eq:posterior}. During an MCMC procedure, an evaluation of the likelihood ${p}({\mathbf{y}}|{\pmb{\rho}})$ in Eq.~\ref{eq:truelike} is needed for each candidate parameter ${\pmb{\rho}}$. The likelihood requires the simulation output ${\pmb{\sigma}}({\pmb{\rho}})$ to be evaluated at a given parameter setting ${\pmb{\rho}}$. MCMC procedures usually require thousands or millions of evaluations of the simulation depending upon the model's complexity and the number of input parameters considered. Therefore, a direct evaluation of a reaction model with MCMC techniques becomes computationally challenging when substantial computational time is required to produce a single model evaluation, as is the case for CDCC calculations (each run described here takes approximately 25 hours on a single core).

Emulators are computationally efficient approximations of the input–output relationships for expensive simulations, and they have been used to address the computational challenges associated with executing the simulation run. To train an emulator, a sample of parameters is generated, and the corresponding simulation outputs are collected. Once an emulator is built, it can be used to efficiently predict the simulation outputs for parameter values that are not in the training set. In this study, we use GP emulators since, in addition to a predictive mean, they can produce a predictive variance, which allows us to quantify the uncertainty on emulation predictions. For an untried parameter ${\pmb{\rho}}$, the GP provides a probabilistic representation of the simulation output with mean ${\pmb{\mu}}^{\rm GP}({\pmb{\rho}})$ and covariance ${\pmb{\Sigma}}^{\rm GP}({\pmb{\rho}})$ based on training data from a simulation experiment.

Once the emulator is built,
for any parameter ${\pmb{\rho}}$
one can
apply Eqs.~\ref{eq:posterior}~and~\ref{eq:truelike} to
approximate the posterior ${p}({\pmb{\rho}}|{\mathbf{y}})$ via
\begin{align}
    p({\pmb{\rho}}|{\mathbf{y}}) \propto & |\mathbf{V}({\pmb{\rho}})|^{-1/2}\exp\left(-\frac{1}{2} \mathbf{m}({\pmb{\rho}})^\top \mathbf{V}({\pmb{\rho}}) ^{-1} \mathbf{m}({\pmb{\rho}})\right) p({\pmb{\rho}}),
\label{eq:emu_posterior}
\end{align}
where, for our GP emulator, the quantities $\mathbf{V}({\pmb{\rho}})$ and $\mathbf{m}({\pmb{\rho}})$ are defined as
\begin{align}
\mathbf{V}({\pmb{\rho}}) = {\pmb{\Sigma}} + {\pmb{\Sigma}}^{\rm GP}({\pmb{\rho}}) \text{ and } \mathbf{m}({\pmb{\rho}}) = {\mathbf{y}} - {\pmb{\mu}}^{\rm GP}({\pmb{\rho}}).
\label{eq:v-m}
\end{align}
In this study, we employ an MCMC method to draw approximate posterior samples from Eq.~\ref{eq:emu_posterior}. Such sampling employs the cheap GP emulator in place of the expensive CDCC simulation, which would require many years of computation. The details of the sampler are explained in Sec.~\ref{sec:results}.

\subsubsection{Emulating breakup reactions}
\label{sec:emulators}

We generated three separate GP emulators for the breakup cross sections, which were then used for calibration. The first emulator considers only breakup angular distributions $\angular$. The second emulator takes into account only the breakup energy distributions $\frac{{\rm d} \sigma}{\rm d E}$. The third emulator includes both the angular and energy distributions. For each of these emulators, we use the same set of uncertain parameters as inputs to the CDCC calculations.
Principal component analysis (PCA) \cite{Ramsay97functionaldata} is used to project the high-dimensional outputs into a low-dimensional space where the projection is a collection of latent outputs. Then, each latent output is modeled using an independent GP model (see details in Apps.~\ref{appendix:emulator}--\ref{appendix:hyper}).

Since computing a single CDCC simulation via {\sc frescox} takes 25 hours on average, we use a GP-based emulator as a cheaper proxy to the CDCC simulation, as discussed before. We collect data via Latin hypercube sampling (LHS) (see the texts \cite{santner2018design, gramacy2020surrogates}) based on the ranges in Table~\ref{tab:parameters}. For evaluation purposes, we randomly split $75 \%$ of these data as training data and the remaining as test data; we then fit an emulator using the training data and evaluate the predictive quality of the trained emulator using the test data. Both the training and test data are illustrated in Fig.~\ref{fig-sample}; diagnostic plots are provided in Sec.~\ref{sec:predictions}.

Once the quality of the emulator has been established, the final emulators are trained with a total of $n=500$ samples, collected via LHS. There were no additional difficulties found for the angular distribution emulator; however, for the energy distribution, we obtained outputs that differed tremendously from the mock data in some regions of the parameter space (see App.~\ref{appendix:pca} for more details). This is due to the existence of resonances in the core-valence system that greatly augment the sensitivity to parameters and reduce the precision in the resulting emulator. To tackle this problem, we developed a filter based on the deviation from the mock data of the predicted energy distributions: we discard from the training set any case in which $\frac{{\rm d} \sigma}{\rm d E}$ is larger than the threshold of $1000$ mb/MeV (corresponding to roughly five times the peak of the ``real'' energy distribution).
Discarding parameter settings based on simulation outputs \cite{Surer2021} has been shown to provide emulators with better accuracy in the reduced space. As a consequence of filtering, the number of samples ultimately included in the training is $n=393$. For the third emulator, where the cross sections for angular and energy distributions are jointly considered, we also use $n = 393$ samples.
Simulation outputs before and after filtering are provided in App.~\ref{appendix:pca}.

\begin{figure}[t]
    \begin{center}
        \includegraphics[width=0.48\textwidth]{sample.png}
    \end{center}
    \caption{Space-filling Latin hypercube sample (LHS) used to specify the set of breakup simulation runs to train the GP emulator (black) and to test the emulator (red). The 2-d projections of this design are shown for each pair of parameters.}
    \label{fig-sample}
\end{figure}

For the rest of this section, let $\{{\pmb{\rho}}_1^{\rm tr}, \ldots, {\pmb{\rho}}_n^{\rm tr}\}$ denote the ${n}$ parameter samples used to train each of the emulators. For each ${\pmb{\rho}}^{\rm tr}_i$, $i=1, \ldots, n$, we obtain the simulation outputs (CDCC cross sections) via {\sc frescox} and use them as training data. The simulation outputs are represented in a $d \times n$ matrix $\Xi = [{\pmb{\sigma}}({\pmb{\rho}}_1^{\rm tr}), \ldots, {\pmb{\sigma}}({\pmb{\rho}}_{n}^{\rm tr})]$. We standardize each row of $\Xi$ so that it has zero mean and unit variance. After standardization, we employ PCA to reduce the dimension of the data from $d$ to $q \leq d$. We let $q_\Omega$, $q_E$, and $q_{\rm tot}$ denote the number of principal components used to fit an emulator for $\angular$, $\frac{{\rm d} \sigma}{\rm d E}$, and both $\angular$ and $\frac{{\rm d} \sigma}{\rm d E}$ as observables, respectively. PCA reduces the ($d_\Omega = 401$)- and ($d_E = 20$)-dimensional cross sections to $q_{\Omega} = 39$ and $q_{E} = 19$, respectively; this reduction is able to capture $99 \%$ of the variance. Note that for energy distribution there is no significant dimensionality reduction since all dimensions contribute to the variance of data. Figs.~\ref{fig-PCAexplained_angular}--\ref{fig-PCAexplained_energy} illustrate how much variance each principal component explains.

To combine $\angular$ and $\frac{{\rm d} \sigma}{\rm d E}$ cross sections, we first reduce the dimension of $\angular$ (i.e., $d_\Omega = 401$) to 25 to treat $\angular$ and $\frac{{\rm d} \sigma}{\rm d E}$ in a fair manner when computing the likelihood. To do that, we consider the cross sections up to the scattering angle of $3\degree$, which reduces the dimension from 401 to 50. This is a reasonable assumption because these breakup simulations are forward focused, and typically this is the angular range measured by detectors. We also find that including the larger angles does not modify the results.
We then take $d_\Omega=25$ evenly spaced angles within these $3\degree$. The final set contains $d_{\rm tot} = 45$ (i.e., 25 for $\angular$ and 20 for $\frac{{\rm d} \sigma}{\rm d E}$) dimensional observables. During emulation, PCA further reduces this set $d_{\rm tot} = 45$ to $q_{\rm tot} = 32$.
\begin{figure}[ht]
    \begin{center}
        \includegraphics[width=0.44\textwidth]{PCA_explained_angular.png}
    \end{center}
    \caption{The individual and cumulative contributions of principal components to the complete data set (angular distribution).}
    \label{fig-PCAexplained_angular}
\end{figure}
\begin{figure}[ht]
    \begin{center}
        \includegraphics[width=0.44\textwidth]{PCA_explained_energy.png}
    \end{center}
    \caption{The individual and cumulative contributions of principal components to the complete data set (energy distribution).}
    \label{fig-PCAexplained_energy}
\end{figure}

The $d \times n$ CDCC cross section matrix $\Xi$ is represented in the $q$-dimensional space via the $q \times n$ matrix ${\mathbf{W}} = {\mathbf{B}}^\top \Xi$ where the $d \times q$ matrix ${\mathbf{B}} = [{\mathbf{b}}_1, \ldots, {\mathbf{b}}_q]$ stores the orthogonal basis vectors $\{{\mathbf{b}}_1, \ldots, {\mathbf{b}}_q\}$. We obtain ${\mathbf{W}}$ with PCA as described in App.~\ref{appendix:emulator}, and fit a GP for each of $q$ latent outputs (i.e., each row of ${\mathbf{W}}$) such that
\begin{equation}
    w_i(\cdot) \sim \text{GP}(\gamma_i, C_i(\cdot, \cdot)),
    \label{eq:gp_prior}
\end{equation}
where $\gamma_i$ is the mean and $C_i(\cdot, \cdot)$ is the covariance function. The details of GP emulator are given in Apps.~\ref{appendix:emulator}--\ref{appendix:hyper}. The full implementation of this GP emulator for multidimensional simulation outputs is provided in the Python package \texttt{surmise} \cite{surmise2021}.

\subsubsection{Predicting simulation outputs with emulators}
\label{sec:predictions}

As mentioned above, we fit three different emulators using training data. The next step is to test each of the emulators at any test parameter ${\pmb{\rho}}$ to check whether they are able to predict the CDCC calculations. Let the $n$-dimensional vector ${\mathbf{w}}_i$ be ${\mathbf{w}}_i = (w_i({\pmb{\rho}}_1^{\rm tr}), \ldots, w_i({\pmb{\rho}}_{n}^{\rm tr}))$ for $i = 1, \ldots, q$, and let the $n \times n$ matrix ${\mathbf{C}}_{i}$ be the covariance matrix resulting from applying Eq.~\ref{eq:covariance} to each pair of the parameter settings ${\pmb{\rho}}_1^{\rm tr}, \ldots, {\pmb{\rho}}_{n}^{\rm tr}$. Following the results on normal distribution presented in \cite{Rasmussen2005}, one can derive that
\begin{equation}
    w_i({\pmb{\rho}})|{\mathbf{w}}_i \sim \text{N}(\nu_i({\pmb{\rho}}), \varsigma^2_i({\pmb{\rho}},{\pmb{\rho}})),
\end{equation}
with mean $\nu_i({\pmb{\rho}})$ and variance $\varsigma^2_i({\pmb{\rho}},{\pmb{\rho}})$ given by
\begin{equation}
    \begin{aligned}
        & \nu_i({\pmb{\rho}}) = \hat{\gamma}_i + {\mathbf{c}}_i({\pmb{\rho}})^\top {\mathbf{C}}_{i}^{-1} ({\mathbf{w}}_i - \hat{\gamma}_{i} {\mathbf{1}}_n) \\
        & \varsigma^2_i({\pmb{\rho}},{\pmb{\rho}}) = C_i({\pmb{\rho}}, {\pmb{\rho}}) - {\mathbf{c}}_i({\pmb{\rho}})^\top {\mathbf{C}}_{i}^{-1} {\mathbf{c}}_i({\pmb{\rho}}). \label{eq:meanvar_latent}
    \end{aligned}
\end{equation}
Here, the $n$-dimensional vector ${\mathbf{c}}_i({\pmb{\rho}})$ is obtained by evaluation of the covariance in Eq.~\ref{eq:covariance} between ${\pmb{\rho}}$ and the training points $\{{\pmb{\rho}}_1^{\rm tr}, \ldots, {\pmb{\rho}}_n^{\rm tr}\}$ and $\hat{\gamma}_i$ is the estimated parameter for the mean function of a GP (see App.~\ref{appendix:hyper}).

For prediction at any test point ${\pmb{\rho}}$, first, we obtain the mean $\nu_i({\pmb{\rho}})$ and variance $\varsigma^2_i({\pmb{\rho}},{\pmb{\rho}})$ in Eq.~\ref{eq:meanvar_latent} for each of the corresponding latent outputs for $i = 1, \ldots, q$. Then, these are transformed back to the original high-dimensional space through the inverse PCA transformation as follows.
Define the $q$-dimensional vector ${\pmb{\mu}}({\pmb{\rho}}) = (\nu_1({\pmb{\rho}}), \ldots, \nu_q({\pmb{\rho}}))$ and $q \times q$ diagonal matrix $ \mathcal{S}({\pmb{\rho}})$ with diagonal elements $\varsigma^2_i({\pmb{\rho}},{\pmb{\rho}})$ for $i = 1, \ldots, q$. Then, due to the inverse PCA transformation, we obtain
\begin{equation}
    {\pmb{\sigma}}({\pmb{\rho}}) \sim \text{MVN}({\pmb{\mu}}^{\rm GP}({\pmb{\rho}}), {\pmb{\Sigma}}^{\rm GP}({\pmb{\rho}}))
    \label{emu_final}
\end{equation}
where ${\pmb{\mu}}^{\rm GP}({\pmb{\rho}}) \coloneqq {\mathbf{B}} {\pmb{\mu}}({\pmb{\rho}})$ is the emulator predictive mean and ${\pmb{\Sigma}}^{\rm GP}({\pmb{\rho}}) \coloneqq {\mathbf{B}}\mathcal{S}({\pmb{\rho}}){\mathbf{B}}^\top$ is the covariance matrix. By plugging ${\pmb{\mu}}^{\rm GP}({\pmb{\rho}})$ and ${\pmb{\Sigma}}^{\rm GP}({\pmb{\rho}})$ in Eq.~\ref{eq:v-m} we obtain $\mathbf{m}({\pmb{\rho}})$ and $\mathbf{V}({\pmb{\rho}})$. Then from Eq.~\ref{eq:emu_posterior}, we compute the posterior $p({\pmb{\rho}}|{\mathbf{y}})$.
In this way, by employing the MCMC method on the emulators, we extract posterior distributions for the parameters of the effective interaction between the halo nucleon and the core without running the computationally expensive CDCC calculations when doing the MCMC sampling.

In order to test the prediction quality of the trained emulator, we use $m$ parameter sets, which correspond to the $25 \%$ of the parameter sets, namely ${\pmb{\rho}}_1^{\rm test}, \ldots, {\pmb{\rho}}_m^{\rm test}$, and their corresponding CDCC calculations as test data. We evaluate the performance by computing the test $r^2$ value, standardized error, and the relative error since they are standard ways to evaluate the emulator quality (see App.~\ref{appendix:diag} for the definition of the metrics). The diagnostics plots are illustrated in Fig.~\ref{fig:diagnostics_combined} for fitting the emulator based on the combined angular and energy distributions. For the sake of brevity, we provide the diagnostics plots for the remaining emulators in App.~\ref{appendix:diag}. For all the emulators, the test $r^2$ value is very close to one, and both relative and standardized errors are centered around zero.
\begin{figure}
     \begin{subfigure}{0.48\columnwidth}
         \includegraphics[width=1\textwidth]{emu_ang_energy_diag2.png}
         \caption{Standardized errors.}
     \end{subfigure}
     \begin{subfigure}{0.48\columnwidth}
         \includegraphics[width=1\textwidth]{emu_ang_energy_diag3.png}
         \caption{Relative errors.}
     \end{subfigure}
        \caption{Diagnostic plots to check the quality of the emulator using test parameters ${\pmb{\rho}}_1^{\rm test}, \ldots, {\pmb{\rho}}_{l}^{\rm test}$ and their corresponding simulation outputs ${\pmb{\sigma}}({\pmb{\rho}}_1^{\rm test}), \ldots, {\pmb{\sigma}}({\pmb{\rho}}_{l}^{\rm test})$ for the emulator based on combined angular and energy distributions with $r^2=0.99$: (a) standardized errors and (b) relative errors. Blue line shows the density of a standard normal random variable.}
        \label{fig:diagnostics_combined}
\end{figure}

\section{Results and Discussion}
\label{sec:results}

\begin{figure}[t]
    \begin{center}
        \includegraphics[width=0.4\textwidth]{Posterior_combined.png}
    \end{center}
    \caption{Univariate marginal estimates of the posterior distribution for the 4-dimensional parameter vector of the breakup reactions for the emulator based on both angular and energy distributions:  prior distributions (light grey) and posterior distributions (dark grey) for the Coulomb radius parameter (a), the Woods-Saxon radius (b), the Woods-Saxon diffuseness (c) and the spin-orbit depth (d). The red dashed line corresponds to the true parameter values.}
    \label{fig-post-combined}
\end{figure}
With trained GP emulators of the computationally expensive CDCC reaction model in hand, we now discuss the posterior distributions obtained for the parameters of the effective interaction between the halo nucleon and the core, as well as the resulting uncertainties on the breakup cross sections.

As mentioned in Sec.~\ref{sec:model}, mock data is generated with the parameters of \cite{capel2007}. Based on the typical statistical errors in breakup experiments of halo nuclei, we include a standard deviation of $10 \%$ on the mock data. Bayesian analysis requires the prior distributions for the core-valence interaction parameters. We use a beta distribution with the shape parameters $\alpha = \beta = 3.5$ as prior for each core-valence parameter ${\pmb{\rho}}=({\rho}_1,\ldots,{\rho}_4)$, meaning
\[
    p({\pmb{\rho}}) = \prod_{i=1}^4 \frac{\Gamma(\alpha + \beta) \bigg(\frac{{\rho}_i - \underline{{\rho}_i}}{\overline{{\rho}_i} - \underline{{\rho}_i}}\bigg)^{\alpha -1}\bigg(1-\frac{{\rho}_i - \underline{{\rho}_i}}{\overline{{\rho}_i} - \underline{{\rho}_i}}\bigg)^{\beta-1}}{\Gamma(\alpha)\Gamma(\beta)},
\]
where $\underline{{\rho}_i}$ and $\overline{{\rho}_i}$ denote, respectively, the lower and upper bounds of the range of parameter ${\rho}_i$ given in Table~\ref{tab:parameters}. The beta distribution is a widely used prior since it is a bounded continuous distribution with a decreasing density at the edges of the parameter range. We implement the parallel-tempering ensemble MCMC algorithm \cite[Chapter 10]{liu2008monte} through the \texttt{surmise} package in order to generate a set of samples from the posterior distribution.

The posteriors obtained by combining both angular and energy distributions are illustrated in Fig.~\ref{fig-post-combined}. For the sake of brevity, the results obtained for angular and energy distributions as separate observables are provided in App.~\ref{appendix:cal}. As compared to the prior distributions represented by the light grey bins, the posterior distributions for $R_{ws}$ and $a_{ws}$ are tightly concentrated around the true parameter values (represented by the red line); this indicates that $R_{ws}$ and $a_{ws}$ are well constrained by the breakup cross section distributions. On the other hand, the posterior distributions of $R_C$ and $V_{so}$ are not very well constrained by this breakup data, implying that those parameters are not very influential in predicting both angular and energy distributions. We did not find significant differences in the results obtained using the other two emulators, which suggests that the information on $V_{cv}$ contained in the breakup energy distribution is similar to that contained in the angular distribution.

Moreover, for the emulator based on the combined angular and energy distributions, we obtain the pair plots using posterior samples of four parameters as in Fig.~\ref{pairplot_combined}. The pair plots indicate a negative correlation between the two most influential parameters $R_{ws}$ and $a_{ws}$.
\begin{figure}[t]
    \begin{center}
        \includegraphics[width=0.48\textwidth]{Corr_Combined.png}
    \end{center}
    \caption{Pair plots for the posterior samples of four parameters obtained using MCMC with angular and energy distribution data.}
    \label{pairplot_combined}
\end{figure}

\begin{figure}[t]
     \begin{subfigure}{0.9\columnwidth}
         \includegraphics[width=1\textwidth]{CombinedIntervalA.png}
         \caption{Angular distribution.}
     \end{subfigure}
     \begin{subfigure}{0.9\columnwidth}
         \includegraphics[width=1\textwidth]{CombinedIntervalE.png}
         \caption{Energy distribution.}
     \end{subfigure}
     \caption{Cross section for $^8$B$ + ^{208}$Pb$ \rightarrow ^7$Be$ + p + ^{208}$Pb at $80$ MeVA (a) angular distribution and (b) relative $^7$Be$+p$ energy distribution: prediction mean (dashed black line) and the 95\% credible interval (shaded gray area) obtained from the Bayesian analysis, compared with mock data. Results correspond to the emulator which uses combined angular and energy distribution data.}
    \label{pred_interval_combined}
\end{figure}
One of the main advantages of using Bayesian techniques is for a rigorous propagation of the quantified uncertainties in the parameters to the predictions. That is, once parameters are sampled from the posterior distribution, they can be used to produce predictions with quantified uncertainties on the breakup cross section distributions. Figure~\ref{pred_interval_combined} displays the $95\%$ credible intervals (grey band), along with the prediction means (black dashed line) and the mock data (red circles) for (a) the breakup angular and (b) the energy distribution. These intervals do not include the uncertainty on the emulation itself. We find the credible interval is roughly within two standard deviations of the predicted mean (corresponding to a $20$\% uncertainty).

Finally, from the  posterior distributions shown in Fig.~\ref{fig-post-combined}, we compute the corresponding ANC parameters. Figure~\ref{ANCdistribution} shows the ANC posterior distribution compared to the prior used in the training set. In the model used for the core-valence interaction, the radius and diffuseness are correlated in such a way that the combined breakup data is able to tightly constrain the ANC. This result is consistent with a study on another halo nucleus \cite{capel2006}. From Fig.~\ref{ANCdistribution}, we can extract the value for the $^8$B ANC squared $C^2=0.49 \pm 0.01$ (fm$^{-1}$) which agrees with the value extracted from a different breakup reaction \cite{b8anc}.
\begin{figure}[ht]
    \begin{center}
        \includegraphics[width=0.3\textwidth]{ANC.png}
    \end{center}
\caption{The estimate of the posterior distribution (dark grey bins) of ANC parameters as compared to its prior distribution obtained from training parameters (light grey bins).}
\label{ANCdistribution}
\end{figure}

\section{Conclusions and Vision} \label{sec:conclusions}

This work consists of the first Bayesian analysis for nuclear breakup reactions.
Here, we focus on the breakup reaction of $^8$B$ + ^{208}$Pb$ \rightarrow ^7$Be$ + p + ^{208}$Pb at $80$ MeVA. The reaction model used is the continuum discretized coupled channel method, inherently a three-body non-perturbative approach to breakup, widely used in the field. Even though there are almost 30  parameters in the model, in this study we consider only the uncertainty from the core-valence effective interaction. The free parameters then are the Coulomb radius $R_C$, the Woods-Saxon radius $R_{ws}$ and diffuseness $a_{ws}$ and the spin-orbit strength $V_{so}$. The depth of the central Woods-Saxon term is fixed by imposing the experimental proton separation energy in the core-valence system.

In order to perform the Bayesian calibration, we constructed three emulators, one based on the breakup angular distribution, another based on the breakup energy distribution, and a third one based on the combined data set (angular and energy distributions). We used full CDCC calculations to train the emulators and found that approximately 400 training sets were sufficient to obtain an uncertainty on the emulation smaller than the typical error bar on the experimental data ($\approx 10$\%). We used PCA to further reduce the dimensionality of the observable space. We found no significant differences between the predictions obtained with the three emulators. In all three cases, posterior distributions obtained for $R_{ws}$ and $a_{ws}$ were significantly constrained, as opposed to the posterior distributions for $R_C$ and $V_{so}$. Most important, from the Bayesian calibration we obtained $95 \%$ credible intervals for the breakup observables. The resulting uncertainties on the breakup cross sections (both for the angular distribution and the energy distribution) were around $20$\%.

From this study, we also obtained a posterior distribution for the ANC which allows for an extraction of the ANC squared consistent with previous studies. We should emphasize that, although the uncertainty on the extracted ANC is less than $10$\%, this value only includes the uncertainty arising from the core-valence interaction in the three-body model. A larger study, including the uncertainty from the optical potentials, need to be done before the full parametric uncertainty in the model is accounted for.

One can expect further challenges when augmenting the parameter set to include all parameters in the core-target and valence-target optical potential. Obviously, increasing the number of parameters, the parameter space to be explored becomes much larger and the percentage of the parameter space corresponding to a region of high posterior becomes much smaller. In such cases, the training data sets need to be increased considerably to better explore the parameter space \cite{Loeppky2009}.
Moreover, the computational cost to fit an accurate emulator and run the MCMC sampler grows when the number of uncertain parameters increases. Considering the amount of computational time for a single run of {\sc frescox} (e.g., 25 hours), the one-shot, emulate-then-calibrate principle can lead to an ineffective calibration procedure. In such a case, a sequential Bayesian inference approach depending on active learning procedures may be a better strategy to perform the uncertainty quantification \cite{Jarvenpa2019, Joseph2019}.

\newpage


\section{Introduction}
Astrophysical and cosmological evidences propose that non-baryonic cold  dark matter (DM) accounts for $84\%$ of the matter density of the Universe~\cite{Aghanim:2018eyx}. Nevertheless, the nature of DM is still an open question for modern physics, as this elusive kind of matter can not be made of any of the Standard Model (SM) particles. Brane-world theory has been put forward as a prospective framework for DM candidates~\cite{2003PhRvL..90x1301C}. The characteristics of the suggested massive brane fluctuations (branons) in this theory match the ones of weakly interacting massive particles (WIMPs), which are one of the most favored candidates of cold DM~\cite{2012PhRvD..86b3506S}.

Dwarf spheroidal galaxies (dSphs) are preferred targets for indirect DM searches, due to their close proximities and high mass-to-light ratios. In this work, we are focusing on the ultra-faint dSph Segue~1 located in the Northern Hemisphere and outside of the Galactic plane ($ \text{RA} = \unit[10.12]{h} $, $ \text{DEC} = 16.08^{\circ}$) only $ \unit[23\pm2]{kpc}$ away from us~\cite{2007ApJ...654..897B}. The Segue~1 observational campaign was carried out by the ground-based gamma-ray telescope MAGIC between 2011 and 2013 and is with 157.9 hours the deepest observation of any dSph by a Cherenkov telescope to date~\cite{2014JCAP...02..008A}.

\section{Expected branon dark matter gamma-ray flux}

The expected gamma-ray flux produced by branon DM annihilation is composed of the \textit{astrophysical} factor (\textit{J}-factor), which depends on both the distance $ l $ and the DM distribution at the source region $ \rho_{\text{DM}} $, and the branon DM annihilation photon yield. It can be expressed in a given region of the sky, $\Delta\Omega$, and observed at Earth as
\begin{equation}
    \label{eq:Branon_Flux}
    \frac{\text{d}\Phi}{\text{d}E}\left( \langle\sigma v\rangle \right) = \underbrace{\int_{\Delta\Omega} d\Omega' \int_{\text{l.o.s.}} dl \, \rho_{\text{DM}}^{2} (l,\Omega')}_{\text{Astrophysics}} \cdot \underbrace{\frac{1}{4\pi} \frac{\langle\sigma v\rangle}{2m^{2}_{\chi}} \sum_{i=1}^{n} \text{Br}_{i} \frac{\text{d}N_{i}}{\text{d}E}}_{\text{Particle physics}},
\end{equation}
\noindent
where $ \langle\sigma v\rangle $ is the thermally-averaged annihilation cross section, $ m_{\chi} $ is the mass of the branon DM particle and $ \text{l.o.s.} $ stands for line-of-sight. The left panel of Fig.~\ref{fig:Branching_ratios_and_spectra} shows the branon branching ratios $ \text{Br}_{i} $ as a function of $ m_{\chi} $. The differential photon yield per branon annihilation $ \text{d}N/\text{d}E $ is depicted for a set of branon DM masses (from light to dark: 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50 and $ \unit[100]{TeV} $) in the right panel of Fig.~\ref{fig:Branching_ratios_and_spectra}.

\begin{figure}[h]
    \centering
    \begin{subfigure}{.49\textwidth}
        \includegraphics[width=\textwidth]{branon_br_plot.pdf}
    \end{subfigure}
    \begin{subfigure}{.49\textwidth}
        \includegraphics[width=\textwidth]{branon_spectra_plot.pdf}
    \end{subfigure}
\caption{Taken from JCAP05(2022)005. See text for more details.}
\label{fig:Branching_ratios_and_spectra}
\end{figure}

In our analysis, we are using the \textit{J}-factor and its statistical uncertainty for Segue~1 from~\cite{2015ApJ...801...74G} with the value of $ \mathrm{log}_{10} \left( J \left( \theta \right) [\unit{GeV^{2} cm^{-5}}] \right) = 19.02_{-0.35}^{+0.32} $. The  differential gamma-ray yields per annihilation $ \text{d}N_{i}/\text{d}E $ (see Eq. \ref{eq:Branon_Flux}) are taken from the PPPC 4 DM ID distribution~\cite{2011JCAP...03..051C}.

\section{Likelihood analysis method}
The data reduction of the Segue~1 observation have been published by the MAGIC Collaboration~\cite{2014JCAP...02..008A,2016JCAP...02..039M}. We re-analyse those high-level data in the context of brane-world extra-dimensional theories using the likelihood analysis described in Aleksi\'{c}, Rico and Martinez~\cite{2012JCAP...10..032A}. The binned ($ N_{\text{bins}} = 30 $) likelihood function of the dataset $ \bm{\mathcal{D}} $ with nuisance parameters $ \bm{\nu} $ reads as:
\begin{equation}
    \label{eq:Binned_lkl}
    \begin{split}
    \mathcal{L}_{\text{bin}} \left( \langle \sigma v \rangle ; \bm{\nu} \mid \bm{\mathcal{D}} \right) &= \mathcal{L}_{\text{bin}}( \langle \sigma v \rangle; \{ b_{i} \}_{i=1,\ldots,N_{\text{bins}}}, \tau \mid \{ N_{\text{ON},i}, N_{\text{OFF},i} \}_{i=1,\ldots,N_{\text{bins}}}) \\
    &= \prod_{i=1}^{N_{\text{bins}}} \Big[ \mathcal{P} (s_{i}(\langle \sigma v \rangle) + b_{i} \mid N_{\text{ON},i}) \cdot \mathcal{P} (\tau b_{i} \mid N_{\text{OFF},i}) \Big] \times \mathcal{T} \left( \tau \mid \tau_{\text{o}}, \sigma_{\tau} \right)
    \end{split}
\end{equation}
\noindent
where $ \mathcal{P} (x | N) $ is the Poisson distribution of mean $x $ and measured value $ N $ and $ s_{i}(\langle \sigma v \rangle) $ and $ b_{i} $ are the expected numbers of signal and background events in the $ i $-th energy bin. Besides $ b_{i} $, the normalization between background and signal regions $ \tau $, described by the likelihood function $ \mathcal{T} $, is also a nuisance parameter in the analysis~\cite{2022PDU....3500912A}. The total number of observed events in a given energy bin $ i $ in the signal (ON) and background (OFF) regions are $ N_{\text{ON},i} $, $ N_{\text{OFF},i} $, respectively.

The Segue~1 observational campaign $ \bm{\mathcal{D}}_{\mathrm{Segue1}} $ results in $ N = 8 $ distinct datasets with an individual set of instrument response functions (IRFs). The joint likelihood function
\begin{equation}
    \mathcal{L}\left( \langle \sigma v \rangle ; J, \bm{\nu} \mid \bm{\mathcal{D}}_{\text{Segue1}} \right) = \prod_{k=1}^{N} \Big[ \mathcal{L}_{\text{bin},k} \left( \langle \sigma v \rangle, \bm{\nu_{k}} \mid \bm{\mathcal{D}_{k}} \right) \Big] \times \mathcal{J} \left( J \mid J_{\text{o}}, \sigma_{\log_{10}J} \right)
\end{equation}
\noindent
is the product of the likelihood function of each dataset. We treat the \textit{J}-factor as a nuisance parameter and include the likelihood $ \mathcal{J} $ for the \textit{J}-factor~\cite{2015PhRvL.115w1301A}. $ \bm{\nu_{k}} $ represents the set of nuisance parameters different from the \textit{J}-factor affecting the analysis of the $ k $-th dataset $ \bm{\mathcal{D}_{k}} $.

The same analysis technique, implemented in the open-source DM analysis software packages~\cite{2021arXiv211201818M} \texttt{gLike}~\cite{Rico:gLike} and \texttt{LklCom}~\cite{Miener:LklCom}, has been used by the MAGIC Collaboration~\cite{2022PDU....3500912A,2018JCAP...03..009A} and other gamma-ray observatories in several DM searches~\cite{Oakes:2019,Armand:2021}.

\section{Results}

As reported in~\cite{2022JCAP...05..05}, we present the observational $ 95 \% $ confidence level (CL) upper limits on the thermally-averaged cross-section $ \langle \sigma v \rangle $ (see left panel of Fig.~\ref{fig:Limits}) and on the brane tension $ f $ (see right panel of Fig. \ref{fig:Limits}) for branon DM annihilation obtained with almost 160 hours of Segue~1 observation with the MAGIC telescopes. We perform a model dependent search for branon DM particles of masses between $ \unit[100]{GeV} $ and $ \unit[100]{TeV} $. Previous branon DM limits from CMS~\cite{2016physletb.755.102} (blue) and AMS-02~\cite{2019physletb.790.345} (orange) as well as the prospects of the future CTA~\cite{2020JCAP...10..041A} (purple) and SKA~\cite{2020PDU....2700448C} (yellow) are depicted in Fig.~\ref{fig:Limits}.

As expected from the no significant gamma-ray excess in the Segue~1 data~\cite{2014JCAP...02..008A}, our constraints for branon DM annihilation are located within the $ 68 \% $ containment band, which is consistent with the no-detection scenario. We obtain our strongest limit $ \langle \sigma v \rangle \simeq \unit[1.4 \times 10^{-23}]{cm^{3}s^{-1}} $ for a $ \sim \unit[0.7]{TeV} $ branon DM particle mass.

\begin{figure}[h]
    \centering
    \begin{subfigure}{.49\textwidth}
        \includegraphics[width=\textwidth]{branon_final_plot.pdf}
    \end{subfigure}
    \begin{subfigure}{.49\textwidth}
        \includegraphics[width=\textwidth]{brane_tension_final_plot.pdf}
    \end{subfigure}
\caption{Taken from JCAP05(2022)005. See text for more details.}
\label{fig:Limits}
\end{figure}

\section{Conclusion and outlook}

This work leads to the most constraining branon DM limits in the multi-TeV mass range, superseding previous constraints by CMS~\cite{2016physletb.755.102} and AMS-02~\cite{2019physletb.790.345} above $ \sim \unit[1]{TeV}$. We can achieve even more stringent exclusion limits by adding further dSph observations of the MAGIC telescopes~\cite{2022PDU....3500912A} or other gamma-ray telescopes~\cite{Oakes:2019,Armand:2021} to this analysis.

\section*{Acknowledgments}
The authors thank the MAGIC Collaboration for providing private data.\newline
TM acknowledges support from PID2019-104114RB-C32. DN and TM acknowledge support from the former {\em Spanish Ministry of Economy, Industry, and Competitiveness / European Regional Development Fund} grant FPA2015-73913-JIN.\newline
VG’s contribution to this work has been supported by Juan de la Cierva-Incorporaci\'on IJC2019-040315-I grant, and by the PGC2018-095161-B-I00 and CEX2020-001007-S projects, both funded by MCIN/AEI/10.13039/501100011033 and by "ERDF A way of making Europe". VG thanks J.A.R. Cembranos for useful discussions. \newline
DK is supported by the European Union's Horizon 2020 research and innovation programme under the Marie Sk\l{}odowska-Curie grant agreement No. 754510. DK and JR acknowledge the support from MCIN/AEI/ 10.13039/501100011033 under grants PID2019-107847RB-C41 and SEV-2016-0588 ("Centro de Excelencia Severo Ochoa"), and from the CERCA institution of the Generalitat de Catalunya.

\section*{References}