2021 JEL Classification]A14,C11, C53,O52 [1,3]Krzysztof Rusek 1] AGH
University of Science and Technology, Krakow, Poland 2] Jan Kochanowski
University of Kielce, Kielce, Poland 3] Barcelona Neural Networking Center,
Universitat Politècnica de Catalunya, Barcelona, Spain

11footnotetext: This work has been submitted to Scientometrics and is under
review process

# Bayesian inference of spatial and temporal relations in AI patents for EU
countries

<EMAIL_ADDRESS>Agnieszka Kleszcz<EMAIL_ADDRESS>Albert Cabellos-
Aparicio<EMAIL_ADDRESS>[  [  [

###### Abstract

In the paper, we propose two models of Artificial Intelligence (AI) patents in
European Union (EU) countries addressing spatial and temporal behaviour. In
particular, the models can quantitatively describe the interaction between
countries or explain the rapidly growing trends in AI patents. For spatial
analysis Poisson regression is used to explain collaboration between a pair of
countries measured by the number of common patents. Through Bayesian
inference, we estimated the strengths of interactions between countries in the
EU and the rest of the world. In particular, a significant lack of cooperation
has been identified for some pairs of countries. Alternatively, an
inhomogeneous Poisson process combined with the logistic curve growth
accurately models the temporal behaviour by an accurate trend line. Bayesian
analysis in the time domain revealed an upcoming slowdown in patenting
intensity.

###### pacs:

[

###### keywords:

Artificial intelligence, Patent cooperation network, European Union

##  1 Introduction

Cooperation between the member states on issues of common interest was one of
the reasons to establish the European Union. The European Union (EU) is an
international organization of contrasting geo-political concepts. On the one
hand it resembles a federation, yet on the other the member states have large
independence and freedom of foreign cooperation. Patent interactions are a
form of cooperation. In this study, we attempt to quantify internal and
external EU patent interactions as well as the European trend in the field of
Artificial Intelligence. Although there is no internationally agreed
definition of Artificial Intelligence it could be generally defined as
intelligence demonstrated by machines. A more formal definition from Copeland
(2020) states that it is “the ability of a digital computer or computer-
controlled robot to perform tasks commonly associated with intelligent
beings”. AI as a concept was first created in the 1950s but its market
relevance has been recognised worldwide in the last decade, mainly due to the
development of high-performance parallel computing chips and the availability
of large datasets that have extended this technology’s applicability Leusin .
(2020).

Modern AI requires two enablers: data and computational power. While the
technology is available to everyone, the groundbreaking work is usually done
by big tech companies that have all the required resources and skilled
employees at their disposal. Those employees usually come from academia and
often work in remote branch offices. This makes the top-tier AI application a
group effort involving multiple countries. Such technologies are an important
driver of economic growth at national and regional level of competitiveness
Klinger . (2018). The technology will significantly contribute to improvements
in human welfare across a wide range of sectors, including among others
healthcare, education, transportation, robotics, public safety, employment or
entertainment (among others). AI could also lead to groundbreaking discoveries
as e.g. predicting protein structure. This development could help supercharge
the discovery of new drugs to treat disease, alongside other applications
Senior . (2020).

The growth in recent years in AI is a topic of interest in multiple
publications e.g. Cioffi . (2020); Goodell . (2021); Niu . (2016). Having said
that we point out that AI can also adversely affect sustainable development
but this does not happen frequently. Vinuesa . (2020) analysed how AI might
impact all aspects of sustainable development in terms of 17 Sustainable
Development Goals (SDGs) and its 169 targets internationally agreed in the
2030 Agenda for Sustainable Development. For that study, the authors divided
the SDGs into the following three main groups: Society, Economy, and
Environment. The results show that 67 targets (i.e., 82%) in the Society main
group could potentially benefit from AI-based technologies. In Economy 42
targets benefit from AI ( i.e., positive impact in 70% of targets). For the
Environment group, results showed that AI could act as an enabler in 93% of
targets.

In an increasingly knowledge‐driven economy, society invariably needs creative
or inventive ideas or concepts to improve existing features and add/develop
useful new features to products. All the positive aspects of AI are accessible
for the customers as a service or product developed by a company and usually
protected by a patent. Patents are an effective method of measuring the
potential for specific technology areas. Patents are also a key measure of
innovation output, as patent indicators reflect the inventive performance of
countries, technologies and firms. Furthermore, patents are a long-term
investment. The patent applicant can commercialize the invention at any point
during that time, either through developing products or services incorporating
the patented technology or by licensing it to others. For scientometrics,
patent literature and patent information are important sources of information
and critical guidance for scientific research, business operations and
technological innovation Tsay Liu (2020); H. Wu . (2021).

AI has been one of the key drivers of the massive increase in industrial
revolution-related patenting over the past decade. Many of the biggest
technology companies have invested heavily into AI related research and
development. For example, Google, Microsoft, IBM, and Samsung have each
submitted thousands of patent applications Liu (2020). The scale of AI patents
is further confirmed by _OECD.Stat_ OECD.Stat (2021) where patent counts are
provided for selected technology areas, and technologies related to AI are
distinguished as separate technology domains and IPC classes. In general the
identification of all the AI patents could be demanding, as the topic is broad
and evolves.

In 2019, World Intellectual Property Organization (WIPO) emphasized that
despite the availability of information in patent documents, it can be
difficult to identify exactly which patent families relate to AI because of
the lack of a standardized definition WIPO (2019). The literature proposes a
variety of strategies for identifying AI-related documents, including the use
of predefined classes based on patent classification schemes e.g. Fujii Managi
(2018), the use of specific keywords e.g. Aristodemou Tietze (2018); Niu .
(2016), or even both. Both strategies have pros and cons, and outcomes from
their use could bring slightly different results Leusin . (2020). Since the
field of AI is dynamically developing the proposed keywords are quickly
becoming outdated e.g set of keywords by IPO (2019) is missing many keywords
characteristic of contemporary AI technologies e.g. TensorFlow, JAX,
generative model or attention model. Thus we relied on predefined patent
properties from IPO (2019). These are more formalized and detailed thus should
better represent the true aspect of the invention. Another important property
of each patent is the country of its inventors (one or more). So a patent with
multiple countries can be interpreted as a relation between countries.

In a globalized world, international cooperation plays a vital role in
tackling global issues including research and science. Hervás-Oliver . (2021)
showed a positive effect on SMEs innovation and R&D due to collaboration with
others. Similarly the European Commission recognises the importance of global
level cooperation to grand societal challenges, hence wants Horizon Europe to
open up eligibility to strong science countries (e.g., Canada, Australia,
etc.). The problem concerning innovation in Europe was highlighted in the
literature. One example is that frequently, the EU-13 (countries that joined
the EU in and after 2004) were found at the lower end of participation
rankings concerning Horizon 2020. Abbott Schiermeier (2019) identified that
the majority of nations with smaller research budgets are former communist
countries in central and eastern Europe, which together with Cyprus and Malta
joined the EU after 2004. Further research into the collaboration structure
between countries indicated the strongest collaboration research network was
between EU-15 (countries that entered the EU before 2004) particularly between
Germany, France, United Kingdom, Italy and Spain.

Besides grouping countries, the network of relations can also reveal their
strengths as reported by B. Wu (2020). This particular analysis focuses on
collaborations on AI-related papers in journals tracked by the Nature Index.
Cooperation based on publications can be biased toward academic interaction.
The best universities in Europe are based in United Kingdom and the findings
of B. Wu (2020) are that three United Kingdom institutions are among the
biggest collaborators in AI research in Europe. Companies may not be willing
to publish their results in scientific journals though they definitely will
benefit from patents. We argue that a patent cooperation network is a better
tool for the job when it comes to cooperation in the industry. We further
argue that analysis of such a cooperation can give a better picture of
technology development trends You . (2017).

Patent cooperation is a common business activity and patents are often owned
by multiple assignees. The patents that have multiple different assignees are
hereafter referred to as _cooperation patents_. The cooperative relationship
formed between the assignee is a channel for the transmission of knowledge,
information, and resources. The network relationship structure also determines
the action opportunities and results of the assignee Tsay Liu (2020).

Considering all aforementioned aspects and the importance of the role of
patents, it is crucial to analyse cooperation between assignees, especially in
the advanced, developing field of AI technologies. It is also worth noting the
importance of AI for EU competitiveness is highlighted in many actions i.e.,
in strategies, regulations, plans, investment and funding programs. Although
there are some publications about patent cooperation analysis (including AI
patents between assignees) there is a lack of reliable cooperation analysis
between EU countries concerning AI patents - including investigation on
cooperation between the Member States and the influence of cooperation on
patents performance.

The closest related work is by Tsay Liu (2020) who utilised the Derwent World
Patents Index, by using the patent metric and basic graph properties like
centrality to investigate the global cooperative network structure of the
assignee’s in AI patents. In this paper, we build upon that work and extend
the research in the following areas:  1. The simple graph theory based
quantities like centrality was replaced with a structure model based state-of-
the-art research in graph neural networks Battaglia . (2018).  2. Model
parameters and their uncertainties are inferred from the observations by the
Bayesian approach for which the prior distributions were proposed.  3. The
relation type or its lack was inferred and interpreted.  4. A parametric trend
model based on point process Goulding . (2016) and logistic growth
Balakrishnan (1991) was proposed for the total amounts of patents. The prior
and posterior distribution of the parameters are discussed.

To the best of the authors’ knowledge there is no other work relating to the
discovery of hidden interactions of AI patent data in EU countries. The same
applies to trend analysis and parametric forecasts. The topic discussed mostly
in technical reports like WIPO (2019) are simple visual plots which are used
to express the rapid growth of AI. We noticed that the majority of trend
analyses presented in scientific journals are based on year to year percentage
growth or simple linear trends estimated from data without uncertainty
estimation or thorough statistical analysis for discrete numbers of patents
e.g. Klinger . (2018); Tsay Liu (2020). In this paper we seek to address this
gap by providing a coherent spatial and temporal description of AI patents in
terms of the Poisson distribution.

The next section describes methods and data acquisition and processing steps
used in the research. Section 3 presents our findings in numerical experiments
and section 4 concludes the paper.

##  2 Materials and methods

Patent information could be obtained from different sources. Starting from
publicly available databases like EPO, OECD or Google Patents to the more
advanced commercial ones like Derwent, PATSTAT or Global Patent Index. In the
following research, we utilized the Global Patent Index database. There is
also a wide variety of possibilities for searching patents on a given topic,
which includes the use of keywords and the use of patent classification
schemes. In our analysis we utilised Cooperative Patent Classification. Figure
1 presents hierarchical structure of some exemplary AI patents which could be
considered as a tree graph. As a example an interpretation of hierarchical
structure of subgroup “H04L25/0254” is as follows:“H” - electricity; “H04” -
electric communication technique; “H04L” transmission of digital information,
e.g. telegraphic communication; “H04L 25/00” - adapted for orthogonal
signalling; “H04L25/0254” using neural network algorithms.

Figure 1: Scheme of hierarchical structure of Cooperative Patent
Classification

###  2.1 Data description

The research examines the cooperation network in European Union countries
referred to their code names: CY- Cyprus, MT- Malta, LV- Latvia, LT-
Lithuania, EE- Estonia, HR- Croatia, BG- Bulgaria, SK- Slovak Republic, RO-
Romania, SI- Slovenia, CZ- Czech Republic, HU- Hungary, PL- Poland, namely the
so-called EU-13 and LU- Luxembourg, PT- Portugal, EL- Greece, IE- Ireland, DK-
Denmark, ES- Spain, BE- Belgium, FI- Finland, AT-Austria, SE- Sweden, NL-
Netherlands, IT- Italy, GB- United Kingdom, FR- France, DE- Germany,
collectively named EU-15 (because the period considered was from 1990 the
United Kingdom was also included in our analysis despite it having left the EU
in 2020). To identify relevant and comparable techniques, we used 76 subgroups
(see Level 5 on figure 1) of the Cooperative Patent Classification (CPC)
scheme proposed in the Intellectual Property Office report, see p.32 in IPO
(2019), to identify AI. To extract AI patents for EU countries in the time
period from 1990-2021 (18th of September) we formulate the query in the Global
Patent Index (GPI 2021/28) database as detailed in appendix A. As a result, we
obtained 10,759 AI patents families (a family is a set of all patents in
different countries that protect the same invention). As a cooperation
indicator, we used _Applicant country of residence_ content. The presence of
authors from at least 2 different countries indicates cooperation. An example
patent having six applicants from [DE, LU, DE, DE, BE, BE] accounts for
cooperation between all pairs of 3 countries [LU, DE, BE]. The strength of the
cooperation is measured by counting the number of relations. From 10,759
patents families, we identified 1,644 (15.3 %) patents that have cooperation
with at least one additional country - we will refer to them as _cooperation
patents_.

The cooperation network was created as a graph with nodes representing
countries and edges representing relations. All the non-EU countries are
represented as a single entity _Others_ and the graph contains an edge only if
there was at least one relation (cooperation patent). Patents assigned to a
single country are discarded in the relation network – no self-connection is
present in the graph (however they contribute to the node feature). Having
said that, the information about total patents issued by a given country is
retained in node attributes. We utilized fractional counts applied to patents
with multiple applicants. When a patent was invented by several inventors from
different countries, the respective contributions of each country were
considered. This is done to eliminate multiple counting of such patents. For
example, a patent co-invented by 2 French, 1 Belgian and 1 German resident
will be counted as 1/2 patent for France; 1/4 for Belgium; and 1/4 patent for
Germany. Total (fractional) number of patents is a node attribute in the
relation network, while edge attribute is the total number of relations.

###  2.2 Relational regression mixture models

Let $C_{ij}\in\mathbb{N}_{0}$ be the total number of patents shared between
countries $i$ and $j$ (_cooperation patents_). Furthermore the amount of
patent of the $i$-th country will be denoted by $C_{i}\in[0,\infty)$.
Following Battaglia . (2018) we propose the relational inductive bias based on
a Poisson regression model to relate the amount of patents and the number of
interactions:

$C_{ij}\sim\text{Pois}\left(\lambda_{ij}=e^{\beta}(C_{i}C_{j})^{\alpha}\right).$
(1)

Such a model has multiple advantages: (i) The model is simple and easy for
interpretation. (ii) Permutation equivariance doesn’t require special
ordering. (iii) Setting $\alpha=1$ yields the popular gravity model (without
distance term). (iv) The canonical link function of the Poisson distribution
is logarithmic so in the $\log$ domain the model can be expressed in its
simple linear form:

$\log\lambda_{ij}=\alpha\log(C_{i}C_{j})+\beta.$ (2)

Conditional distribution of $C_{ij}$ is multimodal and classical Poisson
regression gives poor results. The problem can be overcome by using a mixture
of Poisson distributions as a conditional distribution Lamont . (2016). Such a
model is known as the regression mixture model Leisch (2004).

In this representation every pair of countries $i,j$ gets a hidden categorical
variable $Z_{ij}\in\\{0,1,2\\}$. Denoting $\lambda_{ijk}=\mathsf{E}(C_{ij}\mid
Z_{ij}=k)$ and $x_{ij}=\log(C_{i}C_{j})$ we propose the following model of the
relation intensities $\lambda$:

$\log\lambda_{ijk}=\alpha_{k}x_{ij}+\beta_{k}\quad k\in\\{0,1,2\\}.$ (3)

Equation (3) is a system of three equations parameterized by
$\alpha_{k},\beta_{k}$ defining three possible behaviours of the cooperation.
Since $Z_{ij}$ is an unknown hidden variable we model it by a categorical
probability distribution with parameters
$\pi_{k}=\mathsf{P}(Z_{ij}=k),\quad\sum_{k}\pi_{k}=1$. Hence, the conditional
distribution of $C_{ij}$ is given by:

$\mathsf{P}(C_{ij}=c\mid
x_{ij})=\sum_{k}\pi_{k}p_{\text{Pois}}(\lambda_{ijk},c),$ (4)

where $p_{\text{Pois}}(\lambda,\cdot)$ is the probability mass function of the
Poisson distribution of rate $\lambda$. Equation (4) can be used to estimate
parameters of the model
($\theta=(\alpha_{0},\alpha_{1},\alpha_{2},\beta_{0},\beta_{1},\beta_{2},\pi_{0},\pi_{1},\pi_{2})$)
by using the conditional maximum likelihood. Alternatively one can use it to
construct a Bayesian posterior distribution of the parameters. Given the
parameters $\theta$, the true counts of relational patents $y_{ij}$ and
features $x_{ij}$ one can follow Murphy (2012) and compute the posterior
distribution of the hidden variable $Z_{ij}$ as :

$\ \mathsf{P}(Z_{ij}=c\mid
x_{ij},y_{ij},\theta)=\frac{\pi_{c}p_{\text{Pois}}(\lambda_{ijc},y_{ij})}{\sum_{k}\pi_{k}p_{\text{Pois}}(\lambda_{ijk},y_{ij})}.$
(5)

This way the regression model performs also clusterisation, and allows to
infer the relation type from the patent counts.

###  2.3 Logistic growth of inhomogeneous Poisson process

Patent application and their publication can occur at any time (to within a
granularity of days and working hours). Thus, any aggregation to month or
year, results in some information being lost. A natural way of modelling such
a phenomenon is a point process Daley Vere-Jones (2006). In particular the
Poisson process is a common choice. We follow this approach with the exception
that the process cannot be stationary because a huge change in patenting
dynamics can be observed over the last decade. A realistic model of the patent
publication process is an inhomogeneous Poisson process Daley Vere-Jones
(2006). Such a process is an extension of the stationary Poisson process whose
rate changes over time.

We argue that patents are similar to other goods thus their production can be
described by the logistic curve, commonly used in economics, ecology and other
areas Twiss (1992); Balakrishnan (1991). Since patent production is random in
nature we propose to model cumulative rate function (6) (expected number of
patents produced until time $t$) instead of the cumulative patent counts.

$\Lambda(t)=\int_{-\infty}^{t}\lambda(t)\mathrm{d}t=\frac{L}{1+e^{-\frac{t-t_{0}}{s}}}$
(6)

Such a curve is parameterized by the capacity $L>0$, the midpoint $t_{0}$ and
the scale $s>0$ (an alternative and also commonly used parameterisation
involves the rate equal to $1/s$). The capacity represents the asymptotic
maximum at infinity, while the maximum rate is at the midpoint. This is also
the point in time where the growth starts slowing down. The scale controls the
width of the growth period ($10s$ contains almost 99% of the events). Those
three parameters describe the dynamics of the process totally and can be
estimated from the observations. Let $\bm{t}$ be a vector of $n$ random points
in time representing patent publication date as a real number. The log
likelihood function of those observations is given by Daley Vere-Jones (2006)
as:

$\ell(\bm{t})=\sum_{i=1}^{n}\log\lambda(t_{i})-\int_{t_{1}}^{t_{n}}\lambda(t)\mathrm{d}t.$
(7)

One can estimate the unknown parameters $\theta=(L,t_{0},s)$ by maximizing (7)
or by using $\ell(\bm{t})$ to construct a Bayesian posterior distribution of
the parameters. In this paper we followed the approach used for spatial
relation discovery and estimated parameters via Bayesian inference. Since the
likelihood function involves calculations on thousands of patents, we used
approximate yet much faster variational inference van de Schoot . (2021)

AI is not a static area of research. Over the past few decades we observed
raises and falls of various trends and methods. There are periods of low
interests in AI known as AI winter Schuchmann (2019). Starting around 2012 we
could observe the beginning of a new method known as deep learning that
quickly caught the attention of the researchers and the industry. The
beginning of a new method would definitely break the simple logistic growth
model. Indeed, in the case of patent data, we observed this in numerical
experiments. A single logistic curve doesn’t fit well to the observations.
However a superposition of two logistic curves can explain the data with much
lower error. The final cumulative patent rate can be expressed as:

$\Lambda(t)=\sum_{i=1}^{2}\frac{L_{i}}{1+e^{-\frac{t-t_{0,i}}{s_{i}}}}=L\left(\frac{p_{1}}{1+e^{-\frac{t-t_{0,1}}{s_{1}}}}+\frac{1-p_{1}}{1+e^{-\frac{t-t_{0,2}}{s_{2}}}}\right).$
(8)

The capacity of this new rate is $L=L_{1}+L_{2}$. By setting
$p_{1}=\frac{L_{1}}{L_{1}+L_{2}}$ and
$p_{2}=\frac{L_{2}}{L_{1}+L_{2}}=1-p_{1}$ we can observe that (8) corresponds
to the cumulative density function of a mixture of logistic distributions
Balakrishnan (1991) scaled by $C$ where $p_{1},p_{2}$ are mixing
probabilities. This observation allows us to use numerical procedures for the
mixture distribution in a totally different context. Furthermore $\lambda(t)$
corresponds to the scaled probability density function (pdf). This is
important for numerical estimation as $\log\lambda(t)$ can be calculated from
log pdf which is quite often provided by numerical libraries with an optimized
implementation.

###  2.4 Bayesian inference

To fit the relational model to the graph data we adopt a Bayesian approach. In
this way the uncertainty of the model is captured via a posterior distribution
of parameters van de Schoot . (2021); Murphy (2012).

For Bayesian inference, a joint prior distribution is specified for all the
unknown model parameters. Since all the quantities in the relational model are
real numbers the prior distribution is assumed to be normal. Furthermore, we
assumed a priori parameters independence. The hyperparameters for this
distribution are set using a polynomial fit in the log transform domain and
widened appropriately to avoid biasing the model towards maximum likelihood
estimates. When the experiment is repeated for time series, the Bayesian
approach has the advantage that today’s posterior may play the role of
tomorrow’s prior.

For temporal analysis, we have less informative priors. In particular, the
capacity was modelled using a log-normal distribution with a heavy tail.
Midpoint’s priors were weakly informative independent normal distributions and
the scale priors were exponential distributions - the least entropy prior.
Furthermore, the posterior was estimated using stochastic variational
inference with independent normals used as a surrogate posterior Murphy
(2012).

##  3 Results

As a result of a formulated query (see A) into our data set we observed 9
sections, 99 classes, 328 sub-classes, 460 main groups and 8,767 subgroups (it
is worth noting that single patents could be simultaneously classified as a
few sub-classes or main groups or many different subgroups). Analysing only
subgroups, classified as AI, the most popular subgroups were: G06N20/00
(computer systems based on specific computational models - Machine learning),
G06N 7/005 (computer systems based on specific mathematical models -
Probabilistic networks), G05D1/0088 (systems for controlling or regulating
non-electric variables - control of position, course or altitude of land,
water, air, or space vehicles, e.g. automatic pilot - characterized by the
autonomous decision making process, e.g. artificial intelligence, predefined
behaviours).

###  3.1 Spatial relation

We represent the cooperation network of the AI patents as a graph where nodes
are countries, and edges represent relations based on AI collaboration – see
figure. 2.

Figure 2: Cooperation network of the AI patents.

The size and colour of nodes represents the total amount of patents ($C_{i}$)
in log scale, similarly the width and colour of edges expresses the number of
collaborations also in log scale ($C_{ij}$). As can be observed, the number of
AI patents (as a sum over the years 1990-2021) among the EU member countries
varies greatly between countries. From the analysis a manifestation of the
Pareto principle is observed i.e. most patents are issued by only a few
countries. Accordingly, the greatest number of patents was assigned in
Germany, followed by the United Kingdom, France and Netherlands. While for
Slovenia, Slovakia, Lithuania and Romania there was only a single patent.
Results from the cooperation graph shows that in EU countries the greatest
number of collaboration patents was with United Kingdom. In total 586
cooperation patents include 509 with non-EU countries and 22 with France -
United Kingdom’s top partner in the EU. Germany scored 514 cooperation
patents. Out of those 411 were with non-EU countries and only 22 with France -
again the top partner in the EU. France itself has 373 cooperation patents.
Still, most of them (230) are with non-EU countries, however, France has
stronger cooperation within EU and scored 83 patents with the Netherlands -
its top partner. At the other end of the scale, a lack of cooperation with
other countries was observed in Estonia, Latvia, and Slovenia. For the
majority of the countries with cooperation patents most of their partners were
non-EU countries. In terms of percentage share with non-EU countries, the data
shows the following: Finland - 87%, United Kingdom - 86%, Germany - 80%,
Ireland - 75% and France - 61%. Since “Others” is a meta node representing
multiple countries, its cooperation power is expected to be stronger than for
a single country. Less obvious is the relation between the number of
cooperation patents to the total number of patents for each EU country.
Considering the ratio of former and later, the highest share is observed in
Bulgaria (37%) followed by Romania ( 20%), Hungary (20%) and Czech Republic
(18%) – the countries with a relatively low number of patents. We don’t
consider the Slovak Republic because there was only 1.75 patents in total- all
with different countries. On the other hand, the lowest percentage shares were
observed in Finland (5%), Ireland (4.8%), Luxembourg(4%) and Denmark (4%)
(Croatia, Estonia, Latvia and Lithuania don’t have any cooperation patents and
in Slovenia there are no AI patents registered).

These basic statistics and the cooperation network (figure 2) are definitely
informative about the EU patent network yet their interpretation can be
misleading. We observe a very high dynamic range of data: some countries have
only a few patents while others have hundreds or thousands - hence the use of
$\log$ transformation. Without rigorous statistical analysis we have no means
of assessing uncertainty in the statistics regarding the relations, especially
for the countries at the lower end of patent ranking. Furthermore, the lack of
cooperation patents between two countries does not necessarily imply the lack
of cooperation between them. There might be simply not enough patents for the
relation to be observed.

The numerical results suggest that regression mixture models allow for the
inference of both the presence and the strength of patent cooperation between
EU countries. Let us consider the network (figure 2) as a full graph with with
29 nodes (28 EU countries and “Others”) attributed with the fractional number
of patents for a given country. Edges in the graph are labeled by the total
number of interaction patents. The graph model gives the distribution of edge
features (of total number of cooperation patents) from the node features
(fractional total number of patents). We specifically propose to use a mixture
of three Poisson distributions as described in 2.2. In our Bayesian approach
to statistical inference we work with simple, slightly weakly informative
priors with one exception for $\alpha_{0}$ where we propose a deterministic
distribution at 0 so $\alpha_{0}$ de facto is not a parameter and component 0
doesn’t scale with node features.

For the numerical experiment we propose the following priors for the remaining
parameters. The slopes $\alpha_{1}$ and $\alpha_{2}$ are normal with mean 0.5
and standard deviation 0.5. This prior allows for both positive and negative
slopes with higher probability mass on the positive side. All biases $\beta$
are independent normal distributions with mean -8 and standard deviation 3.
Finally, mixing probabilities $\pi_{k}$ are parameterized in the unconstrained
domain by logits $l_{k}=\log\left(\frac{\pi_{k}}{1-\pi_{k}}\right)$. The
logits of mixing probabilities are independent centered normals with standard
deviation 2. Here we use an overparametrized yet simpler prior where all three
logits are specified, and the normalization is encoded in the model. Since the
posterior distribution is analytically intractable the inference works with
4,000 MCMC (No U-turn sampler) samples from 16 independent chains resulting in
64,000 posterior samples whose marginal distributions are presented in figure
3. Numerical values of means and standard deviations are collected in the
first two rows of table 1.

The independent chains allowed us to assess the convergence to the stationary
distribution. For all the variables the potential scale reduction (R-hat) is
very close to one, so we consider the samples to come from the stationary
distribution.

Figure 3: Marginal posterior distributions Table 1: Marginal moments of posterior distribution | $\alpha_{1}$ | $\alpha_{2}$ | $\beta_{0}$ | $\beta_{1}$ | $\beta_{2}$ | $l_{0}$ | $l_{1}$ | $l_{2}$
---|---|---|---|---|---|---|---|---
mean | 0.719 | 0.955 | -8.146 | -8.044 | -7.966 | -0.186 | 0.818 | -0.613
std | 0.046 | 0.021 | 2.873 | 0.648 | 0.291 | 1.193 | 1.173 | 1.166
2020
mean | 0.726 | 0.906 | -12.099 | -7.910 | -7.147 | -0.225 | 0.951 | -0.579
std | 0.020 | 0.010 | 2.834 | 0.275 | 0.140 | 0.510 | 0.473 | 0.463
2019
mean | 0.745 | 0.909 | -12.031 | -7.863 | -6.962 | -0.396 | 0.994 | -0.491
std | 0.032 | 0.017 | 2.909 | 0.424 | 0.217 | 0.803 | 0.745 | 0.738
2018
mean | 0.766 | 0.911 | -11.990 | -7.859 | -6.744 | -0.734 | 1.142 | -0.360
std | 0.058 | 0.030 | 2.996 | 0.742 | 0.382 | 1.328 | 1.215 | 1.206

Each posterior sample represents a single relational model consisting of three
Poisson regression models being the components of the mixture. Their
collective behaviour is presented in figure 4, where the mean of each
component is visualized independently. The confidence ribbons show 95%
credible intervals for the means obtained from the posterior samples as
empirical quantiles.

Figure 4: Relational mixture regression with the inferred relation type.

By design, the means of the last two components depend on nodes features,
while the first one is independent of the features ($\alpha_{0}=0$). Since
slope (exponent) describes how fast the average number of cooperation patents
grows with nodes features we argue that its value quantifies the strength of
the relation. Both slopes ($\alpha_{1},\alpha_{2}$) are positive and the
posterior probability of having negative slope is
$\mathsf{P}(\alpha_{i}<0)<0.00025$. This means that the increase in total
patents yields an increase in cooperation patents as well (rate is an
increasing function of features). The third component has $\alpha$ almost one
and this relation corresponds to cooperation patents with “Others”.
Interestingly this relation has almost the form of the gravitational model.
The within-EU cooperation is mostly covered by the second component. However,
the most interesting is the first component- the one whose slope was a priori
set to 0. Hence, its average does not depend on node features, and we
associate this component with the lack of long-lasting cooperation between
given countries. The reasoning is that if there were cooperation, an increase
of total patents would increase the number of cooperation patents. The
statistical model allows us to infer which pairs of countries do not show any
evidence of cooperation.

Given the model and its parameters, every pair of node and edge features can
be assigned a distribution of the hidden variable $Z_{ij}$. This way it is
possible to estimate the most probable type of relation between two countries.
This clusterisation is visualized in figure 5.

Figure 5: The most probable type of interaction between countries labeled with
$\log$ edge features. Type 0:FR-AT, NL- AT, IE-BE, DE-CY, FI- DK, FR-DK, GB-
DK, DE-EE, FR-FI, IE-FI, NL-FI, SE-FI, GB-FI, LU-DE, PT-DE, DK-IE, IT-IE. Type
2: mainly concerned EU countries with non-EU countries (Others), as well as
within the EU between AT-BE, BE-IT, BE-NL, RO -ES, FR-NL, HU-SE

For 28 EU countries and “Others”,there are 406 possible connection
combinations (everyone with everyone). The results showed that 17 pairs of
countries have low cooperation (Type 0); a cooperation level of 332 was
average (Type 1) and in 29 cases the cooperation was high (Type 2) (Slovenia
was removed because of no patents). As previously mentioned, type-2 relations
(high cooperation) mainly involved EU countries with non-EU countries
(“Others”), as well as within the EU between AT-BE, BE-IT,BE -NL, RO -ES, FR-
NL, HU-SE.

Type-1 represents most of the pairs with cooperation patents plus some pairs
without cooperation patents. In this case we argue the total amount of patents
is not large enough, and the resulting rate is so low, that a zero
interactions result has high probability. On the other hand, the lack of
cooperation patents between countries with high numbers of total patents is an
indication of a type-0 relation. This interpretation is consistent with the
results from the models at different times.

In table 1 we provide additional statistics of parameters based on time-
limited subsets of the main dataset. For example, the data in the rows for
2018, show parameters for the posterior distribution for patents published by
the end of 2018. Priors were mostly the same except for $\beta_{0}$, where we
were forced to use a more informative prior of average -10 due to greater
noise in the observations. Priors for the subsequent years were just
posteriors from the previous years. The results are stable and estimated
posterior distributions are quite similar over the years - this shows our
current estimates can be useful for further analysis in the future.

Furthermore, the posteriors are also close to the main posteriors estimated
for the whole dataset. A minimal increase of uncertainty can be observed for
all $\alpha$ this is probably because of less informative priors and noisy
observations from the pandemic year 2021. The value of $\beta_{0}$ seems to be
highly determined by the choice of the prior distribution and it makes sense
as the only evidence for it is the lack of cooperation patents. Having said
that, the case of no cooperation patents still can be grouped into their
types. Furthermore, over the three years as the number of patents grow we
observed an increase in the number of type-0 relations. As expected, the
additional observations improve the discovery of a significant lack of
cooperation.

The Bayesian analysis of interactions according to the rule that today’s
posterior is tomorrow’s prior suggests that the proposed model holds over time
and the only change we can expect is the discovery of more non-interacting
pairs of countries.

###  3.2 Temporal dynamics

Analysing the total number of AI patents over the period 1990-2021 from our
data set, we could observe that the growth rate was very dynamic over the last
5 years. However, from a simple quantitative assessment, we can’t infer a
trend and expected future changes. Although in the literature it is hard to
find reliable analysis based on statistical methods. Majority trends in
scientific publications and reports are based on percentage growth. These are
simple trends estimated from data without model regularization, or uncertainty
based on the quantitative patent data or analysing trends in most prominent
areas in AI. However, in periods, when many companies are investing in AI
development, (including investment in infrastructure, data and staff training)
it is important to have a long-term perspective on those technologies.

Patent information is a valuable resource for assessing technology trends. In
the history of AI, few technology booms were observed. The previous section
described spatial relations between countries, however, the patent landscape
is dynamic and changes rapidly over time. In this section, we discuss the
temporal behaviour of the results. Since modelling interactions between 29
elements over time would require a complicated spatio-temporal model, here we
focus on finding the global trends in AI patents.

As the patent can be published at any time we use the homogeneous Poisson
process to describe this process. Furthermore, we propose the use of a
logistic curve to model the patent production process. Such a model is simple
and more technically correct for point data (patent counts are integer numbers
– they must be modelled by a discrete probability distribution) compared to
simple trend lines estimated by the least-squares method. Also, no aggregation
means less noise in the data. As in the case of spatial relations, we use a
Bayesian posterior distribution to capture the uncertainty of the forecasts.

For the numerical experiment, we propose the following priors. The log-normal
capacity $L$ is derived from a normal distribution with an average of 26 and a
standard deviation of 4. The weakly informative prior was selected such that
most of the probability is after 2021 as $L$ is the asymptotic number of
patents at infinity. Midpoints $t_{0}$ are independent normal distributions
with their mean in January 2011 and with standard deviation 27 years(which
results in almost a century above the mean for probable locations of the
midpoints by 6 sigma rule). The scale parameters are independent exponentials
with an average of 137 years (50,000 days). Finally, the fraction parameter
$p_{1}$ (specified as logit $l_{1}$) is normal of unit average and standard
deviation. Since we used, variational inference van de Schoot . (2021); Murphy
(2012) to approximate Bayesian inference. The moments of posterior
distributions (mean-field- an independent normal approximation of the true
posterior) are collected in table 2. The resulting trend is visualized in
figure 6.

Figure 6: Cumulative rate Figure 7: Rate components Table 2: Posterior moments of trend model | $t_{0,1}$ | $t_{0,2}$ | $s_{1}$ | $s_{2}$ | $p_{1}$ | $L$
---|---|---|---|---|---|---
mean | 2015-02-12 | 2024-07-26 | 3342 days | 666 days | 0.123 | 50288
std | 187 days | 24 days | 142 days | 9 days | 0.003 | 542

The trend was estimated on all patents published before ’2020-09-01’ and the
last year of observation was used for model evaluation during the training.
The training involved maximization of the evidence lower bound (ELBO) if we
observed an increase of the expected likelihood on the validation patents
(after 2020-09-01). Both the expected likelihood and Kullback–Leibler
divergence were estimated by the Monte Carlo method with 16,384 samples from
the surrogate posterior.

The estimated trend (together with its 95% credible interval) is compared to
the cumulative number of patents up to a given time. Except for the period in
the early 90’s we observe high accuracy of the trend line. This trend line
extrapolates patenting cumulative rate in the future until 2030. We can expect
the cumulative rate to begin saturating after 2025. Further analysis is
difficult from such a plot and the detailed rate functions depicted in figure
7 reveal the forecasting behaviour. The rates correspond to the two components
of the logistic growth curves we used in the model. Given the times and the
behaviour of the rate, we can speculate about their nature. We can observe
that $\lambda_{2}$ took off after 2012 - a time we consider as the beginning
of the deep learning era, see the dotted line (alexnet) marking the
publication of Krizhevsky . (2012). This suggests that $\lambda_{1}$ must
correspond to the remaining AI technologies - shallow learning. The rate of
this part of patents has already begun slowing down around 2014/2015. We
cannot be more precise as the credible interval for the corresponding midpoint
(see $t_{0,1}$ in table 2 ) is over two years wide ($\pm$ two standard
deviations). Our trend predicts that deep learning will also begin slowing
down around the third quarter of 2024 after reaching the peak rate of over 17
patents per day. The credible interval, in this case, is much narrower (of the
order of months). The estimated capacity is above 50,000 patents so the
current number of patents is expected to increase 5 times. Over 87 %
(1-$p_{1}$) of them can be attributed to the deep learning that will last
around 20 years from its beginning. The shallow learning part will last a lot
longer and its overall time span can be estimated to be 100 years. This again
shows a Pareto rule: 87 % of all patents are produced during the period
spanning only 20% of the AI time.

##  4 Conclusion

AI technologies are currently one of the most rapidly developing areas of ICT,
which could be easily used in all areas of technology, solving many problems
in healthcare, automotive, industry and many others. This article analyzes AI
patents produced in the years 1990-2021 in the EU countries using the EPO
patent database.

It is not a surprise that the strongest economies in the EU have produced most
of the AI patents. The country with the highest number was Germany, followed
by United Kingdom, France, the Netherlands, Italy, Sweden and Finland. The
numbers also support the intuition that the countries having a large number of
patents will most likely have a large number of common patents. In particular,
the number of cooperation patents depends on the product of patents of
individual countries by a power law. When non-EU countries are considered as a
single meta country “Others” we can identify significantly different types of
relationships with “Others”. With regard to most EU countries, it can be
generalized that a significant amount of patent cooperation on AI took place
with countries outside the EU.

The proposed Bayesian model additionally enables the potential to distinguish
between types of interactions between different countries, ordered by the
strength of these relations. The key insight from the paper is that the
absence of common patents doesn’t imply that there is no cooperation between
countries especially those with a low total number of patents. Bayesian
inference allows us to distinguish lack of cooperation versus the alternative
of low cooperation where the probability of interaction patent is low but not
zero. We showed that 17 pairs of countries had no cooperation and the results
hold for four years of analysis. As time goes by and the total number of
patents grows, more pairs are assigned to the none cooperative group. It’s
worth noting that future studies can become our posterior distributions
according to the principle that today’s posterior is tomorrow’s prior. Having
said that we point out that the relationship model is static and does not
account for the dynamics of change.

When it comes to temporal dynamics the trend in AI patent application over
three decades can be accurately explained by the sum of two logistic curves.
By utilising the historical analysis of machine learning publications, the two
components of the trend were associated with two trends in the field. The
first is a family of classic algorithms used mainly in the 90s, the second
corresponds to deep learning methods. Since the emergence of the growing trend
of deep learning, the trend of traditional methods has also been gradually
diminishing. We suspect that further detailed analysis of patents
categorization would confirm our interpretation. One of the curves has already
reached its peak rate in the 2010s and the second one will begin slowing down
in the mid-2020s. Having said that the total curve has many decades till
saturation.

We would like to emphasize that such a forecast can be biased by the pandemic
years 2020-2021 so could be seen as a limitation of our model. Having said
that the model could help to assess what would have happened to AI patents had
the pandemic not happened.

Another limitation is that we ignore the possibility that another technology
(e.g. efficient probabilistic modelling) could emerge just like deep learning
did and start a new growth curve. Historical research on AI is characterized
by cycles of growth and decline over several years. The emergence of new
technologies causes a sharp increase in interest in the new topic; therefore,
the emergence of a new method may result in a further increase in AI patents
based on the new technology, beyond the prediction from our developed model.
Such a paradigm shift is impossible to predict without strong assumptions
about breakthrough occurrence in time. One aspect where the model could be
improved is the forecast uncertainty. Since in its current form the trend
model does not account for increased uncertainty of the long term forecast,
the credible intervals of the forecast are expected to be too narrow. Future
work could solve this issue by modelling logistic growth using a differential
equation with additional random components e.g. Wiener process.

Acknowledgments The research was supported in part by PL-Grid Infrastructure,
POWER 2014-2020 program and the Polish Ministry of Science and Higher
Education with the subvention funds of the Faculty of Computer Science,
Electronics and Telecommunications of AGH University.

## References

  *   * Abbott Schiermeier (2019) Abbott2019Abbott, A. Schiermeier, Q.  2019may.  How European scientists will spend €100 billion How European scientists will spend €100 billion.  Nature5697757472–475.
  * Aristodemou Tietze (2018) ARISTODEMOU201837Aristodemou, L. Tietze, F.  2018\.  The state-of-the-art on Intellectual Property Analytics (IPA): A literature review on artificial intelligence, machine learning and deep learning methods for analysing intellectual property (IP) data The state-of-the-art on intellectual property analytics (ipa): A literature review on artificial intelligence, machine learning and deep learning methods for analysing intellectual property (ip) data.  World Patent Information5537-51.  https://www.sciencedirect.com/science/article/pii/S0172219018300103  Advanced Analytics of Intellectual Property Information for TechMining
  * Balakrishnan (1991) balakrishnan1991handbookBalakrishnan, N.  1991\.  Handbook of the Logistic Distribution Handbook of the logistic distribution.  CRC Press.
  * Battaglia . (2018) gnnBattaglia, P., Hamrick, J.B.C., Bapst, V., Sanchez, A., Zambaldi, V., Malinowski, M.Pascanu, R.  2018\.  Relational inductive biases, deep learning, and graph networks Relational inductive biases, deep learning, and graph networks.  arXiv.  https://arxiv.org/pdf/1806.01261.pdf
  * Cioffi . (2020) cioffi2020artificialCioffi, R., Travaglioni, M., Piscitelli, G., Petrillo, A. De Felice, F.  2020\.  Artificial intelligence and machine learning applications in smart production: Progress, trends, and directions Artificial intelligence and machine learning applications in smart production: Progress, trends, and directions.  Sustainability122492.
  * Copeland (2020) CopelandCopeland, B.  2020\.  artificial intelligence artificial intelligence.  Encyclopedia Britannica..
  * Daley Vere-Jones (2006) daley2006introductionDaley, D. Vere-Jones, D.  2006\.  An Introduction to the Theory of Point Processes: Volume I: Elementary Theory and Methods An introduction to the theory of point processes: Volume i: Elementary theory and methods.  Springer New York.  https://books.google.pl/books?id=6Sv4BwAAQBAJ
  * Fujii Managi (2018) FUJII201860Fujii, H. Managi, S.  2018\.  Trends and priority shifts in artificial intelligence technology invention: A global patent analysis Trends and priority shifts in artificial intelligence technology invention: A global patent analysis.  Economic Analysis and Policy5860-69.  https://www.sciencedirect.com/science/article/pii/S0313592617302539
  * Goodell . (2021) GOODELL2021100577Goodell, J.W., Kumar, S., Lim, W.M. Pattnaik, D.  2021\.  Artificial intelligence and machine learning in finance: Identifying foundations, themes, and research clusters from bibliometric analysis Artificial intelligence and machine learning in finance: Identifying foundations, themes, and research clusters from bibliometric analysis.  Journal of Behavioral and Experimental Finance32100577.  https://www.sciencedirect.com/science/article/pii/S2214635021001210
  * Goulding . (2016) Goulding2016PoissonGoulding, J., Preston, S. Smith, G.  2016\.  Event Series Prediction via Non-Homogeneous Poisson Process Modelling Event series prediction via non-homogeneous poisson process modelling.  2016 IEEE 16th International Conference on Data Mining (ICDM) 2016 ieee 16th international conference on data mining (icdm) ( 161-170).  10.1109/ICDM.2016.0027
  * Hervás-Oliver . (2021) hervas2021driversHervás-Oliver, JL., Parrilli, M.D., Rodríguez-Pose, A. Sempere-Ripoll, F.  2021\.  The drivers of SME innovation in the regions of the EU The drivers of sme innovation in the regions of the eu.  Research Policy509104316.
  * IPO (2019) IntellectualPropertyOffice2019IPO  2019\.  Artificial Intelligence A worldwide overview of AI patents and patenting by the UK AI sector Artificial Intelligence A worldwide overview of AI patents and patenting by the UK AI sector .  Intellectual Property Office.
  * Klinger . (2018) klinger2018deepKlinger, J., Mateos-Garcia, J.C. Stathoulopoulos, K.  2018\.  Deep learning, deep change? Mapping the development of the Artificial Intelligence General Purpose Technology Deep learning, deep change? mapping the development of the artificial intelligence general purpose technology.  Mapping the Development of the Artificial Intelligence General Purpose Technology (August 17, 2018).
  * Krizhevsky . (2012) krizhevsky2012imagenetKrizhevsky, A., Sutskever, I. Hinton, G.E.  2012\.  Imagenet classification with deep convolutional neural networks Imagenet classification with deep convolutional neural networks.  Advances in neural information processing systems251097–1105.
  * Lamont . (2016) mixregLamont, A.E., Vermunt, J.K. Van Horn, M.L.  2016\.  Regression Mixture Models: Does Modeling the Covariance Between Independent Variables and Latent Classes Improve the Results? Regression mixture models: Does modeling the covariance between independent variables and latent classes improve the results?  Multivariate behavioral research51135–52.  https://pubmed.ncbi.nlm.nih.gov/26881956
  * Leisch (2004) JSSv011i08Leisch, F.  2004\.  FlexMix: A General Framework for Finite Mixture Models and Latent Class Regression in R Flexmix: A general framework for finite mixture models and latent class regression in r.  Journal of Statistical Software1181–18.  https://www.jstatsoft.org/index.php/jss/article/view/v011i08
  * Leusin . (2020) LEUSIN2020101988Leusin, M.E., Günther, J., Jindra, B. Moehrle, M.G.  2020\.  Patenting patterns in Artificial Intelligence: Identifying national and international breeding grounds Patenting patterns in artificial intelligence: Identifying national and international breeding grounds.  World Patent Information62101988.  https://www.sciencedirect.com/science/article/pii/S0172219020300806
  * Liu (2020) liu_2020Liu, S.  2020\.  Global AI software market revenue by region 2018-2025. Global AI software market revenue by region 2018-2025.  [2021-12-06]https://www.statista.com/statistics/721747/worldwide-artificial-intelligence-market-by-region/
  * Murphy (2012) murphy2012machineMurphy, K.  2012\.  Machine Learning: A Probabilistic Perspective Machine learning: A probabilistic perspective.  MIT Press.
  * Niu . (2016) ijgi5050066Niu, J., Tang, W., Xu, F., Zhou, X. Song, Y.  2016\.  Global Research on Artificial Intelligence from 1990–2014: Spatially-Explicit Bibliometric Analysis Global research on artificial intelligence from 1990–2014: Spatially-explicit bibliometric analysis.  ISPRS International Journal of Geo-Information55.  https://www.mdpi.com/2220-9964/5/5/66
  * OECD.Stat (2021) oecdstatOECD.Stat  2021\.  Patents by technology. Patents by technology.  https://stats.oecd.org/viewhtml.aspx?datasetcode=PATS_IPC.  [2021-08-27]https://stats.oecd.org/viewhtml.aspx?datasetcode=PATS_IPC&lang=en#
  * Schuchmann (2019) phdthesisaiwinterSchuchmann, S.  2019\.  Analyzing the Prospect of an Approaching AI Winter Analyzing the prospect of an approaching ai winter .
  * Senior . (2020) senior2020improvedSenior, A.W., Evans, R., Jumper, J., Kirkpatrick, J., Sifre, L., Green, T.others  2020\.  Improved protein structure prediction using potentials from deep learning Improved protein structure prediction using potentials from deep learning.  Nature5777792706–710.
  * Tsay Liu (2020) TSAY2020102000Tsay, MY. Liu, ZW.  2020\.  Analysis of the patent cooperation network in global artificial intelligence technologies based on the assignees Analysis of the patent cooperation network in global artificial intelligence technologies based on the assignees.  World Patent Information63102000.  https://www.sciencedirect.com/science/article/pii/S0172219020300922
  * Twiss (1992) twiss1992forecastingTwiss, B.  1992\.  Forecasting for Technologists and Engineers: A Practical Guide for Better Decisions Forecasting for technologists and engineers: A practical guide for better decisions.  P. Peregrinus.  https://books.google.pl/books?id=AINp6drLXdUC
  * van de Schoot . (2021) bayes:naturevan de Schoot, R., Depaoli, S., King, R., Kramer, B., Märtens, K., Tadesse, M.G.Yau, C.  202112\.  Bayesian statistics and modelling Bayesian statistics and modelling.  Nature Reviews Methods Primers11.  http://www.nature.com/articles/s43586-020-00001-2
  * Vinuesa . (2020) vinuesa2020roleVinuesa, R., Azizpour, H., Leite, I., Balaam, M., Dignum, V., Domisch, S.Nerini, F.F.  2020\.  The role of artificial intelligence in achieving the Sustainable Development Goals The role of artificial intelligence in achieving the sustainable development goals.  Nature communications1111–10.
  * WIPO (2019) WorldIntellectualPropertyOrganization2019WIPO  2019\.  WIPO Technology Trends 2019. Artificial Intelligence WIPO Technology Trends 2019. Artificial Intelligence .  Geneva.  World Intellectual Property Organization.
  * B. Wu (2020) 2020Natur.588S.112.Wu, B.  202001\.  Network effect: visualizing AI connections in the natural sciences Network effect: visualizing ai connections in the natural sciences.  Nature5887837S112-S113.
  * H. Wu . (2021) wu2021integratedWu, H., Yi, H. Li, C.  2021\.  An integrated approach for detecting and quantifying the topic evolutions of patent technology: a case study on graphene field An integrated approach for detecting and quantifying the topic evolutions of patent technology: a case study on graphene field.  Scientometrics1–21.
  * You . (2017) Sciento2017developmentYou, H., Li, M., Hipel, K.W., Jiang, J., Ge, B. Duan, H.  2017\.  Development trend forecasting for coherent light generator technology based on patent citation network analysis Development trend forecasting for coherent light generator technology based on patent citation network analysis.  Scientometrics1111297–315.

##  Appendix A Query

⬇

APPC = AT OR BE OR BG OR HR OR CY OR CZ OR DK OR EE OR FI OR FR OR DE OR GR OR
HU OR IE OR IT OR LV OR LT OR LU OR MT OR NL OR PL OR PT OR RO OR SK OR SI OR
ES OR SE OR GB AND CPC=(A61B5/7267 OR G01N33/0034 OR G06F19/24 OR G10H2250/151
OR H04L2025/03464 OR B29C66/965 OR G01N2201/1296 OR G06F19/707 OR G10H2250/311
OR H04N21/4662 OR B29C2945/76979 OR G01S7/417 OR G06F2207/4824 OR
G10K2210/3024 OR H04N21/4663 OR B60G2600/1876 OR G05B13/027 OR G06K7/1482 OR
G10K2210/3038 OR H04N21/4665 OR B60G2600/1878 OR G05B13/0275 OR G06N3/004 OR
G10L25/30 OR H04N21/4666 OR B60G2600/1879 OR G05B13/028 OR G06N3/02 OR
G11B20/10518 OR H04Q2213/054 OR E21B2041/0028 OR G05B13/0285 OR G06N3/12 OR
H01J2237/30427 OR H04Q2213/13343 OR F02D41/1405 OR G05B13/029 OR G06N5 OR
H02P21/0014 OR H04Q2213/343 OR F03D7/046 OR G05B13/0295 OR G06N7 OR
H02P23/0018 OR H04R25/507 OR F05B2270/707 OR G05B2219/33002 OR G06N20 OR
H03H2017/0208 OR Y10S128/924 OR F05B2270/709 OR G05D1/0088 OR G06N99/005 OR
H03H2222/04 OR Y10S128/925 OR F05D2270/707 OR G06F11/1476 OR G06T3/4046 OR
H04L25/0254 OR Y10S706 OR F05D2270/709 OR G06F11/2257 OR G06T9/002 OR
H04L25/03165 OR F16H2061/0081 OR G06F11/2263 OR G06T2207/20081 OR H04L41/16 OR
F16H2061/0084 OR G06F15/18 OR G06T2207/20084 OR H04L45/08 OR G01N29/4481 OR
G06F17/16 OR G08B29/186 OR H04L2012/5686) AND PUD >=19900101


# Constraining light dark matter and mediator with $B^{+}\rightarrow
K^{+}\nu\bar{\nu}$ data

Murat Abdughani<EMAIL_ADDRESS>Yakefu Reyimuaji<EMAIL_ADDRESS>School
of Physical Science and Technology, Xinjiang University,Urumqi 830017, China

###### Abstract

We study the decay of $B^{+}$ meson into $K^{+}$ plus a light mediator $\phi$
followed by production of dark matter (DM) pair $\bar{\chi}\chi$, which can
mimic the flavor changing neural current (FCNC) process of $B^{+}\rightarrow
K^{+}\nu\bar{\nu}$ in the standard model (SM). Constraints from DM relic
density, direct detection, collider and $B-$physics data combined with the
recently reported $B^{+}\rightarrow K^{+}\nu\bar{\nu}$ branching ratio
performed at the Belle II experiment, we analyse couplings between the
mediator and the SM fermions as well as the dark matter particle. We obtain
that: if the process $\phi\rightarrow\bar{\chi}\chi$ is kinematically allowed,
i.e. $m_{\phi}>2m_{\chi}$, then $m_{\phi}$ should be larger than mass of the
muon pair; Otherwise, $\phi$ should play the role of missing energy, and only
allowed decaying channel is $\phi\rightarrow e^{+}e^{-}$.

##  I Introduction

Rare $B$ meson decay with a kaon and neutrino pair in the final state, i.e.
$B^{+}\rightarrow K^{+}\nu\bar{\nu}$, has a clear theoretical prediction among
the flavour changing neutral current (FCNC) processes, and it is severely
suppressed due to the loop effect and the Glashow-Iliopoulos-Maiani mechanism
Glashow  _et al._ (1970). Thus it an plays important role in searching for new
physics and testing the Standard Model (SM). While its SM prediction is
$(5.06\pm 0.14\pm 0.28)\times 10^{-6}$ Bečirević  _et al._ (2023), recent
experiment at the SuperKEKB asymmetric energy electron-positron collider
reported the value of $(2.4\pm 0.7)\times 10^{-5}$ Glazov  with the
significance of about $3\sigma$ above the SM. This indeed calls for the new
physics beyond the SM (BSM).

The weak interaction of neutrinos makes it invisible at the detector, and any
process with missing energy $\cancel{E}$ at the final state, i.e.
$B^{+}\rightarrow K^{+}+\cancel{E}$, contributes to measured branching ratio
of the decay $B^{+}\rightarrow K^{+}\nu\bar{\nu}$ . At the quark level it is
just the $b\rightarrow s+\cancel{E}$ transition. Dark matter (DM) or any other
light particle with sufficiently weak interaction strength involved in this
transition can be considered as the missing energy. Therefore, precision
experiments of $B$ meson decay are crucial in the searching or constraining
the DM or other light BSM particles.

DM on the other side, is an indispensable ingredient of the evolution of the
universe. The thermal DM known to be in equilibrium with SM particles in the
early times and freeze-out later when its annihilation rate could not catch up
with expansion rate of the universe Srednicki  _et al._ (1988). After the
thermal freeze-out, the DM comoving number density remains constant and leads
to the observed DM relic abundance. This, on the other hand, determines the DM
annihilation rate at the freeze-out. Interestingly, annihilation cross-section
of this process is of the same order of magnitude as the weak interaction.
Thus, this weakly-interacting massive particle (WIMP) Ellis  _et al._ (1984)
becomes the mostly well studied and promising DM candidate. The lower limit
for mass of the WIMP is few GeV Lee and Weinberg (1977), which is model
dependent. However, the null result from various DM search experiments have
pushed the lower limit to MeV scale with proper parameter or model selection
Pospelov  _et al._ (2008).

Combined study of DM and new physics hint from B meson decay is an intriguing
window towards the construction of BSM physics. In this work, an extension of
the SM with a light scalar mediator and a Majorana fermion is studies. Relying
the constraints from various experiments, we find optimal ranges of the model
parameters. The paper is structured as follows. In Sec. II a new physics model
is presented. Sec. III gives discussions about the B meson decay into kaon and
SM fermions as well as extra fermionic dark matter pairs. Experimental
constraints and parameter spaces are studied in Sec. IV, followed by
conclusion of the paper drawn in Sec. V.

##  II A simple DM model

In this paper, we study a generic model with light Majorana DM $\chi$ and a
light scalar $\phi$ couples to the SM fermions, the dimension-4 interaction
Lagrangian can be simply written as 111Pseudo-scalar current
$\bar{f}\gamma^{5}f$ vanish at the loop order for $0^{-}\rightarrow 0^{-}$
meson decay processes He  _et al._ (2023); Li  _et al._ (2021).

$\mathcal{L_{\rm
int}}=-\frac{ym_{f}}{v}\phi\bar{f}f-\frac{1}{2}\kappa\phi\bar{\chi}\chi\ ,$
(1)

where $m_{f}$ represent masses of the SM fermions, $y$ is weight of the Yukawa
coupling, $v\simeq 246$ GeV is the electroweak vacuum expectation value (vev).
The Feynman diagrams at the leading-order for DM and fermion pair productions
in the final state are shown in FIG. 1.

We have to note that, simple Lagrangian Eq.(1) is not started from gauge
symmetry, but phenomenologically viable after the electroweak symmetry
breaking of higher-dimensional effective operators Kamenik and Smith (2012).
The constant $\frac{m_{f}}{v}$ in the coupling indicate that the this operator
is induced from integrating out the Higgs portal Patt and Wilczek (2006).
Effective filed theory approach to the decay processes of $B$ with a single
invisible scalar or fermion pair at the final state in dimension-5 or
dimension-6 operators were studied at Refs. Li  _et al._ (2021); Winkler
(2019a); Krnjaic (2016); Filimonova  _et al._ (2020); Kachanovich  _et al._
(2020); He  _et al._ (2023)

Figure 1: The leading-order Feynman diagrams for $b\to s\chi\bar{\chi}$ and
$b\to sf\bar{f}$ processes. The $q$ stands for $u,c,t$ quarks, $f$ include all
quarks and massive leptons.

##  III $B$ meson decay

When the scalar mediator mass $m_{\phi}$ is smaller than $B$ and $K$ meson
mass difference, $m_{\phi}<m_{B}-m_{K}$, the effective Lagrangian for
$b\rightarrow s\phi$ process can be obtained by integrating out heavy
particles running in the loop in the FIG. 1, which is Schmidt-Hoberg  _et al._
(2013a); Batell  _et al._ (2011)

$\mathcal{L}_{b\rightarrow
s\phi}=\frac{ym_{b}}{v}\frac{3\sqrt{2}G_{F}m_{q}^{2}V^{*}_{qs}V_{qb}}{16\pi^{2}}\times\phi\bar{s}_{L}b_{R}+{\rm
h.c.}\ ,$ (2)

where $q$ are the $u,c,t$ quarks in the loop, $G_{F}$ is the Fermi constant
and $V_{qs,qb}$ are CKM matrix elements. It is obvious that contribution from
the quarks other than top are negligible.

The decay width of process $B\rightarrow K+\phi$ is in the form of

$\displaystyle\Gamma_{B\rightarrow K\phi}=$
$\displaystyle\left(\frac{ym_{b}}{v}\frac{3\sqrt{2}G_{F}m_{t}^{2}V^{*}_{ts}V_{tb}}{16\pi^{2}}\right)^{2}\frac{\sqrt{(m_{B}^{2}-(m_{K}+m_{\phi})^{2})(m_{B}^{2}-(m_{K}-m_{\phi})^{2})}}{16\pi
m_{B}^{3}}$ (3) $\displaystyle\times$ $\displaystyle|\langle
K|\bar{s}_{L}b_{R}|B\rangle|^{2}\ ,$

where hadronic transition matrix element is

$\langle
K|\bar{s}_{L}b_{R}|B\rangle=\frac{m_{B}^{2}-m_{K}^{2}}{m_{b}-m_{s}}f_{0}(q^{2})\
,$ (4)

$m_{b,s}$ are bottom and strange quark masses. The form factor is Ball and
Zwicky (2005); Gubernari  _et al._ (2019)

$f_{0}(q^{2})=\frac{0.33}{1-q^{2}/37.46}\ ,$ (5)

with $q^{2}=m_{\phi}^{2}$ for the final state $\phi$ particle.

Subsequently, if kinematically allowable, scalar $\phi$ decay to DM or SM
fermion pairs. The decay width of $\phi$ to fermion pair (includes DM) is

$\Gamma_{\phi\rightarrow\bar{f}f}=\frac{C_{\phi\bar{f}f}^{2}}{8\pi}m_{\phi}\left(1-\frac{4m_{f}^{2}}{m_{\phi}^{2}}\right)^{3/2}\Theta(m_{\phi}-2m_{f})\
,$ (6)

where $C_{\phi\bar{f}f}^{2}$ is scalar-fermion-fermion interaction term
coefficient, $\Theta$ is the Heaviside step function. For the decay width to
the hadronic final states, we consulted to the Refs. Schmidt-Hoberg  _et al._
(2013a); McKeen (2009); Winkler (2019a).

Figure 2: Left: The branching ratios of $B\rightarrow K+\phi$ process as a
function of $m_{\phi}$ and $y$. The red solid line is correspond to central
value from Glazov , blue dashed and green dotted lines are 1$\sigma$ lower and
upper bounds respectively. Right: The branching ratios of process
$B\rightarrow K+\phi$ followed by $\phi\rightarrow\bar{\chi}\chi$ meet with
central value from Glazov  for different $\kappa$ values, where DM masses are
fixed to 0.1 GeV.

In the left panel of FIG. 2, we show the branching ratio of $B\rightarrow
K+\phi$ process satisfying the result of Belle experiment Glazov  in
($m_{\phi},y$) space (assuming that branching ratio of $\phi$ decays into
~~E~~ is 100%). For $m_{\phi}\lesssim 4.5$ GeV, $y$ is around 0.003 to meet
with observed $B\rightarrow K+\nu\bar{\nu}$ value. When $m_{\phi}$ approaches
$m_{B}-m_{K}$, we need larger coupling $y$ due to decreasing phase space. On
the right panel, fixed $m_{\chi}$ to 0.1 GeV, for different $\kappa$ values,
we draw the curves which satisfy same condition as the red line in the left
panel. For the fixed $m_{\chi}$ and $m_{\phi}$, as the value of $\kappa$
decrease, the branching fraction $Br(\phi\rightarrow\bar{\chi}\chi)$ also
decrease. To obtain the value $(2.4\pm 0.7)\times 10^{-5}$, one needs larger
$y$ to increase the $Br(B\rightarrow K+\phi)$, and it can be seen from the
right panel of FIG. 2. Three bumps near $m_{\phi}\sim 1,3.4,3.9$ GeV are due
to the peaks in the decay width to SM particles shown as in FIG. 4 of Ref.
Winkler (2019a), which are caused by Pion pair and Charmonium resonance final
decay products.

##  IV DM relic density and constraints

When the decay process $\phi\rightarrow\bar{\chi}\chi$ is kinematically
allowed, i.e. $m_{\chi}<m_{\phi}/2$, in the vicinity of freeze-out, DM only
can annihilate to the lighter SM fermion pairs final states through
$s$-channel heavy mediator $\phi$. The cross section which determines the DM
relic density for the process $\bar{\chi}\chi\rightarrow\bar{f}f$ is given as

$\sigma_{\bar{\chi}\chi\rightarrow\bar{f}f}=\frac{1}{512\pi}\left(\frac{\kappa
ym_{f}}{v}\right)^{2}\sqrt{1-\frac{4m_{f}^{2}}{s}}\left(1-\frac{2m_{f}^{2}}{s}\right)\left(1-\frac{2m_{\chi}^{2}}{s}\right)\frac{s}{(s-m_{\phi}^{2})^{2}}\,$
(7)

where $s\equiv E^{2}_{\rm cm}$ is the center-of-mass energy squared. Note that
this equation holds for not resonant annihilation. The thermally averaged
cross section is defined as

$\langle\sigma
v\rangle_{\bar{\chi}\chi\rightarrow\bar{f}f}=\int_{4m^{2}_{\chi}}^{\infty}ds\
\sigma_{\bar{\chi}\chi\rightarrow\bar{f}f}v_{\rm
lab}\frac{\sqrt{s-4m_{\chi}^{2}}(s-2m_{\chi}^{2})K_{1}(\sqrt{s}/T)}{8m_{\chi}^{4}TK_{2}^{2}(m_{\chi}/T)}\
,$ (8)

where $K_{i}$ are the modified Bessel functions of order $i$. If DM does not
resonantly annihilate, even for the light fermions $m_{f}\ll m_{\chi}$,
$\langle\sigma v\rangle$ is at least three order less than standard thermal
cross section of $10^{-26}$ cm${}^{3}/s$. Therefore, we have to resort to
resonantly enhanced cross section of DM annihilation. We use Feynrules Alloul
_et al._ (2014) to implement the models, and then import them to
MicroMEGAS-5.3.41 Belanger  _et al._ (2011) for DM relic density and DM-
nucleus cross section calculations.

Parameters | Minimum | Maximum | Prior
---|---|---|---
$m_{\phi}$ | 0.01 GeV | $m_{B}-m_{K}$ | flat/log
$m_{\chi}$ | 0.01 GeV | $m_{\phi}/2$ | flat/log
$\kappa$ | $10^{-6}$ | 4 | log
y | $10^{-6}$ | 0.1 | log
Table 1: Ranges and priors for input parameters adopted in the scans.

To explore full parametric phase space, we scan over the range in the TABLE 1,
where upper limit of 0.1 for $y$ is taken from LEP Acciarri  _et al._ (1996).
We consider likelihoods of DM relic density from PLANCK and direct detection
for various experiments. The total $\chi^{2}_{\rm tot}$ defined as the sum of
individual $\chi^{2}$ values of DM relic density and DD cross section as

$\chi^{2}_{\rm tot}=\chi^{2}_{\Omega h^{2}}+\chi^{2}_{\rm DD}\ .$ (9)

We hire emcee Foreman-Mackey  _et al._ (2013) based on Markov Chain Monte
Carlo (MCMC) method to undertake the task of sampling the parameter space with
the likelihood $\propto{\rm exp}(-\chi^{2}_{\rm tot}/2)$.

The $\chi^{2}$ of DM relic density is described as a Gaussian distribution

$\chi^{2}_{\Omega h^{2}}=\left(\frac{\mu_{t}-\mu_{0}}{\sqrt{\sigma^{2}_{\rm
theo}+\sigma^{2}_{\rm exp}}}\right)^{2}\,$ (10)

where $\mu_{t}$ is the predicted theoretical value we obtained, $\mu_{0}$ is
the experimental central value, and introduce 10% theoretical uncertainty
$\sigma^{2}_{\rm theo}=0.1\mu_{t}$ due to uncertainties from the Boltzmann
equation solver and the entropy table in the early universe. We use PLANCK
2018 data Aghanim  _et al._ (2020) to constrain our predicted relic density
$\Omega h^{2}$. Their reported central value with statistical error is $\Omega
h^{2}=0.118\pm 0.002$.

The estimation of $\chi^{2}$ for the DM-nucleus spin-independent and spin-
dependent cross section is

$\chi^{2}_{\rm DD}=\left(\frac{\sigma_{\rm DD}}{\sigma_{\rm
DD}^{0}/1.64}\right)^{2}\,$ (11)

where $\sigma_{\rm DD}$ and $\sigma_{\rm DD}^{0}$ are predicted theoretical
value obtained from MicroMEGAS-5.3.41 and upper limits of the cross sections
for a given DM mass at 90% confidence level from the experiments DarkSide-50
Agnes  _et al._ (2018), CRESST-III Abdelhameed  _et al._ (2019), PICO-60 Amole
_et al._ (2019), PandaX-4T Meng  _et al._ (2021), and Xenon1T Aprile  _et al._
(2018).

The salient point of this paper over the other previous works is that
branching fraction of $\phi$ decay to all possible SM particles
$Br(\phi\rightarrow{\rm SM\ SM})$ is not unity. Thus, we have to factorize the
limits from the experiments LEP Acciarri  _et al._ (1996), Belle Jia  _et al._
(2022), BESIII Ablikim  _et al._ (2022), LHCb Aaij  _et al._ (2015, 2017),
NA62 Cortina Gil  _et al._ (2021a, b, c), KTeV Winkler (2019b); Alavi-Harati
_et al._ (2004); Abouzaid  _et al._ (2008), CHARM Winkler (2019b); Clarke  _et
al._ (2014); Bezrukov and Gorbunov (2010); Schmidt-Hoberg  _et al._ (2013b);
Alekhin  _et al._ (2016) and PS191 Gorbunov  _et al._ (2021). If we assume
$Br(\phi\rightarrow{\rm SM\ SM})=f_{\rm SM}$ for a sample point $y=y_{0}$ and
other parameters are fixed, $y$ should equivalent to

$y^{\prime}=(y_{0}^{4}/f_{\rm SM})^{1/4}$ (12)

in the case of SM only decay products. Therefore, the exclusion limit become
more robust if branching ratio to the SM final states are less than one. Note
that pure SM final state cannot explain $B$ meson decay anomaly unless the
mediator is long lived compared to the detector size.

Figure 3: The sample points with $\chi^{2}_{\rm tot}<6$ and $Br(B\rightarrow
K+\bar{\chi}\chi)<2\sigma$. The colorbar represents the branching fraction of
$\phi\rightarrow\bar{\chi}\chi$ process. The red and orange solid line are the
upper limits of experimental constraints from LEP Acciarri  _et al._ (1996)
and LHCb Aaij  _et al._ (2015, 2017), respectively.

In the FIG. 3, we show the sample plots with $\chi^{2}_{\rm tot}<6$ (minimum
value of $\chi^{2}_{\rm tot}\simeq 0$) and require that the $Br(B\rightarrow
K+\cancel{E})$ are in the $2\sigma$ range. The relevant constraints from LEP
Acciarri  _et al._ (1996) and LHCb Aaij  _et al._ (2015, 2017) are also shown.
From the figure we know that, in the area above experimental upper limits, the
branching ratio to the DM pair should be equal to survive. The SM pair final
states can be part of $\phi$ decay channels below the experimental exclusion
lines, the closer to the upper limit, the bigger $Br(\phi\rightarrow SM\ SM)$
gets, and it is known from Eq. (12). Generally, in order to obtain definite
value of $Br(B\rightarrow K+\cancel{E})$, $Br(B\rightarrow K+\phi)$ should be
increases as the $\phi$ decay ratio to the SM particles become larger,
consequently we need larger $y$ as shown in the FIG.3. The survived samples
are located at Charmonium resonant area or large $m_{\phi}$ near the
$m_{B}-m_{K}$. There in not sample points survived when $m_{\phi}\lesssim
2m_{\mu}$, where only $\bar{\chi}\chi\rightarrow e^{+}e^{-}$ channel allowed
and corresponding annihilation cross section is too small even if one
introduces the resonance enhancement. When annihilation process
$\bar{\chi}\chi\rightarrow\bar{c}c$ is closed, DM relic density is on the high
side, and that is the reason that there is a gap near $m_{\phi}\simeq 2.5$
GeV.

Note that if the lifetime of $\phi$ particle is long enough to escape the
detector before its decay, $\phi$ is the missing energy regardless of the
decay products. However, to match with DM relic density and $Br(B\rightarrow
K+\cancel{E})$, both $\kappa$ and $y$ can not be too small. In our samples in
FIG. 3, maximum lifetime of the $\phi$ particle is at the order of $\sim
10^{-14}$s, thus it decays before the escaping from the detector.

On the other case, if $m_{\chi}>m_{\phi}/2$, although the decay products of
$\phi$ are only SM particles, long lived $\phi$ relative to the detector size
behaves like missing energy. When $y\simeq 3\times 10^{-3}$ and
$m_{\phi}<2m_{\mu}$, $\Gamma(\phi\rightarrow e^{+}e^{-})$ can be less than
$10^{-15}$ GeV to be invisible, and $m_{\phi}$ should be heavier than electron
pair mass to decay before BBN. Heavy DM in this case easily satisfy relic
density Schmidt-Hoberg  _et al._ (2013b); Winkler (2019a).

##  V conclusion

In this paper, we study a light mediator $\phi$ and a light DM $\chi$ in a
simple model which lead to a signal in the recent Belle $B\rightarrow
K+\bar{\nu}\nu$ measurement. Combined with DM relic density and direct
detection constraints, when $2m_{\chi}<m_{\phi}$, i.e.
$\phi\rightarrow\bar{\chi}\chi$ is available, we obtain that $m_{\phi}$ should
be larger than muon pair mass. In addition, branching fraction of $\phi$ to DM
pair equal to one, except for three narrow areas where $\phi$ can decay to SM
particles partially. When $m_{\chi}>m_{\phi}/2$, DM easily have correct relic
density, $\phi$ should play the role of missing energy in Belle detector, and
small width of $\phi$ require $m_{\phi}<2m_{\mu}$, in other word, $\phi$ only
can decay to electron pair.

## acknowledgment

MA is supported by Tianchi Yingcai Young Doctoral Project. YR is supported by
the Natural Science Foundation of Xinjiang Uygur Autonomous Region of China
under grant No. 2022D01C52.

## References

  * Glashow  _et al._ (1970) S. L. Glashow, J. Iliopoulos, and L. Maiani, Phys. Rev. D 2, 1285 (1970).
  * Bečirević  _et al._ (2023) D. Bečirević, G. Piazza, and O. Sumensari, Eur. Phys. J. C 83, 252 (2023), arXiv:2301.06990 [hep-ph] .
  * (3) A. Glazov (Belle II), “News from Belle II,” EPS conference 2023\.
  * Srednicki  _et al._ (1988) M. Srednicki, R. Watkins, and K. A. Olive, Nucl. Phys. B 310, 693 (1988).
  * Ellis  _et al._ (1984) J. R. Ellis, J. S. Hagelin, D. V. Nanopoulos, K. A. Olive, and M. Srednicki, Nucl. Phys. B 238, 453 (1984).
  * Lee and Weinberg (1977) B. W. Lee and S. Weinberg, Phys. Rev. Lett. 39, 165 (1977).
  * Pospelov  _et al._ (2008) M. Pospelov, A. Ritz, and M. B. Voloshin, Phys. Lett. B 662, 53 (2008), arXiv:0711.4866 [hep-ph] .
  * He  _et al._ (2023) X.-G. He, X.-D. Ma, and G. Valencia, JHEP 03, 037 (2023), arXiv:2209.05223 [hep-ph] .
  * Li  _et al._ (2021) G. Li, T. Wang, J.-B. Zhang, and G.-L. Wang, Eur. Phys. J. C 81, 564 (2021), arXiv:2103.12921 [hep-ph] .
  * Kamenik and Smith (2012) J. F. Kamenik and C. Smith, JHEP 03, 090 (2012), arXiv:1111.6402 [hep-ph] .
  * Patt and Wilczek (2006) B. Patt and F. Wilczek, (2006), arXiv:hep-ph/0605188 .
  * Winkler (2019a) M. W. Winkler, Phys. Rev. D 99, 015018 (2019a), arXiv:1809.01876 [hep-ph] .
  * Krnjaic (2016) G. Krnjaic, Phys. Rev. D 94, 073009 (2016), arXiv:1512.04119 [hep-ph] .
  * Filimonova  _et al._ (2020) A. Filimonova, R. Schäfer, and S. Westhoff, Phys. Rev. D 101, 095006 (2020), arXiv:1911.03490 [hep-ph] .
  * Kachanovich  _et al._ (2020) A. Kachanovich, U. Nierste, and I. Nišandžić, Eur. Phys. J. C 80, 669 (2020), arXiv:2003.01788 [hep-ph] .
  * Schmidt-Hoberg  _et al._ (2013a) K. Schmidt-Hoberg, F. Staub, and M. W. Winkler, Phys. Lett. B 727, 506 (2013a), arXiv:1310.6752 [hep-ph] .
  * Batell  _et al._ (2011) B. Batell, M. Pospelov, and A. Ritz, Phys. Rev. D 83, 054005 (2011), arXiv:0911.4938 [hep-ph] .
  * Ball and Zwicky (2005) P. Ball and R. Zwicky, Phys. Rev. D 71, 014015 (2005), arXiv:hep-ph/0406232 .
  * Gubernari  _et al._ (2019) N. Gubernari, A. Kokulu, and D. van Dyk, JHEP 01, 150 (2019), arXiv:1811.00983 [hep-ph] .
  * McKeen (2009) D. McKeen, Phys. Rev. D 79, 015007 (2009), arXiv:0809.4787 [hep-ph] .
  * Alloul  _et al._ (2014) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, and B. Fuks, Comput. Phys. Commun. 185, 2250 (2014), arXiv:1310.1921 [hep-ph] .
  * Belanger  _et al._ (2011) G. Belanger, F. Boudjema, P. Brun, A. Pukhov, S. Rosier-Lees, P. Salati, and A. Semenov, Comput. Phys. Commun. 182, 842 (2011), arXiv:1004.1092 [hep-ph] .
  * Acciarri  _et al._ (1996) M. Acciarri _et al._ (L3), Phys. Lett. B 385, 454 (1996).
  * Foreman-Mackey  _et al._ (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Publ. Astron. Soc. Pac. 125, 306 (2013), arXiv:1202.3665 [astro-ph.IM] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Agnes  _et al._ (2018) P. Agnes _et al._ (DarkSide), Phys. Rev. Lett. 121, 081307 (2018), arXiv:1802.06994 [astro-ph.HE] .
  * Abdelhameed  _et al._ (2019) A. H. Abdelhameed _et al._ (CRESST), Phys. Rev. D 100, 102002 (2019), arXiv:1904.00498 [astro-ph.CO] .
  * Amole  _et al._ (2019) C. Amole _et al._ (PICO), Phys. Rev. D 100, 022001 (2019), arXiv:1902.04031 [astro-ph.CO] .
  * Meng  _et al._ (2021) Y. Meng _et al._ (PandaX-4T), Phys. Rev. Lett. 127, 261802 (2021), arXiv:2107.13438 [hep-ex] .
  * Aprile  _et al._ (2018) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 121, 111302 (2018), arXiv:1805.12562 [astro-ph.CO] .
  * Jia  _et al._ (2022) S. Jia _et al._ (Belle), Phys. Rev. Lett. 128, 081804 (2022), arXiv:2112.11852 [hep-ex] .
  * Ablikim  _et al._ (2022) M. Ablikim _et al._ (BESIII), Phys. Rev. D 105, 012008 (2022), arXiv:2109.12625 [hep-ex] .
  * Aaij  _et al._ (2015) R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 115, 161802 (2015), arXiv:1508.04094 [hep-ex] .
  * Aaij  _et al._ (2017) R. Aaij _et al._ (LHCb), Phys. Rev. D 95, 071101 (2017), arXiv:1612.07818 [hep-ex] .
  * Cortina Gil  _et al._ (2021a) E. Cortina Gil _et al._ (NA62), JHEP 03, 058 (2021a), arXiv:2011.11329 [hep-ex] .
  * Cortina Gil  _et al._ (2021b) E. Cortina Gil _et al._ (NA62), JHEP 06, 093 (2021b), arXiv:2103.15389 [hep-ex] .
  * Cortina Gil  _et al._ (2021c) E. Cortina Gil _et al._ (NA62), JHEP 02, 201 (2021c), arXiv:2010.07644 [hep-ex] .
  * Winkler (2019b) M. W. Winkler, Phys. Rev. D 99, 015018 (2019b), arXiv:1809.01876 [hep-ph] .
  * Alavi-Harati  _et al._ (2004) A. Alavi-Harati _et al._ (KTeV), Phys. Rev. Lett. 93, 021805 (2004), arXiv:hep-ex/0309072 .
  * Abouzaid  _et al._ (2008) E. Abouzaid _et al._ (KTeV), Phys. Rev. D 77, 112004 (2008), arXiv:0805.0031 [hep-ex] .
  * Clarke  _et al._ (2014) J. D. Clarke, R. Foot, and R. R. Volkas, JHEP 02, 123 (2014), arXiv:1310.8042 [hep-ph] .
  * Bezrukov and Gorbunov (2010) F. Bezrukov and D. Gorbunov, JHEP 05, 010 (2010), arXiv:0912.0390 [hep-ph] .
  * Schmidt-Hoberg  _et al._ (2013b) K. Schmidt-Hoberg, F. Staub, and M. W. Winkler, Phys. Lett. B 727, 506 (2013b), arXiv:1310.6752 [hep-ph] .
  * Alekhin  _et al._ (2016) S. Alekhin _et al._ , Rept. Prog. Phys. 79, 124201 (2016), arXiv:1504.04855 [hep-ph] .
  * Gorbunov  _et al._ (2021) D. Gorbunov, I. Krasnov, and S. Suvorov, Phys. Lett. B 820, 136524 (2021), arXiv:2105.11102 [hep-ph] .