$\hat{\partial}^{1}(\rho)-\hat{\partial}^{1}\left(\tau_{\mathbb{Q}_{p}}(\mathrm{Ind}_{K/{\mathbb{Q}_{p}}}(\chi))\right)_{\chi\in\mathrm{Irr}(G)}=0.$

We recall that these proofs crucially use a fundamental result of M. Taylor
(see [Frö83, Thm. 31]) which computes the quotient of norm resolvents and
Galois Gauss sums. We also note that the so-called non-ramified characteristic
which occurs in these results is an integral unit itself. ∎

###### Proof of Theorem 1.3.

The crucial input in our proof is a result of Picket and Vinatier in [PV13].
In [BC16, Sec. 5.1] we used this result to construct an integral normal basis
generator $b=p^{2}\alpha_{M}\theta_{2}$ of
$\mathcal{L}=\mathfrak{p}_{N}^{p+1}$. Hence the expression in (29) simplifies
to

$\begin{split}\tilde{R}_{N/K}=&r\hat{\partial}^{1}(\rho)-\chi_{\mathbb{Z}_{p}[G],{B_{\mathrm{dR}}}[G]}(M^{\bullet}(\mathcal{L}),0)\\\
&+\hat{\partial}^{1}({}^{*}(\det(1-Fp^{-d_{K}}u^{-d_{K}})e_{I}))-\hat{\partial}^{1}({}^{*}(\det(1-u^{d_{K}}F^{-1})e_{I}))\\\
&+\hat{\partial}^{1}\left(\det(u)^{-d_{K}(s_{K}\chi(1)+m_{\chi})}\right)_{\chi\in\mathrm{Irr}(G)}-r\hat{\partial}^{1}\left(\tau_{\mathbb{Q}_{p}}(\mathrm{Ind}_{K/{\mathbb{Q}_{p}}}(\chi))\right)_{\chi\in\mathrm{Irr}(G)}\end{split}$
(30)

By [BC16, Prop. 5.2.1] the element
$\hat{\partial}^{1}(\rho)-\hat{\partial}^{1}\left(\tau_{\mathbb{Q}_{p}}(\mathrm{Ind}_{K/{\mathbb{Q}_{p}}}(\chi))\right)_{\chi\in\mathrm{Irr}(G)}$
is represented by $\hat{\partial}^{1}(\eta^{-1})$ with $\eta$ as in [BC16,
Prop. 5.2.1]. Furthermore,

$\hat{\partial}^{1}\left(\det(u)^{-d_{K}(s_{K}\chi(1)+m_{\chi})}\right)_{\chi\in\mathrm{Irr}(G)}=\hat{\partial}^{1}\left(\det(u)^{-mm_{\chi})}\right)_{\chi\in\mathrm{Irr}(G)}.$

From now on, we have to distinguish the three conditions in the statement of
the Theorem. Let us start with (a), i.e. Hypothesis (I). In this case
$\tilde{R}_{N/K}$ is represented by

$\begin{split}\tilde{r}_{\chi\phi}&=\begin{cases}\frac{\mathfrak{d}_{K}^{r/2}\mathcal{N}_{K/{\mathbb{Q}_{p}}}(\theta_{2}|\phi)^{r}p^{2rm}\det(1-p^{-m}u^{-m}\phi(b^{-1}))}{p^{rm}\det(1-u^{m}\phi(b))}&\text{if
$\chi=\chi_{0}$}\\\
\frac{\mathfrak{d}_{K}^{r/2}\mathcal{N}_{K/{\mathbb{Q}_{p}}}(\theta_{2}|\phi)^{r}p^{rm}\chi(4)^{r}\phi(b)^{-2r}\det(u)^{-2m}}{(-1)^{r(m-1)}(u^{-m}\phi(b)^{-1}-1)(\chi(a)-1)^{rm(p-1)}}&\text{if
$\chi\neq\chi_{0}$}\end{cases}\\\
&=\begin{cases}(-1)^{r}\frac{\mathfrak{d}_{K}^{r/2}\mathcal{N}_{K/{\mathbb{Q}_{p}}}(\theta_{2}|\phi)^{r}}{\det(1-u^{m}\phi(b))}\det(u)^{-m}\phi(b)^{-r}\det(1-p^{m}u^{m}\phi(b))&\text{if
$\chi=\chi_{0}$}\\\
(-1)^{r}\frac{\mathfrak{d}_{K}^{r/2}\mathcal{N}_{K/{\mathbb{Q}_{p}}}(\theta_{2}|\phi)^{r}}{\det(1-u^{m}\phi(b))}\det(u)^{-m}\phi(b)^{-r}\left(\frac{-p}{(\chi(a)-1)^{p-1}}\right)^{rm}\chi(4)^{r}&\text{if
$\chi\neq\chi_{0}$.}\end{cases}\end{split}$

Let
$W_{\theta_{2}}\in\mathcal{O}_{p}^{t}[G]\subset\overline{\mathbb{Z}_{p}^{\mathrm{nr}}}[G]^{\times}$
be defined by

$\chi\phi(W_{\theta_{2}})=(-1)^{r}\mathfrak{d}_{K}^{r/2}\mathcal{N}_{K/{\mathbb{Q}_{p}}}(\theta_{2}|\phi)^{r}.$

So

$\tilde{r}=\frac{W_{\theta_{2}}\det(u)^{-m}b^{-r}}{\det(1-u^{m}b)}\left(\det(1-p^{m}u^{m}b)e_{a}+(-\tilde{u})^{m}\sigma_{4}(1-e_{a})\right),$

where $\tilde{u}\in\mathbb{Z}_{p}[a]$ is a unit whose augmentation is
congruent to $1$ modulo $p$, as in [BC17, Sect. 8]. Then the proof works as in
[BC17].

Let us now consider case (b), i.e. Hypothesis (T). In this case
$\tilde{R}_{N/K}$ is represented by

$\begin{split}\tilde{r}_{\chi\phi}&=\begin{cases}(-1)^{r}\frac{\mathfrak{d}_{K}^{r/2}\mathcal{N}_{K/{\mathbb{Q}_{p}}}(\theta_{2}|\phi)^{r}p^{2rm}\det(u^{m}\phi(b)-1)\det(1-p^{-m}u^{-m}\phi(b^{-1}))}{\det(u)^{m\tilde{m}}\phi(b)^{r\tilde{m}}p^{rm}\det(1-u^{m}\phi(b))}&\text{if
$\chi=\chi_{0}$}\\\
(-1)^{r(m+1)}\frac{\mathfrak{d}_{K}^{r/2}\mathcal{N}_{K/{\mathbb{Q}_{p}}}(\theta_{2}|\phi)^{r}p^{rm}\chi(4)^{r}\phi(b)^{-2r}\det(u)^{-2m}}{(\det(u)^{m}\phi(b)^{r})^{-\tilde{m}(m-1)}(\chi(a)-1)^{rm(p-1)}}&\text{if
$\chi\neq\chi_{0}$}\end{cases}\\\
&=\begin{cases}\chi\phi(W_{\theta_{2}})\phi(b)^{-r-r\tilde{m}}\det(u)^{-m-m\tilde{m}}\det(1-p^{m}u^{m}\phi(b))&\\!\text{if
$\chi=\chi_{0}$}\\\
\chi\phi(W_{\theta_{2}})\phi(b)^{-r-r\tilde{m}}\det(u)^{-m-m\tilde{m}}\left(\frac{-p}{(\chi(a)-1)^{p-1}}\right)^{rm}\chi(4)\det(u)^{m^{2}\tilde{m}-m}&\\!\text{if
$\chi\neq\chi_{0}$,}\end{cases}\end{split}$

So

$\tilde{r}=W_{\theta_{2}}b^{r-r\tilde{m}}\det(u)^{-m-m\tilde{m}}\left(\det(1-p^{m}u^{m}b)e_{a}+(-\tilde{u})^{m}\sigma_{4}\det(u)^{m^{2}\tilde{m}-m}(1-e_{a})\right).$

As in [BC17], Hypothesis (T) implies that $\det(u)^{m^{2}\tilde{m}-m}\equiv
1\pmod{1-\zeta_{p}}$. Then the proof works as in [BC17], using the
computations of the case of Hypothesis (I).

It remains to consider the case (c). Let $E/K^{\prime}$ be the unramified
extension of degree $\tilde{d}$. By the functoriality result of Proposition
7.2 (a) it is enough to show $R_{NE/K}=0$, which is true since for $NE/K$ we
can apply part (b) of the theorem. ∎

We finally prove Theorems 1.5 and 1.6. For that purpose we recall the
following lemma.

######  Lemma 8.1.

Let $A$ be an $r$-dimensional abelian variety defined over the $p$-adic number
field $L$ with good ordinary reduction and let $U$ be the twist matrix of $A$.
Then $\det(U-1)\neq 0$.

###### Proof.

This is shown in [Maz72, Cor. 4.38] or in the proof of [LR78, Thm. 2]. ∎

###### Proof of Theorems 1.5 and 1.6.

The results are immediate from Theorems 1.1 and 1.3. Note that Hypothesis (F)
and the condition $\det(\rho^{\mathrm{nr}}(F_{N})^{\tilde{d}}-1)\neq 0$ are
automatically satisfied by Lemma 8.1. ∎

## References

  * [BB05] Manuel Breuning and David Burns.  Additivity of Euler characteristics in relative algebraic $K$-groups.  Homology, Homotopy Appl., 7(3):11–36, 2005.
  * [BB08] Denis Benois and Laurent Berger.  Théorie d’Iwasawa des représentations cristallines. II.  Comment. Math. Helv., 83(3):603–677, 2008.
  * [BC16] Werner Bley and Alessandro Cobbe.  Equivariant epsilon constant conjectures for weakly ramified extensions.  Math. Z., 283(3-4):1217–1244, 2016.
  * [BC17] Werner Bley and Alessandro Cobbe.  The equivariant local $\varepsilon$-constant conjecture for unramified twists of $\mathbb{Z}_{p}(1)$.  Acta Arith., 178(4):313–383, 2017.
  * [BK90] Spencer Bloch and Kazuya Kato.  $L$-functions and Tamagawa numbers of motives.  In The Grothendieck Festschrift, Vol. I, volume 86 of Progr. Math., pages 333–400. Birkhäuser Boston, Boston, MA, 1990.
  * [BN20] David Burns and Andreas Nickel.  Equivariant local epsilon constants and Iwasawa theory.  Preprint, 2020.
  * [Bre04a] Manuel Breuning.  Equivariant Epsilon Constants for Galois Extensions of Number Fields and $p$-adic Fields.  PhD thesis, King’s College London, 2004.
  * [Bre04b] Manuel Breuning.  Equivariant local epsilon constants and étale cohomology.  J. London Math. Soc. (2), 70(2):289–306, 2004.
  * [BW09] Werner Bley and Stephen M. J. Wilson.  Computations in relative algebraic $K$-groups.  LMS J. Comput. Math., 12:166–194, 2009.
  * [Cob18] Alessandro Cobbe.  A representative of $R\Gamma(N,T)$ for some unramified twists of $\mathbb{Z}_{p}^{r}$.  arXiv:1804.06132, April 2018.
  * [EN18] Nils Ellerbrock and Andreas Nickel.  On formal groups and Tate cohomology in local fields.  Acta Arith., 182(3):285–299, 2018.
  * [Frö68] A. Fröhlich.  Formal groups.  Lecture Notes in Mathematics, No. 74. Springer-Verlag, Berlin-New York, 1968.
  * [Frö83] Albrecht Fröhlich.  Galois module structure of algebraic integers, volume 1 of Ergebnisse der Mathematik und ihrer Grenzgebiete (3) [Results in Mathematics and Related Areas (3)].  Springer-Verlag, Berlin, 1983.
  * [Haz78] Michiel Hazewinkel.  Formal groups and applications, volume 78 of Pure and Applied Mathematics.  Academic Press, Inc. [Harcourt Brace Jovanovich, Publishers], New York-London, 1978.
  * [IV16] Dmitriy Izychev and Otmar Venjakob.  Equivariant epsilon conjecture for 1-dimensional Lubin-Tate groups.  J. Théor. Nombres Bordeaux, 28(2):485–521, 2016.
  * [Köc04] Bernhard Köck.  Galois structure of Zariski cohomology for weakly ramified covers of curves.  Amer. J. Math., 126(5):1085–1107, 2004.
  * [LR78] Jonathan Lubin and Michael I. Rosen.  The norm map for ordinary abelian varieties.  J. Algebra, 52(1):236–240, 1978.
  * [Mar77] J. Martinet.  Character theory and Artin $L$-functions.  In Algebraic number fields: $L$-functions and Galois properties (Proc. Sympos., Univ. Durham, Durham, 1975), pages 1–87. Academic Press, London, 1977.
  * [Maz72] Barry Mazur.  Rational points of abelian varieties with values in towers of number fields.  Invent. Math., 18:183–266, 1972.
  * [Neu92] Jürgen Neukirch.  Algebraic number theory. (Algebraische Zahlentheorie.).  Berlin etc.: Springer-Verlag, 1992.
  * [Nic20] Andreas Nickel.  An equivariant Iwasawa main conjecture for local fields.  Preprint, 2020.
  * [PV13] Erik Jarl Pickett and Stéphane Vinatier.  Self-dual integral normal bases and Galois module structure.  Compos. Math., 149(7):1175–1202, 2013.
  * [Sil09] Joseph H. Silverman.  The arithmetic of elliptic curves, volume 106 of Graduate Texts in Mathematics.  Springer, Dordrecht, second edition, 2009.
  * [Swa70] Richard G. Swan.  $K$-theory of finite groups and orders.  Lecture Notes in Mathematics, Vol. 149. Springer-Verlag, Berlin-New York, 1970.
  * [Was97] Lawrence C. Washington.  Introduction to cyclotomic fields, volume 83 of Graduate Texts in Mathematics.  Springer-Verlag, New York, second edition, 1997.


# Active Learning in Symbolic Regression with Physical Constraints

Jorge Medina
Department of Chemical Engineering
University of Rochester
<EMAIL_ADDRESS>
Andrew D. White
Department of Chemical Engineering
University of Rochester
<EMAIL_ADDRESS>
Corresponding author

###### Abstract

Evolutionary symbolic regression (SR) fits a symbolic equation to data, which
gives a concise interpretable model. We explore using SR as a method to
propose which data to gather in an active learning setting with physical
constraints. SR with active learning proposes which experiments to do next.
Active learning is done with query by committee, where the Pareto frontier of
equations is the committee. The physical constraints improve proposed
equations in very low data settings. These approaches reduce the data required
for SR and achieves state of the art results in data required to rediscover
known equations.

_K_ eywords Symbolic Regression $\cdot$ Query by Committee $\cdot$ Physical
Constraints

##  1 Introduction

A variety of established methods exist for modeling data, ranging from
traditional machine learning techniques (linear regression, ridge regression,
polynomial regression) to deep learning approaches (neural networks). However,
these methods suffer from constraints and/or interpretability, such as
limiting the model to a particular shape (e.g., linear), or being too complex
to interpret (black box models). Symbolic Regression (SR) is less constrained
and searches through the mathematical space of equations. SR allows for
discovering a broader range of functional relationships, including those with
nonlinear or intricate interactions between variables. [1].

SR is a machine learning technique to propose equations describing data,
aiming to identify optimal functional forms and parameters [2, 1]. As the
number of features increases, the search space grows exponentially, leading to
various strategies to search in such space through evolutionary algorithms,
Bayesian models, sparse models, and neural networks [3, 4, 5, 6, 7, 8]. SR has
been applied successfully to scientific and engineering problems,
demonstrating their potential for uncovering real-world mathematical
equations. Weng et al. (2009) utilized SR to propose a model for pipe
deterioration, while Baichang et al. (2020) leveraged it to introduce a new
descriptor for identifying novel perovskite catalysts. More recently, Li et
al. (2023) employed SR to determine electron transfer models of minerals under
pressure, among other applications.[9, 1, 10, 11, 12, 13, 14, 15, 16, 17, 18,
19].

The complexity of ground-truth equations and the presence of noise in the data
can hinder the discovery of equations. This occurs in two phases: discoverable
and non-discoverable ground truth [20]. Some ground-truth equations possess
such complexity that rediscovering them becomes unattainable, even in the
absence of noise[20]. This underlies the necessity of devising methods to
manage noise and enhance SR. As the ground-truth equation complexity
escalates, SR methods often demand many labeled samples, which can incur high
costs or require extensive acquisition time. Active learning strategies, such
as Query by Committee (QBC), can tackle these challenges by reducing the
number of data needed through the intelligent selection of the most
informative samples for labeling [21, 22]. In this paper, we introduce two
straightforward approaches to reduce the data needed for SR. (I) By
incorporating the QBC strategy into the SR framework, and (II) by including
soft constraints based on background knowledge in the discovery of equations.

The optimization process may not find equations with physical meaning; any
equation that fits the data could be found. Applying strict restrictions to
search only for physically meaningful equations could make finding suitable
expressions as challenging as the optimization itself [23]. As a result, using
a soft constraint can effectively guide the search in the desired direction
[23]. This reasoning supports our second approach: implementing soft
constraints based on prior knowledge to limit the search space appropriately,
allowing for faster convergence and better resilience against noisy settings.
Various strategies have been employed to address this, such as constraining
equations to maintain unit correctness [24], using physical properties to
simplify the problem [5, 25], guiding the search for predefined shapes or
forms known to be present in the system or deduced from the dataset [26, 25],
or checking correctness after equations are rediscovered [27, 19]. Emgle et
al, added auxiliary first and second derivatives to be used for constraint
optimization [28] Recently, a soft constrained approach was applied to
adsorption equations, revealing that some constraints might not work as
expected for genetic algorithms [27].

We employ the Feynman dataset as a benchmark for comparing the data points
required to rederive equations. AI Feynman[5], at its publication presented
state of the art results and is used commonly for performance comparisons. Our
findings indicate that our approach frequently surpasses this method in SR.
The rest of this paper explores the principles of the active learning
technique used, the incorporation of physical constraints, and presents our
findings and analysis. As previously states, genetic algorithms are commonly
used for SR and its structure can be easily modified.[23]

Genetic Algorithms (GAs) are computational methods inspired by the principles
of natural selection, incorporating processes such as reproduction, mutation,
and crossover. They involve a population of candidate solutions, referred to
as individuals, which consist of components known as genes. Individuals can be
represented using various data structures, such as strings, trees, or stacks.
Common terms used in GA literature include parents, children, genes and
fitness.[23]

The fitness assigns a score to each individual based on how well they perform
in solving the target problem. This score is then used to guide the selection
process, allowing fitter individuals to have a higher probability of
contributing their genetic material to the next generation, thereby promoting
the evolution of increasingly effective solutions over time. [23]

In this study, we focus on a specific application of GAs using binary trees as
the data structure for representing candidate solutions. Here, the nodes of
the tree represent mathematical operators (e.g., $+$, $-$, $/$, $\sin$,
$\cos$), while the leaves correspond to either constants or features. This
binary tree representation not only serves as an implementation formalism but
also aids in visualization for better understanding and interpretation. See
Figure 1. For more details on the workings of genetic algorithms, such as
initialization, reproduction, crossover, and mutation, please refer to the
Supplementary Information (SI).

+*AB-$\sin$CD Figure 1: Binary tree representation of the equation
$(A*B)+(\sin(C)-D)$

###  1.1 Inclusion of Background Knowledge: Physical Constraints

In the optimization process, the search is primarily driven by accuracy.
However, there is no guarantee that the proposed equation will comply with all
the physical constraints of the system, such as divergence, symmetry of
variables, or conservation laws. To take advantage of this type of prior
knowledge and guide the search towards physically meaningful solutions, the
fitness function can be modified to incorporate these constraints as can be
seen in Equation 1 where a penalty term is included to balance both accuracy
and physical consistency. Given a dataset $(\vec{x},y)\textrm{ of size }N$ and
a member of the population $\xi$ the fitness function $f$ is defined as:

$f(\xi,\vec{x},y)=L(\xi,\vec{x},y)+\lambda\phi_{k}(\xi,\vec{x})+p_{l}(\xi)$
(1)

In this equation, $L(\xi,\vec{x},y)$ represents the Root Mean Square Error
(RMSE), quantifying the accuracy of each equation over the dataset. The term
$\phi_{k}$ serves as a penalty term based on the constraint $k$, while
$\lambda$ is an adjustable parameter set to balance the trade-off between
accuracy and physical consistency; for our experiments, we’ve chosen a value
of $\lambda=100$. The complexity of the equations introduces an additional
penalty, denoted by $p_{l}$, which promotes size homogeneity in the
population. Ultimately, the fitness or score of the equation $\xi$ is
signified by $f(\xi)$.

Some common types of constraints that can be relevant to chemical engineering
applications include divergence, symmetry, and sign constraints. Basically, it
‘activates’ the penalty when the physical constraint is not complying with the
equation.

Divergence is the behavior of a function as its input approaches a particular
point. In some cases, it may be necessary to ensure that the function does not
diverge (i.e., does not becomes infinite). For example, in chemical
engineering applications, it might be important to guarantee that the
concentration of a reactant remains finite. A simple way to represent a
divergence constraint is as follows:

$\phi_{\textrm{div}}(\xi,a)=\left\\{\begin{array}[]{ll}0,&|\xi(a)|=\infty,\\\
1,&|\xi(a)|\neq\infty.\end{array}\right.$ (2)

Symmetry, in this context, refers to the property of a function where swapping
two variables does not change its value (e.g.,
$\xi(x_{1},x_{2})=\xi(x_{2},x_{1})$) or any example where commutativity
applies. This can be particularly relevant in cases where the order of
variables should not affect the system’s behavior.

$\phi_{\textrm{symm}}(\xi,\vec{x})=\left\\{\begin{array}[]{lr}0,&\xi(x_{1},x_{2})=\xi(x_{2},x_{1})\\\
1,&\xi(x_{1},x_{2})\neq\xi(x_{2},x_{1})\end{array}\right.$ (3)

Sign Constraints (Positivity/Negativity) In some systems, measured properties
are always positive or negative within a certain domain. For example, this
could be relevant when dealing with concentrations, which should always be
non-negative. A simple way to represent sign constraints is as follows:

$\phi_{\textrm{sign}}(\xi,\vec{x})=\left\\{\begin{array}[]{lr}0,&\xi(\vec{x})\geq
0\\\ 1,&\textrm{otherwise}\end{array}\right.$ (4)

For negative is the same concept, with the opposite inequality.

Monotonicity Constraints (Increasing/Decreasing Functions) In some systems,
the relationship between variables can be strictly increasing or decreasing
within a certain domain. For example, this could be relevant when studying the
relationship between temperature and pressure in an ideal gas under constant
volume, where temperature and pressure are directly proportional. Monotonicity
constraints can be represented based on the sign of the derivative as follows:

$\phi_{\textrm{mono}}(\xi,\vec{x}_{i})=\left\\{\begin{array}[]{lr}0,&\frac{d\xi(\vec{x})}{dx_{i}}\geq
0\\\ 1,&\textrm{otherwise}\end{array}\right.$ (5)

This constraint enforces that the function is either strictly increasing or,
flipping the inequality, strictly decreasing (non-positive derivative) within
the specified domain or feature.

In addition to these constraints, chemical engineering applications may
require the incorporation of other types of physical constraints, such as
conservation laws, boundary conditions, or monotonicity constraints, depending
on the specific problem being addressed.

###  1.2 Active Learning

Active learning is a subfield of machine learning that focuses on selecting
the most informative samples to label, with the goal of improving the learning
performance while minimizing the labeling effort[29, 30, 31, 32, 33]. In SR,
active learning can be particularly useful for exploring the search space more
efficiently and identifying meaningful equations that describe the underlying
system. One well-established active learning method is QBC [34, 35, 36].

QBC works by selecting a sample based on the disagreement among a group of
experts (the committee). Members of this group should be distinct enough that
their disagreement provides valuable information, while being consistently
good at describing the system being studied. The algorithm can be described as
follows:

  1. 1.

Train a model and select a committee.

  2. 2.

Estimate labels of an unsampled pool with each member of the committee.

  3. 3.

Measure and select the label with the highest disagreement.

  4. 4.

Label the selected point and repeat until necessary.

Disagreement is defined in terms of an unlabeled dataset $D_{u}$, and $N$
committee members $K=\\{\Xi_{1},\Xi_{2},....,\Xi_{N}\\}$ Each member $i$ makes
a prediction on every unlabeled point $j$, denoted as
$\Xi_{i}(\vec{x_{j}})=Y_{ij}$. For this study, committee members are selected
from the Pareto frontier, as they represent the models with the best trade-off
between accuracy and complexity. They generally provide valuable insights into
the whole domain of the system. A depiction of the process can is shown in
figure 2

Figure 2: Query By Committee Depiction. Symbolic Regression outputs a Pareto
frontier of equations that act as expert models, and measure disagreement from
unlabeled data to select the most informative point to label and iterate.

Upon examining the literature, the authors discovered studies that utilized
QBC in a SR framework. These studies yielded best results when applied to low
complexity equations, where active learning setups were deemed unnecessary
since rediscovery occurred with minimal data. However, the approach faced
limitations when dealing with high complexity equations.[37]. Bongard et al.
[38] discussed the implementation of QBC to generate and sample informative
training examples in a grammatical inference problem.

##### Measurement of Disagreement

Measure of disagreements $d(K,\vec{x})$ take the dataset $\vec{x}$ and the
committee members $K$ Two different measures of disagreements were tested. The
first measure employed is the coefficient of variation (CV), which is the
ratio of the standard deviation ($\sigma$) to the mean ($\mu$) of the
predicted labels by all members of the committee. This offers a comprehensible
metric for assessing the dispersion of the predictions:

$\mu(K,\vec{x})=\frac{1}{N}\sum_{i=1}^{N}\Xi_{i}(\vec{x})$ (6)

$\sigma(K,\vec{x})=\sqrt{\frac{1}{N}\sum_{i=1}^{N}(\Xi_{i}(\vec{x})-\mu(K,\vec{x}))^{2}}$
(7)

$d_{CV}(K,\vec{x})=\frac{\sigma(K,\vec{x})}{\mu(K,\vec{x})}$ (8)

The coefficient of variation facilitates the comparison of dispersion across
different data points or scales.

The second measure is an information-based measure of disagreement [39], which
takes into account the relative differences between the predictions:

$d_{IBMD}(K,\vec{x})=\frac{1}{|C^{N}_{2}|}\sum\limits_{(i,j)\in
C^{N}_{2}}\log\left(\frac{|\Xi_{i}(\vec{x})-\Xi_{j}(\vec{x})|}{\max(\Xi_{i}(\vec{x}),\Xi_{j}(\vec{x}))}+1\right)$
(9)

where $(\xi_{i},\xi_{j})$ represents a pair of equations from the committee,
and $|C^{N}_{2}|$ is the number of possible pair combinations of $N$ members.
The IBMD measure calculates the logarithm of the ratio between the absolute
difference of the predictions and the maximum of the two predictions,
incremented by 1. These logarithmic values are then averaged for all possible
pair combinations according to Equation 9.

##  2 Methodology

Our research methodology is structured into three well-defined phases: (I)
establishing proof of concept, (II) benchmarking against AI Feynman, and (III)
conducting robustness analysis. This design enables a systematic evaluation of
our proposed SR approach.

SymbolicRegression.jl or PySR (python front end name) [3] is a Julia package
designed for discovering mathematical equations through regularized evolution.
This high-performance package offers an efficient implementation of SR,
enabling the identification of underlying mathematical relationships within
complex data sets. Equation representations are based on binary trees, with
numerous hyperparameters for customization. To ensure fair and straightforward
comparisons across multiple cases, the same parameters were chosen for all
experiments. These primarily consisted of default parameters for the
evolutionary process. The maximum number of iterations was adjusted to 100,
and the selected binary operators included addition ($+$), subtraction ($-$),
multiplication ($*$), and division ($/$). The chosen unary operators were
inverse (inv), square, cube, exponent ($\exp$), square root ($sqrt$), and
cosine ($\cos$). Restrictions were applied to prevent nested operators, such
as $\cos(\cos(x))$.

The software’s simplicity and efficiency, combined with its ability to
customize loss functions, made it a suitable choice for our research. PySR’s
use of genetic algorithms leverages Julia’s high-performance computing
capabilities, offering a straightforward and adaptable solution for performing
SR in our study.

In assessing our methodology, we leveraged the observation that rediscovering
an equation often leads to a marked enhancement in accuracy (more than 10
orders of magnitude for noiseless systems), accompanied by a minimal rise in
equation complexity (1 or 2). This insight facilitates an automated process
for configuring experiments and ascertaining success using loss metrics.
Nevertheless, the presence of noise results in a diminished increase in
accuracy, making it essential to employ a more hands-on approach for success
evaluation, which involves visually examining whether the equation has been
successfully rediscovered.

In the initial stage, our primary goal was to demonstrate the effectiveness of
our method through a comprehensive toy case focused on the gravitational
equation. This aimed to validate the program’s ability to rediscover known
equations within a controlled setting. To examine the impact of incorporating
physical constraints, we explored the one-dimensional gravity
equation:$\frac{Gm_{1}m_{2}}{r^{2}}$. All variables were held constant (number
of iterations, operator set, default evolution parameters), with the exception
of the fitness function. We acknowledge that, for this equation, given
sufficient time (e.g., 100 iterations), the algorithm will inevitably
rediscover the equation. Therefore, we limited the number of iterations to 20,
making the task more challenging for the algorithm.

Secondly, in order to evaluate the performance of Query By Committee, we set
up an active learning scenario involving the following equation:
$f=(1/2\pi)^{1/2}\exp(-(x_{1}-x_{2})/\sigma)^{2}/2$ from the Feynman dataset.
We started with three data points and examined the percentage of rediscovery
as new data points were added based on either random sampling, QBC with IBMD,
and CV. This stepwise approach allowed us to systematically compare both
approaches.

Next, we benchmark the model against AIFeynman, a SR framework using the
Feynman equations datasets [5]. This comparison enables us to assess the
performance of our approach relative to existing methods, highlighting its
strengths and potential weaknesses.

Finally, we subject our model to a series of stress tests, introducing two
distinct levels of noise into the data. This stage aims to assess the
robustness of our findings when faced with real-world scenarios where data has
noise from experiments.

Two notable limitations are that as noise levels increase, fully automated
experiments are not possible. Even in cases where the ground truth equation
was exactly rediscovered in the Pareto-frontier, there was no big jump in
accuracy ( around 2 orders of magnitude or less) which could cause false
convergence. Consequently, to allow time for human evaluation of success, the
experiment was narrowed down to just three representative case examples to
assess the model’s performance in the presence of noise. Second, due to the
penalty, the loss change was not enough to detect convergence. This presented
a challenge in determining the optimal stopping point for the algorithm. To
address this issue, it was decided to only use the physical constraints for a
limited time during the evolutionary process. After this initial period, the
constraints were removed, allowing exploration of a wider range of potential
solutions. This approach aimed to strike a balance between guiding the model
towards physically meaningful solutions and giving the algorithm freedom to
take full advantage of the evolutionary operations with a head start in the
right direction. Figure B.3 shows the effectiveness of this approach.

##  3 Results

###  3.1 Query By Committee

Both measures of disagreement mentioned earlier appear to function similarly.
As depicted in Figure 3a, comparing the disagreement between the equations on
the Pareto frontier using either method (CV or IBMD) yielded comparable
outcomes. A correlation coefficient close to one was anticipated, since both
measurements represent the range of predictions from each equation. Figure 3b
demonstrates that employing QBC leads to the more frequent rediscovery of the
same equation using fewer data points compared to merely adding new data at
random. Consequently, the remaining experiments presented in Figure 4 were
conducted using QBC with IBMD.

The findings presented in Figure 4 demonstrate that our model, with the
implementation of QBC, outperforms AIFeynman in 8 out of 12 test cases in
rediscovering equations using fewer data points. However, it is important to
note that AIFeynman’s results were only reported in orders of magnitude of 10
data points, which limited our ability to make a direct comparison between the
precise number of data points required for each model.

a)  b)

Figure 3: Evaluation of disagreement measures. a) Disagreement analysis of 15
equations from PySR and 10,000 data points, resulting in a correlation
coefficient of 0.91. Similar points exhibit maximum disagreement. b)
Comparison of $f=(1/2\pi)^{1/2}\exp(-(x_{1}-x_{2})/\sigma)^{2}/2$ rediscovery
performance versus the number of data points when adding new labeled points
through active learning and random addition. Figure 4: Difference in points
needed for rediscovery of non-trivial equations of Feynman dataset. From the
twelve tested equations PySR with QBC outperformed in eight of them. Arrows
show direction of improvement (lower the better or vice-versa)

###  3.2 Physical Constraints

To evaluate the effectiveness of incorporating physical constraints and Query
By Committee in SR, a proof of concept study was conducted using the well-
known gravitational force equation $F=\frac{Gm_{1}m_{2}}{r^{2}}$ as a test
case. This equation was chosen due to its simplicity and multiple possible
constraints, providing a clear benchmark for assessing the impact of the
proposed methods.

Four types of constraints were considered for this experiment: divergence,
asymptote, monotonicity, and positivity. The goal was to investigate how each
of these constraints affects the ability of SR to rediscover the gravitational
force equation.

Figure 5: Rediscovery of Gravitational Law with different types of
constraints.

Out of the four constraints examined, SR incorporating divergence constraints
in the loss function demonstrated superior performance compared to non-
constrained regression. This improvement led to an approximately 10% increase
in the rediscovery rate, as illustrated in Figure 5. Conversely, the remaining
three constraints resulted in reduced performance. Interestingly, even with a
relatively simple equation, certain constraints appear to hinder the search in
the optimal direction, which may be an indication that even soft constraints
can deteriorate the optimization process.

###  3.3 Robustness Experiments

As previously stated, three distinct equations, selected for their non-
triviality and inherent physical constraint, were utilized to assess the
method both in noiseless and noisy setups. The noise level was increased to
the point where no rediscovery was obtained. In these cases, convergence was
rare and not measured enough to obtain valuable statistics. therefore results
are either successful (rediscovery) or unsuccessful (No Rediscovery). Table
LABEL:rediscoverywithnoise shows another comparison between un-constrained and
constrained optimization where using constraints is able to rediscover the
ground truth in a ‘high’ noise level.

The first equation, as depicted in Table LABEL:results1 and Figure 6,
demonstrates an enhancement in rediscovery by employing symmetry constraints,
with a p-value of 0.061. Intriguingly and against our expectations, the
performance exhibits improvement as noise is added; however, this could be
attributed to the inherent limitations of conducting a fully automated process
with noise. Owing to false convergence, iterations often converged
prematurely, resulting in a higher number of runs with reduced data.
Nevertheless, a subtle improvement can be observed in the box plot for both
noise experiments.

The other two equations results (SI) show less improvement with constrained
optimization. Comparing the complexity of three equations can be seen that for
more complex problems constraints have positive effects for optimization.

Table 1: Rediscovery of $\sqrt{x_{1}^{2}+x_{2}^{2}-2x_{1}x_{2}\cos(\theta_{1}-\theta_{2})}$ Type of Experiment | Noiseless | 0.01 Noise | 0.1 Noise
---|---|---|---
No Constraint | 213 | 71.20 | No Rediscovery
Symmetry | 139 | 59.19 | No Rediscovery
Figure 6: Illustrating the Effects of Noise on Rediscovery: A Comparison
between Constrained and Unconstrained Optimizations. The figure presents
results for the expression
$\sqrt{x_{1}^{2}+x_{2}^{2}-2x_{1}x_{2}\cos(\theta_{1}-\theta_{2})}$. P-values
are 0.061 for noiseless conditions and 0.443 for a noise level of 0.01. The
joint p-value for comparing unconstrained and symmetry-constrained
optimization is 0.127.

##  4 Conclusion

This study presents a simple approach to improve Evolutionary Symbolic
Regression (ESR) performance by integrating QBC and physical constraints
regularization. The proposed method prioritizes informative samples using the
Pareto frontier of candidate equations and guides the algorithm toward
physically meaningful expressions by incorporating physical constraints as
regularization terms.

Our results demonstrate that the application of QBC facilitates the
rediscovery of equations using fewer data points than AI Feynman. Moreover,
incorporating symmetry as a constraint in symbolic regression showed
improvement against a non-constrained regression, while the other constraints,
like asymptotes, presented less performance. The removal of constraints needs
to be explored more indepth to found the most appropiate level of constraints
without hindering the evolutionary process.

The proposed method exhibits resilience against noise, successfully
rediscovering equations under various noise levels for specific cases. By
integrating QBC and physical constraint regularization, this approach provides
a framework for uncovering meaningful mathematical expressions across a wide
range of applications, even when faced with limited data.

## Acknowledgments

We thank the Center for Integrated Research Computing (CIRC) at the University
of Rochester for providing computational resources and technical support.

This work has been supported by funds from the Robert L. and Mary L. Sproull
Fellowship gift and U.S. Department of Energy, Grant No. DE-SC0023354.

## References

  * [1] Michael Schmidt and Hod Lipson.  Distilling free-form natural laws from experimental data.  Science, 324(5923):81–85, 2009.
  * [2] JohnR. Koza.  Genetic programming as a means for programming computers by natural selection.  Statistics and Computing, 4(2), June 1994.
  * [3] Miles Cranmer.  Pysr: Fast & parallelized symbolic regression in python/julia, September 2020.
  * [4] Runhai Ouyang, Stefano Curtarolo, Emre Ahmetcik, Matthias Scheffler, and Luca M. Ghiringhelli.  SISSO: A compressed-sensing method for identifying the best low-dimensional descriptor in an immensity of offered candidates.  Physical Review Materials, 2(8):083802, August 2018.
  * [5] Silviu-Marian Udrescu and Max Tegmark.  AI Feynman: A physics-inspired method for symbolic regression.  Science Advances, 6(16):eaay2631, April 2020.
  * [6] Roger Guimerà, Ignasi Reichardt, Antoni Aguilar-Mogas, Francesco A. Massucci, Manuel Miranda, Jordi Pallarès, and Marta Sales-Pardo.  A Bayesian machine scientist to aid in the solution of challenging scientific problems.  Science Advances, 6(5):eaav6971, January 2020.
  * [7] Richa Ramesh Naik, Armi Tiihonen, Janak Thapa, Clio Batali, Zhe Liu, Shijing Sun, and Tonio Buonassisi.  Discovering equations that govern experimental materials stability under environmental stress using scientific machine learning.  npj Computational Materials, 8(1):72, April 2022.
  * [8] Martin Schmelzer, Richard P. Dwight, and Paola Cinnella.  Discovery of Algebraic Reynolds-Stress Models Using Sparse Symbolic Regression.  Flow, Turbulence and Combustion, 104(2-3):579–603, March 2020.
  * [9] Ben McKay, Mark Willis, and Geoffrey Barton.  Steady-state modelling of chemical process systems using genetic programming.  Computers & Chemical Engineering, 21(9):981–996, 1997.
  * [10] William La Cava, Kourosh Danai, and Lee Spector.  Inference of compact nonlinear dynamic models by epigenetic local search.  Engineering Applications of Artificial Intelligence, 55:292–306, 2016.
  * [11] E.L. Sanjuán, M.I. Parra, and M.M. Pizarro.  Development of models for surface tension of alcohols through symbolic regression.  Journal of Molecular Liquids, 298:111971, 2020.
  * [12] Pascal Neumann, Liwei Cao, Danilo Russo, Vassilios S. Vassiliadis, and Alexei A. Lapkin.  A new formulation for symbolic regression to identify physico-chemical laws from experimental data.  Chemical Engineering Journal, 387:123412, 2020.
  * [13] Yanzhang Li, Hongyu Wang, Yan Li, Huan Ye, Yanan Zhang, Rongzhang Yin, Haoning Jia, Bingxu Hou, Changqiu Wang, Hongrui Ding, Xiangzhi Bai, and Anhuai Lu.  Electron transfer rules of minerals under pressure informed by machine learning.  Nature Communications, 14(1):1815, March 2023.
  * [14] Weihua Cai, Arturo Pacheco-Vega, Mihir Sen, and K.T. Yang.  Heat transfer correlations by symbolic regression.  International Journal of Heat and Mass Transfer, 49(23):4352–4359, 2006.
  * [15] Patrick A. K. Reinbold, Logan M. Kageorge, Michael F. Schatz, and Roman O. Grigoriev.  Robust learning from noisy, incomplete, high-dimensional experimental data via physically constrained symbolic regression.  Nature Communications, 12(1):3219, May 2021.
  * [16] Baicheng Weng, Zhilong Song, Rilong Zhu, Qingyu Yan, Qingde Sun, Corey G. Grice, Yanfa Yan, and Wan-Jian Yin.  Simple descriptor derived from symbolic regression accelerating the discovery of new perovskite catalysts.  Nature Communications, 11(1):3513, July 2020.
  * [17] Silviu-Marian Udrescu and Max Tegmark.  Symbolic pregression: Discovering physical laws from distorted video.  Physical Review E, 103(4):043307, April 2021.
  * [18] Luigi Berardi, Zoran Kapelan, Orazio Giustolisi, and Dragan Savic.  Development of pipe deterioration models for water distribution systems using epr.  10, 07 2008.
  * [19] Arthur Grundner, Tom Beucler, Pierre Gentine, and Veronika Eyring.  Data-Driven Equation Discovery of a Cloud Cover Parameterization.  2023\.
  * [20] Oscar Fajardo-Fontiveros, Ignasi Reichardt, Harry R. De Los Ríos, Jordi Duch, Marta Sales-Pardo, and Roger Guimerà.  Fundamental limits to learning closed-form mathematical models from data.  Nature Communications, 14(1):1043, February 2023.
  * [21] Burr Settles.  Active learning literature survey.  Computer Sciences Technical Report 1648, University of Wisconsin–Madison, 2009.
  * [22] H. S. Seung, M. Opper, and H. Sompolinsky.  Query by committee.  In Proceedings of the Fifth Annual Workshop on Computational Learning Theory, COLT ’92, page 287–294, New York, NY, USA, 1992. Association for Computing Machinery.
  * [23] David E. Goldberg.  Genetic Algorithms in Search, Optimization & Machine Learning.  Addison-Wesley Publishing Company, 1989.
  * [24] Wassim Tenachi, Rodrigo Ibata, and Foivos I. Diakogiannis.  Deep symbolic regression for physics guided by units constraints: toward the automated discovery of physical laws.  2023\.
  * [25] Qiang Lu, Jun Ren, and Zhiguang Wang.  Using Genetic Programming with Prior Formula Knowledge to Solve Symbolic Regression Problem.  Computational Intelligence and Neuroscience, 2016:1–17, 2016.
  * [26] Arijit Chakraborty, Abhishek Sivaram, and Venkat Venkatasubramanian.  Ai-darwin: A first principles-based model discovery engine using machine learning.  Computers & Chemical Engineering, 154:07470, 2021.
  * [27] Charles Fox, Neil Tran, Nikki Nacion, Samiha Sharlin, and Tyler R. Josephson.  Incorporating Background Knowledge in Symbolic Regression using a Computer Algebra System.  2023\.
  * [28] Marissa R. Engle and Nikolaos V. Sahinidis.  Deterministic symbolic regression with derivative information: General methodology and application to equations of state.  AIChE Journal, 68(6), June 2022.
  * [29] Rainier Barrett and Andrew D. White.  Investigating active learning and meta-learning for iterative peptide design.  Journal of Chemical Information and Modeling, 61(1):95–105, 2021\.  PMID: 33350829.
  * [30] Natalie S. Eyke, William H. Green, and Klavs F. Jensen.  Iterative experimental design based on active machine learning reduces the experimental burden associated with reaction screening.  React. Chem. Eng., 5:1963–1972, 2020.
  * [31] Burr Settles.  Active learning literature survey.  Computer Sciences Technical Report 1648, University of Wisconsin–Madison, 2009.
  * [32] Yifan Fu, Xingquan Zhu, and Bin Li.  A survey on instance selection for active learning.  Knowledge and Information Systems, 35(2):249–283, May 2013.
  * [33] Vu-Linh Nguyen, Mohammad Hossein Shaker, and Eyke Hüllermeier.  How to measure uncertainty in uncertainty sampling for active learning.  Machine Learning, 111(1):89–122, January 2022.
  * [34] Nathan Haut, Wolfgang Banzhaf, and Bill Punch.  Active learning improves performance on symbolic regressiontasks in stackgp, 2022.
  * [35] Yoav Freund, H. Sebastian Seung, Eli Shamir, and Naftali Tishby.  Selective Sampling Using Query By Committee.  Machine Learning, 28(2/3):133–168, 1997.
  * [36] Ran Gilad-bachrach, Amir Navot, and Naftali Tishby.  Query by committee made real.  In Y. Weiss, B. Schölkopf, and J. Platt, editors, Advances in Neural Information Processing Systems, volume 18. MIT Press, 2005.
  * [37] Nathan Haut, Wolfgang Banzhaf, and Bill Punch.  Active Learning Improves Performance on Symbolic RegressionTasks in StackGP.  2022\.
  * [38] Josh Bongard and Hod Lipson.  Active coevolutionary learning of deterministic finite automata.  Journal of Machine Learning Research, 6(56):1651–1678, 2005.
  * [39] Teresa Henriques, Luis Antunes, João Bernardes, Mara Matias, Diogo Sato, and Cristina Costa-Santos.  Information-based measure of disagreement for more than two observers: a useful tool to compare the degree of observer disagreement.  BMC Medical Research Methodology, 13(1):47, December 2013.
  * [40] John H. Holland.  Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control, and Artificial Intelligence.  The MIT Press, 1992.
  * [41] K. A. De Jong.  An analysis of the behavior of a class of genetic adaptive systems.  1975\.
  * [42] Melanie Mitchell.  An introduction to genetic algorithms.  1996\.

##  Appendix A Supp. Material

###  A.1 Additional Theory

####  A.1.1 Initialization, Reproduction, and Crossover in Genetic Algorthms

The more members in the population, the broader the search and the better
chance of success, but at the expense of computation time. A random
combination of unary and binary operators (e.g., $\sin$, $\cos$, inv(x), $+$,
$-$, $*$), variables, and constants are formed to get an initial population.
After this, “evolution" can start.

Once the first generation is initialized, reproduction takes place.
Reproduction is the process where the ‘fittest’ members of the population get
passed towards the next generation. To decide which member is better, a metric
called fitness, considering both accuracy and complexity, is used. The fitness
metric can be the Mean Absolute Error or Root Mean Square Error, and
complexity is the number of tree nodes. Although it can be customized as
needed.

After reproduction, the crossover comes next. Crossover, also known as
recombination, is a primary genetic operator in genetic algorithms that
generates new offspring solutions by combining the features (genes) of two
parents [40, 23]. Inspired by the natural process of sexual reproduction,
crossover helps promote genetic diversity and exploration of the search space,
which is crucial for solving complex optimization problems [41]. This allows
searching outside the original population on the mathematical equation space.
This is represented in figure A.1

BeforeAfter+*AB-CD*+EF-GH+*AB-GH*+EF-CD Figure A.1: Crossover depiction
between two members of the population

Despite being predominantly heuristic-based, genetic algorithms benefit from
the “schema theorem" [23], which formalizes the underlying intuition and
provides a theoretical foundation for their effectiveness in optimization and
search tasks. Introduced by David E. Goldberg [23], the schema theorem asserts
that genetic algorithms inherently operate on multiple schemata concurrently.
These schemata, which serve as building blocks or templates, consist of shared
characteristics or patterns found in various solutions within the population.
According to the schema theorem, the evolutionary process intrinsically
rewards advantageous traits (higher fitness) while penalizing undesirable
traits (lower fitness) among these schemata. Through iterative selection,
crossover, and mutation, genetic algorithms progressively converge towards
optimal or near-optimal solutions. Consequently, the schema theorem elucidates
the principles behind genetic algorithms, demonstrating their ability to
efficiently navigate the search space and identify high-quality solutions to
intricate optimization challenges.

####  A.1.2 Mutation: a safeguard

Mutation is another fundamental genetic operator in genetic algorithms that
helps to maintain diversity within the population and prevent premature
convergence to suboptimal solutions [23, 42] It is inspired by the natural
process of mutation in biological organisms, where random changes occur in an
individual’s genetic material. In the context of genetic algorithms, mutation
introduces small random alterations to the candidate solutions, or
chromosomes, which can lead to the exploration of new points in the search
space and thus facilitate the discovery of better solutions.

The mutation step is typically applied after the crossover operation, and its
probability is often much lower than that of crossover [41] There are various
mutation operators, and their choice depends on the problem’s representation
and domain. For example, in a binary-encoded genetic algorithm, a commonly
used mutation operator is bit-flip mutation, which involves flipping the value
of a randomly selected bit (changing 0 to 1 or vice versa) in a chromosome
([40].

BeforeAfter+*AB-CD+*AB-CX Figure A.2: Mutation depiction of a member of the
population

###  A.2 Inclusion of background knowledge

In the following pseudocode, we present an algorithm that incorporates a
penalty term into the loss evaluation process to account for constraint
violations. This concise representation demonstrates how to combine the
objective function and penalty term using a weight parameter, lambda. By
including the penalty term, the algorithm can optimize solutions while
adhering to background knowledge and constraints.

Algorithm 1 Inclusion of a Penalty in Loss Evaluation

1:procedure EvaluateWithPenalty(population, objectiveFunc, penaltyFunc,
$\lambda$)

2: for each individual in population do

3: loss $\leftarrow$ objectiveFunc(individual)

4: penalty $\leftarrow$ penaltyFunc(individual)

5: individual.fitness $\leftarrow$ loss + $\lambda$ penalty

6: end for

7:end procedure

##  Appendix B Additional Results

We provide evidence of our approach with constraints with the subtree counts
on the rediscovery of equation
$\frac{1}{\sqrt{(x_{2}-x_{1})^{2}+(x_{4}-x_{3})^{2}}}$. As can be seen in the
following figure:

Figure B.3: Counts of substructures $(x_{2}-x_{1})$ and $(x_{4}-x_{3})$ for
equation $\frac{1}{\sqrt{(x_{2}-x_{1})^{2}+(x_{4}-x_{3})^{2}}}$. Results from
explorations without rediscovery. It is meant as an indirect measurement if
search is performed in the right part of the search space. Without
constraints, substructures are very homogeneous but in small amounts. With
constraints, substructures seem to appear higher but in a more chaotic way,
and a mix of constrained and unconstrained evolution seemes to outperformed
both methods individually.

We provide tables and plots that illustrate the impact of noise on the
rediscovery of equations and the effectiveness of implementing physical
constraints. Tables 1 and 2 present the results for the rediscovery of
$\frac{1}{2}m(v^{2}+u^{2}+w^{2})$ and
$\frac{1}{\sqrt{(x_{2}-x_{1})^{2}+(y_{2}-y_{1})^{2}}}$, respectively, under
different noise levels and constraint types. Figures 1 and 2 visually
represent these results, providing a clearer comparison between the use of
physical constraints and no constraints in the equation rediscovery process.

Table B.1: Rediscovery of $\frac{1}{2}m(v^{2}+u^{2}+w^{2})$ Type of Experiment | Noiseless | 0.01 Noise | 0.1 Noise
---|---|---|---
No Constraint | 15.61 | 20.1 | 18.33
Symmetry | 17.4 | 19.6 | 21.16
Figure B.4: Illustrating the Effects of Noise on Rediscovery: A Comparison
between Constrained and Unconstrained Scenarios. The figure presents results
for the expression $\frac{1}{2}m(v^{2}+u^{2}+w^{2})$. P-values are 0.6350 for
noiseless conditions, 0.0007 for a noise level of 0.01, and 0.005 for a noise
level of 0.1. The joint p-value for comparing unconstrained and constrained
optimization under the influence of noise is 0.0002.

The results from the third equation, illustrated in Table
LABEL:rediscoverywithnoise and Figure B.5, do not exhibit any enhancement when
compared to unconstrained optimization. In fact, the p-values indicate that a
significant difference exists between the constrained and unconstrained
methods when no constraints are applied favoring the latter, showing the
opposite behaviour. However, when the noise level is increased to 0.1, the
ground truth is only rediscovered in the constrained scenario.

Table B.2: Rediscovery of $\frac{1}{\sqrt{(x_{2}-x_{1})^{2}+(y_{2}-y_{1})^{2}}}$ Type of Experiment | Noiseless | 0.01 Noise | 0.1 Noise
---|---|---|---
No Constraint | 11 | 18 | No Rediscovery
Symmetry | 16 | 23 | Rediscovery
Divergency | 16 | 23 | Rediscovery
Figure B.5: Illustrating the Effects of Noise on Rediscovery: A Comparison
between Constrained and Unconstrained Scenarios. The figure presents results
for the expression $\frac{1}{\sqrt{(x_{2}-x_{1})^{2}+(y_{2}-y_{1})^{2}}}$.
P-values for Symmetry vs No Constraint are 0.0266 for a noise level of 0.01
and 1.22e-11 for noiseless conditions; for Div vs No Constraint, p-values are
0.3319 for a noise level of 0.01 and 3.98e-10 for noiseless conditions. The
joint p-value for comparing constrained and unconstrained optimization under
the influence of noise is 1.04e-18.


EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

 CERN-EP-2022-258 LHCb-PAPER-2022-040 12 January 2023

Evidence of a ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{K}^{0}_{\mathrm{S}}}$
structure in
${{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}$
decays

LHCb collaboration

An amplitude analysis of
${{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}$
decays is performed using proton-proton collision data corresponding to an
integrated luminosity of 9$\text{\,fb}^{-1}$, collected with the LHCb detector
at centre-of-mass energies of 7, 8 and 13 TeV. Evidence with a significance of
$4.0$ standard deviations of a structure in the
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{K}^{0}_{\mathrm{S}}}$ system, named
$T_{\psi s1}^{\theta}(4000)^{0}$, is seen, with its mass and width measured to
be $3991\,_{-10}^{+12}\,{}^{+\phantom{0}9}_{-17}\text{\,Me\kern-1.00006ptV}$
and $105\,_{-25}^{+29}\,{}^{+17}_{-23}\text{\,Me\kern-1.00006ptV}$,
respectively, where the first uncertainty is statistical and the second
systematic. The $T_{\psi s1}^{\theta}(4000)^{0}$ state is likely to be the
isospin partner of the $T_{\psi s1}^{\theta}(4000)^{+}$ state, previously
observed in the ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{K}^{+}}$ system of the
${{{B}^{+}}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{+}}$ decay.
When isospin symmetry for the charged and neutral $T_{\psi s1}^{\theta}(4000)$
states, is assumed, the signal significance increases to $5.4$ standard
deviations.

Submitted to Phys. Rev. Lett.

© 2024 CERN for the benefit of the LHCb collaboration. CC BY 4.0 licence.

Hadrons formed by more than three quarks, namely exotic hadrons, can have
complex colour and flavour structures, and studies of their internal dynamics
shed light on the non-perturbative behaviour of quantum chromodynamics at low
energy. Since the discovery of the $\chi_{c1}(3872)$ state [1], several states
compatible with a four or five quark composition have been observed [2],
including the fully charmed tetraquark states [3, 4, 5] and the pentaquark
states [6, 7]. In 2020, the BESIII collaboration reported the observation of
the $T_{\psi s}(3985)^{+}$ state in the ${{D}^{+}_{s}}{{\kern
1.79993pt\overline{\kern-1.79993ptD}}{}^{*0}}$ and ${{D}^{*+}_{s}}{{\kern
1.79993pt\overline{\kern-1.79993ptD}}{}^{0}}$ mass spectra [8].111Charge
conjugation is implied throughout the Letter unless specified otherwise. Two
similar states, $T_{\psi s1}^{\theta}(4000)^{+}$ and $T_{\psi s1}(4220)^{+}$,
were observed in
${{{B}^{+}}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{+}}$
decays by the LHCb experiment [9].222The new exotic hadron naming convention
proposed by the LHCb collaboration [10] is applied throughout the Letter. The
$T_{\psi s}(3985)^{+}$, $T_{\psi s1}^{\theta}(4000)^{+}$, and $T_{\psi
s1}(4220)^{+}$ states are named $Z_{cs}(3985)^{+}$, $Z_{cs}(4000)^{+}$, and
$Z_{cs}(4220)^{+}$ in the original papers [8, 9], respectively. The
superscript $\theta$ refers to 1/2 isospin and positive parity. With a minimum
quark content of ${c}{\overline{c}}{u}{\overline{s}}$, these states are
explicitly exotic. The observations of the $T_{\psi s1}^{\theta}(4000)^{+}$,
$T_{\psi s1}(4220)^{+}$, and $T_{\psi s}(3985)^{+}$ states have stimulated
many theoretical studies to interpret their internal structures, such as
hadronic molecules [11] or compact tetraquarks [12] or threshold effects [13].

Searching for the isospin partners of the $T_{\psi s}^{+}$ states, the
$T_{\psi s}^{0}$ states, provides another opportunity to better understand
hidden-charm tetraquarks with strangeness. The quark contents for the $T_{\psi
s}^{+}$ and $T_{\psi s}^{0}$ states are ${c}{\overline{c}}{u}{\overline{s}}$
and ${c}{\overline{c}}{d}{\overline{s}}$, respectively. Recently, the BESIII
experiment reported evidence for the $T_{\psi s}(3985)^{0}$ state
(${c}{\overline{c}}{d}{\overline{s}}$) in the ${{D}^{+}_{s}}{{D}^{*-}}$ and
${{D}^{*+}_{s}}{{D}^{-}}$ mass spectra [14]. The
${{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}$
decay is an ideal candidate to search for such $T_{\psi s}^{0}$ states. The
quark level processes of the
${{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}$
decay are similar to those of the
${{{B}^{+}}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{+}}$ decay,
as shown in Fig. 1, and the two decays are related by isospin symmetry.

Figure 1: One of the Feynman diagrams contributing to
${{{B}^{+}}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{+}}$ or
${{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}$
decays.

In this Letter, evidence of a
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{K}^{0}_{\mathrm{S}}}$ structure is
reported from an amplitude analysis of
${{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}$
decays. The analysis is based on proton-proton ($pp$) collision data,
collected with the LHCb detector at centre-of-mass energies of
$7\text{\,Te\kern-1.00006ptV}$, $8\text{\,Te\kern-1.00006ptV}$, and
$13\text{\,Te\kern-1.00006ptV}$, corresponding to integrated luminosities of
1$\text{\,fb}^{-1}$, 2$\text{\,fb}^{-1}$, and 6$\text{\,fb}^{-1}$,
respectively.

The LHCb detector is a single-arm forward spectrometer covering the
pseudorapidity range $2<\eta<5$, described in detail in Refs. [15, 16].
Simulated samples of the signal decays are used to estimate the effect of
reconstruction and event selections. The samples are generated with Pythia
[17, *Sjostrand:2006za], EvtGen [19], and the Geant4 toolkit [20,
*Agostinelli:2002hh] as described in Ref. [22].

For the
${{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}$
decays, the ${J\mskip-3.0mu/\mskip-2.0mu\psi}$, $\phi$, and
${K}^{0}_{\mathrm{S}}$ candidates are reconstructed by combining two
oppositely charged muons, kaons, and pions, respectively, forming vertices
detached from any primary vertex (PV). The ${B}^{0}$ candidate is
reconstructed by combining the ${J\mskip-3.0mu/\mskip-2.0mu\psi}$, $\phi$, and
${K}^{0}_{\mathrm{S}}$ candidates and the resulting decay vertex must be of
high quality. In order to improve the resolution of the mass of the
${{B}^{0}}$ candidates (referred to as
$m_{{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}}$), a
kinematic fit [23] is applied, with the ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}$
and ${{K}^{0}_{\mathrm{S}}}$ masses constrained to their known values [2] and
the ${{B}^{0}}$ candidates constrained to originate from their associated PVs.
The associated PV is the one with smallest impact parameter $\chi^{2}$,
defined as the difference in the vertex-fit $\chi^{2}$ of a given PV with and
without the particle under consideration.

The ${{B}^{0}}$ candidates are selected by a multivariate classifier based on
a multilayer perceptron [24, 25]. The selection is similar to that used in
Ref. [9], with an additional requirement on the significance of the flight
distance of the ${{K}^{0}_{\mathrm{S}}}$ from the ${{B}^{0}}$ decay vertex.
Only candidates with the mass of the ${{K}^{+}}{{K}^{-}}$ system within $\pm
15\text{\,Me\kern-1.00006ptV}$ around the known $\phi$ mass [2] are
kept.333Natural units with $\hbar=c=1$ are used throughout the Letter. The
multivariate classifier uses eleven variables related to the decay-chain
topology, particle transverse momentum, vertex fit quality, and charged
particle identification information. The selection criterion on the classifier
response is chosen by maximising the figure of merit,
$\mathcal{S}^{2}/(\mathcal{S}+\mathcal{B})^{3/2}$ [26], where $\mathcal{S}$
and $\mathcal{B}$ are the yields of signal and background in the signal
region, respectively. The signal region in
$m_{{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}}$ is defined
to be $\pm 15\text{\,Me\kern-1.00006ptV}$ around the known ${{B}^{0}}$ mass
[2].

An extended maximum-likelihood fit is performed to the
$m_{{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}}$
distribution of the selected candidates, as is shown in Fig. 2a. The
${{B}^{0}}$ signal is described by a Hypatia function [27] and the background
is described by an exponential function. The ${{B}^{0}}$ yield is measured to
be $1866\pm 47$ in the signal region. The fraction of combinatorial background
in the signal region is $6\,\%$. Roughly $4\,\%$ of the ${B}^{0}$ yield
corresponds to peaking-background
${{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{K}^{+}}{{K}^{-}}{{K}^{0}_{\mathrm{S}}}$
decays without an intermediate $\phi$ meson (referred to as non-$\phi$). The
non-$\phi$ contribution is neglected in the default amplitude model and is
considered as a source of systematic uncertainty. A further kinematic fit is
performed to improve the momentum resolution of the final-state particles by
constraining the measured ${{B}^{0}}$ mass to its known value. Figure. 2b
shows the distribution of
$m^{2}_{{{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{K}^{0}_{\mathrm{S}}}}$ versus
$m^{2}_{{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi}$ for candidates in the signal
region.

(a)(b)

Figure 2: (a) Invariant-mass distribution of selected ${B}^{0}$ candidates and
corresponding fit result. (b) The distribution of
$m^{2}_{{{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{K}^{0}_{\mathrm{S}}}}$ versus
$m^{2}_{{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi}$ for candidates in the $\pm$15
MeV region around the known ${{B}^{0}}$ mass.

A simultaneous amplitude fit is performed to the
${{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}$
and ${{{B}^{+}}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{+}}$
samples. The
${{{B}^{+}}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{+}}$ sample
is the same as that used for the observation of the $T_{\psi
s1}^{\theta}(4000)^{+}$ and $T_{\psi s1}(4220)^{+}$ states [9]. The likelihood
for the ${{B}^{0}}$ decay is given by

$\mathcal{L}(\vec{\omega})=\prod_{i}\left[(1-\beta)\,\mathcal{P}_{\rm
sig}(m_{\phi{{K}^{0}_{\mathrm{S}}}}^{i},\vec{\Omega}^{i}|\vec{\omega})+\beta\,\mathcal{P}_{\rm
bkg}(m_{\phi{{K}^{0}_{\mathrm{S}}}}^{i},\vec{\Omega}^{i})\right],$ (1)

where the probability density functions (PDFs) for the signal and background
components are given by $\mathcal{P}_{\rm
sig}(m_{\phi{{K}^{0}_{\mathrm{S}}}}^{i},\vec{\Omega}^{i}|\vec{\omega})$ and
$\mathcal{P}_{\rm bkg}(m_{\phi{{K}^{0}_{\mathrm{S}}}}^{i},\vec{\Omega}^{i})$,
respectively. The superscript $i$ refers to the $i$-th candidate. The fraction
of combinatorial background, $\beta$, is fixed to $6\,\%$. The likelihood of
the ${{{B}^{+}}}$ decay is similar to Eq. 1 and is described in detail in Ref.
[28]. The decay kinematics are described by a mass, chosen to be
$m_{\phi{{K}^{0}_{\mathrm{S}}}}$, and five angular variables $\vec{\Omega}$ as
defined below. The signal PDF is proportional to the incoherent sum of the
squared matrix elements for different final-state muon helicities
$(\lambda_{\mu^{+}},\lambda_{\mu^{-}})$, and depends on the set of fit
parameters, $\vec{\omega}$, which includes masses, widths, and helicity
couplings of intermediate states contributing to
${{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}$
decays,

$\mathcal{P}_{\rm
sig}(m_{\phi{{K}^{0}_{\mathrm{S}}}}^{i},\vec{\Omega}^{i}|\vec{\omega})=\frac{1}{I(\vec{\omega})}\sum_{\lambda_{\mu^{+}}\lambda_{\mu^{-}}}\sum_{f}{\big{|}\mathcal{M}_{f}(m_{\phi{{K}^{0}_{\mathrm{S}}}}^{i},\vec{\Omega}^{i}|\vec{\omega})\big{|}^{2}}\Phi(m_{\phi{{K}^{0}_{\mathrm{S}}}}^{i},\vec{\Omega}^{i})\epsilon(m_{\phi{{K}^{0}_{\mathrm{S}}}}^{i},\vec{\Omega}^{i}).$
(2)

The term $\Phi(m_{\phi{{K}^{0}_{\mathrm{S}}}}^{i},\vec{\Omega}^{i})$
represents the phase-space density, the term
$\epsilon(m_{\phi{{K}^{0}_{\mathrm{S}}}}^{i},\vec{\Omega}^{i})$ represents the
efficiency, and the term $I(\vec{\omega})$ is a normalisation factor. The
efficiency is obtained from simulated samples. Since the
${{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}$
decay is a self-charge-conjugated mode, the data sample contains both
${{B}^{0}}$ and ${{\kern 1.79993pt\overline{\kern-1.79993ptB}}{}^{0}}$ decays
with unknown ${b}$ flavour. Therefore, the signal PDF for the neutral $B$
decay averages ${{B}^{0}}$ and ${{\kern
1.79993pt\overline{\kern-1.79993ptB}}{}^{0}}$ contributions, where the flavour
of $B$ candidates is denoted as $f$ in Eq. 2.

The matrix element
$\mathcal{M}_{f}(m_{\phi{{K}^{0}_{\mathrm{S}}}}^{i},\vec{\Omega}^{i}|\vec{\omega})$
is constructed based on the helicity formalism [29, *Richman:1984gh,
*Jacob:1959at]. Three interfering decay sequences are considered:
${{B}^{0}}\\!\rightarrow(K^{*}\rightarrow\phi{{K}^{0}_{\mathrm{S}}}){{J\mskip-3.0mu/\mskip-2.0mu\psi}}$,
${{B}^{0}}\\!\rightarrow(X\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi){{K}^{0}_{\mathrm{S}}}$,
and ${{B}^{0}}\\!\rightarrow(T_{\psi
s1}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{K}^{0}_{\mathrm{S}}})\phi$.444The
${{K}^{*}}$ denotes any excited $K$ state. The $X$ denotes $\chi_{c0,1}$,
$\eta_{c2}$, and $T_{\psi\phi 1}^{\eta}$ states. The decay sequences can be
described by $m_{\phi{{K}^{0}_{\mathrm{S}}}}$ and
$\vec{\Omega}\equiv(\theta_{K^{*}},\theta_{{{J\mskip-3.0mu/\mskip-2.0mu\psi}}},\theta_{\phi},\Delta\varphi_{K^{*}{{J\mskip-3.0mu/\mskip-2.0mu\psi}}},\Delta\varphi_{K^{*}\phi})$,
where the $\theta_{{{K}^{*}}}$, $\theta_{{{J\mskip-3.0mu/\mskip-2.0mu\psi}}}$,
and $\theta_{\phi}$ are the helicity angles of the $\phi$, $\mu$, and
${{K}^{+}}$ particles in the ${{K}^{*}}$,
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}$, and $\phi$ rest frames, respectively.
The angles between the $K^{*}$ and ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}$ decay
planes, and between the ${{K}^{*}}$ and $\phi$ decay planes in the $B$ rest
frame, are denoted as
$\Delta\varphi_{K^{*}{{J\mskip-3.0mu/\mskip-2.0mu\psi}}}$, and
$\Delta\varphi_{K^{*}\phi}$, respectively. The background PDF
$\mathcal{P}_{\rm bkg}(m_{\phi{{K}^{0}_{\mathrm{S}}}}^{i},\vec{\Omega}^{i})$
is determined from studying the candidates in the sideband of the
$m_{{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}}$
distribution with the methods outlined in Ref. [28].

The default model is taken from Ref. [9] and is composed of nine ${{K}^{*}}$
components, seven $X$ components, and two $T_{\psi s1}$ states ($T_{\psi
s1}^{\theta}(4000)$ and $T_{\psi s1}(4220)$), and one
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi$ nonresonant component (NR).
Relativistic Breit–Wigner functions are used to model the lineshapes of these
resonances. As shown in Fig. 1, the
${{{B}^{+}}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{+}}$ and
${{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}$
decays are mediated by the same ${b}\rightarrow{c}{\overline{c}}{s}$ process
and they differ only in the spectator quark, either ${u}$ or ${d}$, which
hadronises to the ${{K}^{+}}$ or ${{K}^{0}_{\mathrm{S}}}$ meson. The
amplitudes of the ${{B}^{0}}$ and the ${{{B}^{+}}}$ decays are formed by
intermediate states that are either identical (namely $X$) or isospin partners
(${{K}^{*}}$ or $T_{\psi s1}$). Under the assumption of isospin symmetry, the
mass, width, and helicity couplings for all the components except for the
$T_{\psi s1}^{\theta}(4000)$ states are constrained to be identical. In order
to test the significance of the $T_{\psi s1}^{\theta}(4000)^{0}$ signal
without prior knowledge of its properties, a fit is performed with the masses,
widths, and helicity couplings of $T_{\psi s1}^{\theta}(4000)^{0}$ and
$T_{\psi s1}^{\theta}(4000)^{+}$ states allowed to vary independently. In an
alternative model, isospin symmetry is imposed for the $T_{\psi
s1}^{\theta}(4000)^{0}$ and $T_{\psi s1}^{\theta}(4000)^{+}$ components. The
parameters for the $T_{\psi s1}(4220)^{0}$ state are always constrained to be
identical to the $T_{\psi s1}(4220)^{+}$ state due to the limited size of the
${{B}^{0}}$ sample.

Figure 3:  Distributions of (left) $m_{\phi K}$, (middle)
$m_{{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi}$, and (right)
$m_{{{J\mskip-3.0mu/\mskip-2.0mu\psi}}K}$, overlaid with the corresponding
projections of the default fit model. The upper and lower rows correspond to
the ${{{B}^{+}}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{+}}$
and
${{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}$
decays, respectively.

Figure 3 shows the $\phi K$, ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi$, and
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}K$ mass distributions and the corresponding
fit projections of the default model for the two $B$ decay modes. Table 1
summarises measurements of the mass, width, fit fraction of the $T_{\psi
s1}^{\theta}(4000)^{0}$ state, and the mass difference between the $T_{\psi
s1}^{\theta}(4000)^{0}$ and $T_{\psi s1}^{\theta}(4000)^{+}$ states, defined
as $\Delta M\equiv M_{T_{\psi s1}^{\theta}(4000)^{0}}-M_{T_{\psi
s1}^{\theta}(4000)^{+}}$. The fit value of $\Delta M$ is zero within
uncertainties, consistent with the two states being isospin partners. The fit
fraction of each component is defined as the integral of the signal PDF
divided by the $I(\vec{\omega})$ term. All the fit parameters, including mass,
width, and helicity couplings of the intermediate states, of the default model
in this analysis are consistent with the corresponding parameters of the
default model for the
${{{B}^{+}}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{+}}$ decay
in Ref. [9].

Table 1: Results for the $T_{\psi s1}^{\theta}(4000)^{0}$ state from the default model. The first uncertainty is statistical and the second systematic. State | Mass (MeV) | Width (MeV) | Fit fraction (%) |  $\Delta M$ (MeV)
---|---|---|---|---
$T_{\psi s1}^{\theta}(4000)^{0}$ | $3991\,_{-10}^{+12}\,{}^{\phantom{0}+9}_{-17}$ | $105\,_{-25}^{+29}\,{}^{+17}_{-23}$ | $7.9\pm{2.5}\,^{+3.0}_{-2.8}$ | $-12\,_{-10}^{+11}\,{}^{+6}_{-4}$
Table 2: Systematic uncertainties associated to the mass, width, fit fraction of the $T_{\psi s1}^{\theta}(4000)^{0}$ state, and the mass difference between the $T_{\psi s1}^{\theta}(4000)^{0}$ and $T_{\psi s1}^{\theta}(4000)^{+}$ states. Source |  Mass ( MeV) |  Width ( MeV) | Fit fraction (%) |  $\Delta M$ ( MeV)
---|---|---|---|---
Efficiency and background models | $+10.46$ | $+\phantom{0}1.91$ | $-0.06$ | $+1.02$
Hadron radius | ${}_{-\phantom{0}1.92}^{+\phantom{0}0.29}$ | ${}_{-\phantom{0}5.85}^{+\phantom{0}2.33}$ | ${}_{-0.39}^{+0.05}$ | ${}_{-1.68}^{+0.92}$
$X(4140)$ Flatté model | $+\phantom{0}0.61$ | $-\phantom{0}3.47$ | $-0.31$ | $-0.18$
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi$ NR representation | $+\phantom{0}1.72$ | $+\phantom{0}1.46$ | $+0.38$ | $+0.12$
Simplified K-Matrix model | $-\phantom{0}7.30$ | $-17.54$ | $-1.73$ | $-3.38$
Widths of the Breit–Wigner | $+\phantom{0}1.90$ | $-\phantom{0}3.41$ | $-0.56$ | $+1.70$
Spin-parity of the $T_{\psi s1}(4220)$ state | $-16.32$ | $-17.76$ | $-2.33$ | $-1.48$
Non-$\phi$ contribution | $+\phantom{0}5.52$ | $+\phantom{0}7.98$ | $+0.92$ | $+5.49$
Additional $1^{+}$ ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi$ NR | $+\phantom{0}0.80$ | $-\phantom{0}3.66$ | $+0.18$ | $-1.28$
Additional $2^{+}$ ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi$ NR | $+\phantom{0}7.66$ | $+\phantom{0}2.57$ | $+2.52$ | $+0.92$
Extended model | $-\phantom{0}6.36$ | $-\phantom{0}2.54$ | $-0.20$ | $-2.98$
Additional $X$ state | $-\phantom{0}5.58$ | $+\phantom{0}0.27$ | $+0.34$ | $-1.48$
Fit fraction constraint | $+\phantom{0}0.07$ | $+\phantom{0}1.94$ | $-0.48$ | $-0.45$
Impact parameter $\chi^{2}$ modelling | $+\phantom{0}1.00$ | $-\phantom{0}2.49$ | $-0.22$ | $+0.72$
Finite simulation sample size | $+\phantom{0}2.90$ | $+11.20$ | $+1.33$ | $+2.17$
Fixed mass and width of ${{K}^{*}}$ states | $-\phantom{0}0.58$ | $-\phantom{0}1.30$ | $-0.43$ | $-0.58$
Samples to construct the $\mathcal{P}_{\rm bkg}$ | ${}_{-\phantom{0}1.38}^{+\phantom{0}1.06}$ | ${}_{-\phantom{0}7.25}^{+\phantom{0}6.38}$ | ${}_{-0.59}^{+0.54}$ | ${}_{-0.68}^{+0.52}$
Background fraction uncertainty | ${}_{-\phantom{0}0.45}^{+\phantom{0}0.37}$ | ${}_{-\phantom{0}1.55}^{+\phantom{0}1.38}$ | ${}_{-0.07}^{+0.06}$ | ${}_{-0.38}^{+0.22}$
Final | ${}_{-16.71}^{+\phantom{0}8.46}$ | ${}_{-23.33}^{+17.09}$ | ${}_{-2.84}^{+3.00}$ | ${}_{-4.23}^{+6.04}$

The estimated systematic uncertainties on the mass, width, fit fraction of the
$T_{\psi s1}^{\theta}(4000)^{0}$ state, and on $\Delta M$ are summarised in
Table 2. Both the background PDF and efficiency function are described by an
expansion with Legendre polynomials and a spherical harmonic function instead
of interpolation. The effective hadron radius in the Blatt–Weisskopf barrier
factor [32], equal to $3\text{\,Ge\kern-1.00006ptV}^{-1}$ in the default
model, is replaced with two alternatives, 1.5 and 4.5 GeV-1. The Flatté model
[33] including ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi$ and
${{D}^{*+}_{s}}{{D}^{-}_{s}}$ channels is used for the lineshape of the
$X(4140)$ state instead of the relativistic Breit–Wigner model. The
representation of the ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi$ NR contribution
is changed from a constant to a linear function. The ${{K}^{*}}$ states are
described by a simplified one-channel K-matrix model [2] instead of a sum of
relativistic Breit–Wigner functions. The mass-dependent widths in the
relativistic Breit–Wigner function for the ${{K}^{*}}$ resonances are
calculated using the decay with the largest branching fraction,
${{K}^{*}}\rightarrow K\pi$ or ${{K}^{*}}\rightarrow{{K}^{*}}(892)\pi$,
instead of using the ${{K}^{*}}\rightarrow\phi K$ decay. The spin-parity of
the $T_{\psi s1}(4220)$ state is changed from $1^{+}$ to $1^{-}$. The
uncertainty originating from the neglected non-$\phi$ contribution is
evaluated by narrowing the $\phi$ mass window from $\pm
15\text{\,Me\kern-1.00006ptV}$ to $\pm 8\text{\,Me\kern-1.00006ptV}$. An
additional ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi$ NR contribution with the
spin-parity equal to $1^{+}$ or $2^{+}$ is included. An extended model with
more ${{K}^{*}}$ states is also studied, which includes all the ${{K}^{*}}$
resonances that are within the allowed phase space, as predicted in Ref. [34].
Possible additional $X$ states, with spins ranging from 0 to 2, are checked in
the extended model. The total fit fraction of the default model is 165.2 %,
which indicates that the interference between the decay amplitudes is large.
An alternative fit constraining the total fit fraction to be smaller than 140
%, corresponding to a reduction by three times its uncertainty, is performed.
The final positive (negative) systematic uncertainty is taken as the maximal
positive (negative) deviation from the model uncertainties above summed in
quadrature with the uncertainty from other sources. These other sources
include mismodelling of the impact parameter $\chi^{2}$ of $B$ candidates, the
finite size of the simulated samples, the fixed masses and widths of the known
${{K}^{*}}$ states, the choice of the data sample used to construct the
background PDF model, and the uncertainty on the background fraction.

The significance of the $T_{\psi s1}^{\theta}(4000)^{0}$ state is evaluated
with a likelihood ratio test. The test statistic is defined as
$t\equiv-2\ln[\mathcal{L}(H_{0})/\mathcal{L}(H_{1})]$ with $H_{1}$ and $H_{0}$
denoting the default model with and without the $T_{\psi
s1}^{\theta}(4000)^{0}$ state. Here, the term $\mathcal{L}$ represents the
likelihood of the ${{B}^{0}}$ decay. Five thousand pseudoexperiments are
generated with the $H_{0}$ model, where the parameters of the model are fixed
to the values determined from data. The sample size of each pseudoexperiment
follows a Poisson distribution with its mean being equal to the number of
candidates in the data sample. A fit is performed to each pseudoexperiment
with the $H_{0}$ and $H_{1}$ models to determine the test statistic $t$. A
$\chi^{2}$ function with the number of degrees-of-freedom allowed to vary is
used to fit the $t$ distribution from these pseudoexperiments. Using the
resulting $\chi^{2}$ function, the probability of $t>t_{\rm data}$ is taken as
the $p$-value. The corresponding significance is estimated to be $4.9\sigma$,
decreasing to $4.0\sigma$ after including systematic uncertainties, which
provides evidence for the $T_{\psi s1}^{\theta}(4000)^{0}$ state. The
significance for the $T_{\psi s1}^{\theta}(4000)^{0}$ state is taken as the
smallest significance found when varying the sources of systematic
uncertainty.

The significance is also evaluated assuming isospin symmetry for the $T_{\psi
s1}^{\theta}(4000)$ states, where the test statistic is defined as
$t^{\prime}\equiv-2\ln[\mathcal{L}(H_{0})/\mathcal{L}(H_{1}^{\prime})]$ with
$\mathcal{L}$ denoting the total likelihood of the ${{{B}^{+}}}$ and
${{B}^{0}}$ decays. The corresponding fit parameters for the $T_{\psi
s1}^{\theta}(4000)^{+}$ and $T_{\psi s1}^{\theta}(4000)^{0}$ states, including
mass, width, and helicity couplings, are constrained to be equal between the
two states in the $H_{1}^{\prime}$ model. A fit is performed to the
$t^{\prime}$ distribution from the pseudoexperiments with a Gaussian function.
The Gaussian function rather than a $\chi^{2}$ distribution is used here
because the number of degree-of-freedoms of the $H_{0}$ and $H_{1}^{\prime}$
models are equal. The significance is estimated to be $7.2\sigma$ and
decreases to $5.4\sigma$ after accounting for systematic uncertainties.

In conclusion, an amplitude analysis of the
${{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{0}_{\mathrm{S}}}$
decay is performed. Evidence of a
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{K}^{0}_{\mathrm{S}}}$ structure, denoted
the $T_{\psi s1}^{\theta}(4000)^{0}$ state, is obtained with a significance of
$4.0\sigma$. The mass and width of this state are measured to be

$\displaystyle M(T_{\psi s1}^{\theta}(4000)^{0})$
$\displaystyle=3991\,_{-10}^{+12}\,{}^{+\phantom{0}9}_{-17}\text{\,Me\kern-1.00006ptV},$
$\displaystyle\Gamma(T_{\psi s1}^{\theta}(4000)^{0})$
$\displaystyle=\phantom{0}105\,_{-25}^{+29}\,{}^{+17}_{-23}\text{\,Me\kern-1.00006ptV},$

where the first uncertainty is statistical and the second systematic. The mass
difference between the $T_{\psi s1}^{\theta}(4000)^{0}$ and $T_{\psi
s1}^{\theta}(4000)^{+}$ states is measured to be

$\displaystyle\Delta
M=-12\,_{-10}^{+11}\,{}^{+6}_{-4}\text{\,Me\kern-1.00006ptV},$

which is consistent with the two states being isospin partners. With isospin
symmetry imposed for the $T_{\psi s1}^{\theta}(4000)$ states, the significance
of the $T_{\psi s1}^{\theta}(4000)^{0}$ structure is measured to be
$5.4\sigma$.

## References

  * [1] Belle collaboration, S. K. Choi et al., _Observation of a narrow charmonium-like state in exclusive $B^{\pm}\rightarrow K^{\pm}\pi^{+}\pi^{-}J/\psi$ decays_, Phys. Rev. Lett. 91 (2003) 262001, arXiv:hep-ex/0309032
  * [2] Particle Data Group, R. L. Workman et al., _Review of particle physics_ , Prog. Theor. Exp. Phys. 2022 (2022) 083C01
  * [3] LHCb collaboration, R. Aaij et al., _Observation of structure in the ${J\mskip-3.0mu/\mskip-2.0mu\psi}$-pair mass spectrum_, Science Bulletin 65 (2020) 1983, arXiv:2006.16957
  * [4] ATLAS collaboration, _Observation of an excess of di-charmonium events in the four-muon final state with the ATLAS detector ,_ ATLAS-CONF-2022-040
  * [5] CMS collaboration, _Observation of new structures in the ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{J\mskip-3.0mu/\mskip-2.0mu\psi}}$ mass spectrum in $pp$ collisions at $\sqrt{s}=13$ TeV,_ CMS-PAS-BPH-21-003
  * [6] LHCb collaboration, R. Aaij et al., _Observation of a narrow pentaquark state, $P_{c}(4312)^{+}$, and of two-peak structure of the $P_{c}(4450)^{+}$_, Phys. Rev. Lett. 122 (2019) 222001, arXiv:1904.03947
  * [7] LHCb collaboration, R. Aaij et al., _Observation of a ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}\Lambda$ resonance consistent with a strange pentaquark candidate in ${{{B}^{-}}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\Lambda{\overline{p}}$ decays_, arXiv:2210.10346, submitted to Phys. Rev. Lett.
  * [8] BESIII collaboration, M. Ablikim et al., _Observation of a near-threshold structure in the $K^{+}$ recoil-mass spectra in $e^{+}e^{-}\rightarrow K^{+}(D_{s}^{-}D^{*0}+D_{s}^{*-}D^{0})$_, Phys. Rev. Lett. 126 (2021) 102001, arXiv:2011.07855
  * [9] LHCb collaboration, R. Aaij et al., _Observation of new resonances decaying to ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}K^{+}$ and ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi$ _, Phys. Rev. Lett. 127 (2021) 082001, arXiv:2103.01803
  * [10] LHCb collaboration, T. Gershon et al., _Exotic hadron naming convention_ , arXiv:2206.15233
  * [11] L. Meng, B. Wang, and S.-L. Zhu, _$Z_{cs}(3985)^{-}$ as the $U$-spin partner of $Z_{c}(3900)^{-}$ and implication of other states in the $\text{SU(3)}_{F}$ symmetry and heavy quark symmetry_, Phys. Rev. D102 (2020) 111502, arXiv:2011.08656
  * [12] L. Maiani, A. D. Polosa, and V. Riquer, _The new resonances $Z_{cs}(3985)$ and $Z_{cs}(4003)$ (almost) fill two tetraquark nonets of broken $SU(3)_{f}$_, Science Bulletin 66 (2021) 1616, arXiv:2103.08331
  * [13] N. Ikeno, R. Molina, and E. Oset, _$Z_{cs}$ states from the $D^{*}_{s}{{\kern 1.79993pt\overline{\kern-1.79993ptD}}{}^{*}}$ and ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{K}^{*}}$ coupled channels: Signal in ${{{B}^{+}}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{+}}$ decay_, Phys. Rev. D 105 (2022) 014012, arXiv:2111.05024, [Erratum: Phys.Rev.D 106, 099905 (2022)]
  * [14] BESIII collaboration, M. Ablikim et al., _Evidence for a neutral near-threshold structure in the $K^{0}_{\rm S}$ recoil-mass spectra in $e^{+}e^{-}\rightarrow K^{0}_{\rm S}D_{s}^{+}D^{*-}$ and $e^{+}e^{-}\rightarrow K^{0}_{\rm S}D_{s}^{*+}D^{-}$_, Phys. Rev. Lett. 129 (2022) 112003, arXiv:2204.13703
  * [15] LHCb collaboration, A. A. Alves Jr. et al., _The LHCb detector at the LHC_, JINST 3 (2008) S08005
  * [16] LHCb collaboration, R. Aaij et al., _LHCb detector performance_ , Int. J. Mod. Phys. A30 (2015) 1530022, arXiv:1412.6352
  * [17] T. Sjöstrand, S. Mrenna, and P. Skands, _A brief introduction to PYTHIA 8.1_ , Comput. Phys. Commun. 178 (2008) 852, arXiv:0710.3820
  * [18] T. Sjöstrand, S. Mrenna, and P. Skands, _PYTHIA 6.4 physics and manual_ , JHEP 05 (2006) 026, arXiv:hep-ph/0603175
  * [19] D. J. Lange, _The EvtGen particle decay simulation package_ , Nucl. Instrum. Meth. A462 (2001) 152
  * [20] Geant4 collaboration, J. Allison et al., _Geant4 developments and applications_ , IEEE Trans. Nucl. Sci. 53 (2006) 270
  * [21] Geant4 collaboration, S. Agostinelli et al., _Geant4: A simulation toolkit_ , Nucl. Instrum. Meth. A506 (2003) 250
  * [22] M. Clemencic et al., _The LHCb simulation application, Gauss: Design, evolution and experience_, J. Phys. Conf. Ser. 331 (2011) 032023
  * [23] W. D. Hulsbergen, _Decay chain fitting with a Kalman filter_ , Nucl. Instrum. Meth. A: 552 (2005) 566
  * [24] W. S. McCulloch and W. Pitts, _A logical calculus of the ideas immanent in nervous activity_ , Bull. Math. Biol 5 (1943) 115
  * [25] F. Rosenblatt, _The perceptron: A probabilistic model for information storage and organization in the brain_ , Psy. Rev. 65 (1958) 386–408
  * [26] LHCb collaboration, R. Aaij et al., _Amplitude analysis of ${{B}^{0}_{s}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{\pm}}{{\pi}^{\mp}}$ decays_, JHEP 06 (2019) 114, arXiv:1902.07955
  * [27] D. Martínez Santos and F. Dupertuis, _Mass distributions marginalized over per-event errors_ , Nucl. Instrum. Meth. A764 (2014) 150, arXiv:1312.5000
  * [28] LHCb collaboration, R. Aaij et al., _Amplitude analysis of ${{{B}^{+}}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi{{K}^{+}}$ decays_, Phys. Rev. D95 (2017) 012002, arXiv:1606.07898
  * [29] S. U. Chung, _Spin formalisms_ , CERN-71-08 (1971)
  * [30] J. D. Richman, _An experimenter’s guide to the helicity formalism_ ,  CALT-68-1148
  * [31] M. Jacob and G. C. Wick, _On the general theory of collisions for particles with spin_ , Annals Phys. 7 (1959) 404
  * [32] J. M. Blatt and V. F. Weisskopf, Theoretical nuclear physics, John Wiley and Sons, New York, 1952.  p. 361
  * [33] S. M. Flatté, _Coupled-channel analysis of the $\pi\eta$ and $K{{\kern 1.79993pt\overline{\kern-1.79993ptK}}}$ systems near $K{{\kern 1.79993pt\overline{\kern-1.79993ptK}}}$ threshold_, Phys. Lett. B63 (1976) 224
  * [34] S. Godfrey and N. Isgur, _Mesons in a relativized quark model with chromodynamics_ , Phys. Rev. D32 (1985) 189