# Few-shot Non-line-of-sight Imaging
with Signal-surface Collaborative Regularization

Xintong Liu1, Jianyu Wang1, Leping Xiao1, Xing Fu1, Lingyun Qiu1,2, Zuoqiang
Shi1,2

1 Tsinghua University
2 Yanqi Lake Beijing Institute of Mathematical Sciences and Applications

###### Abstract

The non-line-of-sight imaging technique aims to reconstruct targets from
multiply reflected light. For most existing methods, dense points on the relay
surface are raster scanned to obtain high-quality reconstructions, which
requires a long acquisition time. In this work, we propose a signal-surface
collaborative regularization (SSCR) framework that provides noise-robust
reconstructions with a minimal number of measurements. Using Bayesian
inference, we design joint regularizations of the estimated signal, the 3D
voxel-based representation of the objects, and the 2D surface-based
description of the targets. To our best knowledge, this is the first work that
combines regularizations in mixed dimensions for hidden targets. Experiments
on synthetic and experimental datasets illustrated the efficiency and
robustness of the proposed method under both confocal and non-confocal
settings. We report the reconstruction of the hidden targets with complex
geometric structures with only $5\times 5$ confocal measurements from public
datasets, indicating an acceleration of the conventional measurement process
by a factor of 10000. Besides, the proposed method enjoys low time and memory
complexities with sparse measurements. Our approach has great potential in
real-time non-line-of-sight imaging applications such as rescue operations and
autonomous driving.

##  1 Introduction

The non-line-of-sight (NLOS) imaging technique enables reconstructions of
targets out of the direct line of sight, which is attractive in various
applications such as autonomous driving, remote sensing, rescue operations and
medical imaging [33, 35, 40, 4, 14, 10, 39, 38, 18, 15, 26, 5, 34, 20]. A
typical scenario of NLOS imaging is shown in Figure 1. Several points on the
visible surface are illuminated by a laser and the back-scattered light from
the target is detected to reconstruct the target. The time-correlated single-
photon counting (TCSPC) technique is applied in the detection process due to
the extremely low photon intensity after multiple diffuse reflections. In
practice, a single-photon avalanche diode (SPAD) in the Geiger-mode can be
used to record the photon events with time-of-flight (TOF) information [2].
The first experimental demonstration of NLOS imaging dates back to 2012, where
the targets are reconstructed with the back-projection (BP) method [37].
Extensions of this approach include its fast implementation [1], the filtering
technique for reconstruction quality enhancement [16], and weighting factors
for noise reduction [11]. The BP-type methods have high computational
complexity, resulting in slow reconstruction processes.

Figure 1: A typical non-line-of-sight imaging scenario. a) The time resolved
signals are measured at only $3\times 3$ focal points. b) The three views of
the reconstructed target obtained with the proposed signal-surface
collaborative regularization (SSCR) method.

A number of efficient methods have been designed for fast reconstructions. The
light cone transform (LCT) method [30] formulates the physical model as a
convolution operator, so that the reconstructions can be obtained using the
Wiener deconvolution method with fast Fourier transform. The directional light
cone transform (D-LCT) [42] generalizes the LCT and reconstructs the albedo
and surface normal simultaneously. The method of frequency wavenumber
migration (F-K) [19] formulates the propagation of light using the wave
equation, and also provides a fast inversion algorithm with the frequency-
domain interpolation technique. Whereas the LCT, D-LCT and F-K methods only
work directly in confocal measurement scenarios, the phasor field (PF) method
[32, 23, 22] converts the NLOS imaging scenarios to LOS cases and works for
the general non-confocal setting with low computation complexity. For high-
quality and noise-robust reconstructions, the signal-object collaborative
regularization (SOCR) method can be applied, but brings additional
computational cost. In recent years, deep learning-based methods are also
introduced to the field of NLOS imaging [6, 7, 27, 43]. Besides, advances in
hardware enhance the distance of NLOS detection to kilometers [39], or make it
possible to reconstruct targets on the scale of millimeters[38].

Despite these breakthroughs, the trade off between the acquisition time and
the imaging quality is inevitable. In the raster scanning mode, the
acquisition time is proportional to the number of measurement points and the
number of pulses used for each illumination. Due to the intrinsic ill-
posedness of the NLOS reconstruction problem[21] and heavy measurement
noise[11], dense measurements are necessary for high quality
reconstructions[30, 19, 22]. The measurement process may take from seconds to
hours, which poses a great challenge for applications such as autonomous
driving, where real-time reconstruction of the video stream is needed. The
acquisition process can be accelerated by reducing the number of pulses used
for each illumination point. In the work [17], the pulse number that record
the first returning photon is used to reconstruct the target. Another way to
reduce the acquisition time is to design array detectors for non-confocal
measurements. For example, the implementation of the phasor field method with
SPAD arrays realizes low-latency real-time video imaging of the hidden scenes
[28]. A third way to accelerate the NLOS detection process is to reduce the
number of measurement pairs. It is shown that $16\times 16$ confocal
measurements are enough to reconstruct the hidden target by incorporating the
compressed sensing technique [41].

In this paper, we study the randomness in the photon detection process of non-
line-of-sight scenarios and propose an imaging method that deals with a very
limited number of measurements. We design joint regularizations of the
estimated signal, the 3D voxel-based representation of the objects, and the 2D
surface-based description of the targets, which leads to faithful
reconstruction results. The main contributions of our work are as follows.

  * •

We propose a novel signal-surface collaborative regularization (SSCR)
framework for few-shot non-line-of-sight reconstructions, which works under
both confocal and non-confocal settings.

  * •

We report the reconstruction of the hidden targets with complex geometric
structures with only $5\times 5$ confocal measurements from public datasets,
indicating an acceleration of the conventional measurement process by a factor
of 10000.

Figure 2: The least squares solution of the statue with $3\times 3$ confocal
measurements[19]. The solution does not contain useful information of the
target even though its simulated signal matches the measurements well (see
Fig. 3). Strong regularizations are needed to reconstruct the target. See also
Fig. 6 for a comparison. Figure 3: Comparisions of the measured data and the
simulated data of the least squares solution for the instance of the
statue[19]. The measured signals are shown in black. The simulated data of the
least squares solution are shown in red. The shapes of the signals are very
close to each other.

##  2 Related work

From a mathematical point of view, the non-line-of-sight imaging task belongs
to the category of inverse problem. The goal is to reconstruct the surface of
the hidden target with the measured signal. The NLOS inverse problem is ill-
posed due to the intrinsic structure of the physical model[21] and heavy
measurement noise[11]. When the number of measurement points is small, the
lack of data leads to rank-deficiency of the measurement matrix, making the
reconstruction task even harder[41]. In such cases, regularization becomes a
vital tool for high-quality reconstructions. In Fig. 2 and Fig. 3 we show the
least squares solution without regularization for the instance of the statue
with $3\times 3$ confocal measurements [19] and compare the simulated signal
with the raw measurements. The simulated signal matches the input data well,
but the target cannot be identified in the reconstruction. More details are
provided in the supplement. Here, we review two regularization based methods
that provide high-quality reconstructions. The proposed method is closely
related with these methods.

The SOCR method The signal-object collaborative regularization (SOCR) method
[24] represents the hidden target with voxels and considers prior terms for
both the reconstructed target and the signal. The joint regularization term
concerns sparseness of the reconstruction, non-local self-similarity of the
target and the smoothness of the signal. The method employs $L_{1}$
regularization, orthogonal optimization and empirical Wiener filter to provide
noiseless reconstructions with sharp boundaries, but at the expense of
additional computational cost.

The CC-SOCR method The SOCR method only works in cases where the signals are
measured at rectangular grid points. For generalizations to cases with
arbitrary illumination and detection pattern, the confocal-complemented
signal-object collaborative regularazation (CC-SOCR) method [25] considers the
virtual confocal signal measured at a focal plane to overcome the rank-
deficiency. Joint sparse representation of the simulated signal and the
virtual confocal signal is used to increase robustness. It was reported that
the method can reconstruct the targets with 200 confocal measurements.

Despite these considerations, reconstruction of the hidden targets may fail
when only a very limited number of measurements (for example, $5\times 5$
confocal measurements) are available. In such cases, the ill-posedness of the
NLOS inverse problem is so strong that these methods provide biased
estimations of the targets. We propose a novel signal-surface collaborative
regularization (SSCR) method that provides faithful reconstructions under few-
shot scenarios. The proposed approach are closely related with SOCR and CC-
SOCR methods but improves these methods in the following two aspects.

  * •

The data fidelity term These methods use the quadratic data misfit as the data
fidelity terms, which assume noise of the Gaussian type. In this work, we
study the randomness of photon event stamping, and use logarithmic data
fidelity terms derived from Bernoulli photon event assumptions.

  * •

The dimensionality of priors These methods represent the hidden targets with
three-dimensional voxels and the regularizations apply to three-dimensional
tensors. In fact, we aim at reconstructing the surface of the hidden target,
which is a two-dimensioanl geometric object in nature. In this work, we study
the intrinsic two-dimensional regularization for the hidden surface and
incorporate priors in both 2D and 3D representations of the targets.

##  3 The physical model

In general, the reconstruction quality and the computation complexity of an
imaging algorithm rely heavily on the forward model that simulates the
physical measurement process. Fine physical models lead to reconstructions
with clear geometric structures, but at the expense of high computational
cost[12, 36]. Instead of putting forward a novel physical model, we aim at
overcoming the ill-posedness of the NLOS imaging task with an extremely small
number of measurements. We adopt a linear physical model which only considers
the square fall-off of the photon intensity. Let $x^{\prime}_{i}$ and
$x^{\prime}_{d}$ be the illumination and detection points on the visible
surface, the photon intensity detected at time $t$ is modeled as

$\displaystyle\tau(x^{\prime}_{i},x^{\prime}_{d},t)=$
$\displaystyle\int_{\Omega}{\frac{u(x)}{\|x^{\prime}_{i}-x\|^{2}\|x^{\prime}_{d}-x\|^{2}}}\cdot$
(1) $\displaystyle\ \ \ \ \ \ \ \ \ \
\delta(\|x^{\prime}_{i}-x\|+\left\|x^{\prime}_{d}-x\right\|-ct)dx,$

in which $c$ is the speed of light, $x$ is a point in the three-dimensional
reconstruction domain $\Omega$. The albedo value of the point $x$ is
represented by $u(x)$. The $\delta$ function describes the intrinsic domain of
integration as the set of points with optical path length $ct$. When
$x^{\prime}_{i}\neq x^{\prime}_{d}$, the domain of integration is a half
ellipsoid with foci $x^{\prime}_{i}$ and $x^{\prime}_{d}$. When
$x^{\prime}_{i}=x^{\prime}_{d}=x^{\prime}$, the domain of integration is a
half sphere with center $x^{\prime}$ and the model reduces to the one used in
the work [30] given below

$\tau(x^{\prime},t)=\int_{\Omega}{\frac{u(x)}{\|x^{\prime}-x\|^{4}}}\cdot\delta(2\|x^{\prime}-x\|-ct)dx.$
(2)

##  4 The proposed framework

In this section, we study the randomness in the measurement process and
propose a novel signal-surface collaborative regularization (SSCR) framework
for few-shot NLOS imaging.

Notation The $L_{2}$ norm of a vector $\mathbf{x}$ is denoted by
$\|\mathbf{x}\|$. We use $[N]$ as an abbreviation of the set
$\\{1,2,\dots,N\\}$. The three-dimensional reconstruction domain $\Omega$ is
discretized with voxels
$V=\\{v_{ijk}=(x_{i},y_{j},z_{k})|i\in[I],j\in[J],k\in[K]\\}$, where $x_{i}$,
$y_{j}$, and $z_{k}$ are the coordinates of the point $v_{ijk}$ in the
horizontal, vertical and depth directions. Each point in $V$ represents a
cubic voxel of the same size. The grid function $\mathbf{u}$ is used to denote
the discrete albedo values, with its components denoted by
$u_{ijk}=\mathbf{u}(v_{ijk})$. To reconstruct the hidden targets, the signal
is measured at $P$ pairs of points
$\\{(x^{\prime}_{p},y^{\prime}_{p})\\}^{P}_{p=1}$, in which $x^{\prime}_{p}$
and $y^{\prime}_{p}$ are the coordinates of the $p^{th}$ illumination and
detection points, respectively. For each measurement pair
$(x^{\prime}_{p},y^{\prime}_{p})$, the one-dimensional time resolved signal
contains $Q$ time bins. The length of each time bin is a constant, usually at
the scale of picoseconds in real applications. We denote by $\tau_{p,q}$ the
photon intensity detected at the $p^{th}$ measurement pair and the $q^{th}$
time bin. The measured dataset is denoted by $\bm{\tau}$.

The reconstructed surface We use $\mathbf{u}\in\mathbb{R}^{I\times J\times K}$
to represent the albedo of the NLOS scene, which is a three-diminsional
tensor. However, it is only possible to reconstruct the portion of the hidden
surface where photons are bounced back, which is a two dimenional geometric
object. With this observation, we define a subset
$\mathcal{G}\subsetneq\mathbb{R}^{I\times J\times K}$ as follows

$\displaystyle\mathcal{G}=$
$\displaystyle\\{\mathbf{g}=(g_{ijk})\in\mathbb{R}^{I\times J\times K}\ |\
\forall(i,j)\in[I]\times[J],$ (3) $\displaystyle\ \ \ \ \ \ \exists\text{ at
most one }k=k_{ij}\in[K],\text{ s.t. }g_{i,j,k_{ij}}\neq 0\\}.$

Each element $\mathbf{g}\in\mathcal{G}$ yields a trival two-dimensional
parameterization. For each pixel $(i,j)\in[I]\times[J]$, only one of the
following cases holds.

  * •

Case 1: $g_{ijk}=0$ for all $k\in[K]$. In this case, the line
$x=x_{i},y=y_{j}$ does not intersect with the target and we call $(i,j)$ a
_background pixel_.

  * •

Case 2: There exists only one $k=k_{ij}$, such that $g_{i,j,k_{ij}}>0$. In
this case, we call $(i,j)$ a _foreground pixel_ and the corresponding $depth$
is $z_{ij}$.

To express elements of $\mathcal{G}$ with matrices, we show that there is a
bijection from the set $\mathcal{G}$ to the following set

$\displaystyle\mathcal{G}^{\prime}=\\{(\mathbf{e},\mathbf{d},\bm{\alpha})\ |\
\mathbf{e}$ $\displaystyle=(e_{ij})_{I\times J},e_{ij}\in\\{0,1\\},$ (4)
$\displaystyle\mathbf{d}$ $\displaystyle=(d_{ij})_{I\times
J},d_{ij}\in[K]\cup\\{\operatorname{NaN}\\},$ $\displaystyle\bm{\alpha}$
$\displaystyle=(\alpha_{ij})_{I\times
J},\alpha_{ij}\in{\mathbb{R}}\cup\\{\operatorname{NaN}\\},$
$\displaystyle\forall(i,j),\ \ e_{ij}=0$ $\displaystyle\iff
d_{ij}=\operatorname{NaN}\iff\alpha_{ij}=\operatorname{NaN}\\}.$

In the definition of $\mathcal{G}^{\prime}$, the placeholder NaN represents a
background pixel and does not operate with real numbers. To construct a
bijection from $\mathcal{G}$ to $\mathcal{G}^{\prime}$, for each element
$g\in\mathcal{G}$, let $\mathbf{e}$ be the indicator function of the set of
foreground pixels, $\mathbf{d}$ be the depths of the foreground pixels and
$\bm{\alpha}$ be the corresponding albedo values. Then, fill the matrices
$\mathbf{d}$ and $\bm{\alpha}$ with NaNs where necessary. It is easy to check
that this map: $\mathcal{G}\rightarrow\mathcal{G}^{\prime}$ is one to one and
onto. We call $\mathbf{e}$, $\mathbf{d}$ and $\bm{\alpha}$ the indicator
matrix, the depth matrix and the albedo matrix of $\mathbf{g}$, respectively.
For simplicity, we also denote $\mathcal{G}^{\prime}$ as $\mathcal{G}$. The
set $\mathcal{G}$ is still too large to capture the set of real-world
surfaces, as the distribution of the foreground pixels and depths can be
disorganized. In the following sections, we will design regularizations for
elements of $\mathcal{G}$, where the matrix representations bring remarkable
convenience.

The photon event histogram In NLOS scenarios, the intensity of the back
scattered light is extremely weak after multiple diffuse reflections. For each
measurement pair, a total of $N$ laser pulses are emitted to the illumination
point, and the TCSPC device can be used to record the photon events. Commonly
used detectors are single photon avalanche diodes (SPAD). Although heavily
corrupted with background noise, the photon event histogram provides an
estimation of the photon intensity. We use the binary variables $d_{p,q,n}$ to
denote the recorded photon events.

$d_{p,q,n}=\left\\{\begin{array}[]{l}1,\ \mathrm{record}\ \mathrm{a}\
\mathrm{photon}\ \mathrm{event}\\\ 0,\ \mathrm{otherwise}\\\
\end{array}\right.,$ (5)

in which $p\in[P]$, $q\in[Q]$, and $n\in[N]$ are indices of the measurement
pair, the time bin and the pulse number. For each measurement pair, the
detector can record at most one photon event in each pulse, which means that
$\sum_{q}d_{p,q,n}\leq 1$. However, the cases of recording more than one
photon event in a single pulse can be neglected [17]. The collection of photon
event stamping is denoted by $\mathbf{d}$.

The Bayesian framework We propose a unified Bayesian framework that
reconstructs the hidden surface with the data of photon event stamping. For
each measurement pair $p$, time bin $q$, and pulse number $n$, it is assumed
that the detection of a photon event $e_{p,q,n}$ follows the Bernoulli
distribution with probability
$\mathbb{P}\left\\{e_{p,q,n}=1\right\\}=1-e^{-\eta\tau_{p,q}}$, in which
$\eta>0$ is the detection efficiency [30]. The collection of random variables
$e_{p,q,n}$ is denoted by $\mathbf{e}$. In NLOS detection scenarios, the
probability of detecting a photon event is extremely small. The first order
approximation of the exponent is adopted and we assume

$\mathbb{P}\left\\{e_{p,q,n}=1\right\\}=\eta\tau_{p,q}.$ (6)

In Eq. 1, the photon intensity is linear with respect to the albedo, we choose
$\eta=1$ without loss of generality. We also assume that the detections of
different photon events are independent. Let $\mathbf{g}\in\mathcal{G}$ be the
reconstructed target. We view $\mathbf{g}$ and $\bm{\tau}$ as random vectors
and find them simultaneously by maximizing the posterior probability
$\mathbb{P}(\mathbf{g},\bm{\tau}|\mathbf{e}=\mathbf{d})$, where $\mathbf{e}$
is related with $\bm{\tau}$ by Eq. 6. Noting that the set $\mathcal{G}$ is not
a convex subset of $\mathbb{R}^{I\times J\times K}$, the resulting
optimization problem is non-convex and hard to solve. Besides, the forward
model we adopt is an approximation of the ideal one, and contains inevitable
bias. To tackle these problems, we introduce the random vector
$\mathbf{u}\in\mathbb{R}^{I\times J\times K}$ as an approximation of the
surface $\mathbf{g}$ and maximize
$\mathbb{P}(\mathbf{g},\mathbf{u},\bm{\tau}|\mathbf{e}=\mathbf{d})$. Using the
Bayesian formula, we obtain

$\displaystyle\mathop{\arg\max}_{\mathbf{g},\mathbf{u},\bm{\tau}}\mathbb{P}(\mathbf{g},\mathbf{u},\bm{\tau}|\mathbf{e}=\mathbf{d})$
(7) $\displaystyle=$
$\displaystyle\mathop{\arg\max}_{\mathbf{g},\mathbf{u},\bm{\tau}}\mathbb{P}(\mathbf{e}=\mathbf{d}|\mathbf{g},\mathbf{u},\bm{\tau})\mathbb{P}(\mathbf{g},\mathbf{u},\bm{\tau})$
$\displaystyle=$
$\displaystyle\mathop{\arg\max}_{\mathbf{g},\mathbf{u},\bm{\tau}}\mathbb{P}(\mathbf{e}=\mathbf{d}|\bm{\tau})\mathbb{P}(\mathbf{g},\mathbf{u},\bm{\tau})$
$\displaystyle=$
$\displaystyle\mathop{\arg\max}_{\mathbf{g},\mathbf{u},\bm{\tau}}\prod_{p,q,n}{\mathbb{P}(\mathbf{e}_{p,q,n}=\mathbf{d}_{p,q,n}|\tau_{p,q})}\mathbb{P}(\mathbf{g},\mathbf{u},\bm{\tau})$
$\displaystyle=$
$\displaystyle\mathop{\mathrm{arg}\max}_{\mathbf{g},\mathbf{u},\bm{\tau}}\prod_{p,q}{(\tau_{p,q})^{d_{p,q}}(1-\tau_{p,q})^{N-d_{p,q}}}\mathbb{P}(\mathbf{g},\mathbf{u},\bm{\tau})$
$\displaystyle=$
$\displaystyle\mathop{\mathrm{arg}\min}_{\mathbf{g},\mathbf{u},\bm{\tau}}\ \
\sum_{p,q}[{(d_{p,q}-N)\ln(1-\tau_{p,q})}$ $\displaystyle\ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \
-d_{p,q}\ln(\tau_{p,q})]+\Gamma(\mathbf{g},\mathbf{u},\bm{\tau}),$

where $d_{p,q}=\sum_{n=1}^{N}d_{p,q,n}$ is the data of photon event histogram.
In the second equality we assume that the conditional probability of
$\mathbf{e}$ only depends on $\bm{\tau}$.
$\Gamma(\mathbf{g},\mathbf{u},\bm{\tau})$ is the joint prior term of
$\mathbf{g}$, $\mathbf{u}$ and $\bm{\tau}$.

##  5 The joint regularization term

The joint regularization term $\Gamma(\mathbf{g},\mathbf{u},\bm{\tau})$ plays
a crucial role in the process of reconstruction. An ingenious design of this
term not only results in faithful reconstructions, but also admits low-cost
algorithms to solve the optimization problem (7). In this paper, we assume

$\displaystyle\Gamma(\mathbf{g},\mathbf{u},\bm{\tau})=$
$\displaystyle\lambda\|\bm{\tau}-A\mathbf{u}\|^{2}$ (8) $\displaystyle+$
$\displaystyle
J_{1}(\mathbf{u})+J_{2}(\mathbf{u},\bm{\tau})+J_{3}(\mathbf{u},\mathbf{g}),$

in which $A$ is the forward operator defined by Eq. 1. $J_{1}$ is the prior of
$\mathbf{u}$, $J_{2}$ is the joint prior of $\mathbf{u}$ and $\bm{\tau}$,
$J_{3}$ is the joint prior of $\mathbf{u}$ and $\mathbf{g}$. $\lambda$ is a
fixed parameter.

The prior $J_{1}(\mathbf{u})$ The priors of the three-dimensional
representations of the targets have been widely used in existing works[19, 30,
13]. Two efficient priors are the sparseness and non-local self-similarity of
the objects. In the work [24], these two priors are considered for the 4D
tensor of the directional-albedo. Here, we simplify the approach and directly
use $L_{1}$ norm of the albedo to impose the sparseness of the target. For
non-local self-similarity prior, we directly follow the block-matching and
sparse representation method in the work [24]. We set

$\displaystyle J_{1}(\mathbf{u})=$ $\displaystyle s_{u}\|\mathbf{u}\|_{1}$ (9)
$\displaystyle+$
$\displaystyle\lambda_{u}\sum_{i}\left[\|\mathcal{B}\mathbf{u}_{i}-D_{s}C_{i}D_{n}^{T}\|^{2}+\lambda_{pu}|C_{i}|_{0}\right],$

in which $s_{u}$, $\lambda_{u}$, and $\lambda_{pu}$ are fixed parameters.
$|\cdot|_{0}$ denotes the number of nonzero elements. The summation is made
over all possible local 3D sub-tensors indexed by $i$. The matrix
$\mathcal{B}\mathbf{u}_{i}$ is constructed by putting the vectorizations of
the $i^{th}$ sub-tensor $\mathbf{u}_{i}$ and its neighbors column by column.
The orthogonal matrices $D_{s}$ and $D_{n}$ capture the local structure and
non-local self-similarity of the albedo $\mathbf{u}$. $C_{i}$ contains the
transform-domain coefficients of the $i^{th}$ block, whose sparseness is
imposed by the term $|C_{i}|_{0}$. For more details, we refer the readers to
[9, 3, 24].

The prior $J_{2}(\mathbf{u},\bm{\tau})$ We seek for a joint local sparse
representation scheme for the estimated signal $\bm{\tau}$ and simulated
signal $A\mathbf{u}$. It is assumed that $\bm{\tau}$ is a three dimensional
tensor of size $N_{x}\times N_{y}\times Q$, in which $N_{x}$ and $N_{y}$ are
the number of measurement points in the horizontal and vertical directions.
$Q$ is the number of time bins. We call a three-dimensional sub-tensor of
$\bm{\tau}$ a local patch. Consider the set of all possible patches of size
$r_{x}\times r_{y}\times r_{q}$. We generate the patch dataset
$\mathcal{P}(\bm{\tau})$ by generating the vectorization of each patch and put
them together column by column. In the work [25], the patch dataset of the
complete virtual confocal signal is sparsely represented by an adaptively
learned dictionary. However, when $N_{x}$ and $N_{y}$ are small, the number of
patches are not sufficient for effective dictionary learning. Here, we use the
orthogonal dictionary $D=D_{q}\otimes Dy\otimes D_{x}$ as the transform basis,
where $D_{q}$, $D_{y}$ and $D_{x}$ are matrices of the discrete consine
transform of orders $q$, $y$ and $x$, respectively. The joint prior of
$\mathbf{u}$ and $\bm{\tau}$ is given by

$\displaystyle J_{2}(\mathbf{u},\bm{\tau})=$
$\displaystyle\lambda_{t}\|\mathcal{P}(\bm{\tau})-DS\|^{2}$ (10)
$\displaystyle+$
$\displaystyle\lambda_{ut}\|\mathcal{P}(A\mathbf{u})-DS\|^{2}+\lambda_{pt}|S|_{0},$

in which $\lambda_{t}$, $\lambda_{ut}$ and $\lambda_{pt}$ are fixed
parameters. $\mathcal{P}(A\mathbf{u})$ represents the patch dataset generated
by the simulated signal of $\mathbf{u}$. $S$ contains the public transform
domain coefficients of $\mathcal{P}(\bm{\tau})$ and
$\mathcal{P}(A\mathbf{u})$, whose sparseness is imposed by the $L_{0}$ term.

The prior $J_{3}(\mathbf{u},\mathbf{g})$ We express the joint regularization
of $\mathbf{u}$ and $\mathbf{g}$ as

$J_{3}(\mathbf{u},\mathbf{g})=\lambda_{g}\left[\|\mathbf{u}-\mathbf{g}\|^{2}+\Upsilon(\mathbf{g})\right],$
(11)

in which $\Upsilon(\mathbf{g})$ describes the prior distribution of
$\mathbf{g}$ defined on $\mathcal{G}$. The set $\mathcal{G}$ is not convex,
making it difficult to design $\Upsilon(\mathbf{g})$ explicitly. In fact, it
suffcies to update $\mathbf{g}$ with fixed $\mathbf{u}$ and vise vera in the
final optimization problem. With fixed $\mathbf{g}$, $\mathbf{u}$ can be
easily updated with the term $\|\mathbf{u}-\mathbf{g}\|^{2}$. To update
$\mathbf{g}$ with any fixed $\mathbf{u}\in\mathbb{R}^{I\times J\times K}$, we
choose an element from the set $\mathcal{G}$ which not only lies in the
neighborhood of $\mathbf{u}$ in the $L_{2}$ sense, but also acts like the
surface of some real-world object. With this motivation in mind, we construct
a map $\mathcal{S}:\mathbb{R}^{I\times J\times K}\rightarrow\mathcal{G}$ and
view $\mathcal{S}(\mathbf{u})$ as the solution to the following optimization
problem.

$\mathcal{S}(\mathbf{u})=\mathop{\arg\min}_{\mathbf{g}}\|\mathbf{u}-\mathbf{g}\|^{2}+\Upsilon(\mathbf{g}).$
(12)

We call $\mathcal{S}(\mathbf{u})$ a _surfaciation_ of $\mathbf{u}$. To
construct the map $\mathcal{S}$, we assign for each
$\mathbf{u}\in\mathbb{R}^{I\times J\times K}$ an indicator matrix
$\mathbf{e}$, a depth matrix $\mathbf{d}$ and an albedo matrix $\bm{\alpha}$
(recall the definition of $\mathcal{G}^{\prime}$). We determine these three
matrices using the following approach.

For each pixel $(i,j)$ of $\mathbf{u}$, let
$u_{i,j,{k^{1}_{ij}}},\dots,u_{i,j,{k^{n_{ij}}_{ij}}}$ be all $n_{ij}$ non-
zero albedo values in the depth direction. Define

$\tilde{e}_{ij}=\left\\{\begin{array}[]{l}1,\ n_{ij}>0\\\ 0,\
n_{ij}=0\end{array}\right.$ (13)

to be the indicator function of the set of foreground pixels. There could be
many mislabeled pixels due to heavy background noise in $\mathbf{u}$. To
provide a noise-robust estimation of the indicator matrix $\mathbf{e}$, we
consider the global correlation of all pixels and solve

$\displaystyle\mathbf{e}^{*}=(e_{ij}^{*})_{I\times J}$
$\displaystyle=\mathop{\arg\min}_{\\{e_{ij}\\}}\sum^{I}_{i=1}\sum^{J}_{j=1}\gamma_{ij}(e_{ij}-\tilde{e}_{ij})^{2}$
(14) $\displaystyle+$
$\displaystyle\sum^{I}_{p=1}\sum^{J}_{q=1}\sum^{I}_{r=1}\sum^{J}_{s=1}w^{\mathbf{e}}_{pq,rs}(e_{pq}-e_{rs})^{2},$

in which $w^{\mathbf{e}}_{pq,rs}$ describes the weight of the pixels $(p,q)$
and $(r,s)$. The parameter $\gamma_{ij}$ describes the confidence of
indication of the originial pixel $\tilde{e}_{ij}$. This least squares problem
has a unique solution and can be solved using the standard LSQR method. The
indicator function $\mathbf{e}=(e_{ij})$ is determined by

$e_{ij}=\left\\{\begin{array}[]{l}1,\ e^{*}_{ij}\geq 0.5\\\ 0,\
e^{*}_{ij}<0.5\end{array}\right..$ (15)

Algorithm 1 The SSCR algorithm

0: $\mathbf{d}$, $N$.

0: $\mathbf{u},\mathbf{g}$.

$\tau_{p,q}^{0}=d_{p,q}/N$

$\mathbf{u}^{0}=\mathop{\arg\min}_{\mathbf{u}}\lambda\|\bm{\tau}^{0}-A\mathbf{u}\|^{2}+s_{u}\|\mathbf{u}\|_{1}$

for $k=1$ to $K$ do

$\mathbf{g}^{k+1}=\mathop{\arg\min}_{\mathbf{g}}\|\mathbf{u}^{k}-\mathbf{g}\|^{2}+\Upsilon(\mathbf{g})$

$(D^{k+1}_{s},\mathbf{C}^{k+1},D^{k+1}_{n})=\mathop{\arg\min}_{D_{s},\mathbf{C},D_{n}}\|\mathcal{B}\mathbf{u}^{k}_{i}-D_{s}C_{i}D_{n}^{T}\|^{2}+\lambda_{pu}|C_{i}|_{0}$

$S^{k+1}=\mathop{\arg\min}_{S}\lambda_{t}\|\mathcal{P(\bm{\tau}}^{k})-DS\|^{2}+\lambda_{ut}\|\mathcal{P}(A\mathbf{u}^{k})-DS\|^{2}+\lambda_{pt}|S|_{0}$

$\bm{\tau}^{k+1}=\mathop{\arg\min}_{\bm{\tau}}\sum_{p,q}\left[(d_{p,q}-N)\ln(1-\tau_{p,q})-d_{p,q}\ln(\tau_{p,q})\right]+\lambda_{t}\|\mathcal{P(\bm{\tau}})-DS^{k+1}\|^{2}+\lambda\|\bm{\tau}-A\mathbf{u}^{k}\|^{2}$

$\mathbf{u}^{k+1}=\mathop{\arg\min}_{\mathbf{u}}\lambda\|\bm{\tau}^{k+1}-A\mathbf{u}\|^{2}+s_{u}\|\mathbf{u}\|_{1}+\lambda_{ut}\|\mathcal{P}(A\mathbf{u})-DS^{k+1}\|^{2}+\lambda_{u}\sum_{i}\|\mathcal{B}\mathbf{u}_{i}-D^{k+1}_{s}C^{k+1}_{i}(D^{k+1}_{n})^{T}\|^{2}+\lambda_{g}\|\mathbf{u}-\mathbf{g}^{k+1}\|^{2}$

end for

To obtain the depth matrix $\mathbf{d}$, we solve the following least squares
problem

$\displaystyle\mathbf{d}^{*}=$ $\displaystyle(d_{ij}^{*})_{I\times J}$ (16)
$\displaystyle=$
$\displaystyle\mathop{\arg\min}_{\\{d_{ij}\\}}\sum^{I}_{i=1}\sum^{J}_{j=1}\sum^{n_{ij}}_{n=1}\lambda_{ijn}(d_{ij}-z_{k^{n}_{ij}})^{2}$
$\displaystyle+\sum^{I}_{p=1}\sum^{J}_{q=1}\sum^{I}_{r=1}\sum^{J}_{s=1}w^{\mathbf{d}}_{pq,rs}(d_{pq}-d_{rs})^{2},$

in which $z_{k^{n}_{ij}}$ is the depth of the voxel $(i,j,k^{n}_{ij})$.
$\lambda_{ijn}$ and $w^{\mathbf{d}}_{pq,rs}$ are fixed parameters that control
the weight of the corresponding terms. The depth matrix $\mathbf{d}=(d_{ij})$
is then determined by

$d_{ij}=\left\\{\begin{array}[]{l}\mathop{\arg\min}_{k}\|z_{k}-d^{*}_{ij}\|,\
e_{ij}=1\\\ \text{NaN},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
e_{ij}=0\end{array}\right..$ (17)

The albedo matrix $\bm{\alpha}=(\alpha_{ij})$ is obtained by solving the
following optimization problem

$\displaystyle\mathop{\min}_{\\{\alpha_{ij}\\}}$
$\displaystyle\sum^{I}_{i=1}\sum^{J}_{j=1}\sum^{n_{ij}}_{n=1}\lambda_{ijn}(\alpha_{ij}-u_{ijk^{n}_{ij}})^{2}$
(18) $\displaystyle\ \ \ \ \ \ \ \ \
+\sum^{I}_{p=1}\sum^{J}_{q=1}\sum^{I}_{r=1}\sum^{J}_{s=1}w^{\bm{\alpha}}_{pq,rs}(\alpha_{pq}-\alpha_{rs})^{2},$

in which $\lambda_{ijn}$ and $w^{\bm{\alpha}}_{pq,rs}$ are fixed parameters.
Finally, the element $\alpha_{ij}$ is reset as $\operatorname{NaN}$ if
$e_{ij}=0$.

##  6 The SSCR reconstruction algorithm

Finally we obtain the optimization problem of the proposed signal-surface
collaborative regularization (SSCR) framework as follows

$\displaystyle\mathop{\arg\min}_{\begin{subarray}{c}\bm{\tau},\mathbf{u},\mathbf{g},\\\
D_{s},D_{n},\mathbf{C},S\end{subarray}}\sum_{p,q}\left[(d_{p,q}-N)\ln(1-\tau_{p,q})-d_{p,q}\ln(\tau_{p,q})\right]$
(19)
$\displaystyle+\lambda_{t}\|\mathcal{P}(\bm{\tau})-DS\|^{2}+\lambda_{ut}\|\mathcal{P}(A\mathbf{u})-DS\|^{2}+\lambda_{pt}|S|_{0}$
$\displaystyle+\lambda\|\bm{\tau}-A\mathbf{u}\|^{2}+s_{u}\|\mathbf{u}\|_{1}+\lambda_{g}\left[\|\mathbf{u}-\mathbf{g}\|^{2}+\Upsilon(\mathbf{g})\right]$
$\displaystyle+\lambda_{u}\sum_{i}\left[\|\mathcal{B}\mathbf{u}_{i}-D_{s}C_{i}D_{n}^{T}\|^{2}+\lambda_{pu}|C_{i}|_{0}\right]$
$\displaystyle s.t.\ \ g\in{\mathcal{G}},\ \ D_{s}D_{s}^{T}=I_{x},\ \
D_{n}D_{n}^{T}=I_{y}.$

in which $x$ is the product of the block size in three directions, $y$ is the
number of neighbors for each block. At first sight, this problem contains many
terms and seems complicated. In fact, it can be solved using the standard
alternating iteration method and solutions to most sub-problems have already
been well studied[24, 25]. We summarize the main step of the proposed SSCR
algorithm in Algorithm 1. In the supplementary material, we provide a detailed
discussion of the solution to this optimization problem and the choices of
parameters.

##  7 Results

To validate the capability of the proposed method in reconstructing the hidden
targets with sparse measurements, we compare our reconstruction results with
the Laplacian of Gaussian filtered back-projection (LOG-BP)[16], SOCR[24] and
CC-SOCR[25] methods. We also bring the F-K [19], LCT [30], and the D-LCT [42]
methods into comparisions by constructing dense measurements with the commonly
used zero-padding technique. In the supplement, we provide a gallery of
reconstruction results of these methods with other signal interpolation
techniques, where similar reconstruction results are shown.

Synthetic data We use the synthetic signal of the instance of the pyramid [24]
to test the proposed method. Only $3\times 3$ of the original $64\times 64$
synthetic signals are chosen. The physical model used to generate the data
considers consine attenuation of the photon intensity, and is finer than Eq.
2. The base length and height of the pyramid are 1 m and 0.2 m, respectively.
The central axis of the regular quadrangular pyramid is vertical to the planer
relay surface. The pyramid is 0.5 m away from the relay surface and the time
resolution is 32 ps. The experimental setup and the three views of our
reconstruction are shown in Fig. 1. Figure 4 shows the reconstruction results
of different methods. The LOG-BP, F-K, LCT and D-LCT methods fail to locate
the hidden object. The SOCR and CC-SOCR methods locate the target correctly,
but contain artifacts in the background. The proposed SSCR method provides a
faithful estimation of the hidden target. More detailed comparisons are
provided in the supplement.

Figure 4: Comparisions of the reconstruction results of the pyramid with
confocal synthetic signal. a) The signal is measured at only $3\ \times\ 3$
focal points. b) - h) The front view of the reconstructions obtained with the
LOG-BP, F-K, LCT, D-LCT, SOCR, CC-SOCR, and the proposed algorithms Figure 5:
Comparisions of the reconstruction results of the statue with $5\times 5$
confocal measurements[19]. a) A photo of the hidden object. b) - h) The front
view of the reconstructions obtained with the LOG-BP, F-K, LCT, D-LCT, SOCR,
CC-SOCR, and the proposed algorithms Figure 6: Reconstruction results of the
proposed method with different illumination settings. The illumination points
are shown in yellow in the first column. The front view, top view and side
view of the reconstructions are shown in the second, third and fourth columns.

Measured data We use the measured data of the instance of the statue [19] to
test the proposed method in real-world applications. The original dataset
contains $512\times 512$ confocal signals measured over a square region of
$2\times 2$ $\text{m}^{2}$ in 10 mins. The hidden target is 1 m from the
visible surface. The time resolution is 32 ps. Reconstruction results with
$5\times 5$ confocal measurements are compared in Fig. 5. The focal points are
shown in the third row and first column of Fig. 6. It is shown that the LOG-BP
and F-K methods fail to reconstruct the target. The LCT and D-LCT
reconstructions contain heavy noise and the target can hardly be identified.
The SOCR method only reconstructs a portion of the target. The CC-SOCR
reconstruction is blurry and contains background noise. In contrast, the
proposed SSCR method reconstructs the target faithfully. In Fig. 6 we show
three views of the SSCR reconstructions with different number of illumination
points. The target can be clearly reconstructed with $7\times 7$ confocal
measurements, which is only 0.01% of the original dataset. With $3\times 3$
illuminations, the SSCR method still provides a reasonable estimation of the
hidden target, which demonstrates the robustness of the proposed algorithm.
More comparisions with existing methods are provided in the supplement.

Figure 7: Reconstruction results of the figure ‘4’ with non-confocal
measurements[22]. The SSCR method provides an estimation of the target with
only $6\times 6$ measurements.

Measured data We use the measured data for the instance of the figure ‘4’
provided by the work [22] to test the SSCR method under the non-confocal
setting. The target is 1 m away from the visible surface and the time
resolution is 16 ps. The illumination pattern is shown in the fourth role of
Fig. 6 and the three views of our reconstruction are shown in Fig. 7.
Comparisions with the LOG-BP[16], phasor field[22] and other methods[24, 25]
are provided in the supplement, where the proposed method performs the best.

##  8 Discussion

When the signal is detected at $P$ measurement pairs and the reconstruction
domain is discretized with $L\times L\times L$ voxels, the time and memory
complexities of the proposed SSCR method are $\mathcal{O}(PL^{3})$ and
$\mathcal{O}(L^{3})$ respectively. In the supplement, we provide a detailed
analysis of these complexities. Notebly, when $P=\mathcal{O}(1)$, the time
complexity of our method is $\mathcal{O}(L^{3})$. This is smaller than the
F-K, LCT, and D-LCT methods, which cost $\mathcal{O}(L^{3}\log L)$. In this
case, the SOCR method also has a time complexity of $\mathcal{O}(L^{3})$ but
provides biased estimations of the targets. This fact indicates the
significance of the proposed two-dimensional regularization of the hidden
surface. As is shown in Fig. 4 and Fig. 5, the CC-SOCR method performs the
second. This method relys on the dense virtual confocal signal measured over a
sqaure region, which leads to a high cost of $\mathcal{O}(L^{5})$.

We also discuss some potential directions of research that help improve the
proposed method. In Eq. 6, we did not consider the dark count of the detector
and the background noise. Modeling the noise distribution within the Bayesian
framework would lead to better reconstruction results. Besides, for the two-
dimensional regularization of $\mathbf{g}$, we simply use pixel-wise
similarity in the SSCR method. The two-dimensional image priors have been
widely studied in the past decades, including non-local self-similarity[9],
low-rankness[31], low-dimensionality[29], just to name a few. Incorporating
these regularizations would lead to better estimations of the hidden surface.
In recent years, it is shown that the family of tensor network proves to be a
powerful tool for dimensionality reduction and data processing[8]. It would
also be an interesting direction of research to reduce the dimension of the
object and the signal using tensor networks.

##  9 Conclusion

We conclude that the two-dimensional surface regularization plays an important
role in few-shot NLOS imaging tasks. The joint regularization of the estimated
signal, the voxel-based representation of the target and the hidden surface
makes it possible to reconstruct the hidden scene with certain complex
geometric structures with only $5\times 5$ confocal measurements, even in
cases with high measurement noise. We would like to incorporate other powerful
two-dimensional regularizations of the surface for higher noise-robustness of
the SSCR framework in the future.

Acknowledgments This work was supported by the National Natural Science
Foundation of China (61975087, 12071244, 11971258).

## References

  * [1] Victor Arellano, Diego Gutierrez, and Adrian Jarabo.  Fast back-projection for non-line of sight reconstruction.  Opt. Express, 25(10):11574–11583, May 2017.
  * [2] Mauro Buttafava, Jessica Zeman, Alberto Tosi, Kevin Eliceiri, and Andreas Velten.  Non-line-of-sight imaging using a time-gated single photon avalanche diode.  Opt. Express, 23(16):20997–21011, Aug 2015.
  * [3] Jian-Feng Cai, Hui Ji, Zuowei Shen, and Gui-Bo Ye.  Data-driven tight frame construction and image denoising.  Applied and Computational Harmonic Analysis, 37(1):89–105, 2014\.
  * [4] Susan Chan, Ryan E. Warburton, Genevieve Gariepy, Jonathan Leach, and Daniele Faccio.  Non-line-of-sight tracking of people at long range.  Opt. Express, 25(9):10109–10117, May 2017.
  * [5] Sreenithy Chandran and Suren Jayasuriya.  Adaptive lighting for data-driven non-line-of-sight 3d localization and object identification.  arXiv preprint arXiv:1905.11595, 2019.
  * [6] Wenzheng Chen, Simon Daneau, Fahim Mannan, and Felix Heide.  Steady-state non-line-of-sight imaging.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [7] Javier Grau Chopite, Matthias B. Hullin, Michael Wand, and Julian Iseringhausen.  Deep non-line-of-sight reconstruction.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2020.
  * [8] Andrzej Cichocki, Namgil Lee, Ivan V Oseledets, A-H Phan, Qibin Zhao, and D Mandic.  Low-rank tensor networks for dimensionality reduction and large-scale optimization problems: Perspectives and challenges part 1.  arXiv:1609.00893, 2016.
  * [9] Kostadin Dabov, Alessandro Foi, Vladimir Katkovnik, and Karen Egiazarian.  Image denoising with block-matching and 3d filtering.  In Nasser M. Nasrabadi, Syed A. Rizvi, Edward R. Dougherty, Jaakko T. Astola, and Karen O. Egiazarian, editors, Image Processing: Algorithms and Systems, Neural Networks, and Machine Learning, volume 6064, page 606414\. International Society for Optics and Photonics, SPIE, 2006.
  * [10] Daniele Faccio, Andreas Velten, and Gordon Wetzstein.  Non-line-of-sight imaging.  Nature Reviews Physics, 2(6):318–327, Jun 2020.
  * [11] Xiaohua Feng and Liang Gao.  Improving non-line-of-sight image reconstruction with weighting factors.  Opt. Lett., 45(14):3921–3924, Jul 2020.
  * [12] Felix Heide, Matthew O’Toole, Kai Zang, David B Lindell, Steven Diamond, and Gordon Wetzstein.  Non-line-of-sight imaging with partial occluders and surface normals.  ACM Transactions on Graphics (ToG), 38(3):1–10, 2019.
  * [13] Mariko Isogawa, Dorian Chan, Ye Yuan, Kris M. Kitani, and Matthew O’Toole.  Efficient non-line-of-sight imaging from transient sinograms.  In 16th European Conference on Computer Vision (ECCV), 2020.
  * [14] Mariko Isogawa, Ye Yuan, Matthew O’Toole, and Kris M. Kitani.  Optical non-line-of-sight physics-based 3d human pose estimation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2020.
  * [15] Masaki Kaga, Takahiro Kushida, Tsuyoshi Takatani, Kenichiro Tanaka, Takuya Funatomi, and Yasuhiro Mukaigawa.  Thermal non-line-of-sight imaging from specular and diffuse reflections.  IPSJ Transactions on Computer Vision and Applications, 11(1):8, Nov 2019.
  * [16] Martin Laurenzis and Andreas Velten.  Feature selection and back-projection algorithms for nonline-of-sight laser-gated viewing.  Journal of Electronic Imaging, 23(6):063003, 2014.
  * [17] Zhupeng Li, Xintong Liu, Jianyu Wang, Zuoqiang Shi, Lingyun Qiu, and Xing Fu.  Fast non-line-of-sight imaging based on first photon event stamping.  Opt. Lett., 47(8):1928–1931, Apr 2022.
  * [18] David B. Lindell, Gordon Wetzstein, and Vladlen Koltun.  Acoustic non-line-of-sight imaging.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [19] David B. Lindell, Gordon Wetzstein, and Matthew O’Toole.  Wave-based non-line-of-sight imaging using fast f-k migration.  ACM Trans. Graph. (SIGGRAPH), 38(4):116, 2019.
  * [20] Jianjiang Liu, Yijun Zhou, Xin Huang, Zheng-Ping Li, and Feihu Xu.  Photon-efficient non-line-of-sight imaging.  IEEE Transactions on Computational Imaging, 8:639–650, 2022.
  * [21] Xiaochun Liu, Sebastian Bauer, and Andreas Velten.  Analysis of feature visibility in non-line-of-sight measurements.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [22] Xiaochun Liu, Sebastian Bauer, and Andreas Velten.  Phasor field diffraction based reconstruction for fast non-line-of-sight imaging systems.  Nature Communications, 11(1):1645, Apr 2020.
  * [23] Xiaochun Liu, Ibón Guillén, Marco La Manna, Ji Hyun Nam, Syed Azer Reza, Toan Huu Le, Adrian Jarabo, Diego Gutierrez, and Andreas Velten.  Non-line-of-sight imaging using phasor-field virtual waveoptics.  Nature, 572(7771):620–623, Aug 2019.
  * [24] Xintong Liu, Jianyu Wang, Zhupeng Li, Zuoqiang Shi, Xing Fu, and Lingyun Qiu.  Non-line-of-sight reconstruction with signal–object collaborative regularization.  Light: Science & Applications, 10(1):198, Sep 2021.
  * [25] Xintong Liu, Jianyu Wang, Leping Xiao, Zuoqiang Shi, Xing Fu, and Lingyun Qiu.  Non-line-of-sight imaging with arbitrary illumination and detection pattern.  arXiv preprint arXiv:2211.00648, 2022.
  * [26] Tomohiro Maeda, Yiqin Wang, Ramesh Raskar, and Achuta Kadambi.  Thermal non-line-of-sight imaging.  In 2019 IEEE International Conference on Computational Photography (ICCP), pages 1–11, 2019.
  * [27] Christopher A. Metzler, Felix Heide, Prasana Rangarajan, Muralidhar Madabhushi Balaji, Aparna Viswanath, Ashok Veeraraghavan, and Richard G. Baraniuk.  Deep-inverse correlography: towards real-time high-resolution non-line-of-sight imaging.  Optica, 7(1):63–71, Jan 2020.
  * [28] Ji Hyun Nam, Eric Brandt, Sebastian Bauer, Xiaochun Liu, Marco Renna, Alberto Tosi, Eftychios Sifakis, and Andreas Velten.  Low-latency time-of-flight non-line-of-sight imaging at 5 frames per second.  Nature Communications, 12(1):6526, Nov 2021.
  * [29] Stanley Osher, Zuoqiang Shi, and Wei Zhu.  Low dimensional manifold model for image processing.  SIAM Journal on Imaging Sciences, 10(4):1669–1690, 2017.
  * [30] Matthew O’Toole, David B. Lindell, and Gordon Wetzstein.  Confocal non-line-of-sight imaging based on the light-cone transform.  Nature, 555(7696):338–341, Mar 2018.
  * [31] A. Rajwade, A. Rangarajan, and A. Banerjee.  Image denoising using the higher order singular value decomposition.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(4), 2012.
  * [32] Syed Azer Reza, Marco La Manna, Sebastian Bauer, and Andreas Velten.  Phasor field waves: A huygens-like light transport model for non-line-of-sight imaging applications.  Opt. Express, 27(20):29380–29400, Sep 2019.
  * [33] Nicolas Scheiner, Florian Kraus, Fangyin Wei, Buu Phan, Fahim Mannan, Nils Appenrodt, Werner Ritter, Jurgen Dickmann, Klaus Dietmayer, Bernhard Sick, and Felix Heide.  Seeing around street corners: Non-line-of-sight detection and tracking in-the-wild using doppler radar.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2020.
  * [34] Sheila W. Seidel, John Murray-Bruce, Yanting Ma, Christopher Yu, William T. Freeman, and Vivek K Goyal.  Two-dimensional non-line-of-sight scene estimation from a single edge occluder.  IEEE Transactions on Computational Imaging, 7:58–72, 2021.
  * [35] Kenichiro Tanaka, Yasuhiro Mukaigawa, and Achuta Kadambi.  Polarized non-line-of-sight imaging.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2020.
  * [36] Chia-Yin Tsai, Aswin C. Sankaranarayanan, and Ioannis Gkioulekas.  Beyond volumetric albedo – a surface optimization framework for non-line-of-sight imaging.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [37] Andreas Velten, Thomas Willwacher, Otkrist Gupta, Ashok Veeraraghavan, Moungi G. Bawendi, and Ramesh Raskar.  Recovering three-dimensional shape around a corner using ultrafast time-of-flight imaging.  Nature Communications, 3(1):745, Mar 2012.
  * [38] Bin Wang, Ming-Yang Zheng, Jin-Jian Han, Xin Huang, Xiu-Ping Xie, Feihu Xu, Qiang Zhang, and Jian-Wei Pan.  Non-line-of-sight imaging with picosecond temporal resolution.  Phys. Rev. Lett., 127:053602, Jul 2021.
  * [39] Cheng Wu, Jianjiang Liu, Xin Huang, Zheng-Ping Li, Chao Yu, Jun-Tian Ye, Jun Zhang, Qiang Zhang, Xiankang Dou, Vivek K. Goyal, Feihu Xu, and Jian-Wei Pan.  Non–line-of-sight imaging over 1.43 km.  Proceedings of the National Academy of Sciences, 118(10), 2021.
  * [40] Shumian Xin, Sotiris Nousias, Kiriakos N. Kutulakos, Aswin C. Sankaranarayanan, Srinivasa G. Narasimhan, and Ioannis Gkioulekas.  A theory of fermat paths for non-line-of-sight shape reconstruction.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [41] Jun-Tian Ye, Xin Huang, Zheng-Ping Li, and Feihu Xu.  Compressed sensing for active non-line-of-sight imaging.  Opt. Express, 29(2):1749–1763, Jan 2021.
  * [42] Sean I. Young, David B. Lindell, Bernd Girod, David Taubman, and Gordon Wetzstein.  Non-line-of-sight surface reconstruction using the directional light-cone transform.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2020.
  * [43] Dayu Zhu and Wenshan Cai.  Fast non-line-of-sight imaging with two-step deep remapping.  ACS Photonics, 9(6):2046–2055, 2022.


[table ]capposition=bottom capbtabboxtable[][]

#  $\Xi_{c}-\Xi_{c}^{\prime}$ mixing From Lattice QCD

Hang Liu  INPAC, Key Laboratory for Particle Astrophysics and Cosmology (MOE),
Shanghai Key Laboratory for Particle Physics and Cosmology, School of Physics
and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China  Liuming
Liu  Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou,
730000, China  University of Chinese Academy of Sciences, Beijing 100049,
China  Peng Sun  Institute of Modern Physics, Chinese Academy of Sciences,
Lanzhou, 730000, China  University of Chinese Academy of Sciences, Beijing
100049, China  Wei Sun  Institute of High Energy Physics, Chinese Academy of
Sciences, Beijing 100049, China  Jin-Xin Tan  INPAC, Key Laboratory for
Particle Astrophysics and Cosmology (MOE), Shanghai Key Laboratory for
Particle Physics and Cosmology, School of Physics and Astronomy, Shanghai Jiao
Tong University, Shanghai 200240, China  Wei Wang  Corresponding author:
<EMAIL_ADDRESS>INPAC, Key Laboratory for Particle Astrophysics and
Cosmology (MOE), Shanghai Key Laboratory for Particle Physics and Cosmology,
School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai
200240, China  Southern Center for Nuclear-Science Theory (SCNT), Institute of
Modern Physics, Chinese Academy of Sciences, Huizhou 516000, Guangdong
Province, China  Yi-Bo Yang  CAS Key Laboratory of Theoretical Physics,
Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190,
China  School of Fundamental Physics and Mathematical Sciences, Hangzhou
Institute for Advanced Study, UCAS, Hangzhou 310024, China  International
Centre for Theoretical Physics Asia-Pacific, Beijing/Hangzhou, China
University of Chinese Academy of Sciences, Beijing 100049, China  Qi-An Zhang
Corresponding author<EMAIL_ADDRESS>School of Physics, Beihang
University, Beijing 102206, China

###### Abstract

In heavy quark limit, the lowest-lying charmed baryons with two light quarks
can form an SU(3) triplet and sextet. The $\Xi_{c}$ in the SU(3) triplet and
$\Xi_{c}^{\prime}$ in the sextet have the same $J^{PC}$ quantum number and can
mix due to the finite charm quark mass and the fact the strange quark is
heavier than the up/down quark. We explore the $\Xi_{c}$-$\Xi_{c}^{\prime}$
mixing by calculating the two-point correlation functions of the $\Xi_{c}$ and
$\Xi_{c}^{\prime}$ baryons from lattice QCD. Based on the lattice data, we
adopt two independent methods to determine the mixing angle between $\Xi_{c}$
and $\Xi_{c}^{\prime}$. After making the chiral and continuum extrapolation,
it is found that the mixing angle $\theta$ is $1.2^{\circ}\pm 0.1^{\circ}$,
which seems insufficient to account for the large SU(3) symmetry breaking
effects found in weak decays of charmed baryons.

##  I  Introduction

Heavy baryons with a bottom or charm quark provide an ideal place to study the
underlying theory for strong interactions. In heavy quark limit, charmed
baryons can be elegantly classified according to heavy quark spin symmetry
Isgur:1989vq ; Isgur:1991wq . Ground states of baryons with one heavy quark
can be categorized into two sets: if the light quark system is spinless, three
baryons form an SU(3) triplet; if the spin of the light quark system is one,
then six charmed baryons form a sextet. However in reality, since the charm
quark mass is finite and the strange quark is heavier than the up/down quark,
charm baryons in the triplet and sextet are likely to mix with each other. An
explicit realization is the $\Xi_{c}-\Xi_{c}^{\prime}$ mixing, which is a
direct consequence of the effective heavy quark and flavor SU(3) symmetry
breaking.

Weak decays of charmed baryons also contribute to the understanding of the
electroweak sector of the standard model. In particular semileptonic decays of
charmed baryons are valuable to extract the CKM matrix element $|V_{cs}|$ (for
some recent experimental measurements and lattice QCD calculations , please
see Refs. BESIII:2015ysy ; BESIII:2016ffj ; Meinel:2016dqj ; Belle:2021crz ;
Zhang:2021oja  and many references therein). It is interesting to notice that
based on the experimental measurements the flavor SU(3) symmetry, which is a
very useful tool widely applied in the analysis of weak decays of heavy
mesons, is found to be sizably broken He:2021qnc . Various mechanisms to
explain this observation were explored in Ref. He:2021qnc  and a competitive
explanation is through the $\Xi_{c}-\Xi_{c}^{\prime}$ mixing Geng:2022yxb .
Very interestingly, several recent works are devoted to determine the mixing
angle from various methods, but the results are controversial. Based on
experimental data on weak decays of $\Xi_{c}$ and $\Xi_{cc}$, it is found that
a large mixing angle is called for: $\theta=24.66^{\circ}\pm 0.90^{\circ}$ in
Ref.Geng:2022yxb ; Geng:2022xfz  and $16.27^{\circ}\pm 2.30^{\circ}$ in
Ref.Ke:2022gxm . A direct calculation in QCD sum rule gives
$\theta=5.5^{\circ}\pm 1.8^{\circ}$ Aliev:2010ra , while in heavy-quark
effective theory, it is found $\theta=8.12^{\circ}\pm 0.80^{\circ}$
Matsui:2020wcc . Besides, an earlier lattice QCD exploration indicates a small
mixing angle Brown:2014ena .

This paper presents an update-to-date exploration of the
$\Xi_{c}-\Xi_{c}^{\prime}$ mixing from lattice QCD. Using the recently
generated lattice configurations, we calculate the matrix elements $\langle
0|O_{\Xi_{c}^{(\bar{3},6)}}(t)\bar{O}_{\Xi_{c}^{(\bar{3},6)}}(0)|0\rangle$,
where $O_{\Xi_{c}^{\bar{3}}}$ and $O_{\Xi_{c}^{6}}$ denote the interpolators
for the anti-triplet and sextet baryon respectively. With the results for the
correlation function matrix, we use the direct fitting method to fit the
energy levels and the mixing angle. As a comparison, the lattice spectroscopy
method will also be used. After extracting the mixing angles with different
lattice spacings and pion masses, we extrapolate the result to the chiral and
continuum limit. We finally find that the mixing angle is at the order of
$1^{\circ}$, which is likely to indicate the necessity to include other SU(3)
symmetry-breaking effects in $\Xi_{c}$ decays.

The rest of this paper is organized as follows. The theoretical framework for
the $\Xi_{c}-\Xi_{c}^{\prime}$ mixing is collected in Sec. II. We present our
lattice simulations of the two-point correlation functions in Sec.III and a
determination of the mixing angle from a direct fit in Sec. IV. Sec. V gives
an analysis of the mixing angle through solving the generalized eigenvalue
problem. The conclusion of this work is given in Sec. VI.

##  II $\Xi_{c}-\Xi_{c}^{\prime}$ Mixing in Heavy Quark Effective Theory

In quark model, a baryon is made of three quarks. In heavy quark limit, heavy
baryons with one charm quark can be well classified according to the angular
momentum $J$ of light quark system. If the spin wave function of the light
quark system is antisymmetric, i.e. $J_{qq}=0$, the Fermi statistics and
antisymmetric feature in color space require that the flavor wave function is
also antisymmetric in $SU(3)_{F}$. This corresponds the decomposition
$(3\times 3)_{\mathrm{antisymmetric}}=\bar{3}$. If $J_{qq}=1$, the $SU(3)_{F}$
flavor wave function should be symmetric and hence transforms as $(3\times
3)_{\mathrm{symmetric}}=6$.

In heavy quark effective theory (HQET) Isgur:1989vq ; Isgur:1991wq , the
generic interpolating operator of a heavy baryon takes the form
$O=\epsilon^{abc}\left(q_{1}^{Ta}C\Gamma
q_{2}^{b}\right)\Gamma^{\prime}\tilde{Q}^{c}$, where $q_{1,2}$ denotes a light
quark, and $\tilde{Q}$ is the heavy quark field in HQET satisfying
$\gamma^{0}\tilde{Q}=\tilde{Q}$. In this interpolator we have explicitly
specified the color indices $a,b,c$. The transposition $T$ acts on a Dirac
spinor, and $C=\gamma^{0}\gamma^{2}$ is the charge conjugation matrix. Summing
over the color indices with the totally antisymmetric tensor $\epsilon^{abc}$
results in a gauge invariant operator. The Dirac matrix $\Gamma$ and
$\Gamma^{\prime}$ are related to the internal spin structures of the heavy
baryon. For the ground state, the angular momentum of light quark system is
either 0 or 1. If $J=0$ corresponds to the current
$q_{1}^{T}C\gamma_{5}q_{2}$, and $J=1$ corresponds to
$q_{1}^{T}C\vec{\gamma}q_{2}$. Therefore, the baryonic current
$O_{1}=\left(q_{1}^{T}C\gamma_{5}q_{2}\right)\tilde{Q}$ corresponds to the
spin $1/2$ baryon, and the current
$\vec{O}_{2}=\left(q_{1}^{T}C\vec{\gamma}q_{2}\right)\tilde{Q}$ contains both
spin-$1/2$ and spin-$3/2$ components. Using $\vec{O}_{2}$, one can construct
the spin-$3/2$ component
$\vec{O}_{2}^{3/2}=\vec{O}_{2}+\frac{1}{3}\vec{\gamma}\vec{\gamma}\cdot\vec{O}_{2}$
which satisfying the condition $\vec{\gamma}\cdot\vec{O}_{2}^{3/2}=0$, and the
spin-$1/2$ component
${O}_{2}^{1/2}=-\frac{1}{3}\vec{\gamma}\vec{\gamma}\cdot\vec{O}_{2}=\left(q_{1}^{T}C\vec{\gamma}q_{2}\right)\cdot\vec{\gamma}\gamma_{5}\tilde{Q}$.

As a result, one can obtain the baryonic currents of $SU(3)_{F}$ eigenstate as
Grozin:1992td

$\displaystyle
O^{\bar{3}}=\epsilon^{abc}\left(\ell^{Ta}C\gamma_{5}s^{b}\right)P_{+}\tilde{c}^{c},$
(1) $\displaystyle
O^{6}=\epsilon^{abc}\left(\ell^{Ta}C\vec{\gamma}s^{b}\right)\cdot\vec{\gamma}\gamma_{5}P_{+}\tilde{c}^{c},$
(2)

with $\ell=(u,d)$. $P^{+}=(1+\gamma^{0})/2$ is the positive parity projector.

If the finite charm quark mass and the differences between the strange quark
and up/down quark are taken into account, the two $\Xi_{c}$ states can mix and
their mixing effect could be described by a $2\times 2$ mixing matrix:

$\displaystyle\left(\begin{array}[]{c}\left|\Xi_{c}\right\rangle\\\
\left|\Xi_{c}^{\prime}\right\rangle\end{array}\right)=\left(\begin{array}[]{cc}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{array}\right)\left(\begin{array}[]{c}|\Xi_{c}^{\bar{3}}\rangle\\\
|\Xi_{c}^{6}\rangle\end{array}\right),$ (9)

where $\theta$ denotes the mixing angle. The mass eigenstates are orthogonal:

$\displaystyle
H_{\mathrm{QCD}}\left|\Xi_{c}\right\rangle=m_{\Xi_{c}}\left|\Xi_{c}\right\rangle,\quad
H_{\mathrm{QCD}}\left|\Xi_{c}^{\prime}\right\rangle=m_{\Xi_{c}^{\prime}}\left|\Xi_{c}^{\prime}\right\rangle,$
(10)

where $m_{\Xi_{c}}/m_{\Xi_{c}^{\prime}}$ denote the physical baryon masses.

##  III Two-point correlation function From Lattice QCD

Ensemble | $\beta$ | $L^{3}\times T$ |  $a$ (fm) | $m_{l}^{\mathrm{b}}$ | $m_{s}^{\mathrm{b}}$ | $m_{c}^{\mathrm{b}}$ | $m_{\pi}$ | $N_{\mathrm{meas}}$
---|---|---|---|---|---|---|---|---
C11P14L |  | $48^{3}\times 96$ |  | $-0.2825$ | $-0.2310$ | $0.4800$ | $135$ | $203\times 48$
C11P22M | 6.20 | $32^{3}\times 64$ | $0.108$ | $-0.2790$ | $-0.2310$ | $0.4800$ | $222$ | $451\times 20$
C11P29S |  | $24^{3}\times 72$ |  | $-0.2770$ | $-0.2315$ | $0.4780$ | $284$ | $432\times 26$
C08P30S | 6.41 | $32^{3}\times 96$ | $0.080$ | $-0.2295$ | $-0.2010$ | $0.2326$ | $297$ | $653\times 26$
C06P30S | 6.72 | $48^{3}\times 144$ | $0.055$ | $-0.1850$ | $-0.1687$ | $0.0770$ | $312$ | $136\times 80$
Table 1: Parameters of the ensembles used in this work, including the gauge
coupling $\beta=10/g^{2}$, the 4-dimensional volume $L^{3}\times T$, lattice
spacing $a$, bare quark masses $m_{l,s,c}^{\mathrm{b}}$, pion mass $m_{\pi}$
and total measurements $N_{\mathrm{meas}}$. The total measurements are equal
to the number of gauge configurations times the measurements from different
time slices on one configuration.

The simulations in this work will be performed on the gauge configurations
generated by the Chinese Lattice QCD (CLQCD) collaboration with $N_{f}=2+1$
flavor stout smeared clover fermions and Symanzik gauge action, at three
lattice spacings and three pion masses including the physical one. The
detailed parameters are listed in Table 1. Some previous applications of these
configurations can be found in Refs. Zhang:2021oja ; Wang:2021vqy ;
Liu:2022gxf ; Xing:2022ijm .

In the numerical simulation, to improve the signal-to-noise ratio of the
simulation, we generate the $f$-flavored wall source to point sink propagators

$\displaystyle
S_{{w-p}}^{f}(\vec{x},t,t_{0};\vec{p})=\sum_{\vec{x}_{s}}e^{-i\vec{p}\cdot(\vec{x}-\vec{x}_{s})}S^{f}(\vec{x},t;\vec{x}_{s},t_{0})$
(11)

on the Coulomb gauge fixed configurations. $S^{f}$ denotes the propagator from
wall source at $(\sum_{\vec{x}_{s}}\vec{x}_{s},t_{0})$ to the point sink
$(\vec{x},t)$. We consider the $2\times 2$ correlation function matrix

$\displaystyle\mathcal{C}(t,t_{0})=\sum_{\vec{x}}\left(\begin{array}[]{cc}\left\langle
O^{\bar{3}}_{p}(\vec{x},t)\bar{O}^{\bar{3}}_{w}(\vec{0},t_{0})\right\rangle&\left\langle
O^{\bar{3}}_{p}(\vec{x},t)\bar{O}^{6}_{w}(\vec{0},t_{0})\right\rangle\\\
\left\langle
O^{6}_{p}(\vec{x},t)\bar{O}^{\bar{3}}_{w}(\vec{0},t_{0})\right\rangle&\left\langle
O^{6}_{p}(\vec{x},t)\bar{O}^{6}_{w}(\vec{0},t_{0})\right\rangle\end{array}\right)$
(14)

that contains the baryonic currents. The subscript “$w$” and “$p$” denote the
wall source and point sink. The related two-point functions can be constructed
by the light, strange, and charm quark propagators $S_{{w-p}}^{\\{l,s,c\\}}$
with momentum $\vec{p}=\vec{0}$

$\displaystyle C_{11}(t,t_{0})=$
$\displaystyle-\sum_{\vec{x}}\epsilon^{abc}\epsilon^{a^{\prime}b^{\prime}c^{\prime}}\Big{(}S_{{w-p}}^{l}(\vec{x},t,t_{0})\Big{)}_{\alpha\alpha^{\prime}}^{aa^{\prime}}\Big{(}C\gamma_{5}S_{{w-p}}^{s}(\vec{x},t,t_{0})\gamma_{5}C\Big{)}_{\alpha\alpha^{\prime}}^{bb^{\prime}}\Big{(}P_{+}S_{{w-p}}^{c}(\vec{x},t,t_{0})P_{+}T\Big{)}^{cc^{\prime}},$
(15) $\displaystyle C_{12}(t,t_{0})=$
$\displaystyle\sum_{\vec{x}}\epsilon^{abc}\epsilon^{a^{\prime}b^{\prime}c^{\prime}}\sum_{i=1,2,3}\Big{(}S_{{w-p}}^{l}(\vec{x},t,t_{0})\Big{)}_{\alpha\alpha^{\prime}}^{aa^{\prime}}\Big{(}C\gamma_{5}S_{{w-p}}^{s}(\vec{x},t,t_{0}){\gamma}^{i}C\Big{)}_{\alpha\alpha^{\prime}}^{bb^{\prime}}\Big{(}P_{+}S_{{w-p}}^{c}(\vec{x},t,t_{0}){\gamma}^{i}\gamma_{5}P_{+}T\Big{)}^{cc^{\prime}},$
(16) $\displaystyle C_{21}(t,t_{0})=$
$\displaystyle-\sum_{\vec{x}}\epsilon^{abc}\epsilon^{a^{\prime}b^{\prime}c^{\prime}}\sum_{i=1,2,3}\Big{(}S_{{w-p}}^{l}(\vec{x},t,t_{0})\Big{)}_{\alpha\alpha^{\prime}}^{aa^{\prime}}\Big{(}C\gamma^{i}S_{{w-p}}^{s}(\vec{x},t,t_{0})\gamma_{5}C\Big{)}_{\alpha\alpha^{\prime}}^{bb^{\prime}}\Big{(}P_{+}\gamma^{i}\gamma_{5}S_{{w-p}}^{c}(\vec{x},t,t_{0})P_{+}T\Big{)}^{cc^{\prime}},$
(17) $\displaystyle C_{22}(t,t_{0})=$
$\displaystyle\sum_{\vec{x}}\epsilon^{abc}\epsilon^{a^{\prime}b^{\prime}c^{\prime}}\sum_{i.j=1,2,3}\Big{(}S_{{w-p}}^{l}(\vec{x},t,t_{0})\Big{)}_{\alpha\alpha^{\prime}}^{aa^{\prime}}\Big{(}C\gamma^{i}S_{{w-p}}^{s}(\vec{x},t,t_{0})\gamma^{j}C\Big{)}_{\alpha\alpha^{\prime}}^{bb^{\prime}}\Big{(}P_{+}\gamma^{i}\gamma_{5}S_{{w-p}}^{c}(\vec{x},t,t_{0})\gamma^{j}\gamma_{5}P_{+}T\Big{)}^{cc^{\prime}}$
(18)

where the subscript $\alpha,\alpha^{\prime}$ denote the Dirac indices and
superscript $a,b,c,a^{\prime},b^{\prime},c^{\prime}$ denote the color indices.
In the above a trace has been taken and we have adopted a unpolarized
projector $T=(1+\gamma^{0})/2$. One can generate the wall source propagators
at several time slices $t_{0}$ and then average the two-point functions
constructed by these propagators, so that the statistical precisions of
numerical results on different ensembles can reach the same level. The
explicit statistics on each ensemble, including the number of gauge
configurations and the measurements from different time slices on one
configuration, are collected in Table 1.

One can extract the mixing parameters from a joint analysis of both diagonal
and non-diagonal terms of $\mathcal{C}$. In the correlation functions, one can
insert the mass eigenstates $|n\rangle$ and convert them to local matrix
elements. The mismatch of flavor basis baryonic currents and mass eigenstates
will be related to the mixing angle under the mixing pattern in Eq.(9). For
example, the correlation $C_{11}$ is given as:

$\displaystyle C_{11}(t,t_{0})=$ $\displaystyle\sum_{\vec{x}}\langle
O^{\bar{3}}_{p}(\vec{x},t)\bar{O}^{\bar{3}}_{w}(\vec{0},t_{0})\rangle$
$\displaystyle=$
$\displaystyle\sum_{n=\Xi_{c},\Xi_{c}^{\prime},...}\frac{e^{-m_{n}(t-t_{0})}}{(La)^{3}(2m_{n})}\langle
0|O^{\bar{3}}_{p}(\vec{0},0)|n\rangle\langle
n|\bar{O}^{\bar{3}}_{w}(\vec{0},0)|0\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\left(La\right)^{3}}\left[\frac{e^{-m_{\Xi_{c}}(t-t_{0})}}{2m_{\Xi_{c}}}\langle
0|O^{\bar{3}}_{p}(\vec{0},0)|\Xi_{c}\rangle\langle\Xi_{c}|\bar{O}^{\bar{3}}_{w}(\vec{0},0)|0\rangle+\frac{e^{-m_{\Xi_{c}^{\prime}}(t-t_{0})}}{2m_{\Xi_{c}^{\prime}}}\langle
0|O^{\bar{3}}_{p}(\vec{0},0)|\Xi_{c}^{\prime}\rangle\langle\Xi_{c}^{\prime}|\bar{O}^{\bar{3}}_{w}(\vec{0},0)|0\rangle+\cdots\right].$
(19)

From this expression one can in general determine the mass of the ground state
with a large time slice $t$. The excited state contributions are greatly
suppressed in this regime, and the correlation functions become

$\displaystyle C_{11}(t,t_{0})=$
$\displaystyle\frac{1}{\left(La\right)^{3}}\left[\frac{e^{-m_{\Xi_{c}}(t-t_{0})}}{2m_{\Xi_{c}}}\langle
0|O^{\bar{3}}_{p}(\vec{0},0)|\Xi_{c}\rangle\langle\Xi_{c}|\bar{O}^{\bar{3}}_{w}(\vec{0},0)|0\rangle+\frac{e^{-m_{\Xi_{c}^{\prime}}(t-t_{0})}}{2m_{\Xi_{c}^{\prime}}}\langle
0|O^{\bar{3}}_{p}(\vec{0},0)|\Xi_{c}^{\prime}\rangle\langle\Xi_{c}^{\prime}|\bar{O}^{\bar{3}}_{w}(\vec{0},0)|0\rangle\right]$
$\displaystyle=$
$\displaystyle\frac{1}{\left(La\right)^{3}}\left[\frac{e^{-m_{\Xi_{c}}(t-t_{0})}}{2m_{\Xi_{c}}}\left\langle
0\left|O^{\bar{3}}_{p}(\vec{0},0)\right.\right.\left(\left|\Xi_{c}^{\bar{3}}\right\rangle\cos\theta+\left|\Xi_{c}^{6}\right\rangle\sin\theta\right)\left(\left\langle\Xi_{c}^{\bar{3}}\right|\cos\theta+\left\langle\Xi_{c}^{6}\right|\sin\theta\right)\left|\bar{O}^{\bar{3}}_{w}(\vec{0},0)\right|0\right\rangle$
$\displaystyle\quad\left.+\frac{e^{-m_{\Xi_{c}^{\prime}}(t-t_{0})}}{2m_{\Xi_{c}^{\prime}}}\left.\left\langle
0\left|O^{\bar{3}}_{p}(\vec{0},0)\right.\right.\left(-\left|\Xi_{c}^{\bar{3}}\right\rangle\sin\theta+\left|\Xi_{c}^{6}\right\rangle\cos\theta\right)\left(-\left\langle\Xi_{c}^{\bar{3}}\right|\sin\theta+\left\langle\Xi_{c}^{6}\right|\cos\theta\right)\left|\bar{O}^{\bar{3}}_{w}(\vec{0},0)\right|0\right\rangle\right]$
$\displaystyle\equiv$ $\displaystyle
A_{p}A_{w}^{\dagger}\left[\frac{\cos^{2}{\theta}}{2m_{\Xi_{c}}}e^{-m_{\Xi_{c}}(t-t_{0})}+\frac{\sin^{2}{\theta}}{2m_{\Xi_{c}^{\prime}}}e^{-m_{\Xi_{c}^{\prime}}(t-t_{0})}\right],$
(20)

where the $t$-independent parameter $A_{p/w}=(La)^{-{3}/{2}}\langle
0|O^{\bar{3}}_{p/w}(\vec{0},0)|\Xi_{c}^{\bar{3}}\rangle$. In the second step,
we have used the mixing matrix in Eq. (9). Similarly, the other matrix
elements can be expressed as

$\displaystyle C_{12}(t,t_{0})=$ $\displaystyle\sum_{\vec{x}}\langle
O^{\bar{3}}_{p}(\vec{x},t)\bar{O}^{6}_{w}(\vec{0},t_{0})\rangle=A_{p}B_{w}^{\dagger}\left[\frac{\cos{\theta}\sin{\theta}}{2m_{\Xi_{c}}}e^{-m_{\Xi_{c}}(t-t_{0})}-\frac{\cos{\theta}\sin{\theta}}{2m_{\Xi_{c}^{\prime}}}e^{-m_{\Xi_{c}^{\prime}}(t-t_{0})}\right],$
(21) $\displaystyle C_{21}(t,t_{0})=$ $\displaystyle\sum_{\vec{x}}\langle
O^{6}_{p}(\vec{x},t)\bar{O}^{\bar{3}}_{w}(\vec{0},t_{0})\rangle=B_{p}A_{w}^{\dagger}\left[\frac{\cos{\theta}\sin{\theta}}{2m_{\Xi_{c}}}e^{-m_{\Xi_{c}}(t-t_{0})}-\frac{\cos{\theta}\sin{\theta}}{2m_{\Xi_{c}^{\prime}}}e^{-m_{\Xi_{c}^{\prime}}(t-t_{0})}\right],$
(22) $\displaystyle C_{22}(t,t_{0})=$ $\displaystyle\sum_{\vec{x}}\langle
O^{{6}}_{p}(\vec{x},t)\bar{O}^{{6}}_{w}(\vec{0},t_{0})\rangle=B_{p}B_{w}^{\dagger}\left[\frac{\sin^{2}{\theta}}{2m_{\Xi_{c}}}e^{-m_{\Xi_{c}}(t-t_{0})}+\frac{\cos^{2}{\theta}}{2m_{\Xi_{c}^{\prime}}}e^{-m_{\Xi_{c}^{\prime}}(t-t_{0})}\right],$
(23)

with $B_{p/w}={(La)^{-{3}/{2}}}\langle
0|O^{{6}}_{p/w}(\vec{0},0)|\Xi_{c}^{6}\rangle$. It should be noticed that the
lattice artifacts generated by the wall source and point sink have been
collected into the local matrix elements $A_{p/w}$ and $B_{p/w}$, and would
not contaminate the energy eigenstates and their mixing. Illustrated as Fig.1,
the ratios $C_{12}/C_{21}$ are consistent with 1 for all the ensembles we used
in this calculation, regardless of the pion masses and lattice spacings. Such
an observation indicates that the color-spacial factor introduced by the wall
source will contribute to an overall factor, which is independent of the
mixing in the flavor space.

Figure 1: Ratios of $C_{12}$ and $C_{21}$ as function of $t$ on different
ensembles. Their values are consistent with 1 indicate the lattice artifacts
generated by wall source are independent of the mixing in the flavor space.

##  IV Mixing Angle from Joint fit

|  $m_{\Xi_{c}}$ (GeV) |  $m_{\Xi_{c}^{\prime}}$ (GeV) |  $\theta$ (∘) |  $\chi^{2}/$d.o.f | fit range (fm)
---|---|---|---|---|---
C11P14L | $2.4256(19)$ | $2.5196(22)$ | $1.083(30)$ | $0.96$ | $1.19-2.81$
C11P22M | $2.4380(27)$ | $2.5351(30)$ | $0.988(49)$ | $1.0$ | $1.19-2.92$
C11P29S | $2.4587(27)$ | $2.5536(29)$ | $1.002(50)$ | $1.1$ | $1.19-3.24$
C08P30S | $2.4753(21)$ | $2.5809(26)$ | $1.080(42)$ | $0.95$ | $1.20-2.40$
C06P30S | $2.4695(37)$ | $2.5815(48)$ | $1.021(67)$ | $1.2$ | $1.32-2.40$
Extrapolated | $2.4380(68)_{\mathrm{stat}}(403)_{\mathrm{syst}}$ | $2.5562(74)_{\mathrm{stat}}(422)_{\mathrm{syst}}$ | $1.20(9)_{\mathrm{stat}}(2)_{\mathrm{syst}}$ |  |
Exp. data ParticleDataGroup:2022pth  | $2.46794_{-0.00020}^{+0.00017}$ | $2.5784\pm 0.0005$ | — |  |
Table 2:  Results of masses and mixing angles from correlated matrix fits on
different ensembles, together with $\chi^{2}/$d.o.f and fit ranges
$t_{\mathrm{min}}-t_{\mathrm{max}}$ for each fits. The following are
extrapolated results, where the experimental data are also listed for
comparison.  Figure 2: Chiral and continuum extrapolations of $\theta$ based
on the fit results of Eq.(24) as a function of $m_{\pi}$. The extrapolated
results from chiral fits at fixed lattice spacings are shown as the red,
green, and blue dashed lines with error bands. The black line corresponds to
the chiral and continuum extrapolation. The black point denotes the
extrapolated values at the physical pion mass and continuum limit.
Uncertainties in all data points are statistical only.

To extract the masses and mixing angle from the correlation functions, one can
apply a joint fit based on the parametrization formula of the correlation
function matrix in Eq.(20-23) with parameters $A$, $B$, $m_{\Xi_{c}}$,
$m_{\Xi_{c}^{\prime}}$ and $\theta$ at $t_{0}=0$.

The choice of fit $t$-range relies on both fit quality and physical
consideration: the starting time slices $t_{\mathrm{min}}$ after which the
data points are included in the fit must be chosen such that contributions
from higher excited states have been suppressed sufficiently with the smallest
possible statistical uncertainties. The $t_{\mathrm{max}}$ are determined by
avoiding the contaminations from backward-propagating in the periodic boundary
condition and avoiding too many degrees of freedom in the $\chi^{2}$ fit
leading to poorly estimating. The optimal choices of the fit range are
different for the different ensembles, so we determine them independently to
obtain a reasonable description of the joint matrix fit.

Table 2 includes the fitted values of $m_{\Xi_{c}}$, $m_{\Xi_{c}^{\prime}}$,
$\theta$ and their statistical uncertainties, together with the correlated,
reduced $\chi^{2}$-values ($\chi^{2}/$d.o.f) and fit ranges
$t_{\mathrm{min}}-t_{\mathrm{max}}$ on each ensemble. Of all fits, we evaluate
the fit qualities by $\chi^{2}/\mathrm{d.o.f}\lesssim 1$ and the smallest
uncertainties for the parameters. As an example, Fig.3 shows the correlated
joint fit results on C11P14L, which is perfectly consistent with the data of
correlation functions.

Figure 3: Joint fit of the $2\times 2$ correlation function matrix using
Eq.(20-23). The data shown here are from C11P14L, the colored bands indicate
the time range and fit results of each matrix element.

After extracting the mixing angles from the ensembles at different pion masses
and lattice spacings, we can extrapolate these results to the physical values
of the quark masses and the continuum limit. We perform the chiral and
continuum extrapolation through the following ansatz:

$\displaystyle\theta(m_{\pi},a)$
$\displaystyle=\theta_{\mathrm{phy}}+c_{1}\left(m_{\pi}^{2}-m_{\pi,\mathrm{phy}}^{2}\right)+c_{2}a^{2},$
$\displaystyle m_{n}(m_{\pi},a)$
$\displaystyle=m_{n,\mathrm{phy}}+c_{1}\left(m_{\pi}^{2}-m_{\pi,\mathrm{phy}}^{2}\right)+c_{2}a^{2}.$
(24)

These extrapolations are performed at the next-to-leading-order in the chiral
expansion. We also included a quadratic dependence of lattice spacings. In
order to estimate the systematic uncertainties associated with this
truncation, we add the higher-order analytic terms to the fit functions, such
as the quartic term $m_{\pi}^{4}$ in chiral expansion, or $a^{3}$ term in
continuum extrapolation. The latter one, which may arise from heavy-quark
discretization errors, gives sizable systematic error. In practice, we
introduce these terms into the fit functions Eq.(24) separately and calculate
the extrapolated masses and mixing angle from the new higher-order fits. The
extrapolated results with both statistic and systematic uncertainties from the
higher-order analytic terms are collected in Table 2.

With the extrapolation uncertainties taken into account, we find the masses
for the $\Xi_{c}$ and $\Xi_{c}^{\prime}$ are both consistent with the
experimental data, while the mixing angle is:

$\displaystyle\theta=(1.200\pm 0.090\pm 0.020)^{\circ}.$ (25)

##  V Mixing angle from generalized eigenvalue problem

Figure 4: Upper panel: effective masses from eigenvalues $\lambda_{n}(t)$
($n=\Xi_{c},\Xi_{c}^{\prime}$), and model averaging results on C11P14L. The
eigenvalues are from solving the GEVP Eq.26 for the correlation function
matrix at $t_{r}=5a$. Lower panel: the model weight factors (solid lines) and
standard $p$-values (dashed lines) which reflect the fit quality.

Besides the correlated fit, the mixing angle can be also extracted by solving
the generalized eigenvalue problem (GEVP) Michael:1985ne ; Luscher:1990ck ;
Luscher:1990ux ; Dudek:2010wm ; Blossier:2009kd

$\displaystyle\mathcal{C}(t)v_{n}(t)=\lambda_{n}(t)\mathcal{C}(t_{r})v_{n}(t),$
(26)

where $n$ labels the states of $\Xi_{c},\Xi_{c}^{\prime}$, $\lambda_{n}(t)$ is
the eigenvalue and follows the boundary condition $\lambda_{n}(t_{r})=1$, and
there is an orthogonality condition for the eigenvectors of different states
$(n,n^{\prime})$ as
$v_{n^{\prime}}^{\dagger}\mathcal{C}(t_{r})v_{n}=\delta_{nn^{\prime}}$. This
orthogonality condition allows us to extract the spectrum of near degenerate
states $\Xi_{c}$ and $\Xi_{c}^{\prime}$. The eigenvalues $\lambda_{n}$ and
eigenvectors $v_{n}$ usually be solved independently on each time slice $t$
from the source at $t_{0}$, while for the nearby masses, their eigenvalues
might fluctuate time slice by time slice. So we choose a reference time slice
$t_{r}$ on which reference eigenvector as $v_{n,\mathrm{ref}}$, and compare
eigenvectors on other time slices by finding the maximum value of
$v_{n^{\prime},\mathrm{ref}}^{\dagger}\mathcal{C}(t_{r})v_{n}$ which
associates a state $n$ with a reference state $n^{\prime}$.

As discussed above, the two-point correlation function can be decomposed into
a time-independent factor and $e^{-m_{n}t}$ at large times. In the first
factor, the matrix element $\langle 0|O^{i}|n\rangle$ is related to the
eigenvectors which contain the mixing effects between the flavor and energy
eigenstates, and the effective masses related to energy eigenstates can be
extracted from the exponential behavior of the correlators $\lambda_{n}$. The
fit function of masses can be expressed as

$\displaystyle\lambda_{n}(t)=c_{0}e^{-m_{n}(t-t_{r})}\left(1+c_{1}e^{-\Delta
E(t-t_{r})}\right),$ (27)

where the fit parameters $c_{0}$ collect the contributions from time-
independent factors and $m_{n}$ denotes the effective mass of
$n=(\Xi_{c},\Xi_{c}^{\prime})$, $c_{1}$ and $\Delta E$ describe the higher
excited states contributions which can be neglected at large time.

Figure 5:  Reference time $t_{r}$-dependence of $\lambda_{\Xi_{c}}$ on
C11P14L, with $t_{r}=\\{5,6,7,8\\}a$. It indicates that the selection of
$t_{r}$ hardly affects the extraction of effective masses.

The extractions of effective masses are illustrated in the upper panel of
Fig.4. The effective masses
$m_{n}=\log\left(\lambda_{n}(t)/\lambda_{n}(t+1)\right)$ are calculated from
the eigenvalues of GEVP at reference time slice $t_{r}=5a$. We also
investigate the $t_{r}$ dependence of GEVP method, shown as Fig.5, we choose
several reference time slices to extract and fit the effective masses and find
that the selection of $t_{r}$ hardly affects the results of effective masses.
The fit results of masses are collected in Table 3.

In addition, the final results also depend on the fit range. Instead of taking
the full difference between the variations of fit ranges as a systematic error
as a conservative estimate, it has been advocated that a technique of Bayesian
model averaging Jay:2020jkz  can weight the selections of different fit
ranges. This approach allows for a fully rigorous estimation of probability
distributions for parameters of interest by combining results from several
fits. We adopt this method to investigate the dependence on the fit range in
particular the $t_{\mathrm{min}}$. Shown as the lower panel of Fig.4, the
weight factors $\mathrm{Pr}$ and standard $p$-values describe the fit quality
with range start from $t_{\mathrm{min}}$, and the probability-weighted average
of fit results from different $t_{\mathrm{min}}$ will contribute to the final
result, which corresponds to the colored band in the upper panel. More details
are shwon in the Appendix.

|  $m_{\Xi_{c}}$ (GeV) |  $m_{\Xi_{c}^{\prime}}$ (GeV) |  $\theta$ (∘)
---|---|---|---
C11P14L | $2.4274(33)$ | $2.5201(20)$ | $1.095(88)$
C11P22M | $2.4322(39)$ | $2.5321(44)$ | $1.075(76)$
C11P29S | $2.4622(59)$ | $2.5603(79)$ | $1.094(60)$
C08P30S | $2.4763(11)$ | $2.5868(13)$ | $1.225(52)$
C06P30S | $2.4682(38)$ | $2.5905(39)$ | $1.212(55)$
Extrapolated | $2.428(11)_{\mathrm{fit}}(44)_{\mathrm{ext.}}$ | $2.547(12)_{\mathrm{fit}}(36)_{\mathrm{ext.}}$ | $1.22(13)_{\mathrm{fit}}(1)_{\mathrm{ext.}}$
Exp. data ParticleDataGroup:2022pth  | $2.46794_{-0.00020}^{+0.00017}$ | $2.5784\pm 0.0005$ | —
Table 3: Results of masses and mixing angles from fitting GEVP eigenvalues and
eigenvectors with the model averaging approach Jay:2020jkz . The errors from
fit contains both static errors and systematic ones associated with choice of
fit range, and the second errors denote systematic ones from extrapolation.
The extrapolated results and the PDG-averaged ones are also listed here for
comparison.  Figure 6: Chiral and continuum extrapolations of $\theta$
extracted from GEVP method. The extrapolated results from chiral fits at fixed
lattice spacings are shown as the red, green, and blue dashed lines with error
bands. The black line corresponds to the chiral and continuum extrapolation.
The black point denoted the extrapolated values at the physical pion mass and
continuum limit. Errors in all data points are statistical only.

After diagonalizing the correlation function matrix $\mathcal{C}$, we find the
GEVP eigenvalues $\lambda_{n}$ are degenerate to the energy basis, thereby the
mixing effects will be collected into the GEVP eigenvectors $v_{n}$. From the
parametrization form in Eq.(20-23), the normalized eigenvectors can be
expressed as

$\displaystyle v_{1}$
$\displaystyle=\sqrt{1+\frac{A_{p}^{2}\cot^{2}\theta}{B_{p}^{2}}}\left(\begin{array}[]{c}\frac{A_{p}}{B_{p}}\cot\theta\\\
1\end{array}\right),$ (30) $\displaystyle v_{2}$
$\displaystyle=\sqrt{1+\frac{A_{p}^{2}\tan^{2}\theta}{B_{p}^{2}}}\left(\begin{array}[]{c}-\frac{A_{p}}{B_{p}}\tan\theta\\\
1\end{array}\right),$ (33)

Therefore the mixing angle $\theta$ can be extracted from $v_{1,2}$.

The fit results on different ensembles are shown in Tabel 3. As discussed
above, we can also extrapolate the mixing angles and masses from different
ensembles to their physical values, illustrated as Fig.6. The comparison of
the extrapolated results as well as the PDG averaged ones are also listed in
Tabel 3. One can see that the extrapolated results agree with the ones from
the correlated matrix fit. The estimation of systematic uncertainties is
similar to the above discussions.

With the extrapolation uncertainties taken into account, we find the mixing
angle is:

$\displaystyle\theta=(1.22\pm 0.13\pm 0.01)^{\circ},$ (34)

which is consistent with Eq. (25) within the uncertainty. This value is much
smaller than the model result extracted from experimental data on weak decays
Geng:2022yxb ; Geng:2022xfz ; Ke:2022gxm ; Xing:2022phq , and thus
insufficient to explain the large SU(3) symmetry breaking effects found in
charmed baryon decays. Other mechanisms such as QED corrections should be
included.

##  VI Charm quark mass dependence

In HQET, the mixing between $\Xi_{c}$ and $\Xi_{c}^{\prime}$ would vanish in
the heavy-quark limit. In this limit, the angular momentum $J$ of the light
quark pairs becomes a conserved quantum number, hence the flavor eigenstate
which the angular momentum of light degrees of freedom has an unambiguous
value would coincide with the energy states. In reality, the mixing occurs
through a finite quark mass correction and thus is proportional to $1/m_{c}$
Falk:1992wt ; Falk:1992ws .

In this section, we also investigate the charm quark mass dependence of the
mixing angle. In practice, we use the ensemble C11P29S to generate the
correlation function matrices with several bare charm quark masses around the
physical one. By applying the correlated matrix fits for each case, we obtain
the results in Table 4. As one can see, the mixing angle decreases with the
increase of charm quark mass, which is within our expectation.

$m_{c}^{\mathrm{b}}$ | $0.2$ | $0.3$ | $0.4$ | $0.44$ | $\mathbf{0.478}$ | $0.5$ | $0.6$ | $0.7$ | $0.8$
---|---|---|---|---|---|---|---|---|---
$m_{\Xi_{c}}$ (GeV) | $2.0987(25)$ | $2.2380(28)$ | $2.3594(26)$ | $2.4069(26)$ | $\mathbf{2.4587(27)}$ | $2.4793(29)$ | $2.5878(30)$ | $2.6898(30)$ | $2.7859(31)$
$m_{\Xi_{c}^{\prime}}$ (GeV) | $2.1834(24)$ | $2.3249(29)$ | $2.4514(24)$ | $2.4999(24)$ | $\mathbf{2.5536(29)}$ | $2.5718(29)$ | $2.6823(29)$ | $2.7859(30)$ | $2.8835(30)$
$\theta$ (∘) | $1.639(75)$ | $1.349(73)$ | $1.116(49)$ | $1.049(46)$ | $\mathbf{1.002(50)}$ | $0.969(53)$ | $0.847(47)$ | $0.751(42)$ | $0.674(39)$
Table 4:  Results of $m_{\Xi_{c}}$, $m_{\Xi_{c}^{\prime}}$ and $\theta$ at
different bare charm quark masses. The bold ones correspond to a physical
charm quark. These results are extracted from the correlated fit of the
correlation function matrix on C11P29S. The total measurements of non-physical
mass cases are $N_{\mathrm{meas}}=432\times 14$. Uncertainties in all results
are statistical only.  Figure 7: Mixing angle $\theta$ as a function of
$1/m_{\Xi_{c}}$. The red data points denote the results based on unphysical
charm quark masses, and the black data point denotes the physical one.
Uncertainties in all data points are statistical only.

In order to investigate the heavy quark mass dependence of the mixing angle,
and furthermore to predict the behavior at heavy quark limit, we employ a
roughly fit ansatz for the mixing angle $\theta$ as a function of
$m_{\Xi_{c}}$

$\displaystyle\theta=\frac{B_{1}}{m_{\Xi_{c}}}+\frac{B_{2}}{m^{2}_{\Xi_{c}}}.$
(35)

The results obtained from fit are $B_{1}=-2.78(52)$GeV and
$B_{2}=12.9(1.3)\mathrm{GeV^{2}}$ with $\chi^{2}/\mathrm{d.o.f}=0.11$. It
should be noticed that the value of $B_{1,2}$ would suffer from sizeable
discretization effects due to the ensemble we used for this investigation.

##  VII Summary

In summary, we have explored the $\Xi_{c}$-$\Xi_{c}^{\prime}$ mixing by
calculating the two-point correlation functions of the operators to
interpolate the $\Xi_{c}$ and $\Xi_{c}^{\prime}$ baryons from the lattice QCD.
Based on the lattice data, we have adopted two independent methods to
determine the mixing angle between $\Xi_{c}$ and $\Xi_{c}^{\prime}$. Direct
analysis of the correlation function is conducted and a small mixing angle is
found. This observation is also confirmed through the analysis of solving the
generalized eigenvalue problem.

After making the chiral and continuum extrapolation, we found that the mixing
angle $\theta$ is $1.2^{\circ}\pm 0.1^{\circ}$, which is significantly small
than the ones from other methods. This indicates the
$\Xi_{c}-\Xi_{c}^{\prime}$ mixing is insufficient to account for the large
SU(3) symmetry-breaking effects found in weak decays of charmed baryons, and
further mechanisms are required. A combined investigation of $\Lambda_{c}$ and
$\Xi_{c}$ decay form factors from lattice QCD that include the mixing effects
is undergoing.

## Acknowledgement

We are grateful to Fengkun Guo, Xiao-Gang He, Jiajun Wu, and Qiang Zhao for
useful discussions. We thank the anonymous referee for drawing our attention
to the model averaging approach. This work is supported in part by Natural
Science Foundation of China under grant No. U2032102, 12061131006, 12125503,
12293060, 12293061, 12293062. The LQCD calculations were performed using the
Chroma software suite Edwards:2004sx  and QUDA Clark:2009wm ; Babich:2011np ;
Clark:2016rdz  through HIP programming model Bi:2020wpt . The computations in
this paper were run on the Siyuan-1 cluster supported by the Center for High
Performance Computing at Shanghai Jiao Tong University, and Advanced Computing
East China Sub-center.

##  Appendix A Model averaging method

Statistical modeling plays a crucial role in obtaining meaningful results from
lattice field theory calculations. While these models are usually rooted in
physics, multiple model variations may exist for the same lattice data. To
account for the uncertainties associated with model selection, we employ model
averaging Jay:2020jkz  approach to improve the stability of our results. This
approach involves taking a probability-weighted average of all the model
variations, providing a more comprehensive and less conservative approach to
addressing systematic errors.

In practice, we consider a set of models $\\{M\\}$ which fitting all data in
the range $[\\{\mathrm{list~{}of~{}}t_{\mathrm{min}}\\},t_{\mathrm{max}}]$,
and combine the statistical results and relevant parameters obtained from
fitting in models $\\{M\\}$ to estimate the the weight factor

$\displaystyle\mathrm{Pr}(M|D)\approx\rm{exp}\big{[}-\frac{1}{2}(\chi^{2}_{aug}(\mathbf{a}^{*})+2k+2N_{cut})\big{]}$
(36)

where $\chi^{2}_{\rm{aug}}(\mathbf{a}^{*})$ represents the standard best-fit
augmented chi-squared, $k$ corresponds to the number of fit parameters, and
$N_{\rm cut}$ represents the number of removed data points. The normalized
$\mathrm{Pr}(M|D)$ are used to estimate the fit qualities of models $\\{M\\}$.
The model-averaged values $\langle a\rangle$ from the fit parameters $\langle
a\rangle_{M}$ can then be determined using the following procedure:

$\displaystyle\langle a\rangle=\sum_{M}\langle a\rangle_{M}\mathrm{Pr}(M|D),$
(37)

and its error can be estimated by

$\displaystyle\sigma=\sqrt{\sum_{M}\left[\sigma^{2}_{M}+\left(\langle
a\rangle_{M}-\langle a\rangle\right)^{2}\right]\mathrm{Pr}(M|D)},$ (38)

where the $\left(\langle a\rangle_{M}-\langle a\rangle\right)^{2}$ term can be
interpreted as the variability across the different models.

In our analysis in Sec. V, we set $t_{\mathrm{min}}=[6,17]a$ and
$t_{\mathrm{max}}=20a$ for the C11 ensembles, $t_{\mathrm{min}}=[4,20]a$ and
$t_{\mathrm{max}}=24a$ for C08P30S, and $t_{\mathrm{min}}=[9,24]a$ and
$t_{\mathrm{max}}=28a$ for C06P30S. The illustration of fit results and fit
quality annotated by normalized weight factor $\mathrm{Pr}$ and $p$-value are
collected in Fig.4 as well as Fig.8.

(a) C11P22M (b) C11P29S (c) C08P30S (d) C06P30S

Figure 8: Upper panel: effective masses from eigenvalues $\lambda_{n}(t)$
($n=\Xi_{c},\Xi_{c}^{\prime}$), and model averaging results on each ensembles.
Lower panel: the standard $p$-values (dashed lines) and weight factors (solid
lines) which reflect the fit quality.

## References

  * (1) N. Isgur and M. B. Wise, Phys. Lett. B 232, 113-117 (1989) doi:10.1016/0370-2693(89)90566-2
  * (2) N. Isgur and M. B. Wise, Phys. Rev. Lett. 66, 1130-1133 (1991) doi:10.1103/PhysRevLett.66.1130
  * (3) M. Ablikim et al. [BESIII], Phys. Rev. Lett. 115, no.22, 221805 (2015) doi:10.1103/PhysRevLett.115.221805 [arXiv:1510.02610 [hep-ex]].
  * (4) M. Ablikim et al. [BESIII], Phys. Lett. B 767, 42-47 (2017) doi:10.1016/j.physletb.2017.01.047 [arXiv:1611.04382 [hep-ex]].
  * (5) S. Meinel, Phys. Rev. Lett. 118, no.8, 082001 (2017) doi:10.1103/PhysRevLett.118.082001 [arXiv:1611.09696 [hep-lat]].
  * (6) Y. B. Li et al. [Belle], Phys. Rev. Lett. 127, no.12, 121803 (2021) doi:10.1103/PhysRevLett.127.121803 [arXiv:2103.06496 [hep-ex]].
  * (7) Q. A. Zhang, J. Hua, F. Huang, R. Li, Y. Li, C. Lü, C. D. Lu, P. Sun, W. Sun and W. Wang, et al. Chin. Phys. C 46, no.1, 011002 (2022) doi:10.1088/1674-1137/ac2b12 [arXiv:2103.07064 [hep-lat]].
  * (8) X. G. He, F. Huang, W. Wang and Z. P. Xing, Phys. Lett. B 823, 136765 (2021) doi:10.1016/j.physletb.2021.136765 [arXiv:2110.04179 [hep-ph]].
  * (9) C. Q. Geng, X. N. Jin and C. W. Liu, Phys. Lett. B 838, 137736 (2023) doi:10.1016/j.physletb.2023.137736 [arXiv:2210.07211 [hep-ph]].
  * (10) C. Q. Geng, X. N. Jin, C. W. Liu, X. Yu and A. W. Zhou, Phys. Lett. B 839, 137831 (2023) doi:10.1016/j.physletb.2023.137831 [arXiv:2212.02971 [hep-ph]].
  * (11) H. W. Ke and X. Q. Li, Phys. Rev. D 105, no.9, 9 (2022) doi:10.1103/PhysRevD.105.096011 [arXiv:2203.10352 [hep-ph]].
  * (12) T. M. Aliev, A. Ozpineci and V. Zamiralov, Phys. Rev. D 83, 016008 (2011) doi:10.1103/PhysRevD.83.016008 [arXiv:1007.0814 [hep-ph]].
  * (13) Y. Matsui, Nucl. Phys. A 1008, 122139 (2021) doi:10.1016/j.nuclphysa.2021.122139 [arXiv:2011.09653 [hep-ph]].
  * (14) Z. S. Brown, W. Detmold, S. Meinel and K. Orginos, Phys. Rev. D 90, no.9, 094507 (2014) doi:10.1103/PhysRevD.90.094507 [arXiv:1409.0497 [hep-lat]].
  * (15) A. G. Grozin and O. I. Yakovlev, Phys. Lett. B 285, 254-262 (1992) doi:10.1016/0370-2693(92)91462-I [arXiv:hep-ph/9908364 [hep-ph]].
  * (16) G. Wang et al. [$\chi$QCD], Phys. Rev. D 106, no.1, 014512 (2022) doi:10.1103/PhysRevD.106.014512 [arXiv:2111.09329 [hep-lat]].
  * (17) H. Liu, J. He, L. Liu, P. Sun, W. Wang, Y. B. Yang and Q. A. Zhang, [arXiv:2207.00183 [hep-lat]].
  * (18) H. Xing, J. Liang, L. Liu, P. Sun and Y. B. Yang, [arXiv:2210.08555 [hep-lat]].
  * (19) R. L. Workman et al. [Particle Data Group], PTEP 2022, 083C01 (2022) doi:10.1093/ptep/ptac097
  * (20) C. Michael, Nucl. Phys. B 259, 58-76 (1985) doi:10.1016/0550-3213(85)90297-4
  * (21) M. Luscher and U. Wolff, Nucl. Phys. B 339, 222-252 (1990) doi:10.1016/0550-3213(90)90540-T
  * (22) M. Luscher, Nucl. Phys. B 354, 531-578 (1991) doi:10.1016/0550-3213(91)90366-6
  * (23) J. J. Dudek, R. G. Edwards, M. J. Peardon, D. G. Richards and C. E. Thomas, Phys. Rev. D 82, 034508 (2010) doi:10.1103/PhysRevD.82.034508 [arXiv:1004.4930 [hep-ph]].
  * (24) B. Blossier, M. Della Morte, G. von Hippel, T. Mendes and R. Sommer, JHEP 04, 094 (2009) doi:10.1088/1126-6708/2009/04/094 [arXiv:0902.1265 [hep-lat]].
  * (25) W. I. Jay and E. T. Neil, Phys. Rev. D 103, 114502 (2021) doi:10.1103/PhysRevD.103.114502 [arXiv:2008.01069 [stat.ME]].
  * (26) Z. P. Xing and Y. j. Shi, Phys. Rev. D 107, no.7, 074024 (2023) doi:10.1103/PhysRevD.107.074024 [arXiv:2212.09003 [hep-ph]].
  * (27) A. F. Falk and M. Neubert, Phys. Rev. D 47, 2965-2981 (1993) doi:10.1103/PhysRevD.47.2965 [arXiv:hep-ph/9209268 [hep-ph]].
  * (28) A. F. Falk and M. Neubert, Phys. Rev. D 47, 2982-2990 (1993) doi:10.1103/PhysRevD.47.2982 [arXiv:hep-ph/9209269 [hep-ph]].
  * (29) R. G. Edwards et al. [SciDAC, LHPC and UKQCD], Nucl. Phys. B Proc. Suppl. 140, 832 (2005) doi:10.1016/j.nuclphysbps.2004.11.254 [arXiv:hep-lat/0409003 [hep-lat]].
  * (30) M. A. Clark, R. Babich, K. Barros, R. C. Brower and C. Rebbi, Comput. Phys. Commun. 181, 1517-1528 (2010) doi:10.1016/j.cpc.2010.05.002 [arXiv:0911.3191 [hep-lat]].
  * (31) R. Babich, M. A. Clark, B. Joo, G. Shi, R. C. Brower and S. Gottlieb, doi:10.1145/2063384.2063478 [arXiv:1109.2935 [hep-lat]].
  * (32) M. A. Clark, B. Joó, A. Strelchenko, M. Cheng, A. Gambhir and R. Brower, [arXiv:1612.07873 [hep-lat]].
  * (33) Y. J. Bi, Y. Xiao, W. Y. Guo, M. Gong, P. Sun, S. Xu and Y. B. Yang, PoS LATTICE2019, 286 (2020) doi:10.22323/1.363.0286 [arXiv:2001.05706 [hep-lat]].