QuNu Collaboration

# Simulating Parton Fragmentation on Quantum Computers

Tianyin Li  Key Laboratory of Atomic and Subatomic Structure and Quantum
Control (MOE), Guangdong Basic Research Center of Excellence for Structure and
Fundamental Interactions of Matter, Institute of Quantum Matter, South China
Normal University, Guangzhou 510006, China  Guangdong-Hong Kong Joint
Laboratory of Quantum Matter, Guangdong Provincial Key Laboratory of Nuclear
Science, Southern Nuclear Science Computing Center, South China Normal
University, Guangzhou 510006, China  Hongxi Xing<EMAIL_ADDRESS>Key
Laboratory of Atomic and Subatomic Structure and Quantum Control (MOE),
Guangdong Basic Research Center of Excellence for Structure and Fundamental
Interactions of Matter, Institute of Quantum Matter, South China Normal
University, Guangzhou 510006, China  Guangdong-Hong Kong Joint Laboratory of
Quantum Matter, Guangdong Provincial Key Laboratory of Nuclear Science,
Southern Nuclear Science Computing Center, South China Normal University,
Guangzhou 510006, China  Southern Center for Nuclear-Science Theory (SCNT),
Institute of Modern Physics, Chinese Academy of Sciences, Huizhou 516000,
China  Dan-Bo Zhang<EMAIL_ADDRESS>Key Laboratory of Atomic and
Subatomic Structure and Quantum Control (Ministry of Education), Guangdong
Basic Research Center of Excellence for Structure and Fundamental Interactions
of Matter, School of Physics, South China Normal University, Guangzhou 510006,
China  Guangdong Provincial Key Laboratory of Quantum Engineering and Quantum
Materials, Guangdong-Hong Kong Joint Laboratory of Quantum Matter, and
Frontier Research Institute for Physics,
South China Normal University, Guangzhou 510006, China

###### Abstract

Parton fragmentation functions (FFs) are indispensable for understanding
processes of hadron production ubiquitously existing in high-energy
collisions, but their first principle determination has never been realized
due to the insurmountable difficulties in encoding their operator definition
using traditional lattice methodology. We propose a framework that makes a
first step for evaluating FFs utilizing quantum computing methodology. The key
element is to construct a semi-inclusive hadron operator for filtering out
hadrons of desired types in a collection of particles encoded in the quantum
state. We illustrate the framework by elaborating on the Nambu-Jona-Lasinio
model with numeral simulations. Remarkably, We show that the semi-inclusive
hadron operator can be constructed efficiently with a variational quantum
algorithm. Moreover, we develop error mitigation techniques tailed for
accurately calculating the FFs in the presence of quantum noises. Our work
opens a new avenue for investigating QCD hadronization on near-term quantum
computers.

Introduction.- Parton fragmentation functions (FFs) quantify the probability
densities for the transition of colored partons to colorless hadrons
Berman:1971xz ; Feynman:1972 ; Collins:2011zzd . It stands out as an essential
ingredient for describing high energy hadron production in electron-positron,
lepton-hadron, and hadron-hadron collisions Metz:2016swz . Parton FFs, as
dynamical non-perturbative quantities, not only can help us to understand the
fundamental mechanism of QCD hadronization, but also provide us with a crucial
tool to access the nucleon structure in the upcoming era of electron-ion
colliders AbdulKhalek:2021gbh ; Anderle:2021wcy , and to probe the transport
properties of quark gluon plasma created in heavy ion collisions at
Relativistic Heavy Ion Collider and the Large Hadron Collider JET:2013cls ;
Qin:2015srf .

Determination of FFs from the fundamental QCD theory is crucial for our
understanding of color confinement. In the past decades, tremendous efforts
have been devoted for extracting FFs from global analysis of world data based
on the QCD factorization theorem, see e.g. deFlorian:2007aj ; Hirai:2007cx ;
Bertone:2018ecm ; AbdulKhalek:2022laj ; Soleymaninia:2022alt ; Moffat:2021dji
; Li:2024etc ; Gao:2024nkz ; Borsa:2022vvp . However, direct calculation of
FFs in a first-principle fashion is still in lack of exploration, in sharp
contrast with their counter part nonperturbative functions like parton
distribution functions (PDFs) and light cone distribution amplitudes, where
the relevant lattice investigation has formed a main branch in particle and
nuclear physics community Lin:2017snn . The main challenges that hinder the
lattice calculation of FFs are mainly due to the following two intrinsic
difficulties. One reason is that FFs are essentially real-time dynamical
quantities, their direct calculation using lattice QCD encounters sign
problems in the traditional Monte Carlo method Alexandru:2016gsd . In
addition, in the operator definition of FFs, the key point is to construct the
semi-inclusive hadron operator (SIHO)
$\mathcal{H}\equiv\sum_{X}\ket{h,X}\bra{h,X}$, in which one needs to sum over
all possible multi-particle states $\ket{h,X}$ containing the identified
hadron $h$. However, the unidentified $X$ can be any state as long as they
satisfy the QCD conservation laws Collins:2023cuo . Such fickle states
significantly challenge the encoding of $X$ in lattice QCD framework, and the
sum over of all possible $X$ involves an exponentially increasing complexity
in the numeral calculation. Addressing these difficulties lies far beyond the
capabilities of classical computing for FFs in a first principle fashion.

Quantum computing, as a promising method alternative to the classical Monte
Carlo method, has attracted lots of attention to simulate generic quantum
systems Feynman:1981tf  involving QCD. Typical examples including the proof of
polynomial time for simulating quantum field theory using quantum computing
Jordan:2012xnu , followed by tremendous ideas for the evaluation of real time
dynamics of QCD, such as the calculation of hadron partonic structure
Lamm:2019uyc ; Echevarria:2020wct ; Mueller:2019qqj ; Qian:2021jxp ;
Li:2021kcs ; Li:2022lyt ; Grieninger:2024cdl , simulation of parton shower
Bauer:2019qxa ; Bepari:2021kwv ; Bauer:2023ujy , jet physics Florio:2023dke ;
Qian:2023vzn ; Barata:2022wim  and string breading in thermal medium
Lee:2023urk ; Xie:2022jgj , as well as the calculation of scattering amplitude
Briceno:2020rar ; Li:2023kex ; Barata:2024bzk . For more broad applications,
see e.g. Refs. Bauer:2022hpo ; Bauer:2023qgm . Inspired by the great advantage
of using quantum computing for QCD real time problems, it is foreseeable to
extend the quantum computing approach, especially the one for direct
calculation of PDFs from light-cone correlators Li:2022lyt , for evaluating
FFs. Nevertheless, how to formulate a semi-inclusive process and efficiently
encode the associated hadron operator with qubits remains unexplored for
quantum computing.

In this letter, we propose a quantum computing approach for the first direct
evaluation of parton fragmentation functions. We present a general framework
to construct the semi-inclusive hadron operator for filtering out the
identified hadrons in a collection of particles encoded in the quantum state.
Utilizing a variational quantum algorithm, FFs are efficiently converted to
multi-point dynamical correlators. We validate the approach in a prototypical
model of Nambu-Jona Lasinio (NJL) using classical hardware Nambu:1961fr ;
Nambu:1961tp ; Gross:1974jv . Moreover, we develop error mitigation techniques
tailed for accurately calculating the FFs in the presence of quantum noises.
Our work clarifies the operational meaning of the semi-inclusive process as
well as provides a first step for calculating FFs on near-term quantum
processors.

Encoding FFs on quantum computer.- The operator definition of collinear quark
fragmentation function in $d$-dimension spacetime can be written as
Collins:2011zzd

$\displaystyle D^{h}_{q}(z)=$ $\displaystyle
z^{d-3}\int\frac{dy^{-}}{4\pi}e^{-iy^{-}p^{+}/z}\Tr\\{\bra{\Omega}\psi(y^{-})$
$\displaystyle\times\sum_{X}\ket{h,X}\bra{h,X}\bar{\psi}(0)\ket{\Omega}\gamma^{+}\Big{\\}}\,,$
(1)

where $z$ is the momentum fraction of the fragmenting quark $q$ carried by the
identified hadron $h$. In Eq. (Simulating Parton Fragmentation on Quantum
Computers), the Wilson line is implicit,
$\psi(y^{-})=e^{iHt}\psi(-y^{z})e^{-iHt}$ is quark field on light-cone with
light-cone coordinates defined as $y^{\pm}=(t\pm y^{z})/\sqrt{2}$.
$\ket{\Omega}$ is the vacuum state and $p^{+}$ is the light cone momentum of
hadron $h$. $\ket{h,X}$ denotes a state that encodes hadron $h$ with momentum
$p^{+}$ as well as $X$ representing additional unidentified fragments. The
summation of $X$ in Eq. (Simulating Parton Fragmentation on Quantum Computers)
stresses that only the hadron $h$ is observed.

Figure 1:  Figure (a) shows qubit configurations of bare vacuum and one-
particle states of in (1+1)-dimensional one flavor theory. (b) is the quantum
circuit for variational quantum algorithms to prepare the vacuum and one-
particle state. (c) is the quantum circuit for the simulation of Eq.
(Simulating Parton Fragmentation on Quantum Computers).

To address the challenge for direct evaluation of FFs, it is inspiring to
compare it with PDFs. PDFs involve dynamical light-cone correlator
$\bra{p}\bar{\psi}(y^{-})\gamma^{+}\psi(0)\ket{p}$, which can be calculated
using the quantum algorithm as elaborated in our previous work Li:2021kcs . In
addition to this dynamical light-cone correlator, the semi-inclusive hadron
operator $\mathcal{H}\equiv\sum_{X}\ket{h,X}\bra{h,X}$ is inserted in the
definition of FFs. The most challenging part of simulating FFs is the
construction of $\mathcal{H}$ due to the existence of unidentified particles
$X$. For quantum computing, figuring out concisely what are $h$ and $X$ in
$\ket{h,X}$ as well as construction of $\mathcal{H}$ is unexplored. The goal
of this work is to make the first step both in concept and in practice on how
to construct the semi-inclusive hadron operator $\mathcal{H}$ and perform a
direct evaluation of parotn FFs.

We propose a scheme for constructing $\mathcal{H}$ in a general setting. In
particular, we present the methodology in $(1+1)$-dimensional system including
the preparation of multi-particle states and the construction of
$\mathcal{H}$, the extension to $(3+1)$-dimensional case can be pursued in a
similar fashion.

To be general, we introduce a computational basis $\ket{\boldsymbol{I}}$ of
quantum computing for a one-dimensional spacial lattice system with total
space site $M$, and each site contains $K$ qubits, i.e.
$\ket{\boldsymbol{I}}=\ket{\boldsymbol{I}_{0},\boldsymbol{I}_{1},...,\boldsymbol{I}_{M-1}}$.
$\boldsymbol{I}_{i}$ is a bitstring of $K$ bits encoding the internal degrees
of freedom, such as spin, flavor, and color, in addition to the stagger
fermion. The total number of qubits for the lattice system is $N=KM$. For each
$\boldsymbol{I}_{i}$, $2^{K}$ configurations can be mapped to bare vacuum
($\boldsymbol{I}^{\Omega}_{i}$), different types of bare elementary particles
($\boldsymbol{I}^{q}_{i},\boldsymbol{I}^{\bar{q}}_{i},\boldsymbol{I}^{g}_{i},...$
for quark, anti-quark and gluon,…), and different types of bare hadrons
($\boldsymbol{I}^{h}_{i}$ with $h=\pi^{\pm},K^{\pm},...$). The bare particles
are defined in the large quark mass limit $m\to\infty$. For example, as shown
in Fig. 1 (a), the bare vacuum, bare quark $q$, bare anti-quark $\bar{q}$ and
bare hadron $h_{q\bar{q}}$ are $\ket{\boldsymbol{I}^{\Omega}_{i}}=\ket{01}$,
$\ket{\boldsymbol{I}^{q}_{i}}=\ket{11}$,
$\ket{\boldsymbol{I}^{\bar{q}}_{i}}=\ket{00}$ and
$\ket{\boldsymbol{I}^{h_{q\bar{q}}}_{i}}=\ket{10}$ in the $(1+1)$-dimensional
one flavor theory. Then, one can convert them to physical vacuum
$\ket{\Omega}$ and physical single-particle states $\ket{\Psi}$ with a unitary
transformation $U$ that preserve the quantum numbers

$\displaystyle\ket{\Omega}=U\ket{\tilde{\Omega}},~{}~{}~{}\ket{\Psi}=\frac{1}{\sqrt{M}}\sum_{i=0}^{M-1}U\ket{\tilde{\Psi}_{i}},$
(2)

where
$\ket{\tilde{\Omega}}\equiv\ket{\boldsymbol{I}^{\Omega}_{0},\boldsymbol{I}^{\Omega}_{1},...,\boldsymbol{I}^{\Omega}_{M-1}}$,
$\ket{\tilde{\Psi}_{i}}\equiv\ket{\boldsymbol{I}^{\Omega}_{0},\boldsymbol{I}^{\Omega}_{1},...,\boldsymbol{I}^{\Omega}_{i-1},\boldsymbol{I}^{\Psi}_{i},\boldsymbol{I}^{\Omega}_{i+1},...,\boldsymbol{I}^{\Omega}_{M-1}}$
is the bare one-particle state at space point $i$, and the wavepacket of a
physical particle $\ket{\Psi}$ centered in the $i$-th space point can be
expressed as $\ket{\Psi_{i}}=U\ket{\tilde{\Psi}_{i}}$. Notice that $\Psi$
stands for all one-particle configurations $\Psi=q,\bar{q},g,...,h$. In this
setup, a state with two particles $\ket{\Psi}$ and $\ket{\Psi^{\prime}}$
centered in the $i$-th and $j$-th space point can be prepared as

$\displaystyle\ket{\Psi_{i},\Psi^{\prime}_{j}}=U\ket{\tilde{\Psi}_{i},\tilde{\Psi}^{\prime}_{j}}\,,$
(3)

where
$\ket{\tilde{\Psi}_{i},\tilde{\Psi}^{\prime}_{j}}\equiv\ket{\boldsymbol{I}^{\Omega}_{0},...,\boldsymbol{I}^{\Omega}_{i-1},\boldsymbol{I}^{{\Psi}}_{i},\boldsymbol{I}^{\Omega}_{i+1}...,\boldsymbol{I}^{{\Psi}^{\prime}}_{j},...,\boldsymbol{I}^{\Omega}_{M-1}}$
is the bare two-particle state. Likewise, the multi-particle state can also be
prepared by acting $U$ on the bare multi-particle state
$\ket{\tilde{\Psi}_{i},\tilde{\Psi}^{\prime}_{j},\tilde{\Psi}^{\prime\prime}_{k},...}$,
which is also the computational basis. The definition of the bare multi-
particle state is similar to the bare two-particle state.

Now, we are ready to define $X$ on a quantum computer. $X$ contains vacuum
and/or all possible particles, i.e. $X=\Omega,q,\bar{q},g,...,h$, as long as
the quantum number of the fragmenting parton is conserved in transition to
$h+X$ system. The $\ket{h,X}$ state with $h$ centered in space point $i$ and
$X$ distributed in the space point except $i$ can be written as

$\displaystyle\ket{h_{i},X_{\\{j\not=i\\}}}=U\ket{\tilde{h}_{i},\tilde{X}_{\\{j\not=i\\}}}\,,$
(4)

where
$\ket{\tilde{h}_{i},\tilde{X}_{\\{j\not=i\\}}}\equiv\ket{\boldsymbol{I}^{X_{0}}_{0},...,\boldsymbol{I}^{X_{i-1}}_{i-1},\boldsymbol{I}^{h}_{i},\boldsymbol{I}^{X_{i+1}}_{i+1},...,\boldsymbol{I}^{X_{M-1}}_{M-1}}$.
Then the local operator $\mathcal{H}_{i}$ can be constructed as

$\displaystyle\mathcal{H}_{i}$
$\displaystyle=U\Tr_{\\{j\not=i\\}}\ket{\tilde{h}_{i},\tilde{X}_{\\{j\not=i\\}}}\bra{\tilde{h}_{i},\tilde{X}_{\\{j\not=i\\}}}U^{\dagger}$
$\displaystyle=U\ket{\boldsymbol{I}^{h}_{i}}\bra{\boldsymbol{I}^{h}_{i}}\otimes{\rm
Id}_{\\{j\not=i\\}}U^{\dagger}\,,$ (5)

where $\Tr_{\\{j\not=i\\}}=\sum_{X_{0},...,X_{i-1},X_{i+1},...X_{M-1}}$ stands
for the summation over all possible $X$. Thanks to the direct product
structure and the completeness of computational basis
$\ket{\tilde{h}_{i},\tilde{X}_{\\{j\not=i\\}}}$, the trace leads to identity
operator ${\rm Id}_{\\{j\not=i\\}}$ that acts on all qubits except for the
space point $i$. For the sake of convenience and considering that FFs are
Lorentz invariant, we perform our calculation in the hadron rest frame. To
select out the zero-momentum hadron state, we construct the zero-momentum
projector $\frac{1}{M}\sum_{j=0}^{M-1}T^{j}$, where $T$ is the translation
operator that commutes with $U$. Therefore, the semi-inclusive hadron operator
in FFs can be written as

$\displaystyle\mathcal{H}=\left(\sum_{i=0}^{M-1}U\ket{\boldsymbol{I}^{h}_{i}}\bra{\boldsymbol{I}^{h}_{i}}U^{\dagger}\right)\times\left(\frac{1}{M}\sum_{j=0}^{M-1}T^{j}\right)\,,$
(6)

where the operator ${\rm Id}_{\\{j\not=i\\}}$ is implicit, and more details
about $\mathcal{H}$ can be found in Sec. A of supplemental material.

The next key ingredient is to construct the unitary matrix $U$ that transforms
the computational basis of bare particles into states of physical particles.
As shown in Fig. 1 (b), to realize $U$ on a quantum computer, we refer to a
variational quantum computing approach, where the unitary matrix is
parameterized as $U(\theta)$ with a Hamiltonian variational ansatz which can
preserve the quantum numbers Li:2021kcs , see details in Sec. B of
supplemental material. Correspondingly, vacuum trial states
$\ket{\Omega(\theta)}$, one-particle trial states $\ket{\Psi(\theta)}$ and
two-particle trial states $\ket{\Psi(\theta),\Psi^{\prime}(\theta)}$ can be
generated by replacing $U$ as $U(\theta)$ in Eqs. (2) and (3). The generation
of multi-particle trial states can be pursued similarly. To optimize the
parameters $\theta$, we implement the subspace-search VQE algorithm by
minimizing a weighted summation of variational energies of
$\ket{\Omega(\theta)}$ and different multi-particle trial states
$\ket{(\Psi,\Psi^{\prime},\Psi^{\prime\prime},...)(\theta)}$
PhysRevResearch.1.033062 . In particular, considering limited computing
resources, we apply a truncation of the multi-particle states. We argue that
it is a reasonable approximation to only generate the vacuum trial state
$\ket{\Omega(\theta)}$ and one-particle trial state $\ket{\Psi(\theta)}$ if
the hadron-hadron interaction energy $V_{hh^{\prime}}$ is much smaller than
the total mass $m_{h}+m_{h^{\prime}}$ of two hadrons. The error of such an
approximation is of $O(\overline{n}^{2}\Delta)$, where $\overline{n}$ is the
average number of particles in $X$ and
$\Delta^{2}=\frac{V_{hh^{\prime}}}{m_{h}+m_{h^{\prime}}}$, see Sec. C in
supplemental material.

Proof of concept study in NJL model.- As a proof of concept, we resort to the
one-flavor 1+1D NJL model to validate the proposed quantum algorithm. The
Lagrangian of the one-flavor NJL model is given by

$\mathcal{L}=\bar{\psi}(i\gamma^{\mu}\partial_{\mu}-m)\psi+g(\bar{\psi}\psi)^{2}\,,$
(7)

where $\psi$ is fermionic operator, $g$ is the strong coupling constant and
$m$ is the bare quark mass. The NJL model describes Dirac fermion interactions
that can mimic the behavior of QCD to some degree, and capture key
characteristics of FFs Matevosyan:2011ey ; Yang:2016gnd ; Yang:2020wvt . In
our calculation, the NJL model is discretized with staggered fermion scheme
and mapped to qubits using the standard Jordan-Winger transformation
Jordan:1928wi ; Fradkin:1989yrd . The details as well as the qubit Hamiltonian
can be found in Ref. Li:2021kcs . The qubit representation of quark FF in the
hadron rest frame can be written as

$\displaystyle zD^{h}_{i}(z)=\sum_{y}e^{-iym_{h}/z}\tilde{D}^{h}_{i}(y)\,,$
(8)

with

$\displaystyle\tilde{D}^{h}_{q}(y)=$
$\displaystyle\sum_{k,l=0}^{1}\bra{\Omega}e^{iHy}\phi_{-2y+k}e^{-iHy}\mathcal{H}\phi^{\dagger}_{l}\ket{\Omega}\,,$
(9)

where $\phi_{n}=\prod_{i<n}\sigma^{3}_{i}\sigma_{n}^{+}$ is a string of Pauli
operators backens_shnirman_makhlin_2019 . Implementing the proposed quantum
algorithm presented in the previous section, the semi-inclusive hadron
operator in the NJL model can be written as,

$\displaystyle\mathcal{H}=M^{-1}$
$\displaystyle\sum_{ij=0}^{M-1}U(\theta^{*})\tilde{\mathcal{H}}_{i}U^{\dagger}(\theta^{*})T^{j}\,,$
(10)

where
$\tilde{\mathcal{H}}_{i}=\ket{\boldsymbol{I}^{h_{q\bar{q}}}_{i}}\bra{\boldsymbol{I}^{h_{q\bar{q}}}_{i}}=(I-\sigma^{z}_{2i})(I+\sigma^{z}_{2i+1})$.
$U(\theta^{*})$ can be realized by the variational quantum algorithm and it
should map the bare vacuum and one-particle state to the corresponding
physical states. In the end, the discretized FFs in position space can be
written as

$\displaystyle\tilde{D}^{h}_{q}(y)=$ $\displaystyle
M^{-1}\sum_{kl=0}^{1}\sum_{ij=0}^{M-1}\bra{\Omega}e^{-iHy}\phi_{-2y+k}e^{iHy}$
$\displaystyle\times
U(\theta^{*})\tilde{\mathcal{H}}_{i}U^{\dagger}(\theta^{*})\phi^{\dagger}_{2j+l}\ket{\Omega}\,,$
(11)

which resembles multiple-point correlation functions.

We show in Fig. 1(c) the quantum circuit for evaluating Eq. (Simulating Parton
Fragmentation on Quantum Computers). With the quantum alternating operator
ansatz (QAOA) farhi2014quantum ; wiersema-PRXQuantum2020  for constructing
$U(\theta^{*})$, the circuit depth for preparing vacuum $\ket{\Omega}$ is
expected to be $O(N)$. The circuit depth for the controlled $U(\theta^{*})$ is
$O(N^{2})$. According to Ref. Li:2021kcs , the time complexity for Trotter
decomposition of $e^{-iHt}$ with an evolution time $t=O(N)$ is
$O(N^{2}/\varepsilon)$ for a desired accuracy $\varepsilon$. The total circuit
depth for the circuit in Fig. 1 is $O(N^{3}/\varepsilon)$. The polynomial
scaling with the number of qubits shows a quantum advantage of calculating FFs
on a quantum computer.

Figure 2: The real (top panel) and imaginary (bottom panel) part of FF in
position space with difference qubit number $N$ in NJL model, the dressed
quark mass is fixed as $m_{q}a=0.4$, and hadron mass is $m_{h}a=0.6$.

Using open source package QuSpin quspin  and projectQ
Steiger2018projectqopensource , we simulate the quantum circuit in Fig. 1 (c)
on a classical device with one CPU and one GPU for different system size. We
present our results of quark FF using different numbers of qubits, i.e.
$N=10$, $N=14$, $N=18$ and $N=22$, with fixed dressed quark mass $m_{q}a=0.4$
and hadron mass $m_{h}a=0.6$. Shown in Fig. 2 are the real and imaginary parts
of FFs in position space, with $\tilde{D}^{h}_{q}$ normalized to ${\rm
Re}[\tilde{D}^{h}_{q}(0)]=1$. We find out that $\tilde{D}^{h}_{q}(y)$
converges with the increase of qubit number $N$, which suggests that the
contribution from high multiplicity events is suppressed in our considered
setup.

The $N$ dependence of quark FF $D^{h}_{q}(z)$ in momentum space is shown in
Fig. 3. The convergence of $D^{h}_{q}(z)$ in the small $z$ region is better
than that in the large $z$ region, this is because the FF is more
significantly affected by the finite volume effect in large $z$ region than
that in the small $z$ region. Such finite volume effect also causes the non-
zero value $D^{h}_{q}(z=1)$, similar to those for PDF calculation in lattice
QCD Ishikawa:2017iym . The shapes shown in Fig. 3 are in qualitative agreement
with the extracted quark FFs Gao:2024nkz  and the calculated charm
fragmentation functions using NJL model Yang:2020wvt  and 1+1D QCD in large
$N_{c}$ limit Jia:2023kiq , such as the vanishing of $D^{h}_{q}(z)$ in small
$z$ region and the peak in large $z$ region.

Figure 3: FF in NJL model with different quantum number $N$, which are
obtained by Fourier transform of $\tilde{D}^{h}_{q}(y)$ in Fig. 2.

Moreover, we investigate the effect from quantum noises and explore how error
mitigation can improve the results. In the simulation, we adopt a
depolarization model to describe noises on quantum hardware. This noisy
channel maps the one-qubit density matrix $\rho$ to
$\epsilon(\rho)=(1-p)\rho+\sum_{i}\frac{p}{3}\sigma^{i}\rho\sigma^{i}$ after
acting arbitrary quantum gates on the qubit. The noise level is set as
$p=10^{-3}$ which is feasible for near-term quantum processors. We find out
that the effect from quantum noises can not be neglected and could lead to
unreliable results, which is caused by the fact that both the vacuum and
hadron states prepared with noisy quantum circuits could be mixed with non-
physical states of single quark or antiquark. Remarkably, the mixing with non-
physical states can make a direct zero-noise extrapolation method
PhysRevLett.119.180509 ; 9259940  fails to work. We find that it is necessary
to first mitigate the noise causing non-physical states by postselection with
a projective measurement on the auxiliary qubit, see Sec. D in supplemetal
material. We incorporate postselection for FFs in position space with higher-
order Richardson extrapolation method PhysRevLett.119.180509

$\displaystyle\tilde{D}^{h,\lambda}_{q}(y)=\sum_{j=0}^{\lambda}\gamma_{j}\tilde{D}^{h}_{q}(y;c_{j}p)\,,$
(12)

where $\gamma_{j}=\prod_{m\not=j}c_{m}(c_{j}-c_{m})^{-1}$,
$\tilde{D}^{h,\lambda}_{q}(y)$ is the $\lambda$-th order error mitigation
result of noisy FF $\tilde{D}^{h}_{q}(y;c_{j}p)$, $c_{j}$ is free scale
parameters. Here, we fix $p=10^{-3}$ and vary $c_{j}=(1+0.1j)\times 10^{-3},\
j=1,2,...,7$.

Figure 4: Real (top panel) and imaginary (bottom panel) part of FF in position
space with $N=6$, $ma=1$, $g=0.8$ and quantum gate fidelity $0.999$, NN is the
result with the zero-noise quantum circuit, $\lambda=1$ is the noisy result
after postselection and $\lambda=k\,,k>1$ is the $k$-th Richardson
extrapolation result after postselection.

As shown in Fig. 4, both $\real[\tilde{D}^{h}_{q}(y)]$ and
$\imaginary[\tilde{D}^{h}_{q}(y)]$ converge to the zero-noise results with
higher-order Richardson extrapolation method, indicating the effectiveness of
error mitigation in our proposed algorithm for direct calculation of FFs.

Summary.- We proposed a framework for the first direct calculation of FFs
utilizing quantum computing methodology. In particular, we constructed the
hadron operator for filtering out semi-inclusive hadrons of desired types in a
collection of multi-particles encoded in the quantum state. As a proof of
concept, we demonstrated the viability of our proposed quantum computing
framework by calculating quark FF in the 1+1 dimensional NJL model up to $22$
qubits. Moreover, we show explicitly that it is necessary to incorporate the
postselection of physical states into the error mitigation for accurately
determining FFs on noisy quantum processors. Our exploratory study manifests
the quantum advantage in the direct calculation of FFs, which paves the way
for investigating fragmentation on near-term quantum hardwares.

Acknowledgements.- We thank Xiaohui Liu, Xingyu Guo and Wai Kin Lai for
helpful discussions. The work is supported by the Guangdong Major Project of
Basic and Applied Basic Research No. 2020B0301030008 and 2023A1515011460, and
NSFC under Grants Nos. 12035007, 12375013 and 12022512.

## References

  * (1) S. M. Berman, J. D. Bjorken, and J. B. Kogut,  Phys. Rev. D 4, 3388 (1971).
  * (2) R. P. Feynman,  Photon-hadron Interactions (CRC Press, 1972).
  * (3) J. Collins,  Foundations of perturbative QCD (Cambridge University Press, 2013).
  * (4) A. Metz and A. Vossen,  Prog. Part. Nucl. Phys. 91, 136 (2016), arXiv:1607.02521.
  * (5) R. Abdul Khalek et al.,  Nucl. Phys. A 1026, 122447 (2022), arXiv:2103.05419.
  * (6) D. P. Anderle et al.,  Front. Phys. (Beijing) 16, 64701 (2021), arXiv:2102.09222.
  * (7) JET, K. M. Burke et al.,  Phys. Rev. C 90, 014909 (2014), arXiv:1312.5003.
  * (8) G.-Y. Qin and X.-N. Wang,  Int. J. Mod. Phys. E 24, 1530014 (2015), arXiv:1511.00790.
  * (9) D. de Florian, R. Sassot, and M. Stratmann,  Phys. Rev. D 75, 114010 (2007), arXiv:hep-ph/0703242.
  * (10) M. Hirai, S. Kumano, T. H. Nagai, and K. Sudoh,  Phys. Rev. D 75, 094009 (2007), arXiv:hep-ph/0702250.
  * (11) NNPDF, V. Bertone, N. P. Hartland, E. R. Nocera, J. Rojo, and L. Rottoli,  Eur. Phys. J. C 78, 651 (2018), arXiv:1807.03310.
  * (12) MAP (Multi-dimensional Analyses of Partonic distributions), R. Abdul Khalek, V. Bertone, A. Khoudli, and E. R. Nocera,  Phys. Lett. B 834, 137456 (2022), arXiv:2204.10331.
  * (13) M. Soleymaninia, H. Hashamipour, and H. Khanpour,  Phys. Rev. D 105, 114018 (2022), arXiv:2202.10779.
  * (14) Jefferson Lab Angular Momentum (JAM), E. Moffat, W. Melnitchouk, T. C. Rogers, and N. Sato,  Phys. Rev. D 104, 016015 (2021), arXiv:2101.04664.
  * (15) M. Li, D. P. Anderle, H. Xing, and Y. Zhao,  (2024), arXiv:2404.11527.
  * (16) J. Gao, C. Liu, X. Shen, H. Xing, and Y. Zhao,  (2024), arXiv:2401.02781.
  * (17) I. Borsa, R. Sassot, D. de Florian, M. Stratmann, and W. Vogelsang,  Phys. Rev. Lett. 129, 012002 (2022), arXiv:2202.05060.
  * (18) H.-W. Lin et al.,  Prog. Part. Nucl. Phys. 100, 107 (2018), arXiv:1711.07916.
  * (19) A. Alexandru, G. Basar, P. F. Bedaque, S. Vartak, and N. C. Warrington,  Phys. Rev. Lett. 117, 081602 (2016), arXiv:1605.08040.
  * (20) J. Collins and T. C. Rogers,  Phys. Rev. D 109, 016006 (2024), arXiv:2309.03346.
  * (21) R. P. Feynman,  Int. J. Theor. Phys. 21, 467 (1982).
  * (22) S. P. Jordan, K. S. M. Lee, and J. Preskill,  Science 336, 1130 (2012), arXiv:1111.3633.
  * (23) NuQS, H. Lamm, S. Lawrence, and Y. Yamauchi,  Phys. Rev. Res. 2, 013272 (2020), arXiv:1908.10439.
  * (24) M. G. Echevarria, I. L. Egusquiza, E. Rico, and G. Schnell,  Phys. Rev. D 104, 014512 (2021), arXiv:2011.01275.
  * (25) N. Mueller, A. Tarasov, and R. Venugopalan,  Phys. Rev. D 102, 016007 (2020), arXiv:1908.07051.
  * (26) W. Qian, R. Basili, S. Pal, G. Luecke, and J. P. Vary,  Phys. Rev. Res. 4, 043193 (2022), arXiv:2112.01927.
  * (27) QuNu, T. Li et al.,  Phys. Rev. D 105, L111502 (2022), arXiv:2106.03865.
  * (28) QuNu, T. Li et al.,  Sci. China Phys. Mech. Astron. 66, 281011 (2023), arXiv:2207.13258.
  * (29) S. Grieninger, K. Ikeda, and I. Zahed,  (2024), arXiv:2404.05112.
  * (30) C. W. Bauer, W. A. de Jong, B. Nachman, and D. Provasoli,  Phys. Rev. Lett. 126, 062001 (2021), arXiv:1904.03196.
  * (31) K. Bepari, S. Malik, M. Spannowsky, and S. Williams,  Phys. Rev. D 106, 056002 (2022), arXiv:2109.13975.
  * (32) C. W. Bauer, S. Chigusa, and M. Yamazaki,  Phys. Rev. A 109, 032432 (2024), arXiv:2310.19881.
  * (33) A. Florio et al.,  Phys. Rev. Lett. 131, 021902 (2023), arXiv:2301.11991.
  * (34) W. Qian,  (2023), arXiv:2307.09987.
  * (35) J. a. Barata, X. Du, M. Li, W. Qian, and C. A. Salgado,  Phys. Rev. D 106, 074013 (2022), arXiv:2208.06750.
  * (36) K. Lee, J. Mulligan, F. Ringer, and X. Yao,  Phys. Rev. D 108, 094518 (2023), arXiv:2308.03878.
  * (37) QuNu, X.-D. Xie et al.,  Phys. Rev. D 106, 054509 (2022), arXiv:2205.12767.
  * (38) R. A. Briceño, J. V. Guerrero, M. T. Hansen, and A. M. Sturzu,  Phys. Rev. D 103, 014506 (2021), arXiv:2007.01155.
  * (39) T. Li, W. K. Lai, E. Wang, and H. Xing,  (2023), arXiv:2301.04179.
  * (40) J. a. Barata, K. Ikeda, S. Mukherjee, and J. Raghoonanan,  (2024), arXiv:2406.03536.
  * (41) C. W. Bauer et al.,  PRX Quantum 4, 027001 (2023), arXiv:2204.03381.
  * (42) C. W. Bauer, Z. Davoudi, N. Klco, and M. J. Savage,  Nature Rev. Phys. 5, 420 (2023), arXiv:2404.06298.
  * (43) Y. Nambu and G. Jona-Lasinio,  Phys. Rev. 124, 246 (1961).
  * (44) Y. Nambu and G. Jona-Lasinio,  Phys. Rev. 122, 345 (1961).
  * (45) D. J. Gross and A. Neveu,  Phys. Rev. D 10, 3235 (1974).
  * (46) K. M. Nakanishi, K. Mitarai, and K. Fujii,  Phys. Rev. Res. 1, 033062 (2019).
  * (47) H. H. Matevosyan, A. W. Thomas, and W. Bentz,  Phys. Rev. D 83, 114010 (2011), arXiv:1103.3085,  [Erratum: Phys.Rev.D 86, 059904 (2012)].
  * (48) D.-J. Yang and H.-n. Li,  Phys. Rev. D 94, 054041 (2016), arXiv:1606.04213.
  * (49) D.-J. Yang and H.-n. Li,  Phys. Rev. D 102, 036023 (2020), arXiv:2007.05678.
  * (50) P. Jordan and E. P. Wigner,  Z. Phys. 47, 631 (1928).
  * (51) E. Fradkin,  Phys. Rev. Lett. 63, 322 (1989).
  * (52) S. Backens, A. Shnirman, and Y. Makhlin,  Scientific Reports 9 (2019).
  * (53) E. Farhi, J. Goldstone, and S. Gutmann,  A quantum approximate optimization algorithm, 2014, arXiv:1411.4028.
  * (54) R. Wiersema et al.,  PRX Quantum 1, 020319 (2020).
  * (55) P. Weinberg and M. Bukov,  SciPost Phys. 2, 003 (2017).
  * (56) D. S. Steiger, T. Häner, and M. Troyer,  Quantum 2, 49 (2018).
  * (57) T. Ishikawa et al.,  Sci. China Phys. Mech. Astron. 62, 991021 (2019), arXiv:1711.07858.
  * (58) Y. Jia, Z. Mo, and X. Xiong,  Eur. Phys. J. C 83, 1169 (2023), arXiv:2310.17640.
  * (59) K. Temme, S. Bravyi, and J. M. Gambetta,  Phys. Rev. Lett. 119, 180509 (2017).
  * (60) T. Giurgica-Tiron, Y. Hindy, R. LaRose, A. Mari, and W. J. Zeng,  Digital zero noise extrapolation for quantum error mitigation,  in 2020 IEEE International Conference on Quantum Computing and Engineering (QCE), pp. 306–316, 2020.

## Supplemental Material

### A. Construction of semi-inclusive hadron operator

The construction of semi-inclusive hadron operator (SIHO) has been presented
in the main text

$\displaystyle\mathcal{H}$
$\displaystyle=\left(\sum_{i=0}^{M-1}U\ket{\boldsymbol{I}^{h}_{i}}\bra{\boldsymbol{I}^{h}_{i}}U^{\dagger}\right)\times\left(\sum_{j=0}^{M-1}\frac{1}{M}T^{j}\right)$
$\displaystyle\equiv\mathcal{N}_{h}\times P(\mathbf{p}=0)\,.$ (13)

In the following part of this section, we will show that $\mathcal{N}_{h}$ and
$P(\mathbf{p}=0)$ can be regarded as number operators of hadron $h$ and zero-
momentum projector, respectively.

We can act $\mathcal{N}_{h}$ on multi-particle states to explain why
$\mathcal{N}_{h}$ is a number operator of hadrons. A multi-particle state with
$n_{h}$ hadrons of $h$ localized in positions $i_{1},i_{2},..,i_{n_{h}}$ can
be denoted as
$\ket{h_{i_{1}},h_{i_{2}},...,h_{i_{n_{h}}},\Psi_{\\{j\not=i_{1},i_{2},..,i_{n_{h}}\\}}}$,
where other particles $\Psi\not=h$ and they distribute at the positions except
$i_{1},i_{2},...,i_{n_{h}}$. Note that $\mathcal{N}_{h}$ commutes with $U$ as
$U$ is quantum-number preserving, then we have

$\displaystyle\quad\
\mathcal{N}_{h}\ket{h_{i_{1}},h_{i_{2}},...,h_{i_{n_{h}}},\Psi_{\\{j\not=i_{1},..,i_{n_{h}}\\}}}$
$\displaystyle=\sum_{i=0}^{M-1}U\ket{\boldsymbol{I}^{h}_{i}}\bra{\boldsymbol{I}^{h}_{i}}\otimes{\rm
Id}_{\\{k\not=i\\}}\ket{\boldsymbol{I}^{h}_{i_{1}},\boldsymbol{I}^{h}_{i_{2}},...,\boldsymbol{I}^{h}_{i_{n_{h}}},\boldsymbol{I}^{\Psi}_{\\{j\not=i_{1},..,i_{n_{h}}\\}}}$
$\displaystyle=Un_{h}\ket{\boldsymbol{I}^{h}_{i_{1}},\boldsymbol{I}^{h}_{i_{2}},...,\boldsymbol{I}^{h}_{i_{n_{h}}},\boldsymbol{I}^{\Psi}_{\\{j\not=i_{1},..,i_{n_{h}}\\}}}$
$\displaystyle=n_{h}\ket{h_{i_{1}},h_{i_{2}},...,h_{i_{n_{h}}},\Psi_{\\{j\not=i_{1},..,i_{n_{h}}\\}}}\,.$
(14)

The above equation tells us that the operator $\mathcal{N}_{h}$ can count the
number of hadrons $h$ in multi-particle states.

As mentioned in the main text, the FF needs to be calculated in the hadron
rest frame. In this frame, the state $\ket{h,X}$ also has zero momentum
because we focus on collinear FFs, which have $p_{h}+p_{X}=p_{h}/z$. We need a
zero-momentum projector to guarantee the operator $\mathcal{N}_{h}$ counting
the zero-momentum hadron. To verify $\tilde{P}(\mathbf{p}=0)$ is such
operator, we calculate $\tilde{P}({\boldsymbol{p}}=0)^{2}$ at first,

$\displaystyle\tilde{P}({\boldsymbol{p}}=0)^{2}$
$\displaystyle=\frac{1}{M^{2}}\sum_{ij=1}^{M}T^{i}T^{j}$
$\displaystyle=\frac{1}{M^{2}}\sum_{i=1}^{M}MT^{i}=\tilde{P}({\boldsymbol{p}}=0)\,.$
(15)

showing $\tilde{P}({\boldsymbol{p}}=0)$ is a projector. Then we show that
$\tilde{P}({\boldsymbol{p}}=0)$ acts on a zero-momentum state will keep it
unchanged while acts on the non-zero-momentum state will vanish,

$\displaystyle\tilde{P}({\boldsymbol{p}}=0)\ket{\boldsymbol{p}=0}=\frac{1}{M}\sum_{j=1}^{M-1}T^{j}\ket{\boldsymbol{p}=0}=\ket{\boldsymbol{p}=0}\,.$
(16)

The last step in Eq. (16) is because
$T\ket{\boldsymbol{p}=0}=\ket{\boldsymbol{p}=0}$.

$\displaystyle\tilde{P}({\boldsymbol{p}}=0)\ket{\boldsymbol{p}}=\frac{1}{M}\sum_{j=1}^{M-1}T^{j}\ket{\boldsymbol{p}}=\frac{1}{M}\sum_{j=1}^{M-1}e^{ip_{k}j}\ket{\boldsymbol{p}}=0\,.$
(17)

The last step in Eq. (17) is because $\sum_{j=0}^{M-1}e^{ip_{k}j}=0$, where
$p_{k}=\frac{2\pi k}{M}\,,k=1,...,M-1$ on one dimensional lattice. In
conclusion, the $\tilde{P}({\boldsymbol{p}}=0)$ is exactly the zero-momentum
projector.

### B. Variational quantum algorithm to construct $U$

The unitary operator $U$ can be constructed by a variational quantum algorithm
(VQA) Li:2021kcs . For our purpose, we require $U$ to prepare the vacuum and
one-hadron states instead of preparing all multiple hadron states. Such $U$
suffices for transforming the computational basis of multiple bare hadrons to
states of multiple physical hadrons, once energy is dominated by the mass of
physical hadrons while their interaction energy can be comparably very small,
as justified later in Sec. C C. Error estimation of SIHO. For demonstration,
we use the NJL model to show how the VQA works. The qubit Hamiltonian of the
NJL model can be found in Li:2021kcs , which can be written as
$H=\sum_{j=1}^{4}H_{j}$. Every $H_{j}$ commutes with all the conservation
charges of the NJL model hence the time evolution $\exp(-iH_{j}t)$ preserves
the quantum number of the NJL model, so the variational ansatz $U$ can be
written as

$\displaystyle
U(\theta)=\prod_{i=1}^{p}\prod_{j=1}^{4}\exp(i\theta_{ij}H_{j})\,.$ (18)

The desired $U(\theta)$ needs to map the corresponding four computational
basis to vacuum $\ket{\Omega}$, anti-quark state $\ket{\bar{q}}$, quark state
$\ket{q}$ and $q\bar{q}$ bound state $\ket{h_{q\bar{q}}}$. Those four
computational basis are denoted as $\ket{\tilde{\Omega}}$,
$\ket{\tilde{\bar{q}}}$, $\ket{\tilde{q}}$ and $\ket{\tilde{h}_{q\bar{q}}}$

$\displaystyle\ket{\tilde{\Omega}}=\ket{\boldsymbol{I}^{\Omega}_{0},...,\boldsymbol{I}^{\Omega}_{M-1}}\,,$
$\displaystyle\ket{\tilde{\bar{q}}}=\frac{1}{M}\sum_{i=0}^{M-1}\ket{\boldsymbol{I}^{\Omega}_{0},\boldsymbol{I}^{\Omega}_{1},...,\boldsymbol{I}^{\bar{q}}_{i},\boldsymbol{I}^{\Omega}_{i+1},...,\boldsymbol{I}^{\Omega}_{M-1}}\,,$
$\displaystyle\ket{\tilde{q}}=\frac{1}{M}\sum_{i=0}^{M-1}\ket{\boldsymbol{I}^{\Omega}_{0},\boldsymbol{I}^{\Omega}_{1},...,\boldsymbol{I}^{q}_{i},\boldsymbol{I}^{\Omega}_{i+1},...,\boldsymbol{I}^{\Omega}_{M-1}}\,,$
$\displaystyle\ket{\tilde{h}_{q\bar{q}}}=\frac{1}{M}\sum_{i=0}^{M-1}\ket{\boldsymbol{I}^{\Omega}_{0},\boldsymbol{I}^{\Omega}_{1},...,\boldsymbol{I}^{h_{q\bar{q}}}_{i},\boldsymbol{I}^{\Omega}_{i+1},...,\boldsymbol{I}^{\Omega}_{M-1}}\,,$
(19)

where $\ket{\boldsymbol{I}^{\Omega}}=\ket{01}$,
$\ket{\boldsymbol{I}^{\bar{q}}}=\ket{00}$, $\ket{\boldsymbol{I}^{q}}=\ket{11}$
and $\ket{\boldsymbol{I}^{h_{q\bar{q}}}}=\ket{10}$ in the NJL model. All four
states in Eq. (B. Variational quantum algorithm to construct $U$) are Dicke-
like states, which can be prepared by the quantum circuit in Li:2022lyt . To
obtain the desired $U(\theta)$ we need to optimize parameters $\theta$ by
minimizing the loss function

$\displaystyle E(\theta)=$ $\displaystyle
w_{\Omega}\braket{\tilde{\Omega}}{U^{\dagger}(\theta)HU(\theta)}{\tilde{\Omega}}$
$\displaystyle+\sum_{\Psi}w_{\Psi}\braket{\tilde{\Psi}}{U^{\dagger}(\theta)HU(\theta)}{\tilde{\Psi}}$
(20)

where $\Psi=q,\bar{q},h_{q\bar{q}}$ and $w_{i}$ could be chosen freely so long
as $w_{\Omega}>w_{\bar{q}}=w_{q}>w_{h_{q\bar{q}}}>0$. $w_{\bar{q}}=w_{q}$ is
due to energies of $q$ and $\bar{q}$ should be the same. eThe optimized
parameters are denoted as $\theta^{*}$, then we have
$\ket{\Omega}=U(\theta^{*})\ket{\tilde{\Omega}}$,
$\ket{\bar{q}}=U(\theta^{*})\ket{\tilde{\bar{q}}}$,
$\ket{q}=U(\theta^{*})\ket{\tilde{q}}$,
$\ket{h_{q\bar{q}}}=U(\theta^{*})\ket{\tilde{h}_{q\bar{q}}}$. Now, in the NJL
model, the operator $U$ can be approximated by $U(\theta^{*})$, because it
maps the computational basis to vacuum and one-particle states and preserves
all the symmetry of the NJL model.

### C. Error estimation of SIHO

In the VQA method, the operator $U$ can only prepare vacuum and one-hadron
state, the multi-hadron states are prepared approximately. So the error of the
SIHO needs to be estimated. We assume the two-hadron state prepared by $U$ is
$\ket{\bar{h}_{\alpha},\bar{h}_{\beta}}$ and the exact two-hadron state is
$\ket{h_{\alpha},h_{\beta}}$, which is an eigenstate of Hamiltonian $H$.
Obviously, the state $\ket{\bar{h}_{\alpha},\bar{h}_{\beta}}$ can not reflect
the interaction between two hadrons, so we have

$\displaystyle\braket{\bar{h}_{\alpha},\bar{h}_{\beta}}{H}{\bar{h}_{\alpha},\bar{h}_{\beta}}$
$\displaystyle=\braket{h_{\alpha},h_{\beta}}{H}{h_{\alpha},h_{\beta}}+\epsilon_{\alpha\beta}$
$\displaystyle=E_{\alpha\beta}+\epsilon_{\alpha\beta}\,,$ (21)

where $E_{\alpha\beta}$ is the energy eigenvalue of the state
$\ket{h_{\alpha},h_{\beta}}$. The energy deviation $\epsilon_{\alpha\beta}$ is
roughly equal to the interaction potential between two hadrons, which is
denoted as $V_{h_{\alpha}h_{\beta}}$. It should be stressed
$V_{h_{\alpha}h_{\beta}}$ is not relevant to the original interaction of the
Hamiltonian. Rather, it is residual interaction. The state
$\ket{\bar{h}_{\alpha},\bar{h}_{\beta}}$ can be decomposed as

$\displaystyle\ket{\bar{h}_{\alpha},\bar{h}_{\beta}}=\frac{1}{\sqrt{1+\Delta_{\alpha\beta}^{2}}}\ket{h_{\alpha},h_{\beta}}+\frac{\Delta_{\alpha\beta}}{\sqrt{1+\Delta_{\alpha\beta}^{2}}}\ket{\varphi_{\alpha\beta}}\,,$
(22)

where $\ket{\varphi_{\alpha\beta}}$ is orthogonal to
$\ket{h_{\alpha},h_{\beta}}$. Plug Eq. (22) into LHS of Eq. (C. Error
estimation of SIHO) and ignore the $O(\Delta_{\alpha\beta}^{2})$ term,
$\Delta_{\alpha\beta}$ can be written as

$\displaystyle\Delta_{\alpha\beta}=\sqrt{\frac{\absolutevalue{\epsilon_{\alpha\beta}}}{E_{\alpha\beta}}}\approx\sqrt{\frac{\absolutevalue{V_{h_{\alpha}h_{\beta}}}}{m_{h_{\alpha}}+m_{h_{\beta}}}}\,.$
(23)

The SIHO needs to sum over all multi-hadron states. For the $n$-hadron state,
it has approximately $n^{2}$ two-hadron pairs. So the error of the SIHO
$\epsilon_{h}$ can be written as

$\displaystyle\epsilon_{h}=O[\bar{n}^{2}\,{\rm
max}({\Delta_{\alpha\beta}})]\equiv O(\bar{n}^{2}\Delta),$ (24)

where $\bar{n}$ is the average number of particles in $X$ and ${\rm
max}({\Delta_{\alpha\beta}})$ is the maximum of the $\Delta_{\alpha\beta}$ for
all $\alpha$ and $\beta$. It can be expected that the mass of a hadron stems
mainly from the original interaction of quarks, while the residue interaction
between hadrons is negligibly small. Thus, for a fixed $\bar{n}$, the error
$\epsilon_{h}$ of the SIHO is small because the hadron mass
$m_{h_{\alpha}}+m_{h_{\beta}}$ will be several orders of magnitude larger than
the interaction potential between hadrons.

### D. Error mitigation of FFs

The depolarization noise channel is used to simulate the noisy quantum
computer. The noise channel maps the noiseless density matrix $\rho$ to
$\epsilon(\rho)$

$\displaystyle\epsilon(\rho)=(1-p)\rho+\frac{p}{3}(\sigma^{x}\rho\sigma^{x}+\sigma^{y}\rho\sigma^{y}+\sigma^{z}\rho\sigma^{z}),$
(25)

where $1-p$ is the noise parameter that reflects the quantum gate fidelity. We
set $p\approx 10^{-3}$, which is reasonable for Noisy Intermediate-Scale
Quantum (NISQ) devices. For a $N=6$ qubits quantum computer, the depth of the
quantum circuit to simulate FFs is about $O(10^{2})$. Due to the existence of
SIHO, FF needs to be calculated by adding up $(N/2)^{2}\approx 10$ output
results from different quantum circuits. The noise effect may accumulate after
adding up those outputs, so the noise effect can not be neglected. To obtain
reasonable output from quantum computing, an error mitigation algorithm is
needed. In the main text, the error mitigation algorithm is made up of two
parts, where the first part is postselection and the second part is Richardson
zero noise extrapolation.

Figure 5: Quantum circuit for postselection. Figure 6: Definition of gate
$C_{L}X$.

In the main text, we point out that the unitary operator $U$ need to protect
the quantum number. However, the $\sigma^{x}$ and $\sigma^{y}$ errors will
change the charge of a quantum state by flipping the qubit. We chose to select
the measurement results with even numbers of $\sigma^{x}$ and $\sigma^{y}$
errors in the quantum circuit. The above procedure is called postselection.
Still, it is possible that, for instance, a $q$ state will become a $qqq$
state due to noise, which can not be detected in the postselection. However,
we expect that for small $p$ the one-flip process dominates. In this regard,
the postselection can restore the charge conservation in a good approximation.
Fig. 5 shows the postselection quantum circuit, where the definition of
$C_{L}X$ gate is shown in Fig. 6. We select measurement results with
$\sigma^{z}_{L+1}=+1$ to calculate the FF because the odd numbers $\sigma^{x}$
and $\sigma^{y}$ error will flip the qubit $L+1$ from $\ket{0}$ to $\ket{1}$.
In Fig. 7, we show how the measurement result of qubit $L$ depends on the
noise parameter $p$, where NPo is no postselection result, Po is the
postselection result and NN is the no-noise result. The NISQ device can not
reach the shadow area in Fig. 7 because the noise of this area is smaller than
$10^{-3}$. From Fig. 7, it can be seen that extrapolating zero noise based on
the data in the range of $p>10^{-3}$ for NPo can easily lead to incorrect
results. Meanwhile, the extrapolation of Po data can give us correctly no
noise result, which shows the effectiveness of postselection.

Figure 7: Results of postselection and no postselection with different noise
parameters $p$. The NISQ device can only obtain the result with $p>10^{-3}$
(outside of the shadow region).

There are still some errors after postselection. Those errors can be mitigated
by Richardson zero noise extrapolation PhysRevLett.119.180509 . Richardson
extrapolation suggests one can run a quantum circuit under different noise
parameters, then the error can be mitigated by extrapolating the result to
zero noise by

$\displaystyle O^{\lambda}(p)=\sum_{j=0}^{\lambda}\gamma_{j}O(c_{j}p)$ (26)

where $\gamma_{j}=\prod_{m\not=j}c_{m}(c_{j}-c_{m})^{-1}$, $O^{\lambda}(p)$ is
the $\lambda$-th order error mitigation result of observable $O$ and $c_{j}$
is the scale parameters which can be chosen freely.


# A PINN Approach to Symbolic Differential Operator Discovery with Sparse Data

Lena Podina,
Cheriton School of Computer Science
University of Waterloo
&Brydon Eastman
Department of Applied Mathematics
University of Waterloo
&Mohammad Kohandel
Department of Applied Mathematics
University of Waterloo
equal contribution

(July 2022)

###### Abstract

Given ample experimental data from a system governed by differential
equations, it is possible to use deep learning techniques to construct the
underlying differential operators. In this work we perform symbolic discovery
of differential operators in a situation where there is sparse experimental
data. This small data regime in machine learning can be made tractable by
providing our algorithms with prior information about the underlying dynamics.
Physics Informed Neural Networks (PINNs) have been very successful in this
regime (reconstructing entire ODE solutions using only a single point or
entire PDE solutions with very few measurements of the initial condition). We
modify the PINN approach by adding a neural network that learns a
representation of unknown hidden terms in the differential equation. The
algorithm yields both a surrogate solution to the differential equation and a
black-box representation of the hidden terms. These hidden term neural
networks can then be converted into symbolic equations using symbolic
regression techniques like AI Feynman. In order to achieve convergence of
these neural networks, we provide our algorithms with (noisy) measurements of
both the initial condition as well as (synthetic) experimental data obtained
at later times. We demonstrate strong performance of this approach even when
provided with very few measurements of noisy data in both the ODE and PDE
regime.

##  1 Introduction

Mechanistic models, used for modeling real-world processes in biology,
chemistry and physics, are in the form of a system of differential equations
(DE). For example, in systems biology, glycolytic oscillations in yeast is
modelled using a set of coupled DEs RUOFF2003179 . In chemistry, the toluene
dehydrogenation reaction scheme belohlav1997application  is represented using
ordinary DEs. In physics, the behaviour of a pendulum with damping is modelled
using a second-order ordinary DE strogatz2018nonlinear . For a system of DEs
to be useful in modelling a process, initial conditions and parameters
(henceforth ‘DE parameters’) specific to the system of DEs need to be
estimated. If the correct values for the initial conditions and DE parameters
are obtained, the DE can not only be used to interpolate between experimental
datapoints, but to predict the future state of a dynamical system.

Machine learning algorithms, for instance neural networks (NN), are
particularly helpful in representing unknown quantities in a data-driven way
belohlav1997application . NNs with a wide enough hidden layer can be used to
approximate any function pinkus_1999  by tuning its parameters (henceforth ‘NN
parameters’). Hence, NNs have been used to infer DE parameters or even entire
DE models, due to their ability to approximate functions. Many recent
applications use NNs augmented with prior knowledge in order to learn
underlying DE models111The (underlying or true) DE model of a process refers
to the DEs that produced the experimental data.  from data
transferlearningChakraborty ; RAISSI2019686 ; lu2021neural ;
rackauckas2020universal ; meng2020composite ; lu2021deep . However, acquiring
sufficient data to fit these values accurately using NNs is difficult. A
method that can function in low-data regimes by leveraging the known structure
of the model is needed.

Two prominent NN-based methods that learn DE models from data are physics-
informed neural networks (PINN) RAISSI2019686  and universal differential
equations (UDE) rackauckas2020universal . In UDEs, each unknown component of
the DE model is approximated by a NN and a hard DE constraint is employed.
That is, the best-fit DE is satisfied at all times during training. However,
UDEs are not robust to noise, require a lot of data, and SINDy, as employed in
rackauckas2020universal , does not succeed in finding the true mechanistic
model reliably. PINNs assume the form of the true DE and fits its parameters
via a soft constraint (relaxing the requirement that the NN should satisfy the
best-fit DE exactly), which is added to the NN loss function as an additional
loss term referred to as the ‘PINN loss’. A drawback of PINNs is that the
structure of the DE model must be determined in advance, and there is no way
to learn its unknown components using the method as originally proposed.
Additionally, as iterative optimization is computationally expensive, PINN
loss can fail on stiff DEs wang2021understanding .

Our approach bridges the limitations of both PINNs (cannot be used when the
structure of the DE is not fully known) RAISSI2019686  and UDEs (not robust to
noise and requires lots of data) rackauckas2020universal . To address this, we
replace the hard constraint of the UDE with that of PINN loss, which allows
the approach to learn unknown components of the DE model from data. This
approach is robust to noise and performs well in low-data regimes.
Additionally, using the AI Feynman algorithm udrescu2020ai  yielded good
results in identifying the underlying hidden terms of the DE model.

In this paper, we claim the following contributions:

  1. 1.

We propose a novel training methodology which combines PINN loss with the UDE
framework to allow a PINN-based approach to learn unknown parts of the DE
model from data

  2. 2.

Our model is robust to noise and learns the unknown DE model components well
in the Lotka-Volterra model, cell apoptosis ODE, and Viscous Burgers’ equation

  3. 3.

Using a symbolic regression algorithm, we reached better identification of the
tested systems than in the original UDE paper

##  2 Background

As per RAISSI2019686 , suppose that the following DE governs a physical
process. $u(t,x)$ is an unknown real-valued function of time $(t)$ and
position ($x$). Its time derivative is related to its value for each tuple
$(t,x)$ with a known function $\mathcal{N}$, and unknown vector of parameters
$\theta$. Furthermore, there are $N$ potentially noisy DE measurements
$\\{t_{i},x_{i},u_{i}\\}$.

$\displaystyle\frac{\partial u(t,x)}{\partial
t}=\mathcal{N}[u;\theta],x\in\Omega,\Omega\in R^{D},t\in[0,T]$ (1)

Although $\theta$ is required in order to find a numerical or analytical
function $u$ that satisfies 1, $\theta$ is unknown in this setup. Using the
given data, the PINN method from RAISSI2019686  can estimate $\theta$ and $u$
simultaneously. Its key component is a neural network $U$, which predicts $u$
given any tuple $(t,x)$. The following loss function is used to train $U$:

$\displaystyle\mathcal{L}=\frac{1}{N}\sum_{i=1}^{N}|U(t^{i},x^{i})-u^{i}|+\frac{1}{M}\sum_{j=1}^{M}\bigg{|}\frac{dU}{dt_{j}}-\mathcal{N}[u;\theta]\bigg{|}$
(2)

where $\frac{dU}{dt_{j}}$ is auto-differentiated through the neural network
and evaluated at time $t_{j}$ and position $x_{j}$. The first term penalizes
$U$ for making predictions that do not match the DE at a predefined set of $M$
collocation points $\\{t_{j},x_{j}\\}$. The second term penalizes $U$ for
making predictions that do not match the data. Note that the second term
contains $\theta$, which allows parameter estimate $\hat{\theta}$ to be
updated using its gradient with respect to $\mathcal{L}$. At every
optimization iteration, the parameters of $U$ are updated along with
$\hat{\theta}$.

A UDE rackauckas2020universal  is a DE which is defined in part using
universal approximators, e.g. a neural network. These NNs can be used to
approximate unknown components of $\mathcal{N}$. We will assume that possibly
noisy data $\\{t_{i},x_{i},u_{i}\\}$ is available from Eq. 1. Suppose that
$\mathcal{N}$ is a function $g$ composed of $k$ unknown functions $h_{i}$ and
known parameters $\theta$:

$\displaystyle\mathcal{N}[u;\theta]=g(u,h_{1}(u;\theta),\ldots,h_{k}(u;\theta);\theta)$
(3)

In rackauckas2020universal , the $h_{i}$ terms are approximated by a single
neural network $H$ with $k$ outputs and fit using iterative optimization such
as Adam adam  or gradient descent bishop2006pattern . Since Eq. (1) always
holds, the loss function only consists mean squared error between the DE
solution and the data. Note that with any particular approximation $H$, the DE
(1) is defined fully and $u$ can be solved for numerically. Hence, the
training loop for UDEs involves numerically solving Eq. (1), computing the
error between the solution and the data, and updating $H$ to better
approximate the unknown components of $\mathcal{N}$. At the end of training,
$H$ will represent the unknown component of $\mathcal{N}$ and the numerical
solution of Eq. (1) will yield $u$ that matches the solution of the true DE.

Symbolic regression is a general technique of finding a model that fits data
while balancing the simplicity of the model with its accuracy. This problem
has been solved with various genetic algorithms (see, among others,
gaSymbRegression1 ; gaSymbRegression2 ) but since this method is
computationally expensive, newer techniques are becoming more popular. For
example, AI Feynman udrescu2020ai  leverages NNs and symmetry, units,
compositionality, etc. and finally returns a list of potential models ranked
by error and complexity. Some work SymbolicGPT2021 ; kamienny2022end  makes
use of transformers to find the correct functional form. In our work, we only
use them at the final stage after our method has learned an approximate
representation of the missing components.

##  3 Methods

Our proposed method is a modification of the PINN approach to discover the
functional form of an unknown term within a differential equation. Suppose
$\vec{u}(\vec{x},t)\in\mathbb{R}^{m}$ for $\vec{x}\in\mathbb{R}^{d}$. Let
$\mathcal{N}$ be a (potentially non-linear) differential operator, then
consider time $t$ in the domain $[0,T]\subset\mathbb{R}$ along with a $d$
dimensional, bounded spatial domain $\Omega\subset\mathbb{R}^{d}$ where
$\partial\Omega$ denotes the boundary of $\Omega$. Notably, if $\mathcal{N}$
contains any derivatives, we assume that those derivatives are with respect to
the spatial variables only. We then consider problems of the form

$\frac{{\rm d}}{{\rm
d}t}\vec{u}(\vec{x},t)=\mathcal{N}[\vec{u}](\vec{x},t),\quad
t\in[0,T],\quad\vec{x}\in\Omega$

subject to initial condition

$\vec{u}(\vec{x},0)=\vec{u}_{0}(\vec{x}),\quad\vec{x}\in\Omega$

and boundary conditions

$\beta[\vec{u}](\vec{x},t)=0,\quad\vec{x}\in\partial\Omega,\quad t\in[0,T]$

where $\beta$ is a (potentially non-linear) differential operator whose
derivatives are only with respect to the spatial variables.

Further, suppose

$\mathcal{N}[\vec{u}](\vec{x},t)=\mathcal{N}_{K}[\vec{u}](\vec{x},t)+\mathcal{F}[\vec{u}](\vec{x},t)$

where $\mathcal{N}_{K}$ is some differential operator with known functional
form and $\mathcal{F}$ represents some unknown, target differential operator.
Similarly, suppose $\beta=\beta_{K}+\mathcal{B}$ for some known $\beta_{K}$
and some unknown $\mathcal{B}$.

Finally, one can consider $\Omega=\varnothing$, in which case the underlying
differential law is governed by an ordinary differential equation (ODE). In
this situation, there is no boundary condition and so no need for $\beta$ (or,
equivalently, $\beta$ is the empty function).

Suppose we have $n$ data points
$D=\\{(t_{k},\vec{x}_{k},\vec{u}_{k})\\}_{k=0}^{n-1}$ where
$\vec{u}_{k}=\vec{u}(t_{k},\vec{x}_{k})+\epsilon_{k}$ where $\epsilon_{k}$ is
some noise term (potentially $\epsilon_{k}=0$). We will use this measured data
to fit the parameters of (up to) three neural networks. The first network,
$F(\vec{u};\theta_{F})$, will approximate the target differential operator
$\mathcal{F}[\vec{u}]$ by using a neural network with parameters $\theta_{F}$.
The second network, $U(\vec{x},t;\theta_{U})$, will approximate the value of
$\vec{u}(x,t)$ by a neural network with parameters $\theta_{U}$. The third
network, $B(\vec{u};\theta_{B})$, will approximate the value of
$\mathcal{B}[\vec{u}]$, the unknown target for the boundary condition, with a
neural network parameterized by $\theta_{B}$. To fit these networks, we
consider another two sets of collocation points: these sets are
$X_{P}=\\{(\vec{x}_{k},t_{k})\\}_{k=0}^{n_{P}-1}\subset\left(\Omega\setminus\partial\Omega\right)\times(0,T]$
and
$X_{B}=\\{(\vec{x}_{k},t_{k})\\}_{k=0}^{n_{B}-1}\subset(\partial\Omega)\times(0,T]$.
These sets correspond to locations in the space-time domain where we enforce
that our network $U$ satisfies the underlying differential equation (in the
case of $X_{P}$) and the boundary conditions (in the case of $X_{B}$).

To calculate the gradients for fitting these networks, we consider the loss
function

$L(\theta_{U},\theta_{B},\theta_{F})=L_{M}(\theta_{U})+L_{B}(\theta_{U},\theta_{B})+L_{P}(\theta_{U},\theta_{F}).$

The first component of the loss is the MSE loss. This loss is the difference
in MSE between the measurement value of $\vec{u}\approx\vec{u}_{k}$ from the
input data with the neural network approximation of $\vec{u}\approx
U(\vec{x}_{k},t_{k})$, evaluated at the same space-time location and is given
by

$L_{M}(\theta_{U})=\frac{1}{n}\sum_{(\vec{x}_{k},t_{k},\vec{u}_{k})\in
D}(U(\vec{x}_{k},t_{k};\theta_{U})-\vec{u}_{k})^{2}.$

The second component of the loss is the boundary loss. This loss is the mean
squared value of the approximated value of the boundary condition and is given
by

$L_{B}(\theta_{U},\theta_{B})=\frac{1}{n_{B}}\sum_{(\vec{x}_{k},t_{k})\in
X_{B}}(\beta_{K}[U](\vec{x}_{k},t_{k};\theta_{U})+B(U(\vec{x}_{k},t_{k};\theta_{U});\theta_{B}))^{2}.$

The final component of the loss is the PINN loss. This loss is the mean
squared error between the value $U_{t}$, the time derivative of the neural
network approximation of $U$, and the value $\mathcal{N}_{K}[U]+F(U)$.

$L_{P}(\theta_{U},\theta_{F})=\frac{1}{n_{P}}\sum_{(\vec{x}_{k},t_{k})\in
X_{P}}(\mathcal{N_{K}}[U](\vec{x}_{k},t_{k};\theta_{U})+F(U(\vec{x}_{k},t_{k};\theta_{U});\theta_{F})-U_{t}(\vec{x}_{k},t_{k};\theta_{U}))^{2}.$

This loss function is quite similar to the loss function for PINNs given in
RAISSI2019686 , however here we insert two additional neural networks into the
loss function corresponding to the unknown parts of the underlying dynamics in
the boundary conditions and the differential equation. To compensate for these
additional parameters, we extend the first component of the loss to include
more than just initial data (but solution data as well). In this way, $D$
could contain data from the initial condition, data from the boundary, or data
from the interior of the domain.

Practically, one way to select $X_{P}$ is to simply choose $n_{P}$ and use
Latin hypercube sampling to select $n_{P}$ points in
$\left(\Omega\setminus\partial\Omega\right)\times(0,T]$. A similar
construction works for selecting $X_{B}$. In this way, we are sampling the
domain in a space-filling manner.

##  4 Results

To demonstrate our approach we trained the above model on the following three
test-cases: the Lotka-Volterra equations (an ODE model), a model for cell
apoptosis (an ODE model), and the viscous Burgers’ equation (a parabolic PDE
model).

###  4.1 Lotka-Volterra System

We begin our analysis by testing our method on the Lotka-Volterra (LV) model
berryman1992orgins  of predator-prey interactions. The DE is formulated as
follows:

$\displaystyle\frac{{\rm d}x}{{\rm d}t}$ $\displaystyle=\alpha x-\beta xy$
$\displaystyle\frac{{\rm d}y}{{\rm d}t}$ $\displaystyle=-\delta y+\gamma
xy~{}.$

We take the known portion of the differential equation as
$\mathcal{N_{K}}[U]=[\alpha\,x,-\delta\,y]$ for known parameters $\alpha$ and
$\delta$, and seek to learn
$F=[F_{1},F_{2}]\approx[-\beta\,x\,y,\gamma\,x\,y]$ from data only, without
knowing the target form and without knowing $\beta$ and $\gamma$.

To generate the synthetic data, $\alpha,\beta,\gamma,\delta$ were fixed at
$(1.3,0.9,0.8,1.8)$ respectively, with initial conditions at
$(x_{0},y_{0})=(0.44249296,4.6280594)$ just as in rackauckas2020universal .
The time interval was chosen as $[0,3]$ and stayed the same throughout every
LV experiment.

An ODE solver was used to generate data satisfying the LV equations. This
yields a set of points $\\{t_{i},x_{i},y_{i}\\}$. Then, Gaussian noise is
added to each $x_{i}$ and $y_{i}$. Given a particular noise level $\epsilon$,
Gaussian noise was added to the data as follows:

$\displaystyle(x_{i})_{noise}$ $\displaystyle=x_{i}+\epsilon\cdot\bar{x}\cdot
N(0,1)$ $\displaystyle(y_{i})_{noise}$
$\displaystyle=y_{i}+\epsilon\cdot\bar{y}\cdot N(0,1)$

where $\bar{x}$ denotes the element-wise mean of $x_{i}$ over all $i$
(similarly for $y$).

First, we demonstrate our approach on noise-free data (Table 3) and data with
$\epsilon=5\times 10^{-3}$ noise (Table 3) for various values of $n$ (number
of data points) and $n_{P}$ (number of collocation points). We want to show
how the hard-to-acquire data can be augmented by taking more collocation
points which require no experiments/measurements and come at only the cost of
increased computing power. We see that, in contrast to a standard PINN
approach, we need to provide more data than just the initial condition.
However, even with very sparse measurement data, we can acquire a good
discovery by only increasing the number of collocation points. The additional
benefit gained from increasing the collocation points is only realized when
there is already ample enough experimental data for the algorithm to leverage.

$n$ $n_{P}$ | $10^{2}$ | $10^{3}$ | $10^{4}$
---|---|---|---
$1$ | $2\times 10^{1}$ | $2\times 10^{1}$ | $2\times 10^{1}$
$5$ | $9\times 10^{-4}$ | $1\times 10^{-3}$ | $9\times 10^{-4}$
$10$ | $2\times 10^{-4}$ | $4\times 10^{-5}$ | $5\times 10^{-6}$
Table 1: Noise Free Data

$n$ $n_{P}$ | $10^{2}$ | $10^{3}$ | $10^{4}$
---|---|---|---
$1$ | $2\times 10^{1}$ | $2\times 10^{1}$ | $2\times 10^{1}$
$5$ | $6\times 10^{-2}$ | $4\times 10^{-3}$ | $5\times 10^{-3}$
$10$ | $1\times 10^{-3}$ | $6\times 10^{-4}$ | $8\times 10^{-4}$
Table 2: Noisy ($\epsilon=5\times 10^{-3}$) Data

Table 3: Tables demonstrating the MSE between $F$ and the true hidden target
after training for various values of $n$ and $n_{P}$.

Next, we compare our method’s performance to the UDE method. We test the two
methods on noiseless sparse data (1(a)) and on noisy data (Fig 1(b)). The
error is computed as a mean squared error (MSE) taken with respect to the true
interaction. At minimal noise level, the UDE approach and PINN approach
perform similarly and for the densest data UDEs slightly outperform PINNs.
Although increasing either noise or sparsity degrades the performance of both
methods, the PINN method consistently attains a lower MSE compared to the UDE
method as noise or sparsity increases.

(a) Sparse data

(b) Noisy data

Figure 1: Mean squared error (MSE) of the recovery of the true interaction in
comparison between the UDE and the PINN method. The spacing parameter
determines how much time passes between datapoints, but the overall time
interval $[0,3]$ remains the same.

Figures 2 and 3 show the surrogate solution and hidden terms as recovered by
the UDE and PINN methods. The noise level of the noisy data was set at 0.1
and, for the noiseless sparse data, there were 5 points each 0.6 units apart.
It is clear that the PINN approach is quite robust to noise and performs well
in low-data regimes. The UDE approach performs reasonably on sparse data, but
is not robust to noise.

(a) Noisy data

(b) Sparse data

Figure 2: UDE method performance

(a) Noisy data

(b) Sparse data

Figure 3: PINN performance

Finally, AI Feynman symbolic regression is run on the neural network output
from both our approach and the UDE approach. In all cases, the NNs
approximating $F$ were given the training data as input. Then each NN’s output
was subsequently given to AI Feynman to find the best functional form. This
data is presented in Table 4. In cases of both sparse and noisy data, AI
Feynman correctly recovers the hidden interaction terms more often for our
method than it does for the UDE method. If a formula is recovered for both
methods, the one recovered for the PINN method is often more accurate.

spacing | noise level |  $F_{1}$ (UDE) |  $F_{1}$ (PINN) |  $F_{2}$ (UDE) |  $F_{2}$ (PINN)
---|---|---|---|---|---
0.1 | 0 | -0.901 (2.8e-7) | – | 0.802 (1.1e-6) | 0.797 (2.5e-6)
0.2 | 0 | – | – | 0.797 (3.4e-6) | 0.799 (3.8e-7)
0.3 | 0 | – | -0.897 (4e-6) | – | 0.798 (1.9e-6)
0.4 | 0 | – | -0.888 (8.2e-5) | – | 0.797 (5.2e-6)
0.5 | 0 | – | -0.889 (8.9e-5) | 0.760 (1e-3) | –
0.6 | 0 | -0.892 (4e-3) | -0.890 (1e-5) | – | 0.800 (1e-32)
0.1 | 8e-3 | -9.25 (1.8e-3) | -0.906 (1e-5) | – | 0.798 (3e-5)
0.1 | 1e-2 | – | -0.911 (3.45e-5) | 0.791 (2.3e-5) | 0.777 (1.5e-4)
0.1 | 3e-2 | – | -0.960 (1e-3) | – | 0.777 (1.5e-4)
0.1 | 5e-2 | – | – | – | 0.740 (1.1e-3)
0.1 | 8e-2 | – | – | – | –
0.1 | 1e-1 | – | – | 0.887 (2e-3) | –
Table 4: Coefficients (with MSE) recovered by AI Feynman from the
approximations $F_{1}$ and $F_{2}$, comparing over datasets (rows) and method
of finding $F_{1}$ and $F_{2}$ (columns). True coefficients are -0.9 for
$F_{1}$ and 0.8 for $F_{2}$. A dash indicates AI Feynman did not recover the
functional form $Cxy$. The best performance is in bold.

The terms $\gamma\,x\,y$ and $-\beta\,x\,y$ in the LV equations correspond to
the predator’s uptake function in the ecological model. This represents the
predators’ feeding habits as a function of prey population and its resulting
effect on both the prey population and the predator’s population. The actual
form of these functions can take various forms in predator-prey models (see,
for instance, harrison ; bolger2020predator ). While we initially modelled
this as two unknown, decoupled functions $F_{1}$ and $F_{2}$ and learned them
independently, we could also have modeled them by a single function with an
additional learned parameter as a scaling factor. That is, we could take
$F_{1}=-\phi\,F_{2}$ and then only explicitly learn $F_{2}$ and a single
parameter $\phi$. This results in regressions that are near identical to the
ones presented above, but showcases an important modelling methodology that
our method is amenable to and, for more complicated models than LV, may be
necessary in order to achieve a high-quality regression.

###  4.2 Cell Apoptosis Model

We also test the method on the Q1 cell apoptosis model from wee2006akt . This
is an ODE with three variables, serine-threonine kinase $Akt_{s}$ (active
Akt), $Akt$ (inactive Akt) and tumour suppressor protein $p53$. p53 promotes
cell apoptosis, or programmed cell death, and Akt inhibits it. All parameter
values were taken from the paper.

$\displaystyle v_{0}$ $\displaystyle=k_{0}$ $\displaystyle v_{1}$
$\displaystyle=k_{1}\cdot Akt\cdot(j_{1}+Akt_{s})$ $\displaystyle v_{m1}$
$\displaystyle=k_{m1}\cdot Akt_{s}/(j_{m1}+Akt_{s})$ $\displaystyle v_{2}$
$\displaystyle=k_{2}\cdot Akt_{s}\cdot p53/(j_{2}+p53)$ $\displaystyle v_{m3}$
$\displaystyle=k_{m3}\cdot p53\cdot Akt_{s}/(j_{m3}+Akt_{s})$
$\displaystyle\frac{dp53}{dt}$ $\displaystyle=v_{0}-v_{2}-k_{d}\cdot p53$
$\displaystyle\frac{dAkt_{s}}{dt}$ $\displaystyle=v_{1}-v_{m1}-v_{m3}$
$\displaystyle\frac{dAkt}{dt}$ $\displaystyle=\frac{-dAkt_{s}}{dt}$

With the initial condition $(p53,Akt_{s},Akt)=(0.248,0.0973,0.0027)$ and 30
noiseless datapoints, both the $v_{1}$ and $v_{2}$ interactions were learned
to a high degree of accuracy. In Fig 4 it can be seen that although the
general shape does not match the true interaction 100%, the mean squared error
between the true interaction and the learned function is in fact very small
(on the order of $10^{-4}$) and the surrogate solution fits the data very
well. As is done in yazdani2020systems , if this method were augmented to
handle very small and very large values through scaling, more accurate
learning of the interaction would be possible. However, cases arise where the
interaction itself is not identifiable and several different interactions
learned by $F$ can lead to the minimization of the loss.

(a) Learning $v_{1}$

(b) Learning $v_{2}$

Figure 4: Learning cell apoptosis components with PINN approach

###  4.3 Viscous Burger’s Equation

Finally, our method is easily applied to PDEs (as in the original PINN
implementation). Here we present the discovery of both the solution to the PDE
where the underlying hidden dynamics of the operator were partially hidden.
This reconstruction used only noisy ($\epsilon=5\times 10^{-3}$) data obtained
from two time points (the initial condition, $t=0$, and a later time at
$t=0.5$). While this method can be used to discover the form of the boundary
condition as well, here we assume that the homogeneous Dirichlet boundary
conditions are known. The PDE in question is

$\frac{\partial u}{\partial t}=-u\,\frac{\partial u}{\partial
x}+\nu\,\frac{\partial^{2}u}{\partial
x^{2}},\quad\nu=\frac{1}{1000\,\pi},\quad u(x,0)=-\sin(\pi\,x)$

Here we took $\mathcal{N_{K}}=\nu\,u_{xx}$ and let the algorithm learn the
hidden term $-u\,u_{x}$. To do this, we gave the $F$ network $u$, $u_{x}$, and
$u_{t}$ as inputs. This represents an inductive prior where we are assuming
that the hidden term depends on first order and lower derivatives of the
solution. In our approach, such a prior is necessary (that is, the algorithm
cannot learn what order of derivatives to include or not include, it can
merely choose which inputs presented to it to utilize). For collocation data
we used $n_{P}=10^{4}$ and $n_{B}=10^{2}$ points sampled from the appropriate
parts of the domain $[-1,1]\times[0,1]$ via Latin hypercube sampling. The PDE
solution was reconstructed with MSE of $3\times 10^{-4}$ and the hidden term
was discovered with MSE of $2\times 10^{-2}$. The resulting solution is
visualized in Figure 5.

Figure 5: The reconstructed solution of Burgers’ equation. The two vertical
dashed white lines indicate the noisy experimental data that were sampled for
the algorithm.

##  5 Conclusion

In conclusion, our approach is able to recover, with a great degree of
accuracy, the symbolic functional form of hidden terms within a differential
operator using very sparse measurements of noisy data by utilizing a
modification of PINNs. This approach is robust to both noise and sparsity of
the data by increasing the number of collocation points (an operation that
doesn’t require any additional experimentation, just stronger compute
capacities). This approach can be applied to discovering the functional form
of an unknown ordinary differential equation (ODE) as well as both the
functional form of a partial differential operator in a partial differential
equation (PDE) and unknown terms in the boundary condition of a PDE. Although
PINNs have been noted to perform sub-optimally on stiff equations without
modification stiff1 ; stiff2 , we have noted promising results in this
direction. However, more investigation is needed.

## References

  * [1] Zdenek Belohlav, Petr Zamostny, Petr Kluson, and Jiri Volf.  Application of random-search algorithm for regression analysis of catalytic hydrogenations.  The Canadian Journal of Chemical Engineering, 75(4):735–742, 1997\.
  * [2] Alan A Berryman.  The orgins and evolution of predator-prey theory.  Ecology, 73(5):1530–1535, 1992.
  * [3] Christopher M Bishop and Nasser M Nasrabadi.  Pattern recognition and machine learning, volume 4.  Springer, 2006.
  * [4] Tedra Bolger, Brydon Eastman, Madeleine Hill, and Gail Wolkowicz.  A predator-prey model in the chemostat with holling type ii response function.  Mathematics in Applied Sciences and Engineering, 1(4):333–354, 2020\.
  * [5] Souvik Chakraborty.  Transfer learning based multi-fidelity physics informed deep neural network.  CoRR, abs/2005.10614, 2020.
  * [6] Gary W Harrison.  Global stability of predator-prey interactions.  Journal of Mathematical Biology, 8(2):159–171, 1979.
  * [7] Weiqi Ji, Weilun Qiu, Zhiyu Shi, Shaowu Pan, and Sili Deng.  Stiff-pinn: Physics-informed neural network for stiff chemical kinetics.  The Journal of Physical Chemistry A, 125(36):8098–8106, 2021.
  * [8] Pierre-Alexandre Kamienny, Stéphane d’Ascoli, Guillaume Lample, and François Charton.  End-to-end symbolic regression with transformers.  arXiv preprint arXiv:2204.10532, 2022.
  * [9] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [10] James Lu, Brendan Bender, Jin Y Jin, and Yuanfang Guan.  Deep learning prediction of patient response time course from early data via neural-pharmacokinetic/pharmacodynamic modelling.  Nature machine intelligence, 3(8):696–704, 2021.
  * [11] James Lu, Kaiwen Deng, Xinyuan Zhang, Gengbo Liu, and Yuanfang Guan.  Neural-ode for pharmacokinetics modeling and its advantage to alternative machine learning models in predicting new dosing regimens.  Iscience, 24(7):102804, 2021.
  * [12] Ben McKay, Mark J Willis, and Geoffrey W Barton.  Using a tree structured genetic algorithm to perform symbolic regression.  In First international conference on genetic algorithms in engineering systems: innovations and applications, pages 487–492. IET, 1995\.
  * [13] Xuhui Meng and George Em Karniadakis.  A composite neural network that learns from multi-fidelity data: Application to function approximation and inverse pde problems.  Journal of Computational Physics, 401:109020, 2020.
  * [14] Christian Moya and Guang Lin.  Dae-pinn: A physics-informed neural network model for simulating differential-algebraic equations with application to power networks.  arXiv preprint arXiv:2109.04304, 2021.
  * [15] Allan Pinkus.  Approximation theory of the mlp model in neural networks.  Acta Numerica, 8:143–195, 1999.
  * [16] Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Dominic Skinner, Ali Ramadhan, and Alan Edelman.  Universal differential equations for scientific machine learning.  arXiv preprint arXiv:2001.04385, 2020.
  * [17] M. Raissi, P. Perdikaris, and G.E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [18] Peter Ruoff, Melinda K Christensen, Jana Wolf, and Reinhart Heinrich.  Temperature dependency and temperature compensation in a model of yeast glycolytic oscillations.  Biophysical Chemistry, 106(2):179–192, 2003.
  * [19] Michael Schmidt and Hod Lipson.  Distilling free-form natural laws from experimental data.  Science, 324(5923):81–85, 2009.
  * [20] Steven H Strogatz.  Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering.  CRC press, 2018.
  * [21] Silviu-Marian Udrescu and Max Tegmark.  Ai feynman: A physics-inspired method for symbolic regression.  Science Advances, 6(16):eaay2631, 2020.
  * [22] Mojtaba Valipour, Maysum Panju, Bowen You, and Ali Ghodsi.  Symbolicgpt: A generative transformer model for symbolic regression.  In Preprint Arxiv, 2021.  Under Review.
  * [23] Sifan Wang, Yujun Teng, and Paris Perdikaris.  Understanding and mitigating gradient flow pathologies in physics-informed neural networks.  SIAM Journal on Scientific Computing, 43(5):A3055–A3081, 2021.
  * [24] Keng Boon Wee and Baltazar D Aguda.  Akt versus p53 in a network of oncogenes and tumor suppressor genes regulating cell survival and death.  Biophysical journal, 91(3):857–865, 2006.
  * [25] Alireza Yazdani, Lu Lu, Maziar Raissi, and George Em Karniadakis.  Systems biology informed deep learning for inferring parameters and hidden dynamics.  PLoS computational biology, 16(11):e1007575, 2020.