XXXXXX This is an open access article distributed under the terms of the
Creative Commons Attribution license
(http://creativecommons.org/licenses/by/3.0/), which permits unlimited use,
distribution and reproduction in any medium so long as the original work is
properly cited. XXXX

Nayak et. al.Big Computing: Where are we heading?

# Big Computing: Where are we heading?

Sabuzima Nayak  1  Ripon Patgiri1 and Thoudam Doren Singh  1  1
11affiliationmark:  National Institute of Technology Silchar, Assam-788010,
India<EMAIL_ADDRESS>
(XXXX; XXXX)

###### Abstract

This paper presents the overview of the current trends of Big data against the
computing scenario from different aspects. Some of the important aspect
includes the Exascale, the computing power and the kind of applications which
offer the Big data. This starts with the current computing hardware constraint
against the need of the rising Big data applications. We highlight the issues
and challenges of energy requirement, software complexity, hardware failure,
fault tolerant computing, and communication. As the complexity of computation
is going to rise in the future. The paper also highlights the future direction
of Big computing systems for Bioinformatics, social media, hardware and
software requirements, data intensive computation and then towards GPU era.

###### keywords:

Big Computation, Big Computing, Big Data, Exascale, Exaflops, Zettaflops, GPU.

††articletype: Research Article

[1]Corresponding author. Email:

##  1 Introduction

Big Computing is the future trend in two-fold, namely, data-intensive
computation (Big Data) and task-intensive computation. Big Data deals with
data-intensive computation while High Performance Computing (HPC) deals with
task-intensive computation. If the scientific world is unable to produce a
Exascale supercomputer (HPC) in the next few years, then the avalanche of data
produced at an Exascale speed cannot be processed. Currently, the world is
producing data of nearly 2.5 Exabytes per day as reported by the IBM [1].
Hence, the size of Big Data in the near future is going to reach Exascale. The
current Big data technologies are unable to handle such Exascale data. If the
world is not prepared for the future, then it will lead to shutting of data
processing and loss of information. Storage of data is not the aim of the
technology, but to process the stored data and utilize the mined knowledge.
Hence, the Exascale Big data require a Big computing system capable of
processing Exabyte sized data.

Figure 1: Big Computing and its enabling technologies

The future is about Exascale which is depicted in Figure 1. Hence, no super
power country wants to be left behind in the race of producing Exascale
computational system. Europe has launched EuroHPC initiative to develop a
supercomputer based on European technology by 2023 and to become an
independent superpower in the race of Exascale. It has provided a funding of 1
billion euro till 2020 [2]. Similarly, the United States is investing $1.8
billion USD for development of three Exascale Big computing systems till
2021-2022 [3]. In TOP500, United states Summit [4] supercomputer achieved the
first position on June 2018. Summit has 9216 CPU (Central processing unit)
chips and tightly connected 27648 GPUs (Graphics processing Unit). TOP500 [5]
is an event on advanced digital technology that ranks top 500 powerful
computers. For ranking the computer systems, LINPACK Benchmark [6] is followed
which measure the ability to solve dense linear equations problem set using
floating point arithmetic. For every country having its own supercomputer is
not just a competition, but answers to many unsolved and complex problems.
Summit is used to study about human diseases, fusion energy, advanced
materials and many more areas [7].

The article is organized as follows- Section 2 briefs on Big Data systems.
Section 3 explores on High Performance Computer. Section 4 presents rich
insight on issues and challenges to be faced by Big computing system in
future. Section 5 exposes future possibilities of ZettaFLOPS systems. Section
6 elaborates future aspects of Big computing systems. Finally, the article
draws a suitable conclusion in Section 7.

##  2 Big Data

The current trend of computing scenario has introduced data explosion in
different facets of life. Thus, Big Data is one of the most popular game
changer paradigm. It has no untouched area, particularly, government sector,
public sector, and private sector [8]. More specifically, Big Data is applied
in business, economy, healthcare, medicine, bioinformatics, earthquake,
weather, stock market, and many emerging areas. The growing revenue is amazing
in Big Data technology. Therefore, billion dollars are invested in Big Data
technology. Moreover, the data are to grow at an exponential pace. In 2009,
the Big Data started to grow at an exponential pace till 2019. However, the
data growth rate will not be exponential after 2020. The data growth rate will
face declination from 2025 to 2030. The data growth rate will be at a linear
pace rather than exponential pace. However, the growth rate of data is
expected to grow constantly at a linear pace in future. There is a requirement
for Big computing system to compute massive scale dataset. Big Computing
should redefine Big Data such that every company can be able to compute their
Exabytes of dataset at their own farmhouse in the near future. Almost, all IT
industry has to face Exabytes of data after 2025. Hence, Big Data 2.0 will be
face off in Big Computation. Big Data 2.0 is yet be conceived. Thus, Industry
will completely transform from Big Data to Big Data 2.0 completely after 2025.
The Big Data 2.0 has to deal with massive amount which includes exabytes and
beyond. Current state-of-the-art unable to deal with massive amount of data
which scale exabytes or beyond. It is expected that Big Data 2.0 will be
evolved in 2025. As data are growing at an exponential pace, the Internet of
Everything (IoE) will generate enormous amount of data as shown in Figure 1.
From 2020 to 2025, Beyong 5G (B5G) will evolve and hence, the Big Data 2.0 is
required to for processing massive amount of data.

##  3 High Performance Computing

Essentially, the knowledge or information extraction from the Big data lies in
the way how the high performance computation is carried out. ZettaFLOPS is the
next era for HPC which can perform $10^{21}$ floating point operations per
second (FLOPS). Current race is with ExaFLOPS in HPC which can perform
$10^{18}$ FLOPS. Many countries are participating to get a position in TOP500.
ZettaFLOPS is yet to be achieved. It is expected to achieve by around 2025.
Hence, the scientific and research world is focusing on creating a
supercomputer/High-Performance Computing (HPC) which can perform beyond
ExaFLOPS. However, there are many issues that need to be dealt before creating
a ZettaFLOPS system.

##  4 Issues and challenges

###  4.1 Energy requirement

A Big computing system approximately requires processing capability of 50
GigaFLOPS/Watt of power. Currently, we are unable to produce energy efficient
computation. Energy is consumed by processors, cooling fans, main board,
network infrastructure and many other components. Moreover, the voltage supply
is reduced to leverage the energy consumption, but it results in data
corruption [9]. Power consumed by datacenter, and storage infrastructure also
add to the total power consumption. Farmhouse bears more cost for energy
compared to establishing entire farmhouse. In other scenario, Data are oil of
the modern epoch. Therefore, data are being amalgamated which forms a data
silo. This data silo requires huge processing power to mine meaningful
insight. Data is growing and therefore, hardware is augmented to the
datacenter to increase scalability. Thus, datacenter requires more power
supplies. Therefore, energy-aware computation is the next emerging area.

###  4.2 Hardware failure

Many node failure occurs due to hardware failure. Currently, suppose a single
processor of a system takes 25 years for mean time to failure. Then, in a Big
computing system with hundred thousands of the processor has a mean time to
failure of 2 hours [10]. The throughput of each component is increased by
using transistor scaling process. Transistor scaling is embedding more
transistors in each computing device. But decreasing feature size leads to
reduction in critical charge. Critical charge is used to change the logic of a
circuit. So, the change results in data corruption. It is called soft errors
in Big computing systems. In addition, DRAM (dynamic RAM) errors lead to
undetected errors. Those errors results in incorrect answer calculation.

###  4.3 Fault tolerance

Big computing systems have millions of components with frequent failures. In
such cases, the network should provide efficient fault-tolerant mechanisms. In
software layer, the upper level layer should provide fault-tolerance
mechanisms to make robust software infrastructure. To devise such mechanism,
deep learning and streaming will play a big role.

###  4.4 Communication

Communication among the component requires energy. Hence, reduction in the
number of communication reduces power consumption. But, in the case of Big
computing system, the communication is between hundreds thousands of
components. Hence, reducing the number of communication is not an easy task.
Currently, zero-round trip time (0-RTT) is emerging [11]. However, in case of
occurrence of failure, it affects the communication among the components. So
fault-tolerant techniques among the components are the necessity. From another
viewpoint, network communication among millions of nodes in an Exascale
distributed system need to be efficient and fault-tolerant. Moreover,
reduction of latency is also dependant of efficient network communication.
Besides, 5G is emerging. 5G will be able to provide 20 GBPS downlink peak data
rate, and 10 GBPS uplink peak data rate. Furthermore, 5G will revolutionize
almost all major industries. Thus, Big Computing will become a more prominent
player in modern IT industries.

###  4.5 Software Complexity

In Big computing system, the complexity of system and application software
increases to many folds [12]. First, increasing computational efficiency. New
complex software is required to include finer temporal and spatial scales.
This needs heavy data assimilation and very complex physics. Second, adding
algorithmic resilience. It eliminates undetected errors. Third, energy-
awareness. It is the most important necessity because very high energy
consumption is making the existence of the Big computing system impossible.
Fourth, distributed system. Storing Exascale data is not possible by using a
single datacenter. Therefore, the programming model is required to re-
engineered for the multicore CPU (e.g. data parallel MPI (Massage Passing
Interface) [13]) and another for the accelerator (e.g. OpenCL [14]). Fifth,
memory levels. The hierarchy of memories increases the software complexity.

##  5 ZettaFLOPS

Achieving Exascale, Big computing system requires efficient hardware to handle
millions of computing per second. Performance of a Big computing system is
evaluated by the number of floating point operations per second (FLOPS). In
this section, some new approaches are discussed from the hardware aspect.

###  5.1 Processor

The processing unit is the main heart of the Big computing system. To design a
Big computing system, a high-performing processor with efficient architectural
design is essential. Currently, high-performing processor are implementing
multicore and SIMD (single instruction multiple data) techniques to reach the
performance target. Heterogeneous design is also considered. It combines a
general-purpose multicore processor and a many-core coprocessor. Multicore
processor performs irregular computation and many-core coprocessor performs
high throughput computation. Till date, the Exascale performing processor is
not achieved, but stitching the new processor technology [15] will make it
possible to achieve the Exascale standards. However, many issues and
challenges has to be resolved to achieve a Big Computing processor. Some
issues are (a) Exascale level performance, (b) power efficient, (c) satisfying
the need of new applications such as deep learning and Big data analytics.
Hence, it is important to realize that a processor has to balance performance
and energy for Big computing system.

ARM [16] processors are currently gaining popularity because power consumption
is one of the main issues in designing the Big computing systems. Advanced
RISC Machines (ARM) are embedded in many smart phones and tablets. ARM [16]
uses Restricted Instruction Set Computer (RISC) instructions. RISC instruction
sizes are fixed, has simple instruction fetch, and simultaneous access to the
opcode and operand. This makes the control unit design simple and consumes
less power. Moreover, RISC instructions execute in one clock cycle, which
reduces the interrupt latency. Cost of server-oriented ARM SoC (single on-
chip) is less compared to Intel x86 CPUs. In addition, ARM is easy to debug.
However, ARM has some issues [17] such as efficiency of ARM reduces in case of
strong scalable vector usage, and the architectural design is complex. On the
other hand, it is estimated that ARM is capable of replacing x86 and POWER-
based servers with some optimization in Big computing systems and hyperscale
data centers [18].

ExaNest [19] is a computing architecture designed for Exascale-class systems.
It is built upon 64 bit ARM processors. ARM processors [16] utilize 2-3 times
less electricity and performs more computation. ExaNest reduces power
consumption by decreasing:- data transfer distance, footprint area and cost of
installation, and the number of devices required to reach performance
ambition. The storage devices are kept close to compute nodes and connected
using fast custom-made interconnects. A single and unified interconnect is
implemented to manage storage and application traffic. A daughter board
contains the basic compute unit to enhance the performance of computation. It
has 4 Xilinx Ultrascale+ FPGAs, 16 cores (each Xilinx Ultrascale+ FPGA with
quad ARM Cortex-A53 64-bit), every FPGA is attached to 16 gigabytes of DDR4
memory, and an NVM in-node SSD for storage.

###  5.2 Storage memory

Many new memory technologies are proposed to replace DRAM. These technologies
should be able to solve many problems such as less power consumption,
achieving higher parallelism, run-time error handling and many more. In
addition, Big computing system will have a heterogeneous memory and storage
model. For example, along with DRAM non-volatile memory (NVM) can be used to
increase the efficiency of I/O operations [20]. Similarly, Phase change memory
(PCM) [21] uses phase-change material to store data. The phase-change material
is either amorphous or crystalline. It provides superior density compared to
DRAM. Moreover, the PCM can have different degree of crystallization in
different cells. This enables each cell to store more than one bit [22]. PCM
also scale better compared to DRAM [23]. But it has some demerits such as it
is slower compared to DRAM, has more memory latency, and performs less write
operations to reduce latency. Another technology is STT-MRAM (Spin-Transfer
Torque Magnetic Random Access Memory) array [24] based on Magnetic Tunneling
Junction (MTJ). An MTJ is a thin tunneling having a dielectric in between two
ferro-magnetic layers. STT-MRAM are compatible with the conventional DRAM
chips and can be used to construct byte-addressable memory devices. It is
capable of replacing DRAM. However, STT-MRAM require quick reads which causes
problems due to small sense margins. During manufacturing, thermal
fluctuations causes high write errors. In Big computing system, efficiency in
the main memory is mainly required. Hence, an upgraded and highly efficient
main memory will be able to handle and manage millions of operations per
seconds.

SAGE (percipient StorAGe for Exascale data centric computing) [25] is an
Exascale capable hardware. It has a multi-tiered I/O hierarchy with an
intelligent management software. SAGE has a software stack based on Unified
Object-Based Storage Infrastructure. The software stack consists of mainly
three layers, namely, Mero, Clovis and third layer consist of tools and HPC
APIs. Mero is the lowest layer in the software stack. It provides a platform
for distributed object storage. Mero comprise of a core that supports basic
object storage features. It also does resource management (locks, caches,
extents, etc.). In addition, Mero provides hardware reliability, for example,
RAID enabled through Server Network Striping. Clovis is the second layer which
is a transactional storage API. User applications directly use this layer. In
the third layer, SAGE uses some tools for I/O profiling and optimizing data
movement. Apache Flink tools enable data analytic jobs. SAGE takes the support
of POSIX compliant storage access. SAGE uses HSM to control data movement in
the SAGE hierarchies using the data usage information. Moreover, SAGE has
developed a new tool called RTHMS [26]. RTHMS analyses parallel applications
to recommend about data placement of memory objects in the heterogeneous
memory systems. It is designed to handle Exascale operations, however, many
new features can be added to increase the efficiency.

##  6 Future direction

The Big computing systems cannot reach its expected performance level by using
the existing methods and techniques. Making the Big computing system a
reality, new or hybrid approaches are required. In this section, new
approaches that are proposed for Exascale system is discussed.

###  6.1 Communications

As depicted in Figure 1, the enabler of Big Computing is Internet of
Everything (IoE). The IoE is backed by mobile communication. Currently, 5G
mobile communication is being deployed in many countries, and also, Beyond 5G
(B5G) is being developed. Moreover, it is expected that 6G mobile
communication will be deployed from 2030 and onward. The 6G will revolutionize
modern lifestyle, society, business and industry, and thus, we will evidence
the transition from Internet of Things (IoT) to IoE [27]. Moreover, there will
be many transitions, for instance, smart devices to intelligent devices. Also,
6G will be able to deliver many new technologies, for instance, holographic
communication, five sense communication, intelligent vehicles, etc., to
provide rich Quality of Experiences (QoE) [27]. Due to 6G technology, IoE will
create a massive small-sized data. These massive small-sized data will pose a
challenge to Big Computing in storing, processing and managing. Therefore, Big
Computing is to be redefined for massive scale small-sized data.

###  6.2 Bioinfromatics

With the advancement of DNA sequence extraction technology (next-generation
sequencing) and Big Data [28] resulted in decreasing the cost of whole genome
sequencing to a reasonable price [29]. However, the data obtained are
unordered small fragments of the sequence numbered to billions. These data
scale to terabytes. However, with the help of next-generation sequencing,
genomes of all non-human organism will be obtained in the future and the
memory requirement is in Exascale or beyond. The DNA collecting repository
currently has a size of petabyte [28] such as The Cancer Genome Atlas (TCGA)
and Human Microbiome Project. The size will increase and reach Exascale in
near future.

###  6.3 Social Media

Today, social media is another area that uses Big Data technology for storage,
processing and analyzing data. The user stores their messages, photos, music,
video, etc. in the storage provided by the social media service provider.
These data are mostly public so it opened the door for the businesses to
analyze these data and understand their customers. These data provide the
demography, cultural preferences, likes, dislikes, etc. These data are
analyzed to obtain the type of product the user prefers. Then, the businesses
do advertisement showcasing the preferred product to the specific user.
Moreover, it allows businesses to make decisions. Using this information,
sentiment of the user are understood to make it favorable to their product
selling. Social media has also become a common platform to advertise not only
products, but also talent, movies, music, messages. New movie reviews and
comments of the users can be analyzed to know the reaction of people and the
level of excitement which may lead to flop or hit of a movie. Similarly,
social media is used to showcase their talent such as dance, singing, mimicry
etc. performed by common people. Along with data the user using social media
is also increasing with a very high speed. In the near future the data
produced by the social media are also going to reach Exascale. The business
analytic has to be prepared for the Exascale Big data processing. It should be
intelligent enough to differentiate between the useful and not useful data.
Another important point is in Exascale era the social media service provider
has to protect their users from cyber bullies, cyber crime, etc. Protecting
clients from damaging, unethical and illegal data (messages, music, videos
etc.). Analyzing and finding such harmful content along with maintaining the
privacy and integrity of the user is going to be a difficult task for the
social media service provider.

The history is going to repeat itself. When the Big data era came the
traditional methods are unable to process the Big data. Again, when the
Exascale era will come the Big Data technologies will be unable to process the
data. In such situation, only storing data will be meaningless. Many business
fields are gaining profits by using the knowledge obtained after processing of
data. They will suffer huge losses. As explained earlier controlling the data
will also be difficult. In social media, when the Big data technologies will
not work, then their users data will be left vulnerable. Hence, Exascale Big
Data requires a Big computing system before the onset of the Exascale era.

###  6.4 Future Application

Big computing helps in solving many complex problems and simulations. Some
fields that have future applications of Big computing have advanced material
science, new energy solution, urban planning, cosmology, astrophysics and many
more [30]. Discovery of new material with required characteristics in the
field of advanced material science needs solution to many complex
calculations. Known compounds are combined to obtain the desired properties.
These calculations help in simulating the behavior of the obtained material in
the nature. Currently, classical simulation and deep learning are used for the
research. And, with Big Computing more complex design and calculation can be
performed. The challenge of finding an alternative energy sources can also be
solved using Big computing. For example, suppose in a wind farm the speed of
wind is very less. Then, for continuous production of electricity quickly the
fossil fuel generator is started. Predicting such events manually is
impossible. Such actions need heavy and Big computation. Big Computing helps
in understanding the chemical composition of the fossil fuel. This helps in
optimizing the energy production and reducing pollution. Big Computing also
will play an important role in increasing the quality of life. Big computing
processes Exascale data to optimize infrastructure access and usage choices.
Data is collected from several sources such as sensors, databases, survey done
by government. These data scale to Exabytes. Big computing analyze and
simulate thousands of scenarios. It obtains new planned city or restructuring
existing large cities. Another field of science that desperately need Big
computing system is cosmology [31] and astrophysics. This field of science is
interested in the evolution of the universe and creation, and destruction of
stars. Research in such field heavily depends on complex calculations and
simulations. Hence, requires the help of Big Computing to unfold the hidden
mysteries of the universe.

Initially the focus of the scientist was to obtain the whole genome to
visualize the genome structure. Recently scientist has realized that apart
from storage there is a requirement for processing of such Exascale data.
After a genome is assembled, computation is performed to analyze the
correlation of mutation with the disease or search for the history of
evolution of organisms. But such computation requires pairwise comparison
between genomes. This leads to computation to a quadratic scale. Hence, with
moving of Exascale Big Data the science world move towards Big Computing
capable of Exascale computing. The DNA with its complex network of chemical
composition is an example of Big Graph. Such Big Graph with higher dimension
requires Big Computing to estimate the parameter and optimize the metabolic
model. This metabolic model helps with discovering important genes that encode
proteins. It is useful in bio-manufacturing and other industries. Big
Computing will also show the path to personalized medicine. The effect of an
administered drug will be studied using Big Computing. Big Computing will help
with couple drug-induced perturbations of channel behavior and cellular action
potential model. Moreover, in case of a drug behaviour in patient’s heart. The
MRI images are used to obtain a 3D image of the heart. But such analysis needs
the support of Big Computing.

###  6.5 Hardware and Software Requirements

Big computing system at Exascale level requires both upgraded hardware and
software [32]. First such systems have to solve many issues such as
synchronization, load balancing, failure handling, and distributed memory
management. One main requisite of Big computing systems is an efficient design
for expressing task dependencies and inter-task parallelism. The system also
requires concurrent communication among a very large set of tasks. Moreover,
heterogeneous computation components (eg. CPU and GPU both integrated in a
single chip) are also proposed. However, a heterogeneous target computing
infrastructures arises more issues. Big computing system have to partition the
Exascale data to small dataset and perform parallel operations numbered to
thousands or millions. Exascale data processing also requires a new set of
languages that are capable of extensive parallelism. Many languages are
developed based on data-centric methods. Some examples are X10 [33], UPC [34],
and Legion [35]. Correspondingly, Big computing requires large memory,
hardware supporting large distributed system. MPI is capable to handle one
million cores. This makes MPI eligible to be included in Big computing
systems. Most important issue the Big computing systems will face is handling
failure. In future the number of compute nodes is going to increase. And,
failure of these nodes will be very common and within a very small duration,
mean time to failure of 2 hours [10]. To maintain reliability in Big computing
systems, multiple software and hardware techniques need to be deployed to
predict crashes and tries to maintain the stability of the system.
P2P-MapReduce [36] is an adaptive MapReduce framework. It manages Namenode
failures, Datanode churn and does decentralized job recovery. It helps to
provide a reliable MapReduce middleware in a dynamic large-scale distributed
system.

###  6.6 Data-intensive computation

Big computing systems will be more data intensive rather than compute-
intensive. Hence, more support is required from the storage systems [37]. The
storage system has to support intensive metadata operations. Optimized file
writing is required. Fault tolerance should also be handled by some new
techniques of checkpointing. The storage system has to flush the memory for
external storage. And, rolling back of memory for checkpoint on occasion of
failure should be rare [38]. Because, it will lead to increment in I/O cost.
Jin et al. [39] proposed a multi-tiered dynamic data staging that integrates
both DRAM and SSD. Solid state disk (SSD) is a non-volatile memory device
which is becoming popular. It is deployed as storage buffers to cache
temporally the checkpoint data [40] or caching data before writing to disks
[41]. The hybrid data staging is able to store large data volume, which
exceeds the capacity of currently used DRAM. It has other merits such as it
supports code coupling, data sharing and management, and has adopted an
application-aware data placement mechanism.

It is estimated that the file based data exchange method will become
infeasible to implement in the Big computing system [42]. The cost of I/O will
not permit to achieve the performance target required by the Big computing
system. Hence, the hardware has to be efficient to provide a solution to this
issue. Another constraint is data locality. The main reasons are high latency
and insufficient bandwidth. One possible solution is in-situ data analysis
[43]. It means moving the computation where the data resides. Two ways are
possible to implement in-situ data analysis. First, unified execution by
integration of the data processing and data production. The data analysis
component access the produced data by making function calls to their shared
memory address space. It is becoming popular because it helps to decrease the
I/O operation time in data-intensive systems. In addition, overly decreases
the end-to-end workflow runtime. Second, implementing data staging. In data
staging the workflow coupling is done by converting some compute nodes into
dedicated I/O nodes. Some applications can use the memory of these nodes for
storage and exchanging intermediate workflow data. In some framework this
mechanism is implemented such as DataStager [44] and PreDatA [45]. However, it
has many issues. API integrated data staging does recompiling after the
application compilation. Data staging need extra compute resources.

###  6.7 GPU era

Compared to CPU, GPU provides more peak performance and bandwidth using less
power. A Big computing system with GPU will be smaller compared to non-GPU
system. Therefore, GPU will be the building block of the Big computing system
[46]. A GPU based Big computing system will require less compute resources to
meet the performance target. Moreover, in case of more I/O bound operation the
GPU will also provide efficient memory bandwidth. However, in Big computing
system the overall consumption of power has to be constraint. Hence, more
power efficient GPU accelerators are essential. In GPU based Big computing
system, data communication between the host and the GPU is only through GPU
memory. In case of parallel application, GPU requires more data access [47].
It results in less performance and more power consumption. The power
consumption by the GPU is based on many factors. Moreover, designing energy
efficient scheduling for the GPU is difficult. Because the scheduler is
embedded in the device firmware [48]. Separate study is required to make the
GPU more efficient because knowledge used to study CPU cannot be used on GPU
due to differences in the architecture and organization [49].

Currently, the idea of combining both the CPU and GPU is proposed. The best
features of both are combined to further increase the overall performance.
This is termed as heterogeneous Computing (HC). CPUs are efficient in latency-
critical applications and GPUs are efficient in throughput-critical
applications. The aim of the HC unit will be load balancing between CPU and
GPU to keep both in idle state for very less duration. It is planned to
embedded both in the same chip. Many similar processors already exist, such as
gem5-gpu [50] and Ivy Bridge [51]. However, integrating both also have some
issues such as [52]: (a) HC can have either discrete or fused architecture,
(b) power consumption (c) data transfer and memory bandwidth overhead (d)
pipelining data transfer between CPU and GPU (e) efficient algorithm to
include high levels of parallelism (f) partition of workload to keep both for
very less duration in idle state. The HC is also going to get affected by the
issues they individually pose. In addition, the optimization techniques
implemented on CPU and GPU individually will not be applicable to HC.

In Big computing system, GPU is going to play a big role to meet the
performance peak. Similarly, the advancement made in making the CPU efficient
can be utilized to increase the performance. Both GPU-CPU integration also
provides another way in achieving an efficient Big computing system. However,
above discussed issues need to be solved before there are fully implemented.
Finally, they have to be very optimal in power consumption because overall
power consumption has to be less to achieve the Big Computing reality.

###  6.8 Energy-aware computation

Power consumption is the most important issue the Big computing system has to
solve. Because the amount of power the Big computing system consume is
threatening its feasibility [53]. The world is facing power shortage issue,
hence, providing megawatts of energy to a Big computing system is very
impractical. For example, when fully utilized China’s Tianhe-2 uses 17
MegaWatt. Currently, it is estimated that a Big computing system on an average
requires 20 megawatts of energy. A Big computing system has many power
consumption sources. Some examples are (a) cooling of the heat produced by the
Big computing system. Its contribution is very high. (b) data storage systems
(c) networking (d) data warehouse having the Big computing system. They are
also restricted to some specific location due to security and historical
reason. Achieving lower power consumption need optimization of every
components of the Big computing system. Using power efficient computing
components (eg. ARM). Data transfer also requires power. So, increasing the
locality of data and reducing compute components to reduce the traffic can be
considered. Reducing the speed automatically when less data need to be
transferred. Variable data transfer speed should be adopted to reduce energy
wherever possible. In different software components (application, scheduler,
and runtime) energy-conserving mechanisms should be adopted [20].

ECOSCALE [54] is a novel energy efficient heterogeneous hierarchical
architecture. It has a hybrid architecture having many-core+OpenCL programming
environment and a runtime component. It partition the physical system into
many compute nodes. These compute nodes are interconnected in a tree-like
fashion. ECOSCALE defines a contiguous global address space where nodes are
hierarchically interconnected by an MPI protocol. The compute nodes further
reduces power consumption and provides resilience by implementing a dual stage
Memory Management Unit. This unit maps the reconfigurable accelerators to
virtual address space.

Another power consuming factor is the cooling of the heat produced by the
system. In Big computing system using traditional methods (Air conditioner) is
not sufficient. Hence, new approaches are proposed to cool the Big computing
system. Liquid cooling is another solution which provides higher cooling
compared to traditional methods. Some examples are spray cooling [55],
microchannel [56] and jet impingement [57]. Currently, a hybrid approach is
also proposed which combines cooling liquid and traditional method. Galileo
[58] is a cooling approach which combines traditional method and rack-level
RDHXs. Rack-level RDHX is mounting liquid-to-air rear door heat exchangers
(RDHX) at rack level. Chen et al. [59] proposed a novel system that does
efficient cooling by reusing waste heat. The system combines a plug-type spray
cooling system and absorption chiller.

##  7 Conclusion

The world with a very high speed is moving towards the Exascale era. The near
future of Exascale is a horrible future, which is capable of impeding the
world with its avalanche of data. And, currently the world is nowhere near to
handle such an era. The first step in Exascale is producing a Big computing
machine. Before achieving the first step many preparations are required. Just
to reach the first step of supercomputer many issues and challenges need to be
overcome. Some issues and challenges are energy requirement, hardware failure,
fault tolerance, communication and software complexity. In this paper, we have
presented Big Computing which is two-fold, namely, data-intensive and task-
intensive computation. Therefore, Exascale data and task computing will pose
many research challenges and issues which are already discussed in this paper.
Energy-aware computing is next big things to do for research. Also, 5G will
revolutionize entire industry towards Big Computing. Big Data 2.0 is yet to
conceptualize for very large scale data. Similarly, ExaFLOP need to be evolved
to reach ZettaFLOP. Thus, we expect Big Computing will be the next big thing.

## References

  * [1] IBM, What will we make of this moment?2013 ibm annual report, Accessed on 16 April 2019 from https://www.ibm.com/annualreport/2013/bin/assets/2013_ibm_annual.pdf.
  * [2] Gagliardi, F., Moreto, M., Olivieri, M. and Valero, M. (2019) The international race towards exascale in europe. _CCF Transactions on High Performance Computing_ 10.1007/s42514-019-00002-y.
  * [3] Morgan, T.P., Bidders off and running after $1.8 billion doe exascale super deals, Accessed on 15 April 2019 from https://www.nextplatform.com/2018/04/09/bidders-off-and-running-after-1-8-billion-doe-exascale-supercomputer-deals/.
  * [4] Hines, J. (2018) Stepping up to summit. _Computing in Science & Engineering_ 20(2): 78–82.
  * [5] Top 500: The list, https://www.top500.org/.
  * [6] Dongarra, J.J., Luszczek, P. and Petitet, A. (2003) The linpack benchmark: past, present and future. _Concurrency and Computation: practice and experience_ 15(9): 803–820.
  * [7] Palmer, J. (2019) More super supercomputers. _Engineering_ https://doi.org/10.1016/j.eng.2019.04.007.
  * [8] Patgiri, R. and Ahmed, A. (2016) Big data: The v’s of the game changer paradigm. In _2016 IEEE 18th International Conference on High Performance Computing and Communications; IEEE 14th International Conference on Smart City; IEEE 2nd International Conference on Data Science and Systems (HPCC/SmartCity/DSS)_ (Sydney, Australia: IEEE): 17–24. 10.1109/HPCC-SmartCity-DSS.2016.0014.
  * [9] Martino, C.D., Kalbarczyk, Z. and et al., R.K.I. (2014) Lessons learned from the analysis of system failures at petascale: The case of blue waters. In _2014 44th Annual IEEE/IFIP International Conference on Dependable Systems and Networks_ : 610–621. 10.1109/DSN.2014.62.
  * [10] Bautista-Gomez, L., Zyulkyarov, F., Unsal, O. and McIntosh-Smith, S. (2016) Unprotected computing: A large-scale study of dram raw error rate on a supercomputer. In _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , SC ’16 (Piscataway, NJ, USA: IEEE Press): 55:1–55:11.
  * [11] Derler, D., Jager, T., Slamanig, D. and Striecks, C. (2018) Bloom filter encryption and applications to efficient forward-secret 0-rtt key exchange. In Nielsen, J.B. and Rijmen, V. [eds.] _Advances in Cryptology – EUROCRYPT 2018_ (Cham: Springer International Publishing): 425–455.
  * [12] Geist, A. and Reed, D.A. (2017) A survey of high-performance computing scaling challenges. _The International Journal of High Performance Computing Applications_ 31(1): 104–113.
  * [13] Gropp, W. (2009) Mpi at exascale: Challenges for data structures and algorithms. In _PVM/MPI_ : 3.
  * [14] Stone, J.E., Gohara, D. and Shi, G. (2010) Opencl: A parallel programming standard for heterogeneous computing systems. _Computing in science & engineering_ 12(3): 66.
  * [15] Shalf, J.M. and Leland, R. (2015) Computing beyond moore’s law. _Computer_ 48(SAND-2015-8039J).
  * [16] Aroca, R.V. and GonÃ§alves, L.M.G. (2012) Towards green data centers: A comparison of x86 and arm architectures power efficiency. _Journal of Parallel and Distributed Computing_ 72(12): 1770 – 1780\. https://doi.org/10.1016/j.jpdc.2012.08.005.
  * [17] Lobet, M., Haefele, M., Soni, V., Tamain, P., Derouillat, J. and Beck, A. (2018) High-performance computing at exascale: challenges and benefits. In _15ème congrés de la Société Française de Physique division Plasma_.
  * [18] Durand, Y., Carpenter, P.M. and et al., S.A. (2014) Euroserver: Energy efficient node for european micro-servers. In _2014 17th Euromicro Conference on Digital System Design_ : 206–213. 10.1109/DSD.2014.15.
  * [19] Katevenis, M., Chrysos, N., Marazakis, M. _et al._ (2016) The exanest project: Interconnects, storage, and packaging for exascale systems. In _2016 Euromicro Conference on Digital System Design (DSD)_ (IEEE): 60–67.
  * [20] Panda, D., Lu, X.Y. and Subramoni, H. (2018) Networking and communication challenges for post-exascale systems. _Frontiers of Information Technology & Electronic Engineering_ 19(10): 1230–1235.
  * [21] Tyson, S., Wicker, G., Lowrey, T., Hudgens, S. and Hunt, K. (2000) Nonvolatile, high density, high performance phase-change memory. In _2000 IEEE Aerospace Conference. Proceedings (Cat. No.00TH8484)_ , 5: 385–390 vol.5. 10.1109/AERO.2000.878512.
  * [22] Bedeschi, F., Fackenthal, R., Resta, C., Donze _et al._ (2008) A multi-level-cell bipolar-selected phase-change memory. In _2008 IEEE International Solid-State Circuits Conference-Digest of Technical Papers_ (IEEE): 428–625.
  * [23] Freitas, R.F. and Wilcke, W.W. (2008) Storage-class memory: The next storage system technology. _IBM Journal of Research and Development_ 52(4/5): 439.
  * [24] Khvalkovskiy, A., Apalkov, D., Watts, S., Chepulskii, R. _et al._ (2013) Basic principles of stt-mram cell operation in memory arrays. _Journal of Physics D: Applied Physics_ 46(7): 074001.
  * [25] Narasimhamurthy, S., Danilov, N. and et al., S.W. (2018) Sage: Percipient storage for exascale data centric computing. _Parallel Computing_ https://doi.org/10.1016/j.parco.2018.03.002.
  * [26] Peng, I.B., Gioiosa, R. and Kestor, G.e.a. (2017) Rthms: A tool for data placement on hybrid memory system. _SIGPLAN Not._ 52(9): 82–91. 10.1145/3156685.3092273.
  * [27] Nayak, S. and Patgiri, R. (2020) 6G: Envisioning the Key Issues and Challenges. _CoRR_ abs/2004.040244. URL https://arxiv.org/abs/2004.04024. 2004.04024.
  * [28] Stephens, Z.D., Lee, S.Y., Faghri, F., Campbell, R.H., Zhai, C., Efron, M.J., Iyer, R. _et al._ (2015) Big data: astronomical or genomical? _PLoS biology_ 13(7): e1002195.
  * [29] Lee, C.T. and Amaro, R.E. (2018) Exascale computing: A new dawn for computational biology. _Computing in Science & Engineering_ 20(5): 18–25.
  * [30] Messina, P. (2017) The exascale computing project. _Computing in Science & Engineering_ 19(3): 63–67.
  * [31] Goz, D., Tornatore, L., Taffoni, G. and Murante, G. (2017) Cosmological simulations in exascale era. _arXiv preprint arXiv:1712.00252_ .
  * [32] Petcu, D., Iuhasz, G. and Pop, D.e.a. (2015) On processing extreme data. _Scalable Computing: Practice and Experience_ 16. 10.12694/scpe.v16i4.1134.
  * [33] Tardieu, O., Herta, B. and Cunningham, D.e.a. (2014) X10 and apgas at petascale. _ACM SIGPLAN Notices_ 49(8): 53–66.
  * [34] Nishtala, R., Zheng, Y., Hargrove, P.H. and Yelick, K.A. (2011) Tuning collective communication for partitioned global address space programming models. _Parallel Computing_ 37(9): 576–591.
  * [35] Bauer, M., Treichler, S., Slaughter, E. and Aiken, A. (2012) Legion: Expressing locality and independence with logical regions. In _SC’12: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis_ (IEEE): 1–11.
  * [36] Marozzo, F., Talia, D. and Trunfio, P. (2012) P2p-mapreduce: Parallel data processing in dynamic cloud environments. _Journal of Computer and System Sciences_ 78(5): 1382–1402.
  * [37] Freeman, P.A., Crawford, D.L., Kim, S. and Munoz, J.L. (2005) Cyberinfrastructure for science and engineering: Promises and challenges. _Proceedings of the IEEE_ 93(3): 682–691. 10.1109/JPROC.2004.842782.
  * [38] Zhao, D., Liu, N., Kimpe, D., Ross, R., Sun, X. and Raicu, I. (2016) Towards exploring data-intensive scientific applications at extreme scales through systems and simulations. _IEEE Transactions on Parallel and Distributed Systems_ 27(6): 1824–1837. 10.1109/TPDS.2015.2456896.
  * [39] Jin, T., Zhang, F., Sun, Q., Bui, H., Romanus, M., Podhorszki, N., Klasky, S. _et al._ (2015) Exploring data staging across deep memory hierarchies for coupled data intensive simulation workflows. In _2015 IEEE International Parallel and Distributed Processing Symposium_ : 1033–1042. 10.1109/IPDPS.2015.50.
  * [40] Prabhakar, R., Vazhkudai, S.S. and et al., Y.K. (2011) Provisioning a multi-tiered data staging area for extreme-scale machines. In _2011 31st International Conference on Distributed Computing Systems_ : 1–12. 10.1109/ICDCS.2011.33.
  * [41] Li, M., Vazhkudai, S.S., Butt, A.R., Meng, F., Ma, X., Kim, Y., Engelmann, C. _et al._ (2010) Functional partitioning to optimize end-to-end performance on many-core architectures. In _SC ’10: Proceedings of the 2010 ACM/IEEE International Conference for High Performance Computing, Networking, Storage and Analysis_ : 1–12. 10.1109/SC.2010.28.
  * [42] Ma, K.L., Wang, C., Yu, H. and Tikhonova, A. (2007) In-situ processing and visualization for ultrascale simulations. _Journal of Physics: Conference Series_ 78: 012043. 10.1088/1742-6596/78/1/012043.
  * [43] Ovsyannikov, A., Romanus, M., Van Straalen, B., Weber, G.H. and Trebotich, D. (2016) Scientific workflows at datawarp-speed: accelerated data-intensive science using nersc’s burst buffer. In _2016 1st Joint International Workshop on Parallel Data Storage and data Intensive Scalable Computing Systems (PDSW-DISCS)_ (IEEE): 1–6.
  * [44] Abbasi, H., Wolf, M., Eisenhauer, G., Klasky, S., Schwan, K. and Zheng, F. (2010) Datastager: scalable data staging services for petascale applications. _Cluster Computing_ 13(3): 277–290. 10.1007/s10586-010-0135-6.
  * [45] Zheng, F., Abbasi, H., Docan, C., Lofstead, J., Liu, Q., Klasky, S., Parashar, M. _et al._ (2010) Predata – preparatory data analytics on peta-scale machines. In _2010 IEEE International Symposium on Parallel Distributed Processing (IPDPS)_ : 1–12. 10.1109/IPDPS.2010.5470454.
  * [46] Vuduc, R. and Czechowski, K. (2011) What gpu computing means for high-end systems. _IEEE Micro_ 31(4): 74–78. 10.1109/MM.2011.78.
  * [47] Che, S., Sheaffer, J.W., Boyer, M., Szafaryn, L.G., and Skadron, K. (2010) A characterization of the rodinia benchmark suite with comparison to contemporary cmp workloads. In _IEEE International Symposium on Workload Characterization (IISWC’10)_ : 1–11. 10.1109/IISWC.2010.5650274.
  * [48] Keckler, S.W., Dally, W.J., Khailany, B., Garland, M. and Glasco, D. (2011) Gpus and the future of parallel computing. _IEEE Micro_ 31(5): 7–17. 10.1109/MM.2011.89.
  * [49] Allen, T. and Ge, R. (2016) Characterizing power and performance of gpu memory access. In _Proceedings of the 4th International Workshop on Energy Efficient Supercomputing_ , E2SC ’16 (Piscataway, NJ, USA: IEEE Press): 46–53. 10.1109/E2SC.2016.8.
  * [50] Power, J., Hestness, J., Orr, M.S., Hill, M.D. and Wood, D.A. (2015) gem5-gpu: A heterogeneous cpu-gpu simulator. _IEEE Computer Architecture Letters_ 14(1): 34–36. 10.1109/LCA.2014.2299539.
  * [51] Damaraju, S., George, V., Jahagirdar, S., Khondker, T., Milstrey, R., Sarkar, S., Siers, S. _et al._ (2012) A 22nm ia multi-cpu and gpu system-on-chip. In _2012 IEEE International Solid-State Circuits Conference_ : 56–57. 10.1109/ISSCC.2012.6176876.
  * [52] Mittal, S. and Vetter, J.S. (2015) A survey of cpu-gpu heterogeneous computing techniques. _ACM Comput. Surv._ 47(4): 69:1–69:35. 10.1145/2788396.
  * [53] Abu Ahmad, W., Bartolini, A., Beneventi, F., Benini, L., Borghesi, A., Cicala, M., Forestieri, P. _et al._ (2017) Design of an energy aware petaflops class high performance cluster based on power architecture. In _2017 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW)_ : 964–973. 10.1109/IPDPSW.2017.22.
  * [54] Mavroidis, I., Papaefstathiou, I., Lavagno, L., Nikolopoulos, D.S., Koch, D., Goodacre, J., Sourdis, I. _et al._ (2016) Ecoscale: Reconfigurable computing and runtime system for future exascale systems. In _2016 Design, Automation & Test in Europe Conference & Exhibition (DATE)_ (IEEE): 696–701.
  * [55] Silk, E.A., Golliher, E.L. and Selvam, R.P. (2008) Spray cooling heat transfer: Technology overview and assessment of future challenges for micro-gravity application. _Energy Conversion and Management_ 49(3): 453 – 468. https://doi.org/10.1016/j.enconman.2007.07.046.
  * [56] Rezania, A. and Rosendahl, L. (2012) Thermal effect of a thermoelectric generator on parallel microchannel heat sink. _Energy_ 37(1): 220 – 227. https://doi.org/10.1016/j.energy.2011.11.043.
  * [57] Hosain, M.L., Fdhila, R.B. and Daneryd, A. (2016) Heat transfer by liquid jets impinging on a hot flat surface. _Applied Energy_ 164: 934 – 943. https://doi.org/10.1016/j.apenergy.2015.08.038.
  * [58] Conficoni, C., Bartolini, A., Tilli, A., Cavazzoni, C. and Benini, L. (2016) Integrated energy-aware management of supercomputer hybrid cooling systems. _Ieee transactions on industrial informatics_ 12(4): 1299–1311.
  * [59] Chen, H., Cheng, W.l., Zhang, W.w., Peng, Y.h. and Jiang, L.j. (2017) Energy saving evaluation of a novel energy system based on spray cooling for supercomputer center. _Energy_ 141: 304–315.


# A Statistical Analysis for the Neutrinoless Double-Beta Decay Matrix element
of 48Ca

M. Horoi  Department of Physics, Central Michigan University, Mount Pleasant,
MI 48859, USA  A. Neacsu  International Center for Advanced Training and
Research in Physics (CIFRA), Magurele, Romania  S. Stoica  International
Center for Advanced Training and Research in Physics (CIFRA), Magurele,
Romania

(August 29, 2024)

###### Abstract

Neutrinoless double beta decay ($0\nu\beta\beta$) nuclear matrix elements
(NME) are the object of many theoretical calculation methods, and are very
important for analysis and guidance of a large number of experimental efforts.
However, there are large discrepancies between the NME values provided by
different methods. In this paper we propose a statistical analysis of the 48Ca
$0\nu\beta\beta$ NME using the interacting shell model, emphasizing the range
of the NME probable values and its correlations with observables that can be
obtained from the existing nuclear data. Based on this statistical analysis
with three independent effective Hamiltonians we propose a common probability
distribution function for the $0\nu\beta\beta$ NME, which has a range of (0.45
- 0.95) at 90% confidence level of, and a mean value of 0.68.

##  I Introduction

The study of the double-beta decay (DBD) is currently a hot research topic
since it is viewed as one of the most promising approaches to clarify
important, as yet unknown, properties of neutrinos and to explore physics
beyond the Standard Model (SM) Avignone  _et al._ (2008); Vergados  _et al._
(2012). Two scenarios are possible for this process to occur: i) two-neutrino
double-beta ($2\nu\beta\beta$) transitions (with emission of two
electrons/positrons and two anti-neutrinos/neutrinos), which conserve the
lepton number and are allowed by the SM and ii) double-beta decay transitions
without emission of neutrinos ($0\nu\beta\beta$), which violate the lepton
number conservation and are only allowed by theories beyond SM (BSM).

Neutrinoless DBD was not experimentally detected so far, but its measurement
would provide us with important information about lepton number violating
(LNV) processes; neutrino properties (neutrino absolute mass scale and mass
hierarchy, neutrino nature as Dirac or Majorana fermion, number of neutrino
flavours); CP and Lorentz simmetries violation; constraining of different BSM
mechanisms that may contribute to this decay mode, etc. The most common
mechanism investigated is the light left-handed (LH) Majorana neutrinos
exchange between two nucleons, but once a LNV operator is introduced in the
Lagrangian, several other mechanisms are also allowed, such as: the exchange
of light and heavy neutrinos in left-right symmetric models, the exchange of
supersymmetric particles, DBD with the emission of Majorons, etc.

The DBD half-life equations can be expressed, in a good approximation, as a
product of some factors. The $2\nu\beta\beta$ half-life is a product of a
phase space factor (PSF), which depends on the atomic charge and energy
released in the decay, and a nuclear matrix element (NME) related to the
nuclear structure of the parent and daughter nuclei. The $0\nu\beta\beta$
half-life, besides the PSF and NME factors, also contains a LNV factor,
related to the particular BSM mechanism that may contribute to the decay. If
several mechanisms are considered, the inverse half-life can be written as a
sum of all the individual contributions and their interference terms Doi  _et
al._ (1985); Vergados  _et al._ (2012); Rodejohann (2012); Deppisch  _et al._
(2012); Horoi and Neacsu (2016a); Neacsu and Horoi (2016); Ahmed  _et al._
(2017). Using the experimental limits of the $0\nu\beta\beta$ decay half-lives
and the theoretical values of PSF and NME, one can constrain the LNV
parameters and the associated BSM scenarios, usually under the assumption that
only one mechanism contributes at one time.

There is currently significant progress in the DBD experiments (in terms of
the amount of source material, decreasing background and improvement in the
detection techniques), leading to the expectation that the next generation of
experiments will be able to cover the entire region of the neutrino inverted
mass hierarchy Engel and Menéndez (2017). Concurrently, the progress of the
theoretical methods now provides us with accurate PSF values for all the
double-beta decay modes and transitions. Kotila and Iachello (2012); Stoica
and Mirea (2013); Mirea  _et al._ (2015). Thus, at present, the uncertainty in
the DBD calculations remains mostly the NME evaluation.

There are several nuclear structure methods for the NME calculation, the most
used being: Shell model methods Caurier  _et al._ (1990, 1996, 2005); Horoi
_et al._ (2007); Horoi and Stoica (2010); Horoi (2013); Horoi and Brown
(2013); Sen’kov and Horoi (2014); Neacsu and Horoi (2015); Horoi and Neacsu
(2016b, 2018), pnQRPA methods Suhonen and Civitarese (1998); Simkovic  _et
al._ (1999); Stoica and Klapdor-Kleingrothaus (2001); Rodin  _et al._ (2006);
Kortelainen and Suhonen (2007); Faessler  _et al._ (2012); Simkovic  _et al._
(2013), IBA methods Barea and Iachello (2009); Barea  _et al._ (2013), Energy
Density Functional method Rodriguez and Martinez-Pinedo (2010), PHFB Rath  _et
al._ (2013), Coupled-Cluster method (CC) Novario  _et al._ (2021), in-medium
generator coordinate method (IM-GCM) Yao  _et al._ (2020), and valence-space
in-medium similarity renormalization group method (VS-IMSRG) Belley  _et al._
(2021). Each of these methods have their strengths and weakness, largely
discussed over time in the literature, and the current situation is that there
are still significant differences between NME values calculated with different
methods, and sometimes, even between NME values calculated with the same
methods (see for example the review Engel and Menéndez (2017)). For the
$2\nu\beta\beta$ decay the NME are products of two Gamow-Teller (GT)
transition amplitudes, and most of the nuclear methods overestimate them, in
comparison with experiment. This drawback is often treated by introducing a
quenching factor that multiplies the GT operator and reduces its strength.
This is equivalent to using a quenched axial vector constant, instead of its
bare value $g_{A}=1.27$.

For the $0\nu\beta\beta$ decay the NME calculation is more complicated, since
besides the GT transitions, other transitions may contribute as well. Also,
the NME values calculated by different methods may differ by factors of 3-4
for most relevant isotopes, and up to 7-8 in the case if 48Ca (see e.g. Fig. 5
of Ref. Engel and Menéndez (2017), and Refs. Yao  _et al._ (2020); Novario
_et al._ (2021)). Uncertainties in the NME values are further amplified when
predicting half-lives, since they enter at the power of two in the lifetime
formulas. In addition, there is no measured lifetime for this decay mode to
compared with, and these uncertainties in the NME computation reflect in the
interpretation of the DBD data and planning the performances of the DBD
experiments.

The shell model-based methods have some advantages such as the inclusion of
all correlations between nucleons around the Fermi surface, preserving all
symmetries of the nuclear many-body problem, and the use of widely tested
nucleon-nucleon (NN) interactions. For different mass regions of nuclei, one
uses several different effective NN effective Hamiltonians that are
appropriate for the corresponding model spaces. These effective Hamiltonians
are usually obtained by starting with a theoretical Bruekner G-Matrix
Hamiltonian that is further fine-tuned to describe the experimental energy
levels for a large number of nuclei that can be investigated in the
corresponding model spaces. These effective Hamiltonians are described by a
small number of single particle energies and a finite number of two-body
matrix elements. As a by-product, the wave functions produced by these
Hamiltonian can be used to describe and predict observables, such as the
electromagnetic transition probabilities, Gamow-Teller transitions
probabilities, nucleon occupation probabilities, spectroscopic factors, etc,
using relative simple changes of the transition operators in terms of
effective charges and quenching factors. These effective charges and quenching
factors are calibrating to the existing data. For $0\nu\beta\beta$ NME such
calibrations are not yet possible due to the lack of data. However, different
existing effective Hamiltonians for nuclei envolved in a given
$0\nu\beta\beta$ decay produce smaller ranges of the NME. In addition, some
recent ab-initio methods, such as IMSRG Yao  _et al._ (2020); Belley  _et al._
(2021), build on the modern advances in the shell model by providing ab-initio
derived effective Hamiltonians and effective transition operators, and they
can provide some guidance for calibrating the shell model $0\nu\beta\beta$
NME.

It would be thus interesting to study the robustness of the $0\nu\beta\beta$
NME to small changes of the parameters of different effective shell model
Hamiltonians and to examine how the NME changes are correlated with other
observables. In this work, we propose a statistical analysis of
$0\nu\beta\beta$ NME of 48Ca calculated with the interacting shell-model using
three independent effective Hamiltonians (FPD6, GXPF1A, KB3G), emphasizing the
range of the NME probable values and their correlations with several
observables that can be compared to existing nuclear data. Based on this
statistical analysis we propose a common probability distribution function for
the $0\nu\beta\beta$ NME. We apply our analysis to 48Ca, which is the lightest
DBD isotope and thus more accessible to ab-initio calculations. We only
consider in this work the standard light LH neutrino exchange mass mechanism,
which is most likely to contribute to the $0\nu\beta\beta$ decay process.

The paper is organized as follows. In section II the calculation methods of
the observables and the statistical model are presented. Then, in section III
we present the results and discussions on their relevance, and in section IV
we end with conclusions and outlook. Finally we included an Appendix with a
short presentation of the Gram-Charlier A series that we used in our
statistical model.

##  II The Statistical Model

| Exp. | Error | FPD6 | GXPF1A | KB3G
---|---|---|---|---|---
$0\nu\beta\beta$ NME | N/A | N/A | 0.79 | 0.559 | 0.693
$2\nu\beta\beta$ NME | 0.035 Barabash (2020) | 0.003 | 0.062 | 0.050 | 0.045
48Ca B(E2)$\uparrow$ | 0.008Pritychenko  _et al._ (2016) | 0.001 | 0.007 | 0.006 | 0.05
48Ca 2+ | 3.832 Chen (2022) | 0.15 | 3.658 | 3.735 | 4.238
48Ca 4+ | 4.503 Chen (2022) | 0.15 | 4.134 | 4.264 | 4.231
48Ca 6+ | 7.953 Chen (2022) | 0.15 | 7.396 | 7.705 | 7.831
48Ca Occ(Nf5) | 0.032∗ | 0.395# | 0.117 | 0.032 | 0.112
48Ca Occ(Nf7) | 7.892∗ | 0.395# | 7.693 | 7.892 | 7.795
48Ca Occ(Np1) | 0.009∗ | 0.395# | 0.029 | 0.009 | 0.024
48Ca Occ(Np3) | 0.067∗ | 0.395# | 0.161 | 0.067 | 0.070
48Ca $\rightarrow$ 48Sc GT | 1.09 Grewe  _et al._ (2007) | 0.28 | 1.01 | 1.226 | 0.051
48Ti B(E2)$\uparrow$ | 0.063Pritychenko  _et al._ (2016) | 0.003 | 0.064 | 0.052 | 0.052
48Ti 2+ | 0.984 Chen (2022) | 0.150 | 1.118 | 1.010 | 0.985
48Ti 4+ | 2.296 Chen (2022) | 0.150 | 2.492 | 2.168 | 2.214
48Ti 6+ | 3.333 Chen (2022) | 0.150 | 3.425 | 2.922 | 3.046
48Ti Occ(Nf5) | 0.168∗ | 0.277# | 0.310 | 0.168 | 0.263
48Ti Occ(Nf7) | 5.535∗ | 0.277# | 5.253 | 5.535 | 5.416
48Ti Occ(Np1) | 0.048∗ | 0.277# | 0.068 | 0.048 | 0.061
48Ti Occ(Np3) | 0.248∗ | 0.277# | 0.369 | 0.248 | 0.260
48Ti Occ(Pf5) | 0.032∗ | 0.092# | 0.101 | 0.032 | 0.097
48Ti Occ(Pf7) | 1.839∗ | 0.092# | 1.672 | 1.839 | 1.763
48Ti Occ(Pp1) | 0.010∗ | 0.092# | 0.031 | 0.010 | 0.021
48Ti Occ(Pp3) | 0.119∗ | 0.092# | 0.196 | 0.119 | 0.120
48Ti $\rightarrow$ 48Sc GT | 0.014Grewe  _et al._ (2007) | 0.005 | 0.050 | 0.032 | 0.056
Table 1: Experimental data, experimental errors, and the calculated values
using 3 effective Hamiltonians for the observables analyzed. The data for
occupation probabilities (∗) is not available and it was replaced with the
GXPF1A results and errors (#) of 5% of the highest nucleon species occupation.

We plan to investigate the effect of small, random variation of the shell
model effective Hamiltonian on the neutrinoless double beta decay NME of 48Ca,
and the NME correlations with other calculated observables, such as $2^{+}$
energies, B(E2)$\uparrow$ values, $2\nu\beta\beta$ matrix elements, Gamow-
Teller transition probabilities, neutron and proton occupation probabilities,
etc.

To achieve that goal we selected a number of often used effective Hamiltonians
describing nuclei around 48Ca in the $fp$-shell ($0f_{7/2}$, $0f_{5/2}$,
$1p_{3/2}$, and $1p_{1/2}$ orbitals for both protons and neutrons), and added
small random contributions to their two-body matrix elements (TBME). For this
project we only considered the FPD6 Hamiltonian Richter  _et al._ (1991), the
KB3G Hamiltonian Caurier  _et al._ (2005), and GXPF1A Hamiltonian Honma  _et
al._ (2004, 2005) as starting effective Hamiltonians. In order to maintain the
magicity of 48Ca we decided to keep the single particle (s.p.) energies in the
perturbed effective Hamiltonians the same as in the starting Hamiltonians.

One important decision to be made about the random contributions to the
starting Hamiltonians is the choice of their maximum amplitude (range). In
this work we were guided by the analysis of the USDA/USDB effective
Hamiltonians Brown and Richter (2006) where one starts with an underlying
G-matrix and modifies linear combinations of two-body matrix elements in a
fine-tuning procedure Brown and Richter (2006) until the root mean square
(RMS) deviation of the calculated energies vs the experimental ones shows some
signs of convergence. In this fine-tuning process one would not want to change
the TBME too much from the original G-matrix values, because the over-fitted
TBME could result in unitary changes of the s.p. wave functions that may
produce slightly better energies, but incorrect observables. For USDA, for
example, the RMS deviation of the TBME was about 300 keV, while a small
improvement in the overall energies given by USDB resulted in an additional
change of 100 keV if the RMS deviation of the TBME. This analysis suggests
that an additional RMS of about 100 keV would not dramatically change the
quality of the TBME in the sd-shell and we extended this choice to the fp-
shell. An analysis of the TBME for all three starting Hamiltonians listed
above indicates that a $\pm 10\%$ range for the random contributions would
suffice.

In the analysis we included as observables the $0\nu\beta\beta$ NME, the
$2\nu\beta\beta$ NME, the Gamow-Teller probability to reach the first $1^{+}$
state in 48Sc from the ground states (g.s.) of the parent, 48Ca, and of the
daughter, 48Ti, the energies of the $2^{+}$, $4^{+}$ and $6^{+}$ states of the
parent and daughter, the B(E2)$\uparrow$ transition probabilities to the first
$2^{+}$ state of parent and daughter, the neutron occupation probabilities of
the $pf$ states of the parent, and neutron and proton occupation probabilities
of the $pf$ states of the daughter nucleus. The experimental occupation
probabilities for the nuclei relevant for the 48Ca $0\nu\beta\beta$ decay are
not available, but we include synthetic (calculated) values in the analysis
because the corresponding occupation probabilities are available for other
nuclei of interest for $0\nu\beta\beta$ decay Schiffer  _et al._ (2008); Kay
_et al._ (2009, 2013), and it might be interesting to see if they have any
correlations with the $0\nu\beta\beta$ NME. All in all, there are 24
observables included in our statistical analysis.

The main goals are: (i) for each starting effective Hamiltonian find
correlations between $0\nu\beta\beta$ NME and the other observables that are
accessible experimentally; (ii) find theoretical ranges for each observables;
(iii) establish the shape of different distributions for each observables and
starting Hamiltonians; (iv) use this information to find weights of
contributions from different starting Hamiltonians to the ”optimal”
distribution of the $0\nu\beta\beta$ NME; (v) find an ”optimal” value of the
$0\nu\beta\beta$ NME and its predicted probable range (theoretical error). One
should mention that similar studies for other observables were recently
proposed Fox  _et al._ (2020).

The $0\nu\beta\beta$ NME is related to the half-life of the respective process
Horoi and Stoica (2010) by

$\left(T^{0\nu}_{1/2}\right)^{-1}=G_{0\nu}(E_{0},Z)g_{A}^{4}\mid
M_{0\nu}\mid^{2}\mid<\eta_{l}>\mid^{2}\ ,$ (1)

where $G_{0\nu}$ and $M_{0\nu}$ are the PSF and nuclear matrix elements for
the $0\nu$ decay, $g_{A}$ is the axial vector coupling constant, and
$<\eta_{l}>\equiv<m_{\beta\beta}>/m_{e}c^{2}$ is a BSM parameter associated
with the light neutrino exchange mechanisms. Here we only consider the
contribution from the light LH neutrino exchange mechanism, which is likely to
contribute to the $0\nu\beta\beta$ decay. The methodology of calculating the
$0\nu\beta\beta$ NME, $M_{0\nu}$, within the shell model was extensively
described elsewhere Horoi and Stoica (2010); Horoi (2013); Horoi and Neacsu
(2018) and it will not be repeated here. Suffices to say that it includes a
short range correlation function that can be viewed as an effective
modification of the bare operator (see below).

The $2\nu\beta\beta$ NME is related to the half-life of the respective process
Horoi  _et al._ (2007) by

$\left(T^{2\nu}_{1/2}\right)^{-1}=G_{2\nu}(E_{0},Z)g_{A}^{4}\mid
m_{e}c^{2}M_{2\nu}\mid^{2}\ .$ (2)

Here, $G_{2\nu}$ is the appropriate PSF, and $M_{2\nu}$ can be calculated with

$M_{2\nu}=\sum_{k}\frac{q^{2}\left<0^{+}_{f}\mid\sigma\tau^{-}\mid
1^{+}_{k}\right>\left<1^{+}_{k}\mid\sigma\tau^{-}\mid
0^{+}_{i}\right>}{E_{k}-E_{0}}\ ,$ (3)

where the summation is on the $1^{+}_{k}$ states in 48Sc and
$E_{0}=Q_{\beta\beta}/2+\Delta M(^{48}Sc-^{48}Ca)$

Often, the shell model calculations of the $0\nu\beta\beta$ NME are described
as using the ”bare” transition operator. This characterization is unfortunate,
while the transition operator (see e.g. Eq. (7-12) of Ref. Horoi and Stoica
(2010)) contains the bare operator from the underlying theory of
$0\nu\beta\beta$ decay, modified by a phenomenological effective short-range
correlation function, $1+f(r)$, which is quenching the $0\nu\beta\beta$ NME.
Therefore, the short-range modification of the bare operator acts as an
effective operator. In practice the parameters of an effective operator need
to be calibrated to the data. Given that the short-range correlator has radial
dependence, its calibration has been only done relative to some ab-initio
results. The standard Miller-Spencer short-range correlator Miller and Spencer
(1976); Horoi and Stoica (2010) produces the highest quenching of the NME,
while the CD-Bonn parameterization of the short-range correlator Šimkovic  _et
al._ (2009) produces little to no quenching. A direct renormalization of the
$0\nu\beta\beta$ NME by a similarity renormalization group (SRG) evolution of
the NME of the bare operator from 200 to 10 major harmonic oscillator shells
using CDBonn two body wave functions, indicates that using a phenomenological
CDBonn parametrization of the short-range correlator is a reasonable approach
Horoi (2016). More recent ab-initio calculations of the $0\nu\beta\beta$ NME
using the $N^{3}LO$ Hamiltonian provides more quenched values, more consistent
with the shell model results base on Miller-Spencer parametrization of the
short-range correlator. In an effort to calibrate the effective operator used
in shell model calculations to the latest ab-initio results we used the
Miller-Spencer correlator in this study.

The other observables used in this study including the excited state energies,
the GT strengths to the first $1^{+}$ state 48Sc, the B(E2)$\uparrow$ to the
first $2^{+}$ state in the parent and daughter, as well as the s.p. occupation
probabilities, are calculated in the standard way. Here we use in all cases
the same effective charges ($e_{p}=1.5$ and $e_{n}=0.5)$) for the
B(E2)$\uparrow$, and the same quenching factor ($q=0.74$) for the the GT
strengths and $M_{2\nu}$.

| $0\nu\beta\beta$ | $2\nu\beta\beta$ |  48Ca |  48Ca |  48Ca |  48Ca |  48Ca Occ |  48Ca Occ |  48Ca Occ |  48Ca Occ |  48Ca $\rightarrow$ |  48Ti
---|---|---|---|---|---|---|---|---|---|---|---|---
| NME | NME | B(E2)$\uparrow$ | 2+ | 4+ | 6+ | (Nf5) | (Nf7) | (Np1) | (Np3) |  48Sc GT | B(E2)$\uparrow$
$0\nu\beta\beta$ NME | 1.00 | 0.90 | 0.43 | 0.22 | 0.24 | 0.18 | 0.62 | -0.58 | 0.30 | 0.41 | 0.12 | -0.06
$2\nu\beta\beta$ NME | 0.90 | 1.00 | 0.38 | 0.28 | 0.30 | 0.17 | 0.44 | -0.47 | 0.25 | 0.37 | 0.10 | -0.18
48Ca B(E2)$\uparrow$ | 0.43 | 0.38 | 1.00 | -0.35 | -0.27 | -0.31 | 0.55 | -0.91 | 0.53 | 0.92 | -0.17 | 0.32
48Ca 2+ | 0.23 | 0.28 | -0.35 | 1.00 | 0.96 | 0.49 | -0.06 | 0.32 | -0.23 | -0.42 | 0.51 | -0.38
48Ca 4+ | 0.24 | 0.30 | -0.27 | 0.96 | 1.00 | 0.44 | -0.01 | 0.26 | -0.17 | -0.36 | 0.54 | -0.33
48Ca 6+ | 0.18 | 0.17 | -0.31 | 0.49 | 0.44 | 1.00 | -0.22 | 0.25 | -0.19 | -0.20 | 0.26 | -0.22
48Ca Occ(Nf5) | 0.62 | 0.44 | 0.55 | -0.06 | -0.01 | -0.22 | 1.00 | -0.77 | 0.38 | 0.42 | -0.04 | 0.12
48Ca Occ(Nf7) | -0.58 | -0.47 | -0.91 | 0.32 | 0.26 | 0.25 | -0.77 | 1.00 | -0.58 | -0.90 | 0.17 | -0.30
48Ca Occ(Np1) | 0.30 | 0.25 | 0.53 | -0.23 | -0.17 | -0.19 | 0.38 | -0.58 | 1.00 | 0.48 | -0.11 | 0.19
48Ca Occ(Np3) | 0.41 | 0.37 | 0.92 | -0.42 | -0.36 | -0.20 | 0.42 | -0.90 | 0.48 | 1.00 | -0.22 | 0.34
48Ca $\rightarrow$ 48Sc GT | 0.12 | 0.10 | -0.17 | 0.51 | 0.54 | 0.26 | -0.04 | 0.17 | -0.11 | -0.22 | 1.00 | -0.27
48Ti B(E2)$\uparrow$ | -0.06 | -0.18 | 0.32 | -0.38 | -0.34 | -0.22 | 0.12 | -0.30 | 0.19 | 0.34 | -0.27 | 1.00
48Ti 2+ | 0.80 | 0.89 | 0.37 | 0.32 | 0.32 | 0.20 | 0.47 | -0.47 | 0.24 | 0.35 | 0.23 | -0.24
48Ti 4+ | 0.78 | 0.85 | 0.35 | 0.35 | 0.35 | 0.25 | 0.45 | -0.45 | 0.23 | 0.34 | 0.22 | -0.19
48Ti 6+ | 0.75 | 0.80 | 0.42 | 0.28 | 0.28 | 0.24 | 0.48 | -0.52 | 0.26 | 0.42 | 0.03 | 0.13
48Ti Occ(Nf5) | 0.14 | -0.03 | 0.39 | -0.16 | -0.11 | -0.36 | 0.75 | -0.53 | 0.25 | 0.24 | -0.00 | 0.27
48Ti Occ(Nf7) | 0.12 | 0.32 | -0.51 | 0.47 | 0.42 | 0.37 | -0.37 | 0.53 | -0.33 | -0.49 | 0.20 | -0.82
48Ti Occ(Np1) | -0.13 | -0.29 | 0.27 | -0.34 | -0.29 | -0.24 | 0.10 | -0.26 | 0.56 | 0.25 | -0.02 | 0.65
48Ti Occ(Np3) | -0.25 | -0.38 | 0.43 | -0.51 | -0.47 | -0.26 | 0.03 | -0.37 | 0.21 | 0.50 | -0.28 | 0.88
48Ti Occ(Pf5) | 0.22 | 0.01 | 0.20 | 0.00 | 0.03 | -0.08 | 0.56 | -0.35 | 0.15 | 0.12 | 0.21 | 0.27
48Ti Occ(Pf7) | 0.21 | 0.41 | -0.16 | 0.24 | 0.20 | 0.16 | -0.10 | 0.16 | -0.11 | -0.16 | -0.01 | -0.85
48Ti Occ(Pp1) | -0.11 | -0.29 | 0.13 | -0.19 | -0.16 | -0.14 | 0.06 | -0.13 | 0.30 | 0.12 | 0.16 | 0.60
48Ti Occ(Pp3) | -0.33 | -0.49 | 0.11 | -0.28 | -0.25 | -0.16 | -0.09 | -0.06 | 0.05 | 0.14 | -0.09 | 0.89
48Ti $\rightarrow$ 48Sc GT | 0.15 | 0.40 | 0.23 | -0.03 | -0.02 | -0.24 | 0.15 | -0.20 | 0.13 | 0.19 | -0.55 | 0.03
Table 2: Correlation matrix for the 24 observables described in the text using the FPD6 starting Hamiltonian (continues in Table III).  |  48Ti |  48Ti |  48Ti |  48Ti Occ |  48Ti Occ |  48Ti Occ |  48Ti Occ |  48Ti Occ |  48Ti Occ |  48Ti Occ |  48Ti Occ |  48Ti $\rightarrow$
---|---|---|---|---|---|---|---|---|---|---|---|---
| 2+ | 4+ | 6+ | (Nf5) | (Nf7) | (Np1) | (Np3) | (Pf5) | (Pf7) | (Pp1) | (Pp3) |  48Sc GT
$0\nu\beta\beta$ NME | 0.80 | 0.78 | 0.75 | 0.14 | 0.12 | -0.13 | -0.25 | 0.22 | 0.21 | -0.11 | -0.33 | 0.15
$2\nu\beta\beta$ NME | 0.89 | 0.85 | 0.80 | -0.03 | 0.32 | -0.29 | -0.38 | 0.01 | 0.41 | -0.28 | -0.48 | 0.39
48Ca B(E2)$\uparrow$ | 0.37 | 0.35 | 0.42 | 0.39 | -0.51 | 0.27 | 0.43 | 0.19 | -0.16 | 0.13 | 0.11 | 0.23
48Ca 2+ | 0.32 | 0.35 | 0.27 | -0.16 | 0.47 | -0.34 | -0.51 | 0.00 | 0.24 | -0.19 | -0.28 | -0.03
48Ca 4+ | 0.32 | 0.35 | 0.28 | -0.11 | 0.42 | -0.29 | -0.47 | 0.03 | 0.20 | -0.16 | -0.25 | -0.02
48Ca 6+ | 0.20 | 0.25 | 0.24 | -0.36 | 0.37 | -0.24 | -0.26 | -0.07 | 0.16 | -0.14 | -0.16 | -0.24
48Ca Occ(Nf5) | 0.47 | 0.45 | 0.48 | 0.75 | -0.37 | 0.10 | 0.03 | 0.56 | -0.10 | 0.06 | -0.09 | 0.15
48Ca Occ(Nf7) | -0.47 | -0.45 | -0.52 | -0.53 | 0.53 | -0.26 | -0.37 | -0.35 | 0.16 | -0.13 | -0.06 | -0.20
48Ca Occ(Np1) | 0.24 | 0.23 | 0.26 | 0.25 | -0.33 | 0.56 | 0.21 | 0.15 | -0.11 | 0.30 | 0.05 | 0.13
48Ca Occ(Np3) | 0.35 | 0.34 | 0.42 | 0.24 | -0.49 | 0.25 | 0.50 | 0.12 | -0.16 | 0.12 | 0.14 | 0.19
48Ca $\rightarrow$ 48Sc GT | 0.23 | 0.22 | 0.03 | 0.00 | 0.20 | -0.02 | -0.28 | 0.21 | -0.01 | 0.16 | -0.09 | -0.55
48Ti B(E2)$\uparrow$ | -0.24 | -0.19 | 0.13 | 0.27 | -0.82 | 0.65 | 0.88 | 0.27 | -0.85 | 0.60 | 0.89 | 0.03
48Ti 2+ | 1.00 | 0.96 | 0.85 | 0.16 | 0.25 | -0.24 | -0.42 | 0.20 | 0.36 | -0.19 | -0.51 | 0.28
48Ti 4+ | 0.96 | 1.00 | 0.89 | 0.16 | 0.21 | -0.21 | -0.36 | 0.22 | 0.29 | -0.15 | -0.44 | 0.24
48Ti 6+ | 0.85 | 0.89 | 1.00 | 0.20 | -0.01 | -0.09 | -0.09 | 0.22 | 0.08 | -0.07 | -0.18 | 0.35
48Ti Occ(Nf5) | 0.16 | 0.16 | 0.20 | 1.00 | -0.64 | 0.34 | 0.22 | 0.81 | -0.42 | 0.35 | 0.18 | -0.05
48Ti Occ(Nf7) | 0.25 | 0.21 | -0.01 | -0.64 | 1.00 | -0.74 | -0.88 | -0.55 | 0.86 | -0.66 | -0.80 | 0.12
48Ti Occ(Np1) | -0.24 | -0.21 | -0.09 | 0.34 | -0.74 | 1.00 | 0.66 | 0.41 | -0.75 | 0.88 | 0.67 | -0.22
48Ti Occ(Np3) | -0.42 | -0.36 | -0.09 | 0.22 | -0.88 | 0.66 | 1.00 | 0.19 | -0.82 | 0.56 | 0.90 | -0.10
48Ti Occ(Pf5) | 0.20 | 0.22 | 0.22 | 0.81 | -0.55 | 0.41 | 0.19 | 1.00 | -0.57 | 0.56 | 0.26 | -0.30
48Ti Occ(Pf7) | 0.36 | 0.29 | 0.08 | -0.42 | 0.86 | -0.75 | -0.82 | -0.57 | 1.00 | -0.81 | -0.94 | 0.31
48Ti Occ(Pp1) | -0.19 | -0.15 | -0.07 | 0.35 | -0.66 | 0.88 | 0.56 | 0.56 | -0.81 | 1.00 | 0.67 | -0.35
48Ti Occ(Pp3) | -0.51 | -0.44 | -0.18 | 0.18 | -0.80 | 0.67 | 0.90 | 0.26 | -0.94 | 0.67 | 1.00 | -0.23
48Ti $\rightarrow$ 48Sc GT | 0.28 | 0.24 | 0.35 | -0.05 | 0.12 | -0.22 | -0.10 | -0.30 | 0.31 | -0.35 | -0.23 | 1.00
Table 3: Correlation matrix: continuation of Table II.

##  III Results

The experimental data used in this study listed in Table I and in the legends
of the rightmost column in Tables 4/5 \- 8/9 are taken from Ref. Chen (2022)
(excitation energies of the $2^{+}$, $4^{+}$ and $6^{+}$ states of 48Ca and
48Ti in MeV), Ref. Barabash (2020) ($2\nu\beta\beta$ NME in MeV-1), Ref.
Pritychenko  _et al._ (2016) ( B(E2)$\uparrow$ in $e^{2}b^{2}$), and Ref.
Grewe  _et al._ (2007) for the GT transition probabilities to the first
excited $1^{+}$ state in 48Sc. The experimental errors for the excitation
energies are very small, and for the calculation the $\chi^{2}$ we use the
typical theoretical RMS value of 150 keV Honma  _et al._ (2004). The
experimental occupation probabilities are not available, and we took as
reference the GXPF1A results assuming an uniform error that we choose to be 5%
of the highest occupation for a each nucleon species in the $fp$-shell. In the
Tables Occ (Nf7) designates the neutron occupation probability of the
$f_{7/2}$ s.p. orbital, Occ (Pf3) designates the proton occupation probability
of the $p_{3/2}$ s.p. orbital, etc.

Tables 4/5 \- 8/9 show the main results of this study. The leftmost columns
indicates the 24 observables discussed in section II, including the
$0\nu\beta\beta$ NME. The middle column shows the scatter plots of the
correlation of each variable with the $0\nu\beta\beta$ NME, and the last
column shows the distribution of each observable when the random term is added
to the respective effective Hamiltonian. The legends in column with
correlations show the standard Pearson correlator $R$, and in the last columns
the legends include the mean, standard deviations, and the skewness
(normalized 3rd moment) of the distributions, as well as the result for the
starting interactions (FPD6, GXPF1A, and KB3G) and the experimental values
when available. Tables 4/5 present the results for the FPD6 effective
Hamiltonian, Tables 6/7 show the results for the GXPF1A effective Hamiltonian,
and Tables 8/9 present the results for the KB3G effective Hamiltonian. For
each starting effective Hamiltonian we use 20,000 random Hamiltonians produced
by the procedure described in section II.

The results in Tables 4/5 \- 8/9 indicate strong correlations between the
$0\nu\beta\beta$ NME and the $2\nu\beta\beta$ NME. Alternative approaches of
obtaining these NME, e.g. QRPA calculations, are calibrating parts of their
nuclear Hamiltonian, such as the isoscalar particle-particle interaction
$g_{pp}$, to describe the experimental value of the $2\nu\beta\beta$ NME and
to approximately restore the isospin symmetry, thus inducing correlation with
the $2\nu\beta\beta$ NME. In the shell model approach, the Hamiltonian remains
unchanged, and all symmetries are enforced. Therefore, we conclude that the
strong correlations between the $0\nu\beta\beta$ NME and the $2\nu\beta\beta$
NME are genuine.

Interestingly, the correlations between the $0\nu\beta\beta$ NME and the
Gamow-Teller (GT) transitions probabilities to the first $1^{+}$ state in 48Sc
are much reduced. One explanation of this phenomena is based on the fact the
distributions of the GT strength from the parent and daughter (see last column
in Tables 4/5 \- 8/9) are asymmetric in opposite direction, thus diminishing
the correlation effects. A quick look to the full correlation matrix in Table
3 shows that the GT strengths to the first $1^{+}$ state in 48Sc from 48Ca and
48Ti are anti-correlated with a correlation coefficient of about -0.5.

Other observables that have relatively high (anti)correlations with the
$0\nu\beta\beta$ NME are the energies of the $2^{+}$, $4^{+}$ and $6^{+}$
states in 48Ti, and the neutron occupation probabilities in 48Ca. Overall, the
correlators $R$ with the $2\nu\beta\beta$ NME are around 0.9, the ones with
the energies of the $2^{+}$, $4^{+}$ and $6^{+}$ states in 48Ti are about
0.77, and correlators with the $0f_{5/2}$ occupation probability is about 0.6,
while the occupation probability of the $0f_{7/2}$ is anti-correlated with the
$0\nu\beta\beta$ NME, $R\approx-0.6$.

Figure 1: PDF of the $0\nu\beta\beta$ NME distributions for the FPD6, GXPF1A
and KB3G Hamiltonians and their weighted sum (see text for details).

Additional interesting information can be extracted from the full correlation
matrix for all 24 observables. Tables 2 \- III show the full correlation
matrix evaluated for the FPD6 starting Hamiltonian. It is interesting to
analyze which other observables are correlated with those that are directly
correlated with the $0\nu\beta\beta$ NME. We already discussed the
correlations between the GT strengths and the $2\nu\beta\beta$ NME. In
addition, one can observe that the $2^{+}$, $4^{+}$ and $6^{+}$ in 48Ti are
correlated with the neutron $f$ and $p$ states occupancies in 48Ca, which in
their turn are correlated with some of the neutron states occupancies in 48Ti.
Also, some of the neutron occupation probabilities in 48Ti are correlated to
the B(E2)$\uparrow$ values. These observations highlight the importance of a
reliable experimental investigation of the occupation probabilities for these
nuclei.

It would be interesting to extract some information about possible range and
mean value of the $0\nu\beta\beta$ NME based on this statistical analysis.
First, it is clear that the value of all observables are quite stable to
reasonably small changes of the effective Hamiltonian. No hints of any wild
departure from the main values that would indicate some phase transitions are
found. This seems to be a consequence of the preservation of nuclear many-body
symmetries in the shell model. One can further try using the distributions of
all available effective Hamiltonians to draw conclusions on some optimal
values for the $0\nu\beta\beta$ NME and its range (error). One direct approach
would be to superpose the distributions of the NME produced in Tables 4/5 \-
8/9 withe some weighting factors $W_{H}$,

$\begin{split}P(x)&=W_{FPD6}P_{FPD6}(x)+W_{GXPF1A}P_{GXPF1A}(x)\\\
&+W_{KB3G}P_{KB3G}(x)\ ,\end{split}$ (4)

where $x$ is the value of the $0\nu\beta\beta$ NME. The normalized weights
$W_{H}$ can be inferred using, for example, the likelyhood probability
$\propto exp(-\chi^{2}/2)$, or replacing the bare $\chi^{2}$ with its
individual contributions weighted by the corresponding correlators $R$. Based
on the data we show in Table 1, we get the following $\chi^{2}$ for each
starting effective Hamiltonian: 7.9 for FPD6, 4.8 for GXPF1A, and 7.3 for
KB3G. Unfortunately, there is no experimental data for the occupation
probabilities that seem to correlate directly and indirectly with the
$0\nu\beta\beta$ NME. Therefore, here we present the results of a ”democratic”
approach in which all $W_{H}$ are 0.33. Fig. 1 shows the probability
distribution functions (PDF) for the three starting effective Hamiltonians and
their weighted sum. To calculate each PDF we use the Gram-Charlier A series
expansion Cramer (1957) (see Appendix for detail), based on the first four
normalized moments of the distributions presented in Tables 4/5 \- 8/9. Based
on the results of our statistical analysis summarized in Fig. 1 (see ”weighted
sum” curve) one can infer that with 90% confidence the $0\nu\beta\beta$ NME
lies in the range between 0.45 and 0.95, with a mean value of about 0.68.

##  IV Conclusion and Outlook

In conclusion, we developed a statistical model for analyzing the distribution
of the $0\nu\beta\beta$ NME of 48Ca using the interactive shell model in the
$fp$-shell model space. In the analysis we started from three widely used
effective Hamiltonians for the low part of the $fp$-shell, FPD6, GXPF1A and
KB3G, to which we added a random contributions to the TBME of $\pm$10%. Using
sample sizes of 20,000 points we analyzed for each starting effective
Hamiltonian: (i) the correlations between $0\nu\beta\beta$ NME and the other
observables that are accessible experimentally; (ii) the theoretical ranges
for each observables; (iii) the shape of different distributions for each
observables and starting Hamiltonians; (iv) the weighted contributions from
different starting Hamiltonians to the ”optimal” distribution of the
$0\nu\beta\beta$ NME; (v) an ”optimal” value of the $0\nu\beta\beta$ NME and
its predicted probable range (theoretical error).

We found that the $0\nu\beta\beta$ NME correlates strongly with the
$2\nu\beta\beta$ NME, but much less with the Gamow-Teller strengths to the
first $1^{+}$ state in 48Sc. We also found that the $0\nu\beta\beta$ NME
exhibits reasonably strong correlations with the energies of the $2^{+}$,
$4^{+}$ and $6^{+}$ states in 48Ti, and with the neutron occupation
probabilities in 48Ca. We also found that there are additional correlations
between observables, such as the energies of the $2^{+}$, $4^{+}$ and $6^{+}$
states in 48Ti and the neutron occupation probabilities, as well as between
B(E2)$\uparrow$ values in 48Ti and proton and neutron occupation
probabilities, which can indirectly influence the $0\nu\beta\beta$ NME.
Therefore, we conclude that reliable experimental values of the occupation
probabilities in 48Ti and 48Ca would be useful for this analysis, potentially
helpful to reduce the uncertainties of the $0\nu\beta\beta$ NME.

Based on this statistical analysis with three independent effective
Hamiltonians we propose a common probability distribution function for the
$0\nu\beta\beta$ NME, which has a range (theoretical error) of (0.45 - 0.95)
at 90% confidence level, and a mean value of 0.68. We also hope that the
present analysis will help ab-initio studies, such as Yao  _et al._ (2020);
Novario  _et al._ (2021); Belley  _et al._ (2021) to better identify
correlations and further reduce the uncertainties of the $0\nu\beta\beta$ NME

Observable | Correlation | PDF
---|---|---
$0\nu\beta\beta$ NME |  |
$2\nu\beta\beta$ NME |  |
48Ca $\rightarrow$ 48Sc GT |  |
48Ti $\rightarrow$ 48Sc GT |  |
48Ca B(E2)$(\uparrow)$ |  |
48Ca $2^{+}$ |  |
48Ca $4^{+}$ |  |
48Ca $6^{+}$ |  |
48Ca Occ(Nf5) |  |
48Ca Occ(Nf7) |  |
48Ca Occ(Np1) |  |
48Ca Occ(Np3) |  |
Table 4: (Color online) Correlations scattered plots and PDFs for the FPD6 starting Hamiltonian.  Observable | Correlation | PDF
---|---|---
48Ti B(E2)$(\uparrow)$ |  |
48Ti $2^{+}$ |  |
48Ti $4^{+}$ |  |
48Ti $6^{+}$ |  |
48Ti Occ(Nf5) |  |
48Ti Occ(Nf7) |  |
48Ti Occ(Np1) |  |
48Ti Occ(Np3) |  |
48Ti Occ(Pf5) |  |
48Ti Occ(Pf7) |  |
48Ti Occ(Pp1) |  |
48Ti Occ(Pp3) |  |
Table 5: (Color online) Continuation of Table 4.  Observable | Correlation | PDF
---|---|---
$0\nu\beta\beta$ NME |  |
$2\nu\beta\beta$ NME |  |
48Ca $\rightarrow$ 48Sc GT |  |
48Ti $\rightarrow$ 48Sc GT |  |
48Ca B(E2)$(\uparrow)$ |  |
48Ca $2^{+}$ |  |
48Ca $4^{+}$ |  |
48Ca $6^{+}$ |  |
48Ca Occ(Nf5) |  |
48Ca Occ(Nf7) |  |
48Ca Occ(Np1) |  |
48Ca Occ(Np3) |  |
Table 6: (Color online) Correlations scattered plots and PDFs for the GXPF1A starting Hamiltonian.  Observable | Correlation | PDF
---|---|---
48Ti B(E2)$(\uparrow)$ |  |
48Ti $2^{+}$ |  |
48Ti $4^{+}$ |  |
48Ti $6^{+}$ |  |
48Ti Occ(Nf5) |  |
48Ti Occ(Nf7) |  |
48Ti Occ(Np1) |  |
48Ti Occ(Np3) |  |
48Ti Occ(Pf5) |  |
48Ti Occ(Pf7) |  |
48Ti Occ(Pp1) |  |
48Ti Occ(Pp3) |  |
Table 7: (Color online) Continuation of Table 6.  Observable | Correlation | PDF
---|---|---
$0\nu\beta\beta$ NME |  |
$2\nu\beta\beta$ NME |  |
48Ca $\rightarrow$ 48Sc GT |  |
48Ti $\rightarrow$ 48Sc GT |  |
48Ca B(E2)$(\uparrow)$ |  |
48Ca $2^{+}$ |  |
48Ca $4^{+}$ |  |
48Ca $6^{+}$ |  |
48Ca Occ(Nf5) |  |
48Ca Occ(Nf7) |  |
48Ca Occ(Np1) |  |
48Ca Occ(Np3) |  |
Table 8: (Color online) Correlations scattered plots and PDFs for the KB3G starting Hamiltonian. Observable | Correlation | PDF
---|---|---
48Ti B(E2)$(\uparrow)$ |  |
48Ti $2^{+}$ |  |
48Ti $4^{+}$ |  |
48Ti $6^{+}$ |  |
48Ti Occ(Nf5) |  |
48Ti Occ(Nf7) |  |
48Ti Occ(Np1) |  |
48Ti Occ(Np3) |  |
48Ti Occ(Pf5) |  |
48Ti Occ(Pf7) |  |
48Ti Occ(Pp1) |  |
48Ti Occ(Pp3) |  |
Table 9: (Color online) Continuation of Table 8.

Acknowledgements. MH acknowledges support from the US Department of Energy
grant DE-SC0022538 ”Nuclear Astrophysics and Fundamental Symmetries”. SS and
AN acknowledge support by a grant of Romanian Ministry of Research, Innovation
and Digitalization through the project CNCS – UEFISCDI number 99/2021 within
PN-III-P4-ID-PCE-2020-2374.

*

##  Appendix A Gram-Charlier A series

In order to have a good representation for the PDF of the $0\nu\beta\beta$
NME, we consider small deviation from normal distribution via the Gram-
Charlier A series Cramer (1957). This is given by:

$\begin{split}&P(x)\approx\frac{1}{\sqrt{2\pi}\sigma}exp\left[-\frac{\left(x-\mu\right)^{2}}{2\sigma^{2}}\right]\\\
\cdot&\left[1+\frac{\mu_{3}}{3!}He_{3}((x-\mu)/\sigma)+\frac{\mu_{4}-3}{4!}He_{4}((x-\mu)/\sigma)\right]\end{split}$
(5)

where $He_{k}(y)$ are the Chebyshev-Hermite polynomials, $He_{3}(y)=y^{3}-3y$
and $He_{3}(y)=y^{4}-6y^{2}+3$, $\mu$ and $\sigma$ are the mean and variance
of a probability density function (PDF), and $\mu_{k}$, with $k=3,4$ are
normalized moments of the same PDF, $P(x)$:

$\mu_{k}=\int\left(\frac{x-\mu}{\sigma}\right)^{k}P(x)dx,$ (6)

In practice, we use the sample moments, $\mu_{3}$ (skewnes) and $\mu_{4}-3$
(kurtosis), which in the limit of very large sample sizes become very close to
the underlying moments.

## References

  * Avignone  _et al._ (2008) F. T. Avignone, III, S. R. Elliott, and J. Engel, Rev. Mod. Phys. 80, 481 (2008).
  * Vergados  _et al._ (2012) J. D. Vergados, H. Ejiri, and F. Simkovic, Rep. Prog. Phys. 75, 106301 (2012).
  * Doi  _et al._ (1985) M. Doi, T. Kotani, and E. Takasugi, Prog. Theor. Phys. Suppl. 83, 1 (1985).
  * Rodejohann (2012) W. Rodejohann, J. Phys. G 39, 124008 (2012).
  * Deppisch  _et al._ (2012) F. F. Deppisch, M. Hirsch, and H. Pas, J. Phys. G 39, 124007 (2012).
  * Horoi and Neacsu (2016a) M. Horoi and A. Neacsu, Phys. Rev. D 93, 113014 (2016a), arXiv:arXiv:1511.00670 [hep-ph] [hep-ph] .
  * Neacsu and Horoi (2016) A. Neacsu and M. Horoi, Advances in High Energy Physics 2016 (2016).
  * Ahmed  _et al._ (2017) F. Ahmed, A. Neacsu, and M. Horoi, Physics Letters B 769, 299 (2017).
  * Engel and Menéndez (2017) J. Engel and J. Menéndez, Reports on Progress in Physics 80, 046301 (2017).
  * Kotila and Iachello (2012) J. Kotila and F. Iachello, Phys. Rev. C 85, 034316 (2012).
  * Stoica and Mirea (2013) S. Stoica and M. Mirea, Phys. Rev. C 88, 037303 (2013).
  * Mirea  _et al._ (2015) M. Mirea, T. Pahomi, and S. Stoica, Rom. Rep. Phys. 67, 872 (2015).
  * Caurier  _et al._ (1990) E. Caurier, A. Poves, and A. P. Zuker, Phys. Lett. B 252, 13 (1990).
  * Caurier  _et al._ (1996) E. Caurier, F. Nowacki, A. Poves, and J. Retamosa, Phys. Rev. Lett. 77, 1954 (1996).
  * Caurier  _et al._ (2005) E. Caurier, G. Martinez-Pinedo, F. Nowacki, A. Poves, and A. P. Zuker, Rev. Mod. Phys. 77, 427 (2005).
  * Horoi  _et al._ (2007) M. Horoi, S. Stoica, and B. A. Brown, Phys. Rev. C 75, 034303 (2007).
  * Horoi and Stoica (2010) M. Horoi and S. Stoica, Phys. Rev. C 81, 024321 (2010).
  * Horoi (2013) M. Horoi, Phys. Rev. C 87, 014320 (2013).
  * Horoi and Brown (2013) M. Horoi and B. A. Brown, Phys. Rev. Lett. 110, 222502 (2013).
  * Sen’kov and Horoi (2014) R. A. Sen’kov and M. Horoi, Phys. Rev. C 90, 051301(R) (2014).
  * Neacsu and Horoi (2015) A. Neacsu and M. Horoi, Phys. Rev. C 91, 024309 (2015).
  * Horoi and Neacsu (2016b) M. Horoi and A. Neacsu, Phys. Rev. C 93, 024308 (2016b).
  * Horoi and Neacsu (2018) M. Horoi and A. Neacsu, Phys. Rev. C 98, 035502 (2018).
  * Suhonen and Civitarese (1998) J. Suhonen and O. Civitarese, Phys. Rep. 300, 123 (1998).
  * Simkovic  _et al._ (1999) F. Simkovic, G. Pantis, J. D. Vergados, and A. Faessler, Phys. Rev. C 60, 055502 (1999).
  * Stoica and Klapdor-Kleingrothaus (2001) S. Stoica and H. Klapdor-Kleingrothaus, Nucl. Phys. A 694, 269 (2001).
  * Rodin  _et al._ (2006) V. Rodin, A. Faessler, F. Simkovic, and P. Vogel, Nucl. Phys. A 766, 107 (2006).
  * Kortelainen and Suhonen (2007) M. Kortelainen and J. Suhonen, Phys. Rev. C 75, 051303 (2007).
  * Faessler  _et al._ (2012) A. Faessler, V. Rodin, and F. Simkovic, J. Phys. G 39, 124006 (2012).
  * Simkovic  _et al._ (2013) F. Simkovic, V. Rodin, A. Faessler, and P. Vogel, Phys. Rev. C 87, 045501 (2013).
  * Barea and Iachello (2009) J. Barea and F. Iachello, Phys. Rev. C 79, 044301 (2009).
  * Barea  _et al._ (2013) J. Barea, J. Kotila, and F. Iachello, Phys. Rev. C 87, 014315 (2013).
  * Rodriguez and Martinez-Pinedo (2010) T. R. Rodriguez and G. Martinez-Pinedo, Phys. Rev. Lett. 105, 252503 (2010).
  * Rath  _et al._ (2013) P. K. Rath, R. Chandra, K. Chaturvedi, P. Lohani, P. K. Raina, and J. G. Hirsch, Phys. Rev. C 88, 064322 (2013).
  * Novario  _et al._ (2021) S. Novario, P. Gysbers, J. Engel, G. Hagen, G. R. Jansen, T. D. Morris, P. Navrátil, T. Papenbrock, and S. Quaglioni, Phys. Rev. Lett. 126, 182502 (2021).
  * Yao  _et al._ (2020) J. M. Yao, B. Bally, J. Engel, R. Wirth, T. R. Rodríguez, and H. Hergert, Phys. Rev. Lett. 124, 232501 (2020).
  * Belley  _et al._ (2021) A. Belley, C. G. Payne, S. R. Stroberg, T. Miyagi, and J. D. Holt, Phys. Rev. Lett. 126, 042502 (2021).
  * Barabash (2020) A. Barabash, UNIVERSE 6, 159 (2020).
  * Pritychenko  _et al._ (2016) B. Pritychenko, M. Birch, B. Singh, and M. Horoi, ATOMIC DATA AND NUCLEAR DATA TABLES 107, 1 (2016).
  * Chen (2022) J. Chen, Nuclear Data Sheets 179, 1 (2022).
  * Grewe  _et al._ (2007) E.-W. Grewe, D. Frekers, S. Rakers, T. Adachi, C. Bäumer, N. T. Botha, H. Dohmann, H. Fujita, Y. Fujita, K. Hatanaka, K. Nakanishi, A. Negret, R. Neveling, L. Popescu, Y. Sakemi, Y. Shimbara, Y. Shimizu, F. D. Smit, Y. Tameshige, A. Tamii, J. Thies, P. v. Brentano, M. Yosoi, and R. G. T. Zegers, Phys. Rev. C 76, 054307 (2007).
  * Richter  _et al._ (1991) A. Richter, M. G. van der Merwe, R. E. Julies, and B. A. Brown, Nucl. Phys. A 523, 325 (1991).
  * Honma  _et al._ (2004) M. Honma, T. Otsuka, B. A. Brown, and T. Mizusaki, Phys. Rev. C 69, 034335 (2004).
  * Honma  _et al._ (2005) M. Honma, T. Otsuka, B. A. Brown, and T. Mizusaki, Eur. Phys. J. A 25 Suppl. 1, 499 (2005).
  * Brown and Richter (2006) B. A. Brown and W. A. Richter, Phys. Rev. C 74, 034315 (2006).
  * Schiffer  _et al._ (2008) J. P. Schiffer, S. J. Freeman, J. A. Clark, C. Deibel, C. R. Fitzpatrick, S. Gros, A. Heinz, D. Hirata, C. L. Jiang, B. P. Kay, _et al._ , Phys. Rev. Lett. 100, 112501 (2008).
  * Kay  _et al._ (2009) B. P. Kay, J. P. Schiffer, S. J. Freeman, T. Adachi, J. A. Clark, C. M. Deibel, H. Fujita, Y. Fujita, P. Grabmayr, K. Hatanaka, _et al._ , Phys. Rev. C 79, 021301 (2009).
  * Kay  _et al._ (2013) B. P. Kay, T. Bloxham, S. A. McAllister, J. A. Clark, C. M. Deibel, S. J. Freedman, S. J. Freeman, K. Han, A. M. Howard, A. J. Mitchell, _et al._ , Phys. Rev. C 87, 011302 (2013).
  * Fox  _et al._ (2020) J. M. R. Fox, C. W. Johnson, and R. N. Perez, Phys. Rev. C 101, 054308 (2020).
  * Miller and Spencer (1976) G. Miller and J. E. Spencer, Ann. Phys. (NY) 100, 562 (1976).
  * Šimkovic  _et al._ (2009) F. Šimkovic, A. Faessler, H. Müther, V. Rodin, and M. Stauf, Phys. Rev. C 79, 055501 (2009).
  * Horoi (2016) M. Horoi, Bull. Amer. Phys. Soc. 61, 123 (2016).
  * Cramer (1957) H. Cramer,  _Mathematical Methods of Statistics_ (Princeton University Press, Princeton, 1957).