# A Duty-Cycle-Efficient Synchronization
Protocol for Slotted-Aloha in LoRaWAN

Amavi Dossa  El Mehdi Amhoud

###### Abstract

In the current context of massive IoT, the Pure-Aloha scheme used in LoRaWAN
is reaching its limit, and Slotted-Aloha is being considered as an
alternative, as it offers twice Pure-Aloha’s packet success rate. It however
requires synchronization accross the nodes. In this paper, we propose a new
slot structure adapted to devices with low quality clock, and a duty-cycle
efficient synchronization protocol for LoRaWAN class A devices with the lowest
overhead to date. We discuss the conditions of its integration into LoRaWAN.
The experimental results confirm that it succeeds in tracking each device’s
synchronization state, identifying the exact moment they desynchronize and
resynchronizing them. The proposed protocol is also proven to be more duty-
cycle efficient than existing fixed-rate synchronization solutions.

###### Index Terms:

LoRaWAN, scalability, slotted-Aloha, synchronization, massive IoT, duty-cycle
efficiency

##  I Introduction

Since its introduction almost a decade ago, the Long Range Wide Area Network
(LoRaWAN) technology has met an unprecedented success among existing Low Power
Wide Area Networks (LPWAN) technologies. Its low cost, low power and ability
to operate over noisy channels made it reach several markets and application
sectors like smart cities, smart grids and metering, agriculture monitoring,
wireless sensor networks and localization [1, 2, 3, 4, 5]. With the increasing
number of nodes, LoRaWAN enters the so-called Massive IoT era. Therefore,
scalability is currently one of the biggest research challenges [6].

To cope with the previous challenges, propositions include the use of slotted-
Aloha as an alternative to pure-Aloha which LoRaWAN currently uses. It
requires all the nodes to be synchronized. LoRaWAN class B devices regularly
receive beacons from the gateway that can be used to synchronize. Similarly,
class C devices are listening to the channel almost all the time, hence can be
synchronized using broadcast signals. The challenge resides with class A
devices that do not share any synchronization signals. Several works have
targeted synchronization mechanisms in wireless sensor networks over the past
years [7][8], but very few were dedicated to LoRaWAN, [9, 10, 11].
Nonetheless, they didn’t take into account the duty cycle limitation on
LoRaWAN gateways, whilst this is a key factor that requires optimization for a
wide slotted-Aloha adoption in LoRaWAN.

With that in mind, we propose in this paper a protocol that tracks the
synchronization state of end-devices based on their uplink packets arrival
time. This allows the network server to detect the specific moment a device
desynchronizes and transmits, through the acknowledgement packet, the
necessary information for its resynchronization. This minimizes the downlink
overhead and results in a duty cycle gain for the gateway, without requiring
any additional hardware. This paper focuses on the synchronization aspect
rather than a slotted-Aloha implementation. Our contribution is twofold:

  * •

We introduce a two-guard-intervals variant of the existing slot structure.
This design allows us to define the acceptable bounds on end-devices drift,
for a given length of the guard intervals.

  * •

We propose a monitoring and synchronization protocol on top of LoRaWAN that is
able to identify any desynchronized devices and resynchronize them. We show
that it requires only two additional bytes to encode the synchronization-
related information. To the best of our knowledge, this is the LoRaWAN
synchronization protocol with the lowest overhead.

The rest of this paper is organized as follows: Section II states the problem
after a short LoRaWAN overview; Section III introduces the proposed
synchronization protocol and Section IV describes how it integrates into the
existing LoRaWAN stack; Section V presents and discusses the experimental
results of our synchronization solution; Section VI concludes the paper and
set forth our perspectives.

##  II  LoRaWAN Overview and Problem Statement

LoRa (Long Range) is a proprietary modulation technique developed by Semtech
Corporation. It is based on the Chirp Spread Spectrum (CSS) technique, and
operates in the sub-Gigahertz frequency bands. LoRaWAN, on the other hand, is
a communication protocol built on top of LoRa, and designed specifically for
wide-area networks with low-power consumption, low datarate IoT devices. It
was developed by the LoRa Alliance [12]. It employs a star of stars
architecture, as depicted in Fig. 1: End-devices are served by one or more
gateways which, in turn, communicates with a network server. The end-devices
transmit their frames to the gateways using the radio interface; the gateways
forward those frames to the network server through a back-haul link; the
network server is in charge of sending the payload to the intended application
server. Similarly, in case of a downlink, the network server sends the frame
to one gateway that forwards it to the intended end-device through the radio
interface. Additionally, three classes of end-devices are supported,
corresponding to different power saving policies: class B and C are
respectively for deterministic and lowest downlink latency, but consume more
power; class A, on the other hand, has the lowest power consumption profile,
but does has a non-deterministic downlink latency since the only available
downlink windows are after an uplink transmission initiated by an end-device.

Figure 1: LoRaWAN Star of Stars topology

Pure-Aloha is an access scheme that allows devices to transmit anytime without
any form of coordination. Because of its simplicity, it is adapted to
resource-constrained devices and is often adopted by low-power and low-
complexity networks. It is however limited to only 18% packet success rate
when the number of nodes increases, due to collisions [13]. On the other hand,
slotted-Aloha doubles the packet success rate, but requires the nodes to be
synchronized. LoRaWAN was initially based on pure-Aloha but, with the
increasing number of end-devices, slotted-Aloha is being seen as a good
alternative.

LoRaWAN targets low-cost end-devices, that can have poor clock quality,
leading to significant time drift. Hence, devices should be resynchronized on
regular basis to avoid slot violation. Existing synchronization solutions use
the end of uplink transmission as a common reference between an end-device and
the network server [9] [10]: the network server appends the timestamp of the
end of uplink into the downlink sent back to the end-device; the latter will
compare this value to the one it locally stored and use the difference to
correct its clock drift. This however results in a greater downlink air-time,
thus consumes more duty cycle. Consequently, the challenge resides in finding
a synchronization rate that keeps the end-devices drift under a certain
threshold, while minimizing the duty cycle overhead. Assessing each node’s
clock quality to set the synchronization rate accordingly is not practical:
one network could serve thousands of nodes with different clock qualities and
this operation would take a huge amount of time; plus, it will require a non-
negligible software overhead to match each node with its synchronization rate.
On the other hand, choosing a single synchronization rate for all nodes seems
practical, but is not optimal: nodes with better clock quality will be
resynchronized more often that required, adding unnecessary overhead, while
those with worse clock quality will not be resynchronized on time, leading to
significant drift, hence slot violation. Moreover, a node’s drift is also
subject to external conditions such as voltage level and temperature. Hence, a
static synchronization becomes very limited facing the previous challenges.

In what follows, we propose a synchronization protocol that dynamically adapt
to each node by tracking its drift in time.

##  III Dynamic Synchronization Protocol Design

In this section, we describe the layout of the proposed synchronization
protocol, starting by detailing the underlying slot structure.

We consider a LoRaWAN gateway (GW) connected to a network server (NS) and
serving $N$ class A end-devices (ED). We note $t_{ns}(t)$ the local time of NS
at a moment t and $t_{ed}^{k}(t)$ that of the k-th ED, $k=\\{0,1,...,N-1\\}$.
Network servers run on dedicated hardware with much more stable clock compared
to end-devices. So, for the remainder of this paper, we will consider the NS
clock perfect. We also suppose that NS and the EDs share the same time origin.
Thus, at any moment t, the time drift $\Delta t^{(k)}$ of the k-th ED is given
by Eq. (1).

$\Delta t^{(k)}(t)=t_{ns}(t)-t_{ed}^{(k)}(t)$ (1)

$|\Delta t^{(k)}(t)|\leq\Delta_{max}$ (2)

Given a maximum tolerable drift $\Delta_{max}$, an ED is considered in-sync if
the condition in Eq. (2) holds, and out-sync otherwise. If the NS can track a
ED’s drift accross time, it can then detect the exact moment it goes out-sync
and issue a resynchronization.

###  III-A Slot Structure Design

With slotted-Aloha, time is divided into slots of predefined length and nodes
can transmit only at the beginning of each slot. The classic slot consists of
a transmission interval $T_{r}$ and a confidence interval $T_{b}$: the
device’s transmission normally starts at the beginning of the slot and lasts
$T_{r}$; $T_{b}$ ensures that the transmission ends within the slot duration
even in case of clock drift. As depicted by Fig. 2, $T_{b}$ can be seen as
virtually made of two guard intervals: a backward guard interval $T_{b1}$ –
that prevents against drift towards the slot’s start – and a forward guard
interval $T_{b2}$ – that prevents against drift towards the slot’s end.
Moreover, the propagation delay is neglected: uplink air-times range from
hundreds of ms while the propagation delay equals $50\mu s$, assuming a line
of sight and a typical LoRa coverage of 15 km; it can then be neglected even
in case of multipath propagation. The slot is therefore made of the uplink
air-time $T_{tx}$ and the confidence interval $T_{b}$.

Devices – that are microcontroller-based – keep time by counting periodic
pulses generated from the system clock, which frequency $f_{0}$ is known. The
clock however has frequency precision, depending on its quality [14]. The
instantaneous frequency is given by Eq. (3). $\Delta f$ is a real number, so
the frequency can either increase or decrease.

$f(t)=f_{0}+\Delta f(t)$,
---
$\Delta f(t)$ being the instantaneous frequency drift.
(3)

Moreover, relatively to a device with an ideal clock running at $f_{0}$,
another device running at $f_{1}>f_{0}$ will be ahead in time while another
one running at $f_{2}<f_{0}$ will be late, as illustrated by Fig. 3. The
maximum tolerated backward and forward drifts being respectively $T_{b1}$ and
$T_{b2}$, the frame in-sync condition can be derived through Eq. (4).
Moreover, Eq. (5) gives the in-sync condition if the end of uplink
transmission is used as a reference. Therefore, the NS can compute the k-th
ED’s drift upon a frame reception, where $T_{Tx-ns}$ is the end of uplink
transmission expected by the NS and $T_{Tx-ed}$ the moment the uplink
transmission effectively ends. If an out-sync is detected, the device is
resynchronized using the algorithm described in the following subsection.

| $-T_{b2}<\Delta t^{(k)}(t)<T_{b1}$
---|---
$\implies$ | $-T_{b2}<t_{ns}(t)-t_{ed}^{(k)}(t)<T_{b1}$
$\implies$ | $t_{ns}(t)-T_{b1}<t_{ed}^{(k)}(t)<t_{ns}(t)+T_{b2}$
(4)

$T_{Tx-ns}-T_{b1}<T_{Tx-ed}^{(k)}<T_{Tx-ns}+T_{b2}$ (5)

Figure 2: Slotted-Aloha: Virtual double confidence intervals Figure 3:
Apparent drift effect on slot position

The uplink air-time $T_{tx}$ depends on physical layer parameters and can be
calculated from Eq. (6), [15]. Table I describes the symbols involved in the
equations.

$T_{packet}=T_{preamble}+T_{payload}$ (6)

$T_{preamble}=(n_{preamble}+4.25)T_{s}$ (7)

$T_{s}=\frac{2^{SF}}{BW}$ (8)

$T_{payload}=n_{payload}\times T_{s}$ (9)

$n_{bits}$ | $=$ | $ceil\left(\frac{8PL-4SF+28+16CRC-20IH}{4(SF-2DE)}\right)$
---|---|---
$n_{payload}$ | $=$ | $8+max\left(n_{bits}(CR+4),0\right)$
(10)

TABLE I: Air-time equations symbols description Symbols |  Parameters
---|---
Tpacket |  Total packet air-time
Tpreamble |  Preamble duration
Tpayload |  Payload duration
Ts |  Symbol period
npayload |  Number of payload symbols
PL |  Payload size (in bytes)
SF |  Spreading Factor
IH |  Implicit Header flag
DE |  Low data rate optimization flag
CRC |  Cyclic Redundancy Check flag
CR |  Coding Rate

The guard intervals $T_{b1}$ and $T_{b2}$ represent unconsumed time slots, so
they should be calculated as an acceptable ratio of $T_{tx}$, based on the
application.

###  III-B Synchronization and Monitoring Protocol

Algorithm 1 Network Server

1:receive()

2:arrival $\leftarrow$ current_time() - ref

3:arrival $\leftarrow$ arrival mod T_slot

4:if out_sync() then

5: remain $\leftarrow$ T_slot - arrival

6: add_remain_to_ack()

7:wait_for_rx_window()

8:transmit_ack()

Algorithm 2 end-device

1:if !is_first_tx then

2: is_first_tx $\leftarrow$ 1

3: slot_start $\leftarrow$ current_time()

4:else

5: wait_for_next_slot()

6:transmit()

7:beg $\leftarrow$ current_time()

8:wait_for_rx_window()

9:get_ack()

10:end $\leftarrow$ current_time()

11:if remain_in_ack() then

12: remain $\leftarrow$ get_remain_from_ack()

13: elapsed $\leftarrow$ end - beg

14: t $\leftarrow$ remain - elapsed

15: if $t<0$ then

16: t $\leftarrow$ t mod T_slot

17: slot_start $\leftarrow$ current_time() + t

The core idea of our protocol lies in a unique time reference, ref, stored at
the network server and is used to synchronize all the end-devices and later
monitor their synchronization from the incoming packets arrival time. ref is
set only once with the current time when the network server boots up, and it
is considered as the beginning of the very first slot. From there it is easy
to get the future slots start through Eq. (11).

$slot\\_start[n]=ref+n\times T\\_slot$, | $n\in\mathbf{N^{*}}$
---|---
(11)

At any time, Eq. (12) gives the relative position in the ongoing slot – that
is, how much time elapsed since the beginning of this slot. Similarly, Eq.
(13) gives the remaining time before the next slot.

$position=(time-ref)\mbox{ mod }T\\_slot$ (12)

$remaining\\_time=T\\_slot-((time-ref)\mbox{ mod }T\\_slot)$ (13)

The protocol is divided into two algorithms: one for the NS and the other one
for EDs. For an ED, the synchronization consists in simply saving the
beginning of a slot into a variable at the right time, so that it can serve to
identify future slots start using Eq. (11). To keep the algorithms short, we
substituted some portions with self-explanatory function names such as
current_time(), wait_for_next_slot(), or wait(). Those are just steps of the
algorithms that can be simple sets of instructions instead of functions.
Below, we explain the flow of both algorithms.

####  III-B1 Algorithm 1: Network Server

Lines 1-3: The server waits for incoming packets. Once a packet is received,
the server immediately uses Eq. (12) to calculate its position inside the
ongoing slot.
Lines 4-7: The server checks whether the packet is out-sync; in which case, it
gets the time remaining before the next slot from Eq. (13), and piggybacks it
to the ACK packet.
Lines 8-9: Once the ACK packet is ready, the server waits for the next receive
window and sends it.

####  III-B2 Algorithm 2: end-device

Lines 1-6: An ED that transmits for the first time does not have any
synchronization reference on its side. So, it takes the transmission moment as
the reference. If that first frame arrives in-sync, this reference will be
kept, otherwise it will be corrected. The corrected reference will be used to
determine slot start for future transmissions (line 5). The flag is_first_tx,
by default initialized to false, is used to identify the first transmission.
Lines 8-11: The end of the transmission is timestamped into beg. This instant
coincides with the one when the NS takes the packet arrival time. The moment
the ACK packet is received from the NS is also timestamped in end.
Lines 12-20: If the received ACK packet contains remaining-time bytes, it
means the last transmission was received out-sync. The remaining time is
extracted and used to resynchronize the device. The elapsed time between the
moment the network server computed the remaining-time and the moment the
device received the ACK packet – equal to end - beg – is substracted from
remaining-time. Then, The ED resynchronizes by updating its slot start
reference with remaining-time added to its current local time. It can happen
that the time elapsed is greater than the remaining time. This scenario is
handled in lines 16-18 by the mean of a modulo operation. Additionally, the
presence of remaining-time in the ACK packet can be checked by different
means: for instance, the implementation used for the experimental results
section below assumed a downlink with no MAC commands; so, checking the packet
size was enough in that case.

##  IV Protocol Integration into LoRaWAN

With the protocol already established, we now discuss how it integrates into
the LoRaWAN stack. LoRaWAN frame structure includes a Frame Options field for
sharing MAC commands between end-devices and the network server. We use this
field to send the synchronization-related bytes, the remaining-time for
instance. Its value is between 0 – for a frame received exactly at the end of
the slot – and $T_{slot}$ – for a frame received at the exact beginning of a
slot. Table II specifies the radio parameters that gives the longest packet
air-time: LoRa spreading factors ranges from 5 to 12, the greater the
spreading factor the more the signal is spread in time; LoRa defines coding
rates 1, 2, 3 and 4 corresponding resp. to effective rates $\frac{4}{5}$,
$\frac{4}{6}$, $\frac{4}{7}$ and $\frac{4}{8}$; LoRaWAN uses only 125 kHz, 250
kHz and 500 kHz bandwidth, the smallest bandwidth giving the greatest symbol
duration thus the longest air-time; LoRa allows a payload limit of 255 bytes.
The resulting maximum uplink air-time equals 11936 ms, requiring
$\log_{2}(11936)$ = 14 bits for encoding. 2 bytes are therefore enough to
encode the slot length $T\\_slot$ and leave a margin of 53599 ms for
$T_{b1}+T_{b2}$. Those 2 bytes are the overhead induced by our protocol,
against the 8 bytes in [9].

To enable coexistence of devices using this protocol for slotted-Aloha with
the preexisting ones in LoRaWAN, we propose to dedicate special FPort (Frame
Port) values to this protocol in LoRaWAN specifications. This is an efficient
method, as it will allow the NS to distinguish the ED requiring
synchronization from those that don’t, and this without any additional
overhead on the transmitted payload. In the performance evaluation section, we
set the FPort to 198; it is not used for any services in the current
specification, except for proprietary protocols.

TABLE II: Longest air-time LoRa Radio Parameters. Parameter | Value
---|---
Coding Rate | 4
Spreading Factor | 12
Bandwidth (kHz) | 125
Payload | 255
Figure 4: Test bench equipment: A National Instruments USRP-2920 SDR as
Gateway; An Adafruit Feather M0 and a TTGO Esp32 LoRa are used as end-devices;
The network server was running on a 6-cores 12-threads 4.0 GHz laptop.

##  V Experimental Results

Figure 5: Protocol performance: devices synchronization and monitoring

(a)

(b)

Figure 6: Fixed-rate synchronization algorithm performance, 6.5 hours
experiment: (a) 1 hour round duration; (b) 30 mins round duration

In this section, we set up a test bench to evaluate the performances of our
protocol. We set up a custom gateway by providing few additional components to
the Open Source GNU Radio implementation of LoRa physical layer in [16]. We
also developed a LoRaWAN network server with the necessary mechanisms for this
protocol. All the related codes are made available on Github [17]. Fig. 4
shows our experimental test bench: a USRP-2920 serves as a gateway, and an
Adafruit Feather M0 – that has a very unstable clock – and a TTGO Esp32 LoRa –
with a much more stable clock – boards were used as end-devices. The network
server was running locally, on a 6-cores 4.0 GHz laptop.

The performance of our protocol will be compared to that of the algorithm used
in [9], where upchirp-modulated symbols were used for both the uplink and
downlink; as a consequence, the slot was made of the uplink air-time $T_{tx}$,
the downlink air-time $T_{rx}$, the 1s gap, and the guard interval $T_{b1}$
and $T_{b2}$. So, we did the same here for a fair comparison.

We first run an experiment using the proposed synchronization protocol. The
end-devices transmit a packet every 30s. The parameters in Table III together
with Eqs. (6-10) gives the air-time in Table IV. Total slot length refers to
the slot duration and accounts for both the uplink and downlink air-time, the
1000 ms delay, and the slot guards $T_{b1}$ and $T_{b2}$ of 180 ms each. We
run the setup for 6.5 hours, logged the relative position of each frame within
the ongoing slot, and plotted the evolution of the devices synchronization in
Fig. 5. The upper and lower bounds mark the limit of the synchronization
interval, and the ideal arrival line marks the position of frames received
from a perfectly synchronized node. The first frames of both devices were
received out-of-sync; that was expected since they have a chance of only
$\frac{T_{b1}+T_{b2}}{T_{slot}}\approx$ 20% to be received in-sync. They got
synchronized immediately and their next transmission matches the ideal arrival
line. The evolution of the synchronization graphs shows the difference in the
devices’ clocks stability and the success of the protocol in resynchronizing
them: The feather desynchronized three more times – twice during the first two
hours – and get corrected right after, while the TTGO never did because its
drifts are very small compared to the guard intervals value; hence it didn’t
need any resynchronization. In total, both devices were synchronized 5 times
during the 6.5 hours.

We run another two other experiments, this time using the fixed-rate
synchronization algorithm described in [9]: the EDs transmit a packet every
30s; after each round, the NS measures their clock drift, logs it and re-
synchronizes them; each round last 1 hour for the first experiment and 30
minutes for the second one. The resulting drifts are plotted as comparative
bar diagrams between both devices and for both experiments. For the 1 hour
round duration in Fig. 6a, the fixed-rate algorithm encounters two slot
violations by the Feather and none by TTGO, and both were resynchronized in
total 12 times during the whole experiment; this is 2.4 times more overheads
than our adaptive algorithm that does even prevent slot violation . And, for
the 30 mins round duration from Fig. 6b, there has been no slot violation for
any of the devices, but they were however resynchronized 26 times; this is
five times more overheads than our algorithm and for the same result.

In summary, it clearly appears that the fixed-rate synchronization algorithm
is not able to prevent slot violations while keeping a low synchronization
overhead at the same time. Our algorithm, on the other hand is able to
accomplish both by adapting to each device separately, hence succeeds in being
more duty-cyle efficient.

TABLE III: Experiment LoRa radio parameters Parameter | Value
---|---
Frequency (MHz) | 868
Bandwidth (kHz) | 125
Coding Rate | 1
Preamble size | 8
Uplink Payload size | 193
Uplink Spreading Factor | 7
Downlink Payload size | 19
Downlink Spreading Factor | 8
TABLE IV: Experiment Packets air-time Parameter | Value (ms)
---|---
Uplink air-time | 306
Downlink air-time | 91
Total slot length | 1757

##  VI Conclusion

In this paper, we proposed and tested a synchronization protocol for LoRaWAN
class A devices, designed to optimize gateways duty cycle consumption by
reducing the downlink air-time overhead, along with a novel slot structure.
The experimental results confirmed that our protocol was able to track each
device synchronization state separately, and resynchronize it only when it
detects a desynchronization. Furthermore, its low overhead (only two bytes),
duty-cycle efficiency and seamless integration with existing schemes position
it as a powerhouse for slotted-Aloha adoption in LoRaWAN deployments within
the current massive IoT context.

Our solution could be improved by reducing downlink ACK packets, since it
currently requires a downlink acknowledgement for each uplink so that devices
can identify collisions and retransmit if necessary. Furthermore, future works
will investigate the impact of the increasing number of nodes on the
synchronization accuracy.

## Acknowledgment

This work was sponsored by the Junior Faculty Development program under the
UM6P-EPFL Excellence in Africa Initiative.

## References

  * [1] P. Ferrari, E. Sisinni, P. Bellagente, D. F. Carvalho, A. Depari, A. Flammini, M. Pasetti, S. Rinaldi, and I. Silva, “On the use of lorawan and cloud platforms for diversification of mobility-as-a-service infrastructure in smart city scenarios,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 71, pp. 1–9, 2022.
  * [2] M. de Castro Tomé, P. H. Nardelli, and H. Alves, “Long-range low-power wireless networks and sampling strategies in electricity metering,” _IEEE Transactions on Industrial Electronics_ , vol. 66, no. 2, pp. 1629–1637, 2018.
  * [3] H. M. Jawad, R. Nordin, S. K. Gharghan, A. M. Jawad, and M. Ismail, “Energy-efficient wireless sensor networks for precision agriculture: A review,” _Sensors_ , vol. 17, no. 8, p. 1781, 2017.
  * [4] N. Podevijn, D. Plets, J. Trogh, L. Martens, P. Suanet, K. Hendrikse, and W. Joseph, “Tdoa-based outdoor positioning with tracking algorithm in a public lora network,” _Wireless Communications and Mobile Computing_ , vol. 2018, pp. 1–9, 2018.
  * [5] Y. Etiabi, M. Jouhari, A. Burg, and E. M. Amhoud, “Spreading factor assisted lora localization with deep reinforcement learning,” in _2023 IEEE 97th Vehicular Technology Conference (VTC2023-Spring)_ , pp. 1–5, 2023.
  * [6] M. Jouhari, N. Saeed, M.-S. Alouini, and E. M. Amhoud, “A survey on scalable lorawan for massive iot: recent advances, potentials, and challenges,” _IEEE Communications Surveys & Tutorials_, pp. 1841–1876, 2023.
  * [7] L.-A. Phan and T. Kim, “Enabling rapid time synchronization with slow-flooding in wireless sensor networks,” _IEEE Communications Letters_ , vol. 26, no. 4, pp. 947–951, 2022.
  * [8] X. Huan, H. He, T. Wang, Q. Wu, and H. Hu, “A timestamp-free time synchronization scheme based on reverse asymmetric framework for practical resource-constrained wireless sensor networks,” _IEEE Transactions on Communications_ , vol. 70, no. 9, pp. 6109–6121, 2022.
  * [9] T. Polonelli, D. Brunelli, and L. Benini, “Slotted aloha overlay on lorawan-a distributed synchronization approach,” in _IEEE 16th international conference on embedded and ubiquitous computing_ , pp. 129–132, 2018.
  * [10] T. Polonelli, D. Brunelli, A. Marzocchi, and L. Benini, “Slotted aloha on lorawan-design, analysis, and deployment,” _Sensors_ , vol. 19, no. 4, p. 838, 2019.
  * [11] L. Beltramelli, A. Mahmood, P. Österberg, M. Gidlund, P. Ferrari, and E. Sisinni, “Energy efficiency of slotted lorawan communication with out-of-band synchronization,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 70, pp. 1–11, 2021.
  * [12] “Lorawan specification v1.1,” https://resources.lora-alliance.org/technical-specifications/lorawan-specification-v1-1, accessed : 2023-08-29.
  * [13] N. Abramson, “The aloha system: Another alternative for computer communications,” in _Proceedings of the fall joint computer conference_ , pp. 281–285, 1970.
  * [14] R. Tjoa, K. L. Chee, P. Sivaprasad, S. Rao, and J. G. Lim, “Clock drift reduction for relative time slot tdma-based sensor networks,” in _IEEE 15th International Symposium on Personal, Indoor and Mobile Radio Communications (IEEE Cat. No. 04TH8754)_ , vol. 2, pp. 1042–1047, 2004.
  * [15] “Semtech sx1272/73 technical datasheet,” accessed: 2023-07-19.
  * [16] J. Tapparel, O. Afisiadis, P. Mayoraz, A. Balatsoukas-Stimming, and A. Burg, “An open-source lora physical layer prototype on gnu radio,” in _IEEE 21st International Workshop on Signal Processing Advances in Wireless Communications_ , pp. 1–5, 2020.
  * [17] “Mini-lorawan implementation,” https://github.com/dossam/mini-lorawan/, accessed: 2024-04-03.


# Structure and production mechanism of the enigmatic $X(3872)$ in high-energy
hadronic reactions

Anna Cisek<EMAIL_ADDRESS>College of Natural Sciences, Institute of
Physics, University of Rzeszów, ul. Pigonia 1, PL-35-310 Rzeszów, Poland
Wolfgang Schäfer<EMAIL_ADDRESS>Institute of Nuclear Physics,
Polish Academy of Sciences, ul. Radzikowskiego 152, PL-31-342 Kraków, Poland
Antoni Szczurek111also at University of Rzeszów, PL-35-959 Rzeszów, Poland
<EMAIL_ADDRESS>Institute of Nuclear Physics, Polish Academy of
Sciences, ul. Radzikowskiego 152, PL-31-342 Kraków, Poland

###### Abstract

We calculate the total cross section and transverse momentum distributions for
the production of enigmatic $\chi_{c,1}(3872)$ (or X(3872)) assuming different
scenarios: $c\bar{c}$ state and $D^{0*}{\bar{D}}^{0}+D^{0}{\bar{D}}^{0*}$
molecule. The derivative of the $c\bar{c}$ wave function needed in the first
scenario is taken from a potential model calculations. Compared to earlier
calculation of molecular state we include not only single parton scattering
(SPS) but also double parton scattering (DPS) contributions. The latter one
seems to give smaller contribution than the SPS one. The upper limit for the
DPS production of $\chi_{c,1}(3872)$ is much below the CMS data. We compare
results of our calculations with existing experimental data of CMS, ATLAS and
LHCb collaborations. Reasonable cross sections can be obtained in either
$c\bar{c}$ or molecular $D{\bar{D}}^{*}$ scenarios for $X(3872)$. Also a
hybrid scenario is not excluded.

##  I Introduction

The $X(3872)$ state was discovered already some time ago by the Belle
collaboration Belle_X . Since then its existence has been confirmed in several
other processes and numerous theoretical studies have been performed, see for
example the review articles Lebed:2016hpi ; Olsen:2017bmm ; Brambilla:2019esw
. There is at present agreement that the $X(3872)$ has the axial vector
quantum numbers $J^{PC}=1^{++}$, and accordingly the state is named as
$\chi_{c1}(3872)$ PDG .

The internal structure of $X(3872)$ stays rather enigmatic. While its quantum
numbers are not exotic – it could indeed be a quarkonium $c\bar{c}$ state,
e.g. a radial excitation of the $\chi_{c1}$, there are strong arguments for a
non-$c\bar{c}$ scenario, manifested e.g. by the violation of isospin in its
decays PDG .

More importantly, the mass of $X$ is very close to the $D{\bar{D}}^{*}$
threshold. It is therefore rather popular to consider $X(3872)$ as very weakly
bound state of the $D{\bar{D}}^{*}$ system - a hadronic molecule, see the
review Guo:2017jvc . A tetraquark scenario was considered in MPPR2005 . The
$c\bar{c}$ quarkonium scenario, where $X(3872)$ is the $\chi_{c1}(2P)$ state
has been advocated in Achasov:2015oia . Other approaches treat the $X(3872)$
as a $c\bar{c}$ bound state in the meson-meson continuum taking into account
coupling of $c\bar{c}$ and meson-meson channels Cincioglu:2016fkm ; CRB2013 ;
Giacosa:2019zxw . The possible mixture of quarkonium and molecule/virtual
state is considered in Kang:2016jxw  and found to be consistent with current
data.

Recently, the transverse momentum distributions of $X(3872)$ were measured at
the LHC CMS_X ; ATLAS_X ; LHCb_X  by the CMS, ATLAS and LHCb collaborations.

There is a debate in the literature Bignamini:2009sk ; Artoisenet:2009wk ;
Albaladejo:2017blx ; Wang:2017gay , whether the rather large production rate
at large $p_{T}$ allows one to exclude the molecular scenario. With few
exceptions Butenschoen:2013pxa ; MHC2013 ; Ilten:2021jcb  the discussion in
the literature is limited to estimates of orders of magnitude.

In this paper we shall consider two scenarios of prompt $X(3872)$ production,
which are not mutually exclusive and in fact both can contribute depending on
the structure of $X(3872)$. Both scenarios have in common that the production
is initiated by the production of a $c\bar{c}$ pair in a hard process.

In the first scenario we shall consider that the $X(3872)$ is a pure
$c\bar{c}$ state, the first radial excitation of $\chi_{c1}$. The
corresponding wave function and its derivative were calculated in potential
models e.g. in EQ2019 .

In the second scenario, where the $X(3872)$ is treated as a weakly bound
$s$-wave state in the $D\bar{D}^{*}+D^{*}\bar{D}$-system, we exploit the
connection between the low-energy scattering amplitude in the continuum and at
the bound state pole below threshold, well known from effective range theory.
In this work we follow Artoisenet:2009wk  and give an estimate of a
$p_{T}$-dependent upper bound for $X(3872)$ production in the molecule
scenario.

As far as the hard production mechanism is concerned, we employ the
$k_{T}$-factorization framework Catani:1990eg ; Levin:1991ry ; Collins:1991ty
. For the quarkonium scenario, the dominant mechanism of $C=+$1 is probably
color singlet two virtual gluon fusion. The production of
$\chi_{c0},\chi_{c1},\chi_{c2}$ production was considered e.g. in Samara ;
Baranov:2015yea ; CS2018 . Recently we have shown that the transverse momentum
distribution of $\eta_{c}$ measured by the LHCb LHCb_etac  can be nicely
described as $g^{*}g^{*}\to\eta_{c}$ fusion within $k_{T}$-factorization
approach BPSS2020 .

For the molecular scenario, in addition to the single-parton scattering (SPS)
mechanism of fusion of two off-shell gluons $g^{*}g^{*}\to c\bar{c}$, we will
also consider production through the double-parton scattering (DPS) mode.

##  II Formalism

In Fig.1 we show two generic Feynman diagrams for $X(3872)$ quarkonium
production in proton-proton collision via gluon-gluon fusion: for the
quarkonium (left) and molecule (right). These diagrams illustrate the
situation adequate for the $k_{T}$-factorization calculations used in the
present paper.

Figure 1: Generic diagrams for the inclusive process of $X(3872)$ production
in proton-proton scattering via two gluons fusion.

The inclusive cross section for $X(3872)$-production via the $2\to 1$ gluon-
gluon fusion mode is obtained from

$\displaystyle{d\sigma\over dyd^{2}\mbox{$\vec{p}_{T}$}}$ $\displaystyle=$
$\displaystyle\int{d^{2}\mbox{$\vec{q}_{T}$}_{1}\over\pi\mbox{$\vec{q}_{T}$}_{1}^{2}}{\cal{F}}(x_{1},\mbox{$\vec{q}_{T}$}_{1}^{2},\mu_{F}^{2})\int{d^{2}\mbox{$\vec{q}_{T}$}_{2}\over\pi\mbox{$\vec{q}_{T}$}_{2}^{2}}{\cal{F}}(x_{2},\mbox{$\vec{q}_{T}$}_{2}^{2},\mu_{F}^{2})\,\delta^{(2)}(\mbox{$\vec{q}_{T}$}_{1}+\mbox{$\vec{q}_{T}$}_{2}-\mbox{$\vec{p}_{T}$})\,$
(1) $\displaystyle\times$
$\displaystyle{\pi\over(x_{1}x_{2}s)^{2}}\overline{|{\cal{M}}_{g^{*}g^{*}\to
X(3872)}|^{2}}\,.$

Here the matrix element squared for the fusion of two off-shell gluons into
the ${}^{3}P_{1}$ color singlet $c\bar{c}$ charmonium is (see e.g. Samara ;
Cisek:2017ikn  for a derivation):

$\displaystyle\overline{|{\cal{M}}_{g^{*}g^{*}\to X(3872)}|^{2}}$
$\displaystyle=$
$\displaystyle(4\pi\alpha_{S})^{2}\,{4|R^{\prime}(0)|^{2}\over\pi
M_{X}^{3}}{\mbox{$\vec{q}_{T}$}_{1}^{2}\mbox{$\vec{q}_{T}$}_{2}^{2}\over(M_{X}^{2}+\mbox{$\vec{q}_{T}$}_{1}^{2}+\mbox{$\vec{q}_{T}$}_{2}^{2})^{4}}$
(2) $\displaystyle\times$
$\displaystyle\Big{(}(\mbox{$\vec{q}_{T}$}_{1}^{2}+\mbox{$\vec{q}_{T}$}_{2}^{2})^{2}\sin^{2}\phi+M_{X}^{2}(\mbox{$\vec{q}_{T}$}_{1}^{2}+\mbox{$\vec{q}_{T}$}_{2}^{2}-2|\mbox{$\vec{q}_{T}$}_{1}||\mbox{$\vec{q}_{T}$}_{2}|\cos\phi)\Big{)}\,,$

where $\phi$ is the azimuthal angle between
$\mbox{$\vec{q}_{T}$}_{1},\mbox{$\vec{q}_{T}$}_{2}$. The momentum fractions of
gluons are fixed as $x_{1,2}=m_{T}\exp(\pm y)/\sqrt{s}$, where
$m_{T}^{2}=\mbox{$\vec{p}_{T}$}^{2}+M_{X}^{2}$. The derivative of the radial
quarkonium wave function at the origin is taken for the first radial $p$-wave
excitation from Ref.EQ2019 , $|R^{\prime}(0)|^{2}=0.1767\,{\rm GeV}^{5}$.

The unintegrated gluon parton distribution functions (gluon uPDFs) are
normalized such, that the collinear glue is obtained from

$\displaystyle
xg(x,\mu_{F}^{2})=\int^{\mu_{F}^{2}}{d^{2}\mbox{$\vec{k}_{T}$}\over\pi\mbox{$\vec{k}_{T}$}^{2}}{\cal{F}}(x,\mbox{$\vec{k}_{T}$}^{2},\mu_{F}^{2})\,.$
(3)

The hard scale is taken to be always $\mu_{F}=m_{T}$, the transverse mass of
the $X(3872)$. In order to estimate the production cross section for the
molecule we also start from a hard production of a $c\bar{c}$-pair, which we
then hadronize into a $D\bar{D}^{*}+h.c$ system using a prescription given
below.

The parton-level differential cross section for the $c\bar{c}$ production,
formally at leading-order, reads:

$\displaystyle\frac{d\sigma(pp\to
Q\bar{Q}\,X)}{dy_{1}dy_{2}d^{2}\mbox{$\vec{p}_{T}$}_{1}d^{2}\mbox{$\vec{p}_{T}$}_{2}}$
$\displaystyle=$
$\displaystyle\int\frac{d^{2}\mbox{$\vec{k}_{T}$}_{1}}{\pi\mbox{$\vec{k}_{T}$}^{2}_{1}}{\cal
F}(x_{1},\mbox{$\vec{k}_{T}$}_{1}^{2},\mu_{F}^{2})\int\frac{d^{2}\mbox{$\vec{k}_{T}$}_{2}}{\pi\mbox{$\vec{k}_{T}$}^{2}_{2}}{\cal
F}(x_{2},\mbox{$\vec{k}_{T}$}_{2}^{2},\mu_{F}^{2})$ (4)
$\displaystyle\times\delta^{(2)}\left(\mbox{$\vec{k}_{T}$}_{1}+\mbox{$\vec{k}_{T}$}_{2}-\mbox{$\vec{p}_{T}$}_{1}-\mbox{$\vec{p}_{T}$}_{2}\right)\
\frac{1}{16\pi^{2}(x_{1}x_{2}s)^{2}}\;\overline{|{\cal M}^{\mathrm{off-
shell}}_{g^{*}g^{*}\to c\bar{c}}|^{2}}\,.$

where ${\cal M}^{\mathrm{off-shell}}_{g^{*}g^{*}\to Q\bar{Q}}$ is the off-
shell matrix element for the hard subprocess Catani:1990eg , we use its
implementation from maciula .

Here, one keeps exact kinematics from the very beginning and additional hard
dynamics coming from transverse momenta of incident partons. Explicit
treatment of the transverse momenta makes the approach very efficient in
studies of correlation observables. The two-dimensional Dirac delta function
assures momentum conservation. The gluon uPDFs must be evaluated at
longitudinal momentum fractions:

$\displaystyle x_{1}$ $\displaystyle=$
$\displaystyle\frac{m_{T1}}{\sqrt{s}}\exp(+y_{1})+\frac{m_{T2}}{\sqrt{s}}\exp(+y_{2})\;,$
(5) $\displaystyle x_{2}$ $\displaystyle=$
$\displaystyle\frac{m_{T1}}{\sqrt{s}}\exp(-y_{1})+\frac{m_{T2}}{\sqrt{s}}\exp(-y_{2}),$
(6)

where $m_{Ti}=\sqrt{p_{Ti}^{2}+m_{c}^{2}}$ is the quark/antiquark transverse
mass.

In the present analysis we employ the heavy $c$-quark approximation and assume
that three-momenta in the $pp$-cm frame are equal:

$\vec{p}_{D}=\vec{p}_{c}\;.$ (7)

This approximation could be relaxed in future. The hadronization is then
included only via fragmentation branching fractions:

$\displaystyle P(c\to D^{0})=P({\bar{c}}\to{\bar{D}}^{0})=0.56\;,$ (8)
$\displaystyle P(c\to D^{+})=P({\bar{c}}\to D^{-})=0.225\;,$ (9)
$\displaystyle P(c\to D^{*0})=P({\bar{c}}\to{\bar{D}}^{*0})=0.236\;,$ (10)
$\displaystyle P(c\to D^{*+})=P({\bar{c}}\to{\bar{D}}^{*-})=0.236\;.$ (11)

The first number is from Lisovyi:2015uqa  while the other numbers are from
ATLAS_D-mesons .

The cross section for $c\bar{c}$ production are then multiplied by

$\displaystyle\frac{1}{2}[P(c\to D^{0})P({\bar{c}}\to{\bar{D}}^{*0})+P(c\to
D^{*0})P({\bar{c}}\to{\bar{D}}^{0})]=0.1322$ (12)

and, when comparing to experimental data by the relevant branching fraction
i.e. by ${\rm BR}(X\to J/\psi\pi^{+}\pi^{-})$ for the case of CMS and LHCb,
and by ${\rm BR}(X\to J/\psi\pi^{+}\pi^{-}){\rm BR}(J/\psi\to\mu^{+}\mu^{-})$
for the ATLAS data. The branching fractions are taken from PDG . The factor
$\frac{1}{2}$ is related to the factor $\frac{1}{\sqrt{2}}$ in the definition
of the molecular wave function,

$|{\Psi_{mol}}\rangle=\frac{1}{\sqrt{2}}\left(|{D^{0}{\bar{D}}^{*0}}\rangle+|{D^{*0}{\bar{D}}^{0}}\rangle\right)\;.$
(13)

In our calculation, we control the dependence on the relative momentum of
quark and antiquark in the rest frame of the pair:

$k_{rel}={1\over 2}\sqrt{M^{2}_{c\bar{c}}-4m_{c}^{2}}\;,$ (14)

where $M_{c\bar{c}}$ is invariant mass of the $c\bar{c}$ system and $m_{c}$ is
the quark mass. In order to obtain an upper bound for the molecule production
cross section, we should integrate the $D\bar{D}^{*}$ cross section over the
relative momentum $k^{D\bar{D}^{*}}_{rel}$ up to a cutoff $k^{D\bar{D}}_{max}$
Artoisenet:2009wk . We will instead impose a cutoff $k_{max}$ on the relative
momentum $k_{rel}$. Within our kinematics the latter will be similar to
$k_{rel}^{D\bar{D}^{*}}$, but somewhat larger. In reality, for larger
$\mbox{$\vec{p}_{T}$}_{X}\approx\mbox{$\vec{p}_{T}$}_{1}+\mbox{$\vec{p}_{T}$}_{2}$,
we have $k_{DD}<k_{rel}$. We therefore estimate, that $k_{max}=0.2\,{\rm GeV}$
corresponds roughly to $k^{DD}_{max}\approx 0.13\,\rm{GeV}$. A better
approximation would be to add simultaneous $c\to D$, ${\bar{c}}\to{\bar{D}}$
fragmentation to our Monte Carlo code, which however means at least two more
integrations. We do not consider here any model of the $D{\bar{D}}^{*}$ wave
function.

What is the appropriate choice for $k^{D\bar{D}^{*}}_{max}$ was a matter of
discussion in the literature. In Ref. Bignamini:2009sk  it was suggested, that
$k^{D\bar{D}^{*}}_{max}$ should be of the order of the binding momentum
$k_{X}=\sqrt{2\mu\varepsilon_{X}}$, where $\mu$ is the reduced mass of the
$D\bar{D}^{*}$-system, and $\varepsilon_{X}$ is the binding energy of
$X(3872)$. This would lead to a very small value for $k^{D\bar{D}^{*}}_{max}$,
similar to $k^{D\bar{D}^{*}}_{max}=35\,\rm{MeV}$ used in Bignamini:2009sk .
However problems with this estimate have been pointed out in Artoisenet:2009wk
; Albaladejo:2017blx ; Wang:2017gay . As has been argued in these works, the
integral should be extended rather to a scale $k^{D\bar{D}^{*}}_{\max}\sim
m_{\pi}$. In our choice of $k_{max}$, we follow this latter prescription.

In the following for illustration we shall therefore assume
$k_{max}=0.2\,\rm{GeV}$. The calculation for the SPS molecular scenario is
done using the VEGAS algorithm for Monte-Carlo integration Lepage:1977sw .

We also include double parton scattering contributions (see Fig.2).

Figure 2: A generic diagram for the inclusive process of $X(3872)$ production
in proton-proton scattering via the double parton scattering mode.

The corresponding cross section is calculated in the so-called factorized
ansatz as:

$\Delta\sigma=\frac{1}{2\sigma_{\rm
eff}}\int\frac{d\sigma_{c\bar{c}}}{dy_{1}d^{2}\mbox{$\vec{p}_{T}$}_{1}}\frac{d\sigma_{c\bar{c}}}{dy_{2}d^{2}\mbox{$\vec{p}_{T}$}_{2}}\;\;dy_{1}d^{2}\mbox{$\vec{p}_{T}$}_{1}dy_{2}d^{2}\mbox{$\vec{p}_{T}$}_{2}\Big{|}_{k_{rel}<k_{max}}\;.$
(15)

Above the differential distributions of the first and second parton scattering
$\frac{d\sigma}{dy_{i}d^{2}\mbox{$\vec{p}_{T}$}_{i}}$ are calculated in the
$k_{T}$-factorization approach as explained above. In the following we take
$\sigma_{\rm eff}$ = 15 mb as in Maciula:2013kd . The differential
distributions (in $p_{T}$ of the $X(3872)$ or $y_{\rm diff}=y_{1}-y_{2}$,
etc.) are obtained by binning in the appropriate variable. We include all
possible fusion combinations leading to $X(3872)$:

$\displaystyle c_{1}\to D^{0},{\bar{c}}_{2}\to{\bar{D}}^{*0}\;,$ (16)
$\displaystyle c_{1}\to D^{*0},{\bar{c}}_{2}\to{\bar{D}}^{0}\;,$ (17)
$\displaystyle{\bar{c}}_{1}\to{\bar{D}}^{0},c_{2}\to D^{*0}\;,$ (18)
$\displaystyle{\bar{c}}_{1}\to{\bar{D}}^{*0},c_{2}\to D^{0}\;.$ (19)

This leads to the multiplication factor two times bigger than for the SPS
contribution (see Eq.(12)).

##  III Results

In this section we shall show our results for recent CMS CMS_X  ATLAS ATLAS_X
and LHCb LHCb_X  data. The CMS data is for $\sqrt{s}$ = 7 TeV and -1.2
$<y_{X}<$ 1.2, the ATLAS data for $\sqrt{s}=$ 8 TeV, -0.75 $<y_{X}<$ 0.75 and
the LHCb data for $\sqrt{s}=$ 13 TeV, 2 $<y_{X}<$ 4.5. In all cases the
$X(3872)$ was measured in the $J/\psi\pi^{+}\pi^{-}$ channel. We have used a
number of different unintegrated gluon distributions, firstly a distribution
obtained from the solution of a BFKL equation with kinematic constraints by
Kutak and Staśto denoted KS Kutak:2004ym , secondly a gluon uPDF obtained from
a modified Kimber-Martin-Ryskin (KMR) procedure Kimber:2001sc ; Martin:2009ii
based on Durham group collinear PDFs Harland-Lang:2014zoa . Finally we also
employ a gluon uPDF obtained by Hautmann and Jung Hautmann:2013tba  from a
description of precise HERA data on deep inelastic structure function by a
solution of the CCFM evolution equations.

###  III.1 $c\bar{c}$ state

Here we show our results for $\chi_{c,1}(3872)$ production treated as a pure
$c\bar{c}$ state. In Figs.3, 4,5 we show the transverse momentum distribution
of the $X(3872)$ together with the data from the CMS, ATLAS and LHCb
experimental data. A surprisingly good description is obtained with different
gluon uPDFs specified in the figure legend without any free parameters. It is
worth to mention approximately good slope of the $p_{TX}$ distributions which
is due to effective inclusion of higher-order corrections. The corresponding
result within the collinear leading-order approximation would be equal to
zero! A slightly different slope is obtained for the molecular scenario (solid
line) discussed in detail somewhat below.

Figure 3: Transverse momentum distribution of $X(3872)$ for the CMS
experiment. Shown are results for 3 different gluon uPDFs. Here BR = 0.038.
The upper limit for the SPS molecular scenario is shown as the solid line. We
also show corresponding distribution for the DPS mechanism (dotted line).
Figure 4: Transverse momentum distribution of $X(3872)$ for the ATLAS
experiment. Shown are results for 3 different gluon uPDFs. Here ${\rm
BR}=0.038\cdot 0.0596$. The upper limit for the SPS molecular scenario is
shown as the solid line.  Figure 5: Transverse momentum distribution of
$X(3872)$ for the LHCb experiment. Shown are results for 3 different gluon
uPDFs. Here BR = 0.038. The upper limit for the SPS molecular scenario is
shown as the solid line.

###  III.2 Molecular picture

Here we show our predictions for $\chi_{c,1}(3872)$ production treated as the
$\bar{D}^{*0}D^{0}+h.c.$ molecule. Then the $\chi_{c,1}(3872)$ can be produced
if $k_{rel}$ is small, and an estimate (or upper bound) for its production
cross section can be obtained from the continuum cross section at small
$k_{rel}$. In Fig.6 we show the distribution in this variable for different
windows of $|\mbox{$\vec{p}_{T}$}_{c\bar{c}}|$ (0,5 GeV), (5,10 GeV), (10,15)
GeV, (15,20) GeV, (20,25) GeV, (25,30) GeV.

Similar shapes in $k_{rel}$ are obtained for the different windows of
$p_{Tc\bar{c}}$.

Figure 6: Distribution in $k_{rel}$ for different windows of
$p_{t,c\bar{c}}=p_{t,X}$ (left panel) for the CMS kinematics as specified in
the main text. In the right panel we show the cross sections divided by
$k_{rel}^{2}$. In these calculations the KMR UGDF with the MMHT NLO collinear
gluon distribution was used. Branching fractions are included here.

To visualize this better we show in the right panel of Fig.6
$k_{rel}^{-2}d\sigma/dk_{rel}$. As expected, the so-obtained distributions
closely follow phase-space, and are almost flat in a broad range of $k_{rel}$.
Therefore, the cross section has essentially the phase-space behaviour

$\displaystyle d\sigma\propto k_{rel}^{2}dk_{rel}\,,$ (20)

which implies the strong, cubic dependence $\propto k_{max}^{3}$ on the upper
limit $k_{max}$ of the $k_{rel}$-integration.

The calculation in the whole phase space $p_{T1}\in(0,20)\,\rm{GeV}$ and
$p_{T2}\in(0,20)\,\rm{GeV}$ leads to fluctuations at large
$p_{Tc\bar{c}}>20\,\rm{GeV}$. This can be understood as due to steep
dependence of the cross section on $p_{T1},p_{T2},y_{1},y_{2},\phi$. Only a
narrow range in the $p_{T1}\otimes p_{T2}$ space with $p_{T1}\approx p_{T2}$
fulfills the condition $k_{rel}<k_{max}=0.2\,\rm{GeV}$.

We remind the reader that this value of $k_{max}$ is imposed on the $c\bar{c}$
final state, and would correspond to a smaller $k_{max}^{DD^{*}}\sim
0.14\,\rm{GeV}$ for the $D\bar{D}^{*}$ mesons. Due to the behaviour shown in
Eq.(20), the cross section for $k_{max}=0.14\,\rm{GeV}$ or
$k_{max}=0.1\,\rm{GeV}$ imposed on the $c\bar{c}$ state would go down by a
factor three or eight, repsectively.

The distribution in relative azimuthal angle between $c\bar{c}$ in the
$pp$-frame with the cut $k_{rel}<0.2\,\rm{GeV}$ is shown in Fig.7. This is
rather a steep distribution around $\phi$ = 0o. This is not a typical region
of the phase space. Recall, that in the leading-order collinear approach
$\phi$ = 180o. It is obvious that the region of $\phi\approx$ 0 cannot be
obtained easily in the collinear approach. As discussed in maciula  the
$k_{T}$-factorization approach gives a relatively good description of
$D^{0}{\bar{D}}^{0}$ correlations in this region of the phase space.

Figure 7: Azimuthal correlations between $c\bar{c}$ that fullfill the
condition $k_{rel}<$ 0.2 GeV. Here the CMS cuts were imposed. We show
contribution of SPS (solid line) and DPS (dashed line).

Having understood the kinematics of the $X(3872)$ production we can improve
the description of transverse momentum distribution of $X(3872)$. The rather
strong correlation $p_{T1}\approx p_{T2}$ allows to perform VEGAS calculation
simultaneously limiting to $p_{T1},p_{T2}>p_{Tmin}$. We have also verified
that

$p_{T1}+p_{T2}\approx
p_{Tc\bar{c}}=|\mbox{$\vec{p}_{T}$}_{1}+\mbox{$\vec{p}_{T}$}_{2}|\;.$ (21)

Therefore imposing lower limits on $p_{T1}>p_{Tmin}$ and $p_{T2}>p_{Tmin}$
means also lower limit on $p_{Tc\bar{c}}>2p_{Tmin}$. The solid line in Fig.3
shows result of such a calculation. The fluctuations are gone. Combining the
two calculations at say $p_{Tc\bar{c}}=15\,\rm{GeV}$ gives a smooth result
everywhere.

###  III.3 Hybrid approach

Still another option is to assume that $X$ is the combination:

$|{X(3782)}\rangle=\alpha|{c\bar{c}}\rangle+\frac{\beta}{\sqrt{2}}\left(|{D{\bar{D}}^{*}}\rangle+|{{\bar{D}}D^{*}}\rangle\right)\;.$
(22)

Such a state could be called hybrid. Since the cross section for production of
conventional and molecular $X(3872)$ are rather similar, it is difficult to
select one combination of $\alpha$ and $\beta$. For example in Fig.8 we show
also such a result for $\alpha^{2}=\beta^{2}$ = 0.5 (dash-dotted line). The
corresponding result is similar to the result for $\alpha^{2}$ = 1,
$\beta^{2}$ = 0 (pure $c\bar{c}$ state) and $\alpha^{2}$ = 0, $\beta^{2}$ = 1
(pure molecular state) calculated for the KMR UGDF. The fifty-fifty mixture
gives quite good representation of all the experimental data.

Figure 8: Transverse momentum distribution of $X(3872)$ for the CMS, ATLAS and
LHCb experiments. Shown are results for the KMR UGDF. Here BR = 0.038 for CMS
and LHCb, and ${\rm BR}=0.038\cdot 0.0596$ for ATLAS. We show results for
different combinations of $\alpha$ and $\beta$ as specified in the figure
legend.

In summary, it is very difficult at present to definitely conclude what is the
mechanism of the $X(3872)$ production. The same is true for the wave function
of the $X(3872)$ meson (see Eq.(22)). At present we cannot answer the question
whether $X(3872)$ is of conventional, molecular or hybrid type. More work is
definitely needed in future.

##  IV Conclusions

We have performed the calculation of $X(3872)$ production at the LHC energies.
We have performed two independent calculations: one within nonrelativistic QCD
approach assuming pure $c\bar{c}$ state and second assuming a coalescence of
$D$ and ${\bar{D}}^{*}$ or $\bar{D}$ and $D^{*}$, consistent with molecular
state assumption. The first calculation requires usage of derivative of the
$c\bar{c}$ wave function. In the present analysis we have used the wave
function obtained in EQ2019 . The resulting cross section was calculated
within the $k_{T}$-factorization approach with a few unintegrated gluon
distributions. In the second approach first the hard production of a
$c\bar{c}$ pair is calculated. Next a simple hadronization is performed giving
a correlation distribution of $D$ and $D^{*}$ mesons. Imposing limitations
(upper limit) on relative momenta of $D$ and ${\bar{D}}^{*}$ we get a
$p_{T}$-dependent upper limit of the cross section for $D$-${\bar{D}}^{*}$ or
${\bar{D}}$-$D^{*}$ fusion (coalescence). We compare, for the first time, to
all available experimental data on the $p_{T}$-dependent cross section. The
both, quite different, scenarios give the cross section of the right order of
magnitude, very similar to the CMS, ATLAS and LHCb experimental data.

Also a mixture of both mechanism leads to correct description of the CMS,
ATLAS and LHCb data. Our study within $k_{T}$-factorization approach shows
that all the options are in the game.

###### Acknowledgements.

This work is partially supported by the Polish National Science Centre under
Grant No. 2018/31/B/ST2/03537 and by the Center for Innovation and Transfer of
Natural Sciences and Engineering Knowledge in Rzeszów (Poland).

## References

  * (1) S.K. Choi et al. (Belle collaboration), “Observation of a narrow charmonium-like state in exclusive $B^{+}\to K^{+}\pi^{+}\pi^{-}J/\psi$ decays”, Phys. Rev. Lett. 91 (2003) 262001.
  * (2) R. F. Lebed, R. E. Mitchell and E. S. Swanson, “Heavy-Quark QCD Exotica,” Prog. Part. Nucl. Phys. 93 (2017), 143-194 [arXiv:1610.04528 [hep-ph]].
  * (3) S. L. Olsen, T. Skwarnicki and D. Zieminska, “Nonstandard heavy mesons and baryons: Experimental evidence,” Rev. Mod. Phys. 90 (2018) no.1, 015003 [arXiv:1708.04012 [hep-ph]].
  * (4) N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C. P. Shen, C. E. Thomas, A. Vairo and C. Z. Yuan, “The $XYZ$ states: experimental and theoretical status and perspectives,” Phys. Rept. 873 (2020), 1-154 [arXiv:1907.07583 [hep-ex]].
  * (5) P.A. Żyła et al. (Particle Data Group), “The Review of Particle Physics (2021)”, Prog. Theor. Exp. Phys. 2020 (2020) no.8, 083C01 and 2021 update.
  * (6) F. K. Guo, C. Hanhart, U. G. Meißner, Q. Wang, Q. Zhao and B. S. Zou, “Hadronic molecules,” Rev. Mod. Phys. 90 (2018) no.1, 015004 [arXiv:1705.00141 [hep-ph]].
  * (7) L. Maiani, F. Piccini, A. D. Polosa and V. Riquer, “Diquark-antidiquark with hidden or open charm and the nature of $X(3872)$”, Phys. Rev. D71 (2005) 014028.
  * (8) N. N. Achasov and E. V. Rogozina, “$X(3872),I^{G}(J^{PC})=0^{+}(1^{++})$, as the $\chi_{1c}(2P)$ charmonium,” Mod. Phys. Lett. A 30 (2015) no.33, 1550181 [arXiv:1501.03583 [hep-ph]].
  * (9) E. Cincioglu, J. Nieves, A. Ozpineci and A. U. Yilmazer, “Quarkonium Contribution to Meson Molecules,” Eur. Phys. J. C 76 (2016) no.10, 576 [arXiv:1606.03239 [hep-ph]].
  * (10) S. Coito, G. Rupp and E. van Beveren, “$X(3872)$ is not a true molecule”, Eur. Phys. J. C73 (2013) 2351.
  * (11) F. Giacosa, M. Piotrowska and S. Coito, “$X(3872)$ as virtual companion pole of the charm–anticharm state $\chi_{c1}(2P)$,” Int. J. Mod. Phys. A 34 (2019) no.29, 1950173 [arXiv:1903.06926 [hep-ph]].
  * (12) X. W. Kang and J. A. Oller, “Different pole structures in line shapes of the $X(3872)$,” Eur. Phys. J. C 77 (2017) no.6, 399 [arXiv:1612.08420 [hep-ph]].
  * (13) S. Chatrchyan et al. [CMS], “Measurement of the $X$(3872) Production Cross Section Via Decays to $J/\psi\pi^{+}\pi^{-}$ in $pp$ collisions at $\sqrt{s}$ = 7 TeV,” JHEP 04 (2013), 154 [arXiv:1302.3968 [hep-ex]].
  * (14) M. Aaboud et al. (ATLAS collaboration), “Measurement of $\psi(2S)$ and $X(3872)\to J/\psi\pi^{+}\pi^{-}$ production in pp collisions at $\sqrt{s}$ = 8 TeV with the ATLAS detector”, JHEP01(2017) 117.
  * (15) R. Aaij et al. [LHCb], “Measurement of $\chi_{c1}(3872)$ production in proton-proton collisions at $\sqrt{s}=8$ and $13$ TeV,” [arXiv:2109.07360 [hep-ex]].
  * (16) C. Bignamini, B. Grinstein, F. Piccinini, A. D. Polosa and C. Sabelli, “Is the X(3872) Production Cross Section at Tevatron Compatible with a Hadron Molecule Interpretation?,” Phys. Rev. Lett. 103 (2009), 162001 [arXiv:0906.0882 [hep-ph]].
  * (17) P. Artoisenet and E. Braaten, “Production of the X(3872) at the Tevatron and the LHC,” Phys. Rev. D 81 (2010), 114018 [arXiv:0911.2016 [hep-ph]].
  * (18) M. Albaladejo, F. K. Guo, C. Hanhart, U. G. Meißner, J. Nieves, A. Nogga and Z. Yang, “Note on X(3872) production at hadron colliders and its molecular structure,” Chin. Phys. C 41 (2017) no.12, 121001 [arXiv:1709.09101 [hep-ph]].
  * (19) W. Wang, “On the production of hidden-flavored hadronic states at high energy,” Chin. Phys. C 42 (2018) no.4, 043103 [arXiv:1709.10382 [hep-ph]].
  * (20) M. Butenschoen, Z. G. He and B. A. Kniehl, “NLO NRQCD disfavors the interpretation of X(3872) as $\chi_{c1}(2P)$,” Phys. Rev. D 88 (2013), 011501 [arXiv:1303.6524 [hep-ph]].
  * (21) C. Meng, H. Han and K. T. Chao, “X(3872) and its production at hadron colliders,” Phys. Rev. D 96 (2017) no.7, 074014 [arXiv:1304.6710 [hep-ph]].
  * (22) P. Ilten and M. Utheim, “Forming molecular states with hadronic rescattering,” Eur. Phys. J. A 58 (2022) no.1, 1 [arXiv:2108.03479 [hep-ph]].
  * (23) E.J. Eichten and C. Quigg, “Quarkonium wave function at the origin: an update”, arXiv:1904.11542.
  * (24) S. Catani, M. Ciafaloni and F. Hautmann, “High-energy factorization and small x heavy flavor production,” Nucl. Phys. B 366 (1991), 135-188
  * (25) E. M. Levin, M. G. Ryskin, Y. M. Shabelski and A. G. Shuvaev, “Heavy quark production in semihard nucleon interactions,” Sov. J. Nucl. Phys. 53 (1991), 657
  * (26) J. C. Collins and R. K. Ellis, “Heavy quark production in very high-energy hadron collisions,” Nucl. Phys. B 360 (1991), 3-30
  * (27) B.A. Kniehl, D.V. Vasin and V.A. Saleev, “Charmonium production at high energy in the $k_{t}$-factorization approach”, Phys. Rev. D73, 074022 (2006).
  * (28) S. P. Baranov, A. V. Lipatov and N. P. Zotov, “Prompt charmonia production and polarization at LHC in the NRQCD with kt-factorization. Part II: $\chi_{c}$ mesons,” Phys. Rev. D 93 (2016) no.9, 094012 [arXiv:1510.02411 [hep-ph]].
  * (29) A. Cisek and A. Szczurek, “Prompt inclusive production of $J/\psi$ , $\psi^{\prime}$ and $\chi_{c}$ mesons at the LHC in forward directions within the NRQCD $k_{t}$-factorization approach: Search for the onset of gluon saturation”, Phys. Rev. D97 (2018) 034035.
  * (30) R. Aaij et al. [LHCb], “Measurement of the $\eta_{c}(1S)$ production cross-section in proton-proton collisions via the decay $\eta_{c}(1S)\rightarrow p\bar{p}$,” Eur. Phys. J. C 75 (2015) no.7, 311 [arXiv:1409.3612 [hep-ex]].
  * (31) I. Babiarz, R. Pasechnik, W. Schäfer and A. Szczurek, “Hadroproduction of $\eta_{c}(1S,2S)$ in the $k_{T}$-factorization approach”, JHEP 2002 (2020) 037.
  * (32) A. Cisek, W. Schäfer and A. Szczurek, “Production of $\chi_{c}$ pairs with large rapidity separation in $k_{T}$ factorization,” Phys. Rev. D 97 (2018) no.11, 114018 [arXiv:1711.07366 [hep-ph]].
  * (33) R. Maciula and A. Szczurek, “Open charm production at the LHC - $k_{t}$-factorization approach”, Phys. Rev. D87 (2013) 094022.
  * (34) M. Lisovyi, A. Verbytskyi and O. Zenaiev, “Combined analysis of charm-quark fragmentation-fraction measurements,” Eur. Phys. J. C 76, no.7, 397 (2016) [arXiv:1509.01061 [hep-ex]].
  * (35) G Aad et al. (ATLAS collaboration), “Measurement of $D^{*\pm}$, $D^{\pm}$ and $D_{s}^{\pm}$ meson production cross section in $pp$ collisions at $\sqrt{s}$ = 7 TeV with the ATLAS detector”, Nucl. Phys. B907 (2016) 717.
  * (36) G. P. Lepage, “A New Algorithm for Adaptive Multidimensional Integration,” J. Comput. Phys. 27, 192 (1978).
  * (37) R. Maciula and A. Szczurek, “Production of $c\bar{c}c\bar{c}$ in double-parton scattering within $k_{t}$-factorization approach – meson-meson correlations,” Phys. Rev. D 87, no.7, 074039 (2013) [arXiv:1301.4469 [hep-ph]].
  * (38) K. Kutak and A. M. Stasto, “Unintegrated gluon distribution from modified BK equation,” Eur. Phys. J. C 41 (2005), 343-351 [arXiv:hep-ph/0408117 [hep-ph]].
  * (39) M. A. Kimber, A. D. Martin and M. G. Ryskin, “Unintegrated parton distributions,” Phys. Rev. D 63, 114027 (2001) [hep-ph/0101348].
  * (40) A. D. Martin, M. G. Ryskin and G. Watt, “NLO prescription for unintegrated parton distributions,” Eur. Phys. J. C 66, 163 (2010) [arXiv:0909.5529 [hep-ph]]; G. Watt, A. D. Martin and M. G. Ryskin, “Unintegrated parton distributions and electroweak boson production at hadron colliders,” Phys. Rev. D 70, 014012 (2004) Erratum: [Phys. Rev. D 70, 079902 (2004)] [hep-ph/0309096].
  * (41) L. A. Harland-Lang, A. D. Martin, P. Motylinski and R. S. Thorne, “Parton distributions in the LHC era: MMHT 2014 PDFs,” Eur. Phys. J. C 75, no. 5, 204 (2015) [arXiv:1412.3989 [hep-ph]].
  * (42) F. Hautmann and H. Jung, “Transverse momentum dependent gluon density from DIS precision data,” Nucl. Phys. B 883, 1 (2014) [arXiv:1312.7875 [hep-ph]].


# Classifying Graphemes in English Words Through the Application of a Fuzzy
Inference System

Samuel Rose  School of Computer Science
University of Hull
Hull, UK
<EMAIL_ADDRESS>Chandrasekhar Kambhampati  School of Computer
Science
University of Hull
Hull, UK
<EMAIL_ADDRESS>
###### Abstract

In Linguistics, a grapheme is a written unit of a writing system corresponding
to a phonological sound. In Natural Language Processing tasks, written
language is analysed through two different mediums, word analysis, and
character analysis. This paper focuses on a third approach, the analysis of
graphemes. Graphemes have advantages over word and character analysis by being
self-contained representations of phonetic sounds. Due to the nature of
splitting a word into graphemes being based on complex, non-binary rules, the
application of fuzzy logic [11] would provide a suitable medium upon which to
predict the number of graphemes in a word. This paper proposes the application
of a Fuzzy Inference System to split words into their graphemes. This Fuzzy
Inference System results in a correct prediction of the number of graphemes in
a word 50.18% of the time, with 93.51% being within a margin of +- 1 from the
correct classification. Given the variety in language, graphemes are tied with
pronunciation and therefore can change depending on a regional accent/dialect,
the +- 1 accuracy represents the impreciseness of grapheme classification when
regional variances are accounted for. To give a baseline of comparison, a
second method involving a recursive IPA mapping exercise using a pronunciation
dictionary was developed to allow for comparisons to be made.

##  I Introduction

A grapheme is the smallest functional unit of a writing system [2]. In
English, graphemes are the written representation of a phonetic sound,
phonemes. In English, a grapheme varies in length from 1-4 characters. These
work as the building blocks for words in English, giving both structure for
how a word is spelt, and representing a words pronunciation. The decoding of
words into graphemes is a fundamental process in learning English and is
included in the National Curriculum for Key Stage 1 Students in England [3].

While the English language contains 44 phonemes [5], there is no fixed number
of graphemes due to the mixed linguistic base of English words. For this
project, a figure of 284 graphemes is used. This includes 89 graphemes in the
``main system'' which account for the majority of words, and a further 195
graphemes that cover niche cases not covered in the ``main system'' [1].

###  I-A Grapheme Decoding Challenges

One of the challenges in decoding graphemes in English is, due to the vast
number of graphemes compared to phonemes, multiple graphemes represent each
phoneme and the same grapheme can represent multiple phonemic sounds depending
on the context. One example of this is the ⟨ei⟩ grapheme which represents
three different phonemes, /i:/ where it is sounded as a long e as in
``deceive'', /eI/ where it is sounded ay as in ``weigh'', and also
/$\mathcal{E}$/ where it is sounded as a short e as in ``heifer''. Another
challenge in decoding a word's graphemes is that characters in some graphemes
are contained as part of another grapheme. Using the ⟨ei⟩ grapheme as an
example again, ⟨ei⟩ is a grapheme consisting of the characters e and i, ⟨e⟩
and ⟨i⟩ are also single-character graphemes, and the characters ei are part of
the ⟨eir⟩ and ⟨eigh⟩ graphemes.

Because of these challenges in decoding graphemes, the decoding of a word's
graphemes does not fit into discrete categories. Due to this, we hypothesise
that by applying Fuzzy Set Theory, English words can be computationally
decoded into their graphemes by applying fuzzy sets and fuzzy inference
systems.

##  II Development of Fuzzy Sets

###  II-A Grapheme Count Probabilities

The first step in developing fuzzy sets for deciding a word's graphemes is to
work out how many graphemes are in a word. To achieve this task, frequency
analysis can be undertaken on an annotated corpora of English words with
graphemes. For this task the training corpus used is the Google 10,000 [10], a
corpora of the top 10,000 words ever searched in Google, with the
EnglishGrammarHere10000 as a validation corpus [6]. Each word in this list was
then decoded manually using the rules set out for decoding graphemes in the
National Curriculum [3]. The decision to choose the top 10,000 most common
English words is that according to analysis undertaken by the Oxford English
Corpus [9], the 7000 most common words in English account for approximately
90% of word usage, making a corpus of 10,000 a suitable size for computational
analysis.

To obtain frequency distributions for the number of graphemes in a word, the
corpus was split into subsets of words length 1-$n$ where $n$ is the longest
word in the corpus. In the corpus, the longest word was the word
``telecommunications'' at 18 characters long. In each subset, the number of
graphemes for each word, from 1 to $n$, is counted. After this process was
completed, the relative probabilities were calculated. This was calculated by
dividing the word by the number of words of length $n$ which are in the
dataset. These probabilities were then plotted, showing that the number of
characters in each grapheme roughly follows a normal distribution, shown in
Figure 1. Based on this result, the conclusion made is that pursuing this line
of experimentation has significant potential for yielding a positive result.

Figure 1: Probabilities of Grapheme count compared to word length

###  II-B Development of Fuzzy Sets

With the number of graphemes when compared to word length loosely following a
normal distribution, by calculating the mean, $\mu$, and standard deviation,
$\sigma$, for the number of characters, vowels and consonants for words of $x$
Graphemes, these can be seen, rounded to 3 decimal places, in Table I. Having
calculated the mean and standard deviation for words consisting of 1-14
graphemes, the first 12 of these have a standard deviation greater than zero
and can be approximated using a normal distribution. The formula of a normal
probability density distribution, $f(x)$ is as follows:

$f(x)=\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^{2}}$
(1)

TABLE I: Grapheme Means and Standard Deviations Grapheme Count | Characters | Vowels | Consonants
---|---|---|---
Mean |  | Standard
---
Deviation
Mean |  | Standard
---
Deviation
Mean |  | Standard
---
Deviation
1 | 2 | 0.508 | 1.219 | 0.608 | 0.781 | 0.608
2 | 3.062 | 0.680 | 1.245 | 0.593 | 1.817 | 0.671
3 | 4.202 | 0.820 | 1.573 | 0.651 | 2.629 | 0.720
4 | 5.556 | 0.897 | 2.031 | 0.730 | 3.525 | 0.758
5 | 6.853 | 0.966 | 2.550 | 0.752 | 4.303 | 0.848
6 | 8.063 | 1.036 | 3.058 | 0.810 | 5.005 | 0.913
7 | 9.221 | 1.072 | 3.537 | 0.891 | 5.684 | 0.962
8 | 10.290 | 1.099 | 4.020 | 0.869 | 6.270 | 1.014
9 | 11.298 | 1.011 | 4.560 | 0.892 | 6.738 | 0.989
10 | 12.18 | 0.967 | 4.924 | 0.882 | 7.260 | 0.925
11 | 13.396 | 1.067 | 5.438 | 0.873 | 7.958 | 0.898
12 | 14.125 | 0.957 | 5.625 | 0.806 | 8.5 | 0.632
13 | 14 | 0 | 6 | 0 | 8 | 0
14 | 18 |  | Div/0
---
Error
8 |  | Div/0
---
Error
10 |  | Div/0
---
Error

Using the Mean and Standard deviations calculated in Table I, frequency sets
for words with 1-12 graphemes were plotted, shown in Figure 2.

Figure 2: Frequency Sets for Grapheme count compared to word length

###  II-C Fuzzy Inference System

Having plotted the frequency sets for words with 1-12 graphemes, a Mamdani
Fuzzy Inference System was developed [7]. This system is designed to have
three inputs, word length, vowel count, and consonant count. These are
represented by membership functions with the mean and standard deviations
outlined above. Using a set of fuzzy rules a mapping is developed that links
the input functions to output functions. The result of this is an aggregate
set correlating how much the inputs relate to each output rule. This set is
then evaluated to result in a non-fuzzy output, in this case, the output is a
number, with a precision of 2 decimal points in the range [1-14].

####  II-C1 Input Membership Functions

Fuzzy Membership Functions for each feature were created using the mean and
standard deviations outlined above in Table I. These are Gaussian functions
that give a membership value in the range [0-1] depending on the input. These
curves were used for graphemes 2-12 for all three features. The sets for 1 &
14 graphemes were developed as Z and S-shaped curves to catch exceptions at
the beginning and end of all three workspaces. This decision was made based on
the characteristics of the corpus and the frequencies shown in Figure 1 which
show that all words of length one and a quarter of words length two have a
single grapheme. Conversely, all words length 18 in the corpus have 14
graphemes. The membership functions for all three features are shown in
Figures 3, 4, and 5. While Figure 3 is very distinct, there is a great deal of
similarity between Figures 4 and 5. The main difference being that the Vowel
sets have a substantially greater degree of overlap when compared to Consonant
sets.

Figure 3: Membership Functions for Grapheme count compared to word length

Figure 4: Membership Functions for Grapheme count compared to a words vowel
count

Figure 5: Membership Functions for Grapheme count compared to a words
consonant count

####  II-C2 Output Membership Functions

As each input membership function corresponds to a number of graphemes in a
word, the output membership functions are in a workspace between [1-14]
representing the number of graphemes in a word. Because of this, the output
membership functions also correspond to the number of graphemes predicted in a
word. Unlike the input membership functions, where the different functions
overlap, the output membership functions were defined as discrete triangular
functions peaking on the whole number of graphemes that function corresponds
to, with the tails placed on the edge of where a rounding algorithm would
round to that number (0.5 below and 0.499 above). These membership functions
are shown in Figure 6.

Figure 6: Membership Functions for Grapheme count

####  II-C3 Fuzzy Rules

Fuzzy rules follow an IF $x$ AND/OR $y$ THEN $z$ format, with $x$ and $y$
being fuzzy input membership functions and $z$ being the fuzzy output
membership function [4]. In the use case of this project, with the input
membership functions directly corresponding to the same feature as the output
membership functions, the rules in this fuzzy inference system follow a
simplistic methodology. With 13 input sets for each feature, the 13 rules
correspond to output membership functions 1-12, and 14. The rules used can be
seen in Figure 7.

Figure 7: Fuzzy Rules for predicting a words grapheme count

####  II-C4 Grapheme Count Prediction Results

The Fuzzy Inference System outputs a number to 2 decimal places in the range
[1-14]. This is calculated using the rules outlined above in Section II-C3
followed by a centroid defuzzification method [8]. The initial results are
listed below in Table II.

TABLE II: Fuzzy Inference System Results Predicted Result | Count | Percentage of Corpus
---|---|---
Correct Result | 4838 | 50.18%
One Greater Than Expected | 2577 | 26.73%
One Lower Than Expected | 1600 | 16.6%
Prediction Wrong By More Than 1 | 626 | 6.49%

###  II-D Grapheme Count Mapping

After the Fuzzy Inference system has predicted the number of graphemes in a
specified word, this predicted value is passed into a mapping method. This
method iterates through the list of graphemes, from largest to smallest to
find a matching grapheme. This then repeats until the entire word has been
split into graphemes. If at any point the mapping results in a failure to find
a valid grapheme, the iteration restarts until a valid combination of
graphemes has been selected.

##  III Results

To allow for comparison, a second method, IPA mapping, was developed. This
method involved using a pronunciation dictionary to return the IPA characters
for the word. This then followed a similar iterative system to find a matching
IPA Phoneme from the characters, which have all the same problems as grapheme
mapping, and from there, use known correspondences to select the valid
grapheme.

Using both methods, two different accuracy scores were measured. Firstly, did
the method select the correct number of graphemes, secondly, where the correct
number of graphemes were selected, were these the correct mapping for the
word. The results of these experiments, shown in Table III, show two
contrasting conclusions. Firstly, in its current state, IPA mapping results in
a greater proportion of words being successfully split into the correct number
of graphemes. However, the second conclusion drawn from these results is that
the Fuzzy Inference System, followed by mapping graphemes results in a greater
proportion of words being split into the correct, British received
pronunciation, mapping of graphemes.

TABLE III: Results Method | Training | Training | Validation | Validation
---|---|---|---|---
| Dataset | Accuracy | Dataset | Accuracy
IPA | Google | 69.27% | EnglishGrammarHere | 71.20%
Mapping | 10,000 | 36.41% | 10,000 | 38.94%
Fuzzy | Google | 50.18% | EnglishGrammarHere | 55.68%
Inference | 10,000 | 62.16% | 10,000 | 59.93%

##  IV Conclusion

This paper has presented a methods to split words into graphemes using a Fuzzy
Inference System This methods can be used for phonological word structure
analysis in NLP and in NLG applications. The experiments undertaken have shown
that the general trend of grapheme count can be predicted using a fuzzy
inference system. However, the fine tuning is dependant on the specific
dialect of the transcriber/reader and therefore a definitive count of N
graphemes in a word cannot be classified using a fuzzy inference system and
would require manual inference to classify precisely.

Comparing the two methods, one of the limitations of the IPA mapping method is
that it only works with linguistically correct spelling due to relying on the
use of a dictionary for converting a word into IPA. Where words are not in the
dictionary, or are spelled incorrectly due to a typo or lack of phonological
awareness due to another reason such as Specific Learning Disorders. The fuzzy
methods developed counter this issue by not being tied to a dictionary lookup
or other limiting resource. Instead they allow for breaking up words into
graphemes regardless of correctness in spelling.

## References

  * [1] Brooks, G. (2019) Dictionary of the British English Spelling System. Open Book Publishers.
  * [2] Coulmas, F. (1999) The Blackwell Encyclopedia of Writing Systems. Wiley.
  * [3] Department for Education (2014) National curriculum in England: English programmes of study. GOV.UK.
  * [4] Dubois, D. and Prade, H. (1996) What are fuzzy rules and how to use them. Fuzzy Sets and Systems, 84, 169-185.
  * [5] Elam, S. (2007) Phonics Primer. National Right to Read Foundation.
  * [6] EnglishGrammarHere (2020) 10000 Most Common Words in English.
  * [7] H, M.E. (1974) Application of fuzzy algorithms for control of simple dynamic plant. Proceedings of the Institution of Electrical Engineers, 121, 1585–1588.
  * [8] Hellendoorn, H. and Thomas, C. (1993) Defuzzification in Fuzzy Controllers. Journal of Intelligent & Fuzzy Systems, 1, 109–123.
  * [9] Oxford Dictionaries (2023) Oxford English Corpus - Facts About the Language. Oxford English Corpus - Facts About the Language.
  * [10] WorldlyWisdom (2021) Google 10000 English.
  * [11] Zadeh, L.A. (1983) A computational approach to fuzzy quantifiers in natural languages. Computers & Mathematics with Applications, 9, 149–184.