11institutetext: Department of Computer Science, University of New Mexico
11email<EMAIL_ADDRESS>22institutetext: Center for Advanced Research
Computing, University of New Mexico
22email<EMAIL_ADDRESS>
# Beatnik: A Novel Global Communication Mini-Application

Jason R. Stewart 11 Patrick G. Bridges 22 0000-0003-4801-0390

(June 2024)

###### Abstract

Beatnik is a novel open source mini-application that exercises the complex
communication patterns often found in production codes but rarely found in
benchmarks or mini-applications. It simulates 3D Raleigh-Taylor instabilities
based on Pandya and Shkoller’s Z-Model formulation using the Cabana
performance portability framework. This paper presents both the high-level
design and important implementation details about Beatnik along with four
benchmark setups for evaluating different aspects of HPC communication system
performance. Evaluation results demonstrate both Beatnik’s scalability on
modern accelerator-based systems using weak and strong scaling tests up to
1024 GPUs, along with Beatnik’s ability to expose communication challenges in
modern systems and solver libraries.

###### Keywords:

Communication Patterns Benchmarking Scalability

##  1 Introduction

Mini-applications that reproduce and exercise key features of important
numerical methods and production applications have become an essential part of
the toolbox for designing and evaluating high-performance computing (HPC)
systems. Existing mini-applications with significant communication components,
however, generally exercise only simple communication patterns and primitives,
most frequently either:

  * •

Neighbor exchanges for regular meshes or particle distributions [8], irregular
or adaptive meshes [6, 7, 12], or sparse matricies [15]; or

  * •

Simple global collectives such as `MPI_Allreduce` for computing dot products
[4], `MPI_Allgather` for neighbor identification [15, 12], or `MPI_Alltoall`
for matrix transposition.

Few if any HPC mini-applications, however, exercise the varied, complex
communication patterns common in production applications. This is, we
hypothesize, because both the numerical methods that require these complex
communication patterns are inherently complex and beyond the scope of what is
normally considered appropriate for a mini-application. The result, however,
is that the performance of important communication patterns (e.g. those
resulting the from remeshing, particle/mesh migration in ALE codes, or tree
sweeps in fast multipole methods) are generally not considered when optimizing
MPI communication operations.

This paper presents Beatnik, a novel open source [2] mini-application designed
to exercise the complex global communication patterns found in real production
codes on modern GPU-accelerated systems. Beatnik does so in a modestly-sized
code base (less than 5000 source lines of C++ code) by leveraging recent
numerical method innovations and modern performance portability libraries [14,
3, 1, 9, 13]. Specifically, this paper describes the following contributions:

  * •

The overall goals and architecture of Beatnik, a mini-application that
simulates Raleigh-Taylor instabilities in 3D based on Pandya and Shkoller’s
Z-Model formulation [10]. The chosen numerical method and benchmark
architecture support multiple solution approaches that can exercise a wide
range of important global communication patterns not found in many other mini-
applications (Section 2);

  * •

The initial Beatnik software design and implementation built using the Cabana
performance portability framework and associated libraries, including high-,
medium-, and low-order solver implementations that leverage large-scale Fast
Fourier Transforms (FFTs) and multi-dimensional data redistributions (Section
3);

  * •

The identification of Beatnik test cases that exercise different elements of
MPI communication such as network bandwidth, network latency, and load
imbalance (Section 4); and,

  * •

The demonstration of the scalability of the Beatnik implementation on these
test cases and its ability to exercise communication performance on modern
machines up to 1024 GPUs (Section 5).

Following the presentation of these contributions, the paper discusses
directions for future work (Section 6), and concludes (Section 7).

##  2 Benchmark Goals and Architecture

The primary goal of Beatnik is to provide a communication benchmark that
exercises numerical solvers with irregular and/or global communication
patterns while still being relatively straightforward to understand and
implement. As a first step, we first sought to identify a numerical method
that is scalable, straightforward to understand, straightforward to implement
on modern architectures, and can leverage a range of numerical algorithms
which require global communication.

We settled on Pandya and Shkoller’s so-called “Z-Model” [10] for simulating
fluid interface instabilities as a numerical method that met these
requirements. In particular, the Z-Model is an ideal method for our global
communication benchmarking purposes because it includes low-, medium-, and
high-order methods that exercise a range of important numerical algorithms
with global communication requirements. Specifically, the low-order solver
uses Fast Fourier Transforms (FFTs)—and by extension all-to-all
communication—to approximate the Birkhoff-Rott velocity integral. In contrast,
the high-order model requires solving the Birkhoff-Rott integral directly
using any of a number of far-field force solution methods, for example the
fast multipole method (FFM) or particle–particle –particle–mesh (P3M), both of
which require comnplex global communication. Finally, the medium-order solver
_couples_ the FFT solver and the far-field force solver, resulting multiple
complex coupled communication patterns and data redistributions.

The overall architecture of Beatnik consists of the following high-level
modules:

SurfaceMesh:

This module implements the distributed 2D mesh that represents the interface
between two fluids; each point on this mesh stores the x/y/z coordinate of the
point and two vorticity components. The current Z-Model formulation is based
on an open regular rectangular mesh and this formulation is the basis of
Beatnik; an unstructured Z-model formulation that could more easily support
closed surfaces is planned for development but has not yet been finalized.

BRSolver:

This module calculates the Birkhoff-Rott velocity integral for points on the
SurfaceMesh from the vorticity and position of every point on the mesh. This
requires using an appropriate far-field force calculation method, either
internally implemented or provided by an external solver. This module, used
only by high- and medium- order Beatnik solves, is the main source of
irregular global communication in the benchmark;

ZModel:

This module calculates the derivatives of position and vorticity on points on
the SurfaceMesh using either the low-, medium-, or high- order method. To do
so, it uses an external FFT solver, the BRSolver module, and laplacians and
surface normals calculated from the SurfaceMesh as determined by the requested
order of the ZModel solve;

TimeIntegration:

This module uses derivatives calculated using the ZModel module to perform
time integration of points on SurfaceMesh, resulting in the update of the
SurfaceMesh for the next timestep.

The fundamental decomposition of data resulting in MPI communication in
Beatnik is of points on the SurfaceMesh. A 2D block decompostion of this mesh
is most natural decause every derivative calculation in the ZModel module
requires calculating surface normals and laplacians, distributed FFT solvers
generally expect block-decomposed data, and and the current Z-Model
formulation is based on a regular 2D mesh. However, other data decompositions
are possible, for example ones based on the x/y/z location of the mesh point
to optimize far-field force calculations.

##  3 Implementation

Starting from the high-level architecture described in the previous section,
we implemented Beatnik in C++ using the Cabana particle/grid library [14].
Cabana uses the Kokkos performance portability ecosystem [3] to support
portable execution on a wide range of modern hardware architectures, and
provides robust distributed mesh abstractions as well as interfaces to
multiple external solvers valuable for implementing Beatnik. Key among the
external libraries to which Cabana provided access were the heFFTe GPU-
accelerated fast Fourier transform library [1], the ArborX library for
geometric neighbor search [9], and the Silo I/O library [13]. The availability
of robust mesh data structures and key solver interfaces in Cabana
dramatically reduced the size and complexity of the Beatnik code base.

Beatnik’s current C++ implementation consists of approximately 3800 source
lines of code grouped into two set of modules: core modules which implement
fundamental data structures and algorithms common to all solution methods, and
Birchoff-Rott (BR) solver modules which implement different approaches to
estimating the far-field forces required by the Z-Model’s high-order and
medium-order solvers. We detail these modules and their communication behavior
in the remainder of this section.

###  3.1 Core Modules

The core module structure of Beatnik was taken from the Cabana ExaMPM
benchmark [5]. In particular, Beatnik’s core code is arranged as a C++ header-
only library for use by driver programs that specify specific problems to
solve. The following C++ classes providing core Z-Model solution
functionality:

Solver

initializes and invokes other classes based on parameters passed by the driver
program and runs the simulations for the specified number of timesteps.

SurfaceMesh

holds the Cabana structured grid object that encompasses the 2D domain of the
interface between the two fluids. This mesh is decomposed between MPI
processes using a regular 2D block decomposition. Beatnik uses Cabana’s halo
exchange primitives on the SurfaceMesh to perform two-node-deep stencils for
calculating surface normals, finite differences, and Laplacians along the
surface.

BoundaryCondition

applies periodic or non-periodic boundary conditions to the nodes in the
SurfaceMesh; most boundary condition handling is provided by the underlying
Cabana mesh object, but the BoundaryCondition explicitly corrects x/y/z
coordinates in ghost cells for periodic boundary conditions and extrapolates
position and vorticity into boundary cells for non-periodic boundary
conditions. No additional communication is required for these calculations.

ProblemManager

stores the mesh state shared between multiple classes and maintained across
multiple processes, primarily the position and vorticity values on the
SurfaceMesh for the current time step. It also provides interfaces to the
Cabana halo exchange for both SurfaceMesh state stored by the ProblemManager
as well for use by any solvers that store their own SurfaceMesh objects for
internal state.

ZModel

computes derivatives of interface position and velocity using heFFTe and BR
solver modules as necessary for different model orders. The model order used
by an instance of the ZModel class is determined by a template type tag
(Order::Low, Order::Medium, Order::High) that chooses the appropriate
specialized derivative calculation method. The ZModel class does _not_
directly communicate with other MPI processes but instead invokes other
classes which carry out communication specific to the relevant numerical
method, including SurfaceMesh regular halos, heFFTe all-to-all communication,
and any needed global communication by BRSolver modules.

TimeIntegrator

solves for the position and vorticity of each point on the SurfaceMesh using
third-order Runge-Kutta interpolation. Because this is a third-order method,
it calculates three derivatives and hence invokes the ZModel object three
times per timestep.

SiloWriter

uses the Silo library to write surface mesh data for visualization.

###  3.2 Birchoff-Rott Solvers

The Birchoff-Rott solvers for medium and high-order models are generally the
most complex and computationally and communication-intensive portions of
Beatnik. Because of this, the Beatnik implementation supports multiple
Birchoff-Rott solver classes which can be selected by the invoking driver code
and ZModel solver with appropriate template arguments. Each Birchoff-Rott
solver computes the far-field forces between all points on the SurfaceMesh,
and a wide range of numerical methods have been developed for scalability
approximating the far-field forces, all of which require complex communication
and computation.

Beatnik currently implements two relatively simple strategies for calculating
Birchoff-Rott integrals, with implementation of additional strategies planned
for future work as described in Section 6.

####  3.2.1 Exact Integral Calculation

The ExactBRSolver class calculates the exact Birkhoff-Rott velocity integral
using a brute force all-pairs approach. We include an exact solver for the
Birchoff-Rott integral despite its high computational complexity ($O(n^{2})$
in the number of SurfaceMesh points) to enable evaluation of the
accuracy/performance tradeoffs of approximate Birchoff-Rott solvers.
ExactBRSolver uses a standard ring-pass communication algorithm to circulate
SurfaceMesh points between processes in parallel with computing forces between
mesh points. This results in regular communication patterns, but the
ExactBRSolver is compute bound, not communication bound.

####  3.2.2 Cutoff-based Approximate Integral Calculation

The CutoffBRSolver is the first (and currently only) scalable far-field force
solver currently implemented in Beatnik. It approximately solves the Birchoff-
Rott integral by only considering the points within a specified 3D distance of
each point on the SurfaceMesh being evaluated. The cutoff distance determines
the tradeoff between accuracy and performance in this method; small cutoff
distances result in better scalability at the expense of numerical inaccuracy,
while high cutoff distances increase numerical accuracy but reduce performance
by increasing the number of pairs of forces thtat must be calculated. Note
that unlike more sophisticated far-field solution methods like the fast
multipole method, this method does not allow directly setting numerical
accuracy tolerances; instead the programmer choosing the cutoff distance must
set the cutoff distance most appropriate for the problem being simulated.

The cutoff solver produces dynamic and irregular communication patterns
because its computation is based on the _3D spatial location_ of SurfaceMesh
points which are (1) decomposed between MPI processes based on the _2D surface
index_ of the point and (2) which change as the interface surface evolves.
Specifically, the cutoff solver works by, for _each derivative calculation_ :

  1. 1.

Migrating each 2D SurfaceMesh node into a 3D spatial mesh data structure
decomposed by x/y/z location.

  2. 2.

Haloing points between spatial mesh blocks on differnt processors so that each
process has copies of all points within the cutoff distance of each point in
the spatial mesh that a processor owns.

  3. 3.

Using the ArborX solver to compute local neighbor lists for each owned
interface point.

  4. 4.

Calculating the force on each owned interface point by each point on the
neighbor list.

  5. 5.

Migrating each interface point back to its original 2D SurfaceMesh
decomposition for use by the calling ZModel module (which works in the 2D
surface decomposition).

To do so, CutoffSolver uses two supporting C++ classes:

SpatialMesh

holds the Cabana structured grid object that encompasses the 3D domain of the
interface between the two fluids. We currently use a 2D x/y block
decomposition of the 3D space to mirror the initial distribution of 2D surface
points and reduce load imbalance, though load imbalance between these blocks
can develop as surface point spatial locations change.

HaloComm

handles the migration and haloing of nodes between the surface and spatial
mesh based on their x/y/z positions. This redistribution code originated from
the CabanaPD benchmark [11] with slight modifications to fit the design of
Beatnik.

##  4 Mini-Application Driver Problem and Input Decks

Building on the solver implementation described in the previous section, we
implemented a driver for the well-known rocket rig problem as the primary
numerical problem to use for Beatnik benchmarking. This problem simulates the
Rayleigh-Taylor instabilities that develop on the interface between two fluids
of different densities due to their acceleration in the Z direction. The
rocket rig driver program, comprising approximately 700 source lines of code,
initializes and solves variants of this problem with configurable initial
conditions, boundary conditions, and simulation parameters (e.g., cutoff
distance).

For benchmarking MPI communication, we identified two specific rocket rig test
cases: a multi-mode periodic test case shown in Figure 1 and a single-mode
non-periodic test case shown in Figure 2. The multi-mode test case results in
relatively even particle distributions in the X/Y direction, limiting load
imbalance. It is also more amenable to low-order solutions that require FFTs
due to its periodic nature; Beatnik does not currently support low- and
medium-order solutions of non-periodic boundaries because of its reliance on
periodic FFT solvers. In contrast, the non-periodic single-mode test cases
results in significant load imbalance due to surface rollup but requires a
high-order solver to fully resolve these effects. It is important to note that
the single-mode test case is also a useful test case for evaluating the
accuracy of Birchoff-Rott solvers; the outer rollup show in Figure 2 will not
develop without inclusion of distant far-field surface points.

Figure 1: 4-GPU multi-mode test case output using low-order solver at timestep
20 with mesh points colored by vorticity magnitude Figure 2: 36-GPU single-
mode test case output using the cutoff solver at timestep 340 with mesh points
colored by vorticity magnitude

Based on these two test cases, we identified the following benchmark test
cases to stress different MPI communication aspects:

  * •

A multi-mode low-order weak scaling test case that stresses network bandwidth
through all-to-all communication of the growing global SurfaceMesh by the FFT
solver.

  * •

A multi-mode low-order strong scaling test case that stresses network latency
through all-to-all communication of the fixed size SurfaceMesh on increasing
numbers of process by the FFT solver.

  * •

A multi-mode high-order weak scaling test case using the cutoff solver to test
the general scalability of MPI communication with increasing process counts.

  * •

A single-mode high-order strong scaling test case using the cutoff solver to
create load imbalances in communication and computation, and to test MPI
communication performance with dynamic and irregular communication patterns.

This wide range of potential test cases highlights the flexibility of the
Beatnik mini-application design to exercise a wide range of global
communication features.

##  5 Benchmark Evaluation

To demonstrate Beatnik efficacy in exercising communication system, we present
scaling studies for each of the test cases identified in Section 4. We also
demonstrate the use of the Beatnik low-order weak scaling test case to
evaluate the performance of different heFFTe communication configurations.

###  5.1 Experimental Setup

All tests described in this section were run on the Lassen system at Lawrence
Livermore National Laboratories. Lassen is an IBM Power9 system with 40 usable
cores per node, 4 V100 16GB GPUs per node, and an Infiniband EDR interconnect.
All test cases were run with one MPI process and one Power9 core per GPU and
varying the number of processes and GPUs from 4 to 1024 for each test case.

Beatnik on Lassen was compiled using Spack (a Spack configuration for
compiling Beatnik is available in the standard Spack repository). On Lassen,
this results in linking against the current master branch of Cabana, Kokkos
4.2.0, heFFTe 2.4.0, ArborX 1.5, Silo 4.11.1, CUDA 11.1.1, and the rolling
release of IBM Spectrum MPI. Spectrum MPI was configured to use GPU-aware MPI,
preventing us from using newer CUDA versions due to a known bug in Spectrum
MPI.

For both low-order and high-order tests, we set a base problem size to use all
of the memory of each GPU and then either weak or strong scale this problem
size to the appropriate number of processes and GPUs for the test case. For
low-order problems, this results in a (-19, -19, -19) to (19, 19, 19) spatial
domain and 4864 mesh points in the x and y dimensions (roughly 24M mesh
points). For high-order cutoff solver test cases, this results in a (-3, -3,
-3) to (3, 3, 3) spatial domain. For the the multi-mode high-order weak
scaling test we use a single GPU base mesh size of 768x768 (roughly 600k
points) and a cutoff distance of 0.2. For the single-mode high-order strong
scaling test, we use a mesh size of 512x512 (roughly 250k points) and a cutoff
distance of 0.5; tests with smaller cutoff distances resulted in significant
numerical inaccuracy.

Note that the high-order solver is more memory intensive than the low-order
solver because the high-order solver must construct neighbor lists for each
point to do derivative calculations. Similarly, the single mode high-order
test case is more memory intensive that the multi-mode high-order test case
because of longer neighbor lists on some processes due to load imbalance.

###  5.2 Low-Order Solver Weak and Strong Scaling

To test Beatnik’s ability to exercise communication bandwidth and latency
limits on modern systems, we weak and strong scaled the low-order solver,
whose performance is limited by FFT performance, from 4 to 1024 GPUs on
Lassen.

Figure 3: Low-Order Weak Scaling of Beatnik on Lassen Figure 4: Low-Order
Strong Scaling of Beatnik on Lassen

Figures 3 and 4 show the results of these tests. Both cases demonstrate
significant scalability problems due to the intense communication demands of
large-scale FFTs on modern systems. Specifically, low-order solver runtime
increases approximately linearly between 4 and 196 processes and between 256
and 1024 processes but with a smaller slope despite the fact that number of
mesh points per GPU is held constant in this case.

Similarly, strong-scaling performance shows a significant performance increase
when scaling from 4 to 64 GPUs, a parallel efficiency of only 21% (3.5x
speedup when moving from 4 to 64 GPUS). Similarly, performance turns over and
begins to decrease after 64 GPUs due to the small amount of computation and
large number of messages that must be sent in this case. In particular, each
GPU stores and computes on only a 76 by 76 section of the surface mesh in the
64-node case, but must perform all-to-all communication with dozens of other
processes.

###  5.3 High-Order Cutoff Solver Weak Scaling

To test the scalability of the Beatnik benchmark and evaluate the ability of
the underlying communication libraries, we weak scaled the high-order cutoff
solver from 4 to 1024 GPUs on Lassen. Figure 5 shows the results of these
tests.

Figure 5: Total Runtime When Weak-Scaling the Cutoff Solver

As expected, weak scaling Beatnik from 4 to 1024 GPUs results in only modest
(approximately 20%) increases in runtime. This is because the amount of
computation per GPU remains constant and because communication in the weak-
scaling multi-mode case, where little load imbalance develops, is primarily
halo exchanges between neighboring ranks in the 2D SurfaceMesh and 3D
SpatialMesh. We hypothesize that the modest increase in runtime is due to the
overheads of the surface mesh to spatial mesh to surface mesh migration, but
again these overheads are generally small, despite a factor of 256 increase in
problem size, due to the lack of load imbalance in this test case.

###  5.4 Strong Scaling of the Cutoff Solver

As the last scalability evaluation of Beatnik, we tested its ability to create
load imbalances between processes, resulting in dynamic and irregular
communication. We also evaluated the impact of such communication on be
scaling Beatnik from 4 to 256 GPUs on Lassen. As described in Section 4, this
test uses the single mode problem which is prone to load imbalances as the
center of the mesh rolls up as time advances (see Figure 2).

Figure 6: Particles Owned by Each of the 256 Ranks at Timestep 80

Figure 7: Particles Owned by Each of the 256 Ranks at Timestep 340

Figures 7 and 7 show the point distribution amongst 256 processes at timestep
80 and timestep 340. As time advances, the number of particles each processes
owns in 3D space becomes increasingly more imbalanced as the interface surface
rolls up. At timestep 80, the load is evenly distributed, with all processes
owning slightly under 0.4% of all points. At timestep 340, the load is
imbalanced. The processes that own sections of the mesh outside of the rollup
have their load stay the same at about 0.4% of all points. On the other hand,
processes within the rollup own between 0.2% to 0.65% of all points.

Figure 8: Total Runtime of Strong-Scaling the Cutoff Solver

Figure 8 shows the runtime impact of strong scaling on runtime. Scaling from 4
GPUs to 64 GPUs reduces runtime by factor of 3.3 while the number of GPUs is
increased by a factor of 16, a parallel efficiency of 21%. While performance
turns over beyond this point, the performance reduction from additional GPUs
is modest because of the localization of communication provided by the cutoff
solver.

###  5.5 Evaluation of HeFFTe Configurations

Finally, we used the low-order solver in Beatnik to test the performance of
different communication configurations of the HeFFTe Fast Fourier Transform
library [1]. HeFFTe provides three parameters, which can be set to True or
False, to control how Fast Fourier Transform communication and computation is
performed: All-to-All, Penciling, and Reordering.

To do so, we measured low-order solver performanec on all eight permutations
of these configurations to demonstrate Beatnik’s ability to expose how changes
in communication parameters affect benchmark run time; the tested
configurations are shown in Table 1 and numbered for ease of reference in
figures.

Table 1: HeFFTe parameter configurations on the low-order solver Configuration |  AllToAll |  Pencils |  Reorder
---|---|---|---
0 |  False |  False |  False
1 |  False |  False |  True
2 |  False |  True |  False
3 |  False |  True |  True
4 |  True |  False |  False
5 |  True |  False |  True
6 |  True |  True |  False
7 |  True |  True |  True

For these tests, we used the base problem size specified in Section 5.1 weak
scaled from 4 to 1024 processes.

Figure 9 show the Beatnik runtimes of each of these configurations shaded by
the total runtime.

Figure 9: Total Runtime When Weak-Scaling HeFFTe Configurations

Beatnik exposes important performance differences between different heFFTe
configurations. Most notably, heFFTe performance using the _AllToAll_
parameter, which instructs it whether or not to use the builtin `MPI_Alltoall`
primitive, varies significantly between small and large numbers of processes.
On small numbers of processes, heFFTe performance is better when using its
custom communication routines and not using Spectrum MPI’s `MPI_Alltoall`
primitive. In contrast, on large numbers of processes, heFFTe performance
improves if the AllToAll parameter is true.

##  6 Future Work

We have identified multiple areas of future work in Beatnik. First, Beatnik
would be significantly improved by the enhancement of existing and
implementation of additional Birchoff-Rott solvers. For the existing cutoff
solver, implementation of Cabana sparse mesh support would add additional load
balancing communication steps to the benchmark and increase the variety of
communication patterns that could be studied. Similarly, the addition of fast
multipole and P3M far-field force solvers would also increase Beatnik’s
capabilities for examining a wide variety of global communication patterns.
Second, adding remeshing support to Beatnik, a numerical feature still in
development by the Z-model’s authors, would allow it to redistribute or add
points to the surface mesh as the simulation developed. This would both enable
the benchmark to evolve to greater levels of load imbalance and benchmark
additional important global communication patterns. Third, we would like to
examine both the performance and accuracy of the medium-order model when used
with the cutoff solver. Because the medium-order model uses FFTs for
calculating changes in vorticity and supports larger timesteps than the high-
order model, the performance and accuracy tradeoffs between the two models are
potentially interesting. Finally, we would like to enhance Beatnik to support
additional test cases and boundary conditions, including more test cases
examined by the original Z-Model authors, non-periodic boundary conditions for
low- and medium-order solves, and periodic boundary conditions for scalable
high-order solves.

##  7 Conclusion

Beatnik is a novel mini-application designed to support the evaluation of a
wide range of complex communication patterns not frequently tested by other
mini-applications or benchmarks. Beatnik’s implementation is open source and
available for download from github [2], and leverages multiple open-source
solver libraries and performance portability libraries to provide to
significant computational capabilities in a compact, easy-to-understand code
base. As part of Beatnik’s development, we have identified multiple test
cases, each of which exercises a different, important communication feature in
modern HPC systems. Finally, we have demonstrated both the scalability of
Beatnik on modern GPU systems and its ability to evaluate the impact of
communication parameter changes on benchmark performance.

####  7.0.1 Acknowledgements

This work was performed with partial support from the National Science
Foundation under Grants Nos. OAC-2103510, by the U.S. Department of Energy’s
National Nuclear Security Administration (NNSA) under the Predictive Science
Academic Alliance Program (PSAAP-III), Award DE-NA0003966, and by the U.S.
Department of Energy through the Los Alamos National Laboratory (LANL). LANL
is operated by Triad National Security, LLC, for the National Nuclear Security
Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).

The authors would also like to thank Prof. Steve Shkoller for providing
Z-Model Matlab implementations and test cases for use in the development of
the Beatnik benchmark.

####  7.0.2

The authors have no competing interests to declare that are relevant to the
content of this article.

## References

  * [1] Alan Ayala, Stanimire Tomov, Azzam Haidar and Jack Dongarra  “heFFTe: Highly Efficient FFT for Exascale”  In _Computational Science – ICCS 2020_ 12137 Cham: Springer International Publishing, 2020, pp. 262–275  DOI: 10.1007/978-3-030-50371-0_19
  * [2] Patrick G. Bridges et al.  “Beatnik”  URL: https://github.com/CUP-ECS/beatnik
  * [3] H. Carter Edwards, Christian R. Trott and Daniel Sunderland  “Kokkos: Enabling manycore performance portability through polymorphic memory access patterns”  In _Journal of Parallel and Distributed Computing_ 74.12, 2014, pp. 3202–3216  DOI: 10.1016/j.jpdc.2014.07.003
  * [4] Jack Dongarra, Michael A. Heroux and Piotr Luszczek  “A new metric for ranking high-performance computing systems”  In _National Science Review_ 3.1, 2016, pp. 30–35  DOI: 10.1093/nsr/nwv084
  * [5] “ExaMPM”  Co-design center for Particle Applications  URL: https://github.com/ECP-copa/ExaMPM
  * [6] Kenneth Franko, Travis C. Fisher, Paul Lin and Steven W. Bova  “CFD for Next Generation Hardware: Experiences with Proxy Applications”  In _22nd AIAA Computational Fluid Dynamics Conference_ Dallas, TX: American Institute of AeronauticsAstronautics, 2015  DOI: 10.2514/6.2015-3053
  * [7] Paul Henning  “UME”  Los Alamos National Laboratory  URL: https://github.com/lanl/UME
  * [8] Michael Heroux et al.  “The Mantevo Performance Co-design Project”  Sandia National Laboratories  URL: https://github.com/Mantevo
  * [9] D. Lebrun-Grandié, A. Prokopenko, B. Turcksin and S.. Slattery  “ArborX: A Performance Portable Geometric Search Library”  In _ACM Transactions on Mathematical Software_ 47.1, 2021, pp. 1–15  DOI: 10.1145/3412558
  * [10] Gavin Pandya and Steve Shkoller  “3D Interface Models for Rayleigh-Taylor Problems”, 2023  arXiv:2201.04538 [physics.flu-dyn]
  * [11] Sam Reeve and Pablo Seleson  “ORNL/CabanaPD: Version 0.1.0”  [object Object], 2022  DOI: 10.5281/ZENODO.7087781
  * [12] Bob Robey, Neal Davis, David Nicholaeff and Dennis Trujillo  “CLAMR”  Los Alamos National Laboratory  URL: https://github.com/lanl/CLAMR
  * [13] “Silo”  Lawrence Livermore National Laboratory  URL: https://github.com/LLNL/Silo
  * [14] Stuart Slattery et al.  “Cabana: A Performance Portable Library for Particle-Based Simulations”  In _Journal of Open Source Software_ 7.72, 2022, pp. 4115  DOI: 10.21105/joss.04115
  * [15] Ulrike Yang and Rui Peng Li  “AMG2023”  Lawrence Livermore National Laboratory  URL: https://github.com/LLNL/AMG2023


# Spin-1 quarkonia in a rotating frame and their spin contents

HyungJoo Kim<EMAIL_ADDRESS>Department of Physics and Institute of
Physics and Applied Physics, Yonsei University, Seoul 03722, Korea  Sungtae
Cho<EMAIL_ADDRESS>Division of Science Education, Kangwon National
University, Chuncheon 24341, Korea  Su Houng Lee<EMAIL_ADDRESS>Department of Physics and Institute of Physics and Applied Physics, Yonsei
University, Seoul 03722, Korea

###### Abstract

By studying the spin-1 system composed of a heavy quark and its anti-quark in
a rotating frame, we investigate its spin-rotation coupling and thereby reveal
the spin contents of the spin-1 quarkonia. We demonstrate that both the
perturbative and non-perturbative contributions in the operator product
expansion follow a universal relation that identifies the spin-rotation
coupling with unit strength. This allows us to recognize the total spin-1 of
the vector quarkonia in terms of the total angular momentum of quarks and
gluons. We also find that a similar behavior follows for the axial vector case
but identify the different composition of the spin due to the P-wave nature of
the state.

###### pacs:

11.10.Gh,12.38.Bx

The spin polarization of hadrons in heavy-ion collisions has drawn great
interests in recent years. Following the measurements of the global spin
polarization of $\Lambda$ hyperons[1], the global spin alignments of vector
mesons such as $K^{*0}$, $\phi$, and $J/\psi$ have been measured recently by
ALICE and STAR collaborations[2, 3, 4]. While the hyperon polarization
phenomena are consistent with theoretical works[5, 6, 7, 8, 9, 10], the
various aspects of the measured vector meson spin alignments are still not
clearly understood. As the medium in a heavy ion collision will evolve through
different phases, the overall spin polarization will probably depend on both
the total spin and the spin content that varies for different hadrons.

Because hadrons are composite particles, their spin-rotation coupling should
be described by the responses of quarks and gluons in a rotating frame. This
therefore implies that by studying the spin-rotation coupling of a hadron, we
can understand its spin content, which is one of the fundamental questions of
nuclear physics especially in the future EIC(Electron-Ion Collider)
experiment[11]. Therefore, as a first step we study the spin contents of
quarkonia through their spin-rotation coupling.

To address this question, let us first consider two reference frames: an
inertial frame and a (non-inertial) rotating frame which rotates with an
angular velocity $\vec{\Omega}$ with respect to the inertial frame. Then, for
a classical particle, it is easy to verify that the Hamiltonian in the
inertial frame($H_{i}$) and the Hamiltonian in the rotating frame($H_{r}$) are
related by $H_{r}=H_{i}-\vec{L}\cdot\vec{\Omega}$ where $\vec{L}$ is the
orbital angular momentum of the particle. For a particle with intrinsic spin,
it seems natural to generalize this relation to
$H_{r}=H_{i}-(\vec{L}+\vec{S})\cdot\vec{\Omega}$. The appearance of the spin-
rotation coupling was proposed by Mashhoon[12] for particles of any spin but
the idea relied on special relativity with a tentative extension of the
hypothesis of locality and hence warrants a more general derivation based on
general relativity. Indeed, for spin-1/2 Dirac particles, the spin-rotation
coupling can be explicitly derived from the Dirac equation in a rotating
frame,

$\displaystyle[i\gamma^{\mu}D_{\mu}-m+\gamma^{0}(\hat{L}_{q}+\hat{S}_{q})\cdot\vec{\Omega}]\Psi=0,$
(1)

where $D_{\mu}=\partial_{\mu}+igA_{\mu}$ is the covariant derivative and
$\hat{L}_{q}=\vec{x}\times(-i\vec{D})$ and
$\hat{S}_{q}=\frac{1}{2}\gamma^{0}\vec{\gamma}^{\,}\gamma^{5}$ are the orbital
and spin angular momentum operator for Dirac fields, respectively. Expressing
this equation in the form $i\partial_{0}\Psi=H\Psi$, we indeed find
$H_{r}=H_{i}-(\vec{L}+\vec{S})\cdot\vec{\Omega}$ for Dirac particles; see [13,
14, 15, 16, 17, 18, 19, 20] for more details.

On the other hand, for massive spin-1 particles, only limited discussions are
given in [21, 18] and no concrete derivation exists. Recently, Kapusta et al.
investigated the spin-rotation coupling for massive spin-1 vector particles
starting from the Proca equation in a rotating frame[22]. However, they found
that the Hamiltonian in the non-relativistic limit reduces to
$H_{r}=H_{i}-(\vec{L}+\frac{1}{2}\vec{S})\cdot\vec{\Omega}$ at the leading
order in $\vec{\Omega}$(or vorticity), which is not consistent with our
natural expectations. Hence, it is imperative to derive the strength of the
spin-rotation coupling for spin-1 particles especially for vector mesons in a
model independent way based on Quantum Chromodynamics(QCD).

In this work, we introduce a free parameter $g_{\Omega}$, so-called
gravitomagnetic moment in [23], which represents the strength of the spin-
rotation coupling for spin-1 system composed of a heavy quark and its anti-
quark,

$\displaystyle H_{r}=H_{i}-g_{\Omega}\vec{S}\cdot\vec{\Omega}.$ (2)

To derive the value of $g_{\Omega}$ on the basis of quarks and gluons degrees
of freedom, let us first consider a two-point correlation function for the
vector current,

$\displaystyle\Pi^{\mu\nu}(q)=i\int
d^{4}xe^{iqx}\big{\langle}0|\textrm{T}\big{[}j^{\mu}(x)j^{\nu}(0)\big{]}|0\big{\rangle}.$
(3)

Without loss of generality, we put the system at the center of the rotation
and pick out a right circularly polarized state with
$\epsilon_{\mu}^{+}=(0,1,i,0)/\sqrt{2}$ in the rotating frame in which the
polarization axis and the angular velocity are along the same $z$-direction,
i.e. $q_{\mu}=(\omega,0)$ and $\vec{\Omega}=(0,0,\Omega)$. Under this
circumstance, the orbital motion of the system is absent at the center.
Furthermore, the angular velocity is assumed to be small enough so that
rotation effects can be expanded in terms of $\Omega$ by perturbation theory.
Because we are mainly interested in the terms linear in $\Omega$, we define
the relevant component of the correlation function as

$\displaystyle\Pi^{+}(\omega)=\omega^{2}\Pi^{\rm{vac}}(\omega^{2})+\omega\Omega\Pi^{\rm{rot}}(\omega^{2}),$
(4)

where
$\Pi^{+}(\omega)=\Pi^{\mu\nu}(\omega,0)\epsilon_{\mu}^{+}\epsilon_{\nu}^{+*}$.
$\Pi^{\rm{vac}}(\omega^{2})$ is the vacuum invariant function which is used to
study vacuum properties of hadrons and $\Pi^{\rm{rot}}(\omega^{2})$ is a new
function appearing in the rotating frame. Once we obtain
$\Pi^{\rm{rot}}(\omega^{2})$, $g_{\Omega}$ can be extracted by comparing it
with $\Pi^{\rm{vac}}(\omega^{2})$.

To this end, let us first discuss the phenomenological structure of
$\Pi^{+}(\omega)$. In an inertial frame, $\Omega$ dependent terms do not
appear so that $\Pi^{+}(\omega)=\omega^{2}\Pi^{\text{vac}}(\omega^{2})$. Once
we turn on the rotation, the energy of a right circularly polarized state is
shifted by $-g_{\Omega}\Omega$ as we assumed. If this energy shift happens
equally, i.e. $g_{\Omega}$ is universal, for all vector states contributing in
Eq.(3), we can expect that
$\Pi^{+}(\omega)\to\Pi^{+}(\omega+g_{\Omega}\Omega)$ in the rotating frame. By
comparing this with Eq.(4), we can infer a simple relation between
$\Pi^{\text{vac}}(\omega^{2})$ and $\Pi^{\text{rot}}(\omega^{2})$,

$\displaystyle\Pi^{\text{rot}}_{\text{phen}}(\omega^{2})=2g_{\Omega}\Big{\\{}\Pi^{\rm{vac}}(\omega^{2})+\omega^{2}\frac{\partial\Pi^{\rm{vac}}(\omega^{2})}{\partial\omega^{2}}\Big{\\}},$
(5)

where the subscript ‘phen’ indicates that this relation is derived from a
phenomenological point of view. Furthermore, by applying this relation to the
standard vacuum dispersion relation, we obtain the following useful expression

$\displaystyle\Pi^{\rm{rot}}_{\text{phen}}(\omega^{2})=\frac{g_{\Omega}}{\pi}\int^{\infty}_{4m^{2}}ds\frac{2s\mathrm{Im}\Pi^{\rm{vac}}(s)}{(s-\omega^{2})^{2}}.$
(6)

In the following, we will show that the Operator Product Expansion(OPE)
satisfies a universal relation where the ‘rot’ part of the perturbative and
non-perturbative contributions are related to their respective ‘vac’ parts as
in Eq.(5) but with $g_{\Omega}=1$. For that purpose, let us move on to the OPE
side. For a heavy quark system, we need to compute the following function in
the deep Euclidean space ($Q^{2}=-q^{2}\gg 0$),

$\displaystyle\Pi^{\mu\nu}(q)=i\int
d^{4}xe^{iqx}\big{\langle}0|\textrm{Tr}\big{[}\gamma^{\mu}S(x,0)\gamma^{\nu}S(0,x)\big{]}|0\big{\rangle},$
(7)

where $S(x,0)$ and $S(0,x)$ denote quark propagators. To compute this function
in the rotating frame, we should use quark propagators in the rotating frame.
Referring to Eq.(1), we use the following expansion of the quark propagator in
terms of the coupling constant and $\Omega$,

$\displaystyle S(x$
$\displaystyle,0)=S^{(0)}(x)+\sum^{\infty}_{n=1}(-1)^{n}\int dz_{1}\ldots
dz_{n}S^{(0)}(x-z_{1})$
$\displaystyle\times\big{[}\Delta\mathcal{I}(z_{1})\big{]}S^{(0)}(z_{1}-z_{2})\ldots\big{[}\Delta\mathcal{I}(z_{n})\big{]}S^{(0)}(z_{n}),$
(8)

where $S^{(0)}(x)$ is the free quark propagator and
$\Delta\mathcal{I}=g\not{A}+\gamma^{0}((\hat{L}_{q})_{z}+(\hat{S}_{q})_{z})\Omega$
includes all the interaction terms. For later convenience, we also distinguish
the orbital angular momentum operator into two pieces,
$\hat{L}_{q}=\hat{L}_{k}+\hat{L}_{p}$ where $\hat{L}_{k}=\vec{x}\times\vec{p}$
is the kinetic part and $\hat{L}_{p}=\vec{x}\times(-g\vec{A}(x))$ is the
potential part, respectively. By collecting terms linear in $\Omega$ from
Eq.(8), we compute the OPE for $\Pi^{\rm{rot}}$ upto operators of dimension 4,

$\displaystyle\Pi^{\text{rot}}_{\text{OPE}}(Q^{2})=\Pi^{\text{rot}}_{\text{I}}(Q^{2})+\Pi^{\text{rot}}_{G_{0}}(Q^{2}),$
(9)

where $\Pi_{\text{I}}^{\text{rot}}$ denotes the leading perturbative part and
$\Pi^{\text{rot}}_{G_{0}}=C(Q^{2})\cdot G_{0}$ denotes the leading non-
perturbative part with $C(Q^{2})$ being the Wilson coefficients for the scalar
gluon condensates
$G_{0}\equiv\langle{\frac{\alpha_{s}}{\pi}G^{a}_{\mu\nu}G^{a,\mu\nu}}\rangle=(0.35\,\text{GeV})^{4}$.

Figure 1: The leading perturbative diagram

The leading perturbative diagram which contributes to
$\Pi_{\rm{I}}^{\rm{rot}}$ is shown in Fig.1. Here, the propagator with one
star symbol($\bigstar$) indicates the insertion of interaction terms linear in
$\Omega$ on the quark propagator. The perturbative result can be expressed by
a simple dispersion relation,

$\displaystyle\Pi^{\rm{rot}}_{\rm{I}}(Q^{2})$
$\displaystyle=\frac{1}{\pi}\int^{\infty}_{4m^{2}}ds\frac{\mathrm{Im}\Pi^{\rm{rot}}_{\text{I}}(s)}{s+Q^{2}},$
(10)

where
$\Pi^{\text{rot}}_{\text{I}}(s)=\Pi^{\rm{rot}}_{\text{I},S_{q}}(s)+\Pi^{\rm{rot}}_{\text{I},{L_{k}}}(s)$
and

$\displaystyle\mathrm{Im}\Pi^{\rm{rot}}_{\text{I},S_{q}}$
$\displaystyle(s)=\frac{3m^{2}}{2\pi\sqrt{s(s-4m^{2})}},$ (11)
$\displaystyle\mathrm{Im}\Pi^{\rm{rot}}_{\text{I},{L_{k}}}$
$\displaystyle(s)=\frac{(s-m^{2})\sqrt{s(s-4m^{2})}}{2\pi s^{2}}.$ (12)

Here, $\hat{L}_{p}$ has no contribution because it already contains one gluon
field. Now we can extract $g_{\Omega}$ in the perturbative region by comparing
this OPE result with the corresponding phenomenological expression given in
Eq.(6),

$\displaystyle\frac{1}{\pi}\int^{\infty}_{4m^{2}}ds\frac{\mathrm{Im}\Pi^{\rm{rot}}_{\text{I}}(s)}{s+Q^{2}}=\frac{g_{\Omega}}{\pi}\int^{\infty}_{4m^{2}}ds\frac{2s\mathrm{Im}\Pi^{\text{vac}}_{\text{I}}(s)}{(s+Q^{2})^{2}},$
(13)

where $\Pi^{\text{vac}}_{\text{I}}(s)$ is the leading perturbative
contribution for the vacuum OPE which is given by $u(3-u^{2})/(8\pi)$ with
$u=\sqrt{1-4m^{2}/s}$. Here, the left hand side is directly computed by
collecting all responses of the quarks in the rotating frame and the right
hand side is derived from a phenomenological point of view. Interestingly, we
find

$\displaystyle\text{Im}\Pi^{\text{rot}}_{\text{I}}(s)=2\Big{\\{}\text{Im}\Pi^{\rm{vac}}_{\text{I}}(s)+s\frac{\partial\text{Im}\Pi^{\rm{vac}}_{\text{I}}(s)}{\partial
s}\Big{\\}}.$ (14)

Therefore, by integration by parts of the right hand side, we can explicitly
show that $g_{\Omega}=1$ in the leading perturbative region. The leading
perturbative diagram describes the free quark system so this finding indicates
that when free quarks form a spin-1 vector state, their total energy shift is
just given by the sum of the energy shift of each quark in the rotating frame.
It should be noted that this goes beyond the non-relativistic quark model as
both the spin and orbital angular momentum of quarks are taken into account
relativistically.

Figure 2: The leading non-perturbative diagrams

Now let us examine the non-perturbative region of the OPE. The leading non-
perturbative gluon contribution comes from diagrams shown in Fig.2. Here, each
curly line indicates one background gluon field inside the propagator in
Eq.(8). Expressing the Wilson coefficients with
$J_{N}=\int^{1}_{0}dx[1+x(1-x)Q^{2}/m^{2}]^{-N}$, the result is given as

$\displaystyle C_{S_{q}}$
$\displaystyle=-\frac{2-y+12J_{2}-26J_{3}+12J_{4}}{12Q^{4}},$ (15)
$\displaystyle C_{L_{k}}$
$\displaystyle=\frac{11-4y-(8+y)J_{1}+13J_{2}-16J_{3}}{72Q^{4}},$ (16)
$\displaystyle C_{L_{p}}$
$\displaystyle=-\frac{13+2y-(12+5y)J_{1}+3J_{2}-4J_{3}}{72Q^{4}}.$ (17)

By comparing this result with the counterpart of the vacuum OPE using Eq.(5),
we find that $g_{\Omega}$ falls short from being 1. In fact, what we have
computed in the OPE until now is nothing else but the expectation value of
quark’s total angular momentum operator with respect to a right circularly
polarized vector current. But there could be additional contribution from
gluon’s angular momentum. At this moment, it is not clearly understood how its
contribution comes out in the rotating frame. Instead, regardless of how it is
derived in the rotating frame, we can add its contribution directly by
computing the expectation value of the gluon’s total angular momentum
operator, $\hat{J}_{g}=\vec{x}\times(\vec{E}\times\vec{B})$. Following the
method discussed in [24], we compute the Wilson coefficients for
$\hat{J}_{g}$,

$\displaystyle C_{J_{g}}$
$\displaystyle=\frac{13-(2+2y)J_{1}-23J_{2}+12J_{3}}{36Q^{4}}.$ (18)

After including this contribution, we indeed find

$\displaystyle\Pi^{\rm{rot}}_{G_{0}}(Q^{2})=2\Big{\\{}\Pi^{\rm{vac}}_{G_{0}}(Q^{2})+Q^{2}\frac{\partial\Pi^{\rm{vac}}_{G_{0}}(Q^{2})}{\partial
Q^{2}}\Big{\\}},$ (19)

where
$\Pi^{\text{rot}}_{G_{0}}=(C_{S_{q}}+C_{L_{k}}+C_{L_{p}}+C_{J_{g}})G_{0}$ and
$\Pi^{\text{vac}}_{G_{0}}=\frac{1}{12Q^{4}}(-1+3J_{2}-2J_{3})G_{0}$.
Therefore, $g_{\Omega}=1$ even in the non-perturbative region, once the total
angular momentum of quarks and gluons are considered. From this result, we can
conclude that
$g_{\Omega}=\langle{(\vec{J}_{\text{QCD}})_{z}}\rangle/\langle{(\vec{S}\,)_{z}}\rangle$
where $\vec{J}_{\text{QCD}}=\int
d^{3}x[\frac{1}{2}\bar{\psi}\vec{\gamma}\gamma_{5}\psi+\psi^{\dagger}(\vec{x}\times(-i\vec{D}))\psi+\vec{x}\times(\vec{E}\times\vec{B})]$
is the total angular momentum in QCD[25] and $\vec{S}$ is the spin-1 operator.
Then our finding, i.e. $g_{\Omega}=1$, shows that the total spin of the system
is equal to the total angular momentum of quarks and gluons. We also compute
OPE for the axial vector case and observe that it satisfies the same universal
relation as in the vector case. The OPE results for the axial vector are
listed below.

$\displaystyle\mathrm{Im}\Pi^{\rm{rot}}_{\text{I},{S_{q}}}$
$\displaystyle(s)=\frac{3m^{2}\sqrt{s(s-4m^{2})}}{2\pi s^{2}},$ (20)
$\displaystyle\mathrm{Im}\Pi^{\rm{rot}}_{\text{I},{L_{k}}}$
$\displaystyle(s)=\frac{(s-m^{2})\sqrt{s(s-4m^{2})}}{2\pi s^{2}},$ (21)
$\displaystyle C_{S_{q}}$ $\displaystyle=-\frac{6-y-12J_{2}+6J_{3}}{12Q^{4}},$
(22) $\displaystyle C_{L_{k}}$
$\displaystyle=\frac{7-4y-(8+y)J_{1}+25J_{2}-24J_{3}}{72Q^{4}},$ (23)
$\displaystyle C_{L_{p}}$
$\displaystyle=-\frac{9+2y+(4-5y)J_{1}-17J_{2}+4J_{3}}{72Q^{4}},$ (24)
$\displaystyle C_{J_{g}}$
$\displaystyle=\frac{1+(6-2y)J_{1}-3J_{2}-4J_{3}}{36Q^{4}}.$ (25)

Now we can calculate the spin contents of spin-1 quarkonia using QCD sum
rules(QCDSR). Although $g_{\Omega}$ should be 1 even for bound states, we are
particularly interested in the fraction of $g_{\Omega}$ carried by each
angular momentum operator($\hat{S}_{q}$, $\hat{L}_{k}$, $\hat{L}_{p}$,
$\hat{J}_{g}$). To this end, let us recall the phenomenological expression
once again. In Eq.(6), $\text{Im}\Pi^{\text{vac}}(s)$ represents the spectral
density which is modeled to have a single ground state pole and perturbative
continuum,

$\displaystyle\text{Im}\Pi^{\text{vac}}(s)=\pi
f_{0}\delta(s-m_{0}^{2})+\theta(s-s_{0})\text{Im}\Pi^{\text{vac}}_{\text{I}}(s),$
(26)

where $f_{0}$ is the residue, $m_{0}$ is the ground state mass, $s_{0}$ is the
threshold for the leading perturbative continuum. Then the phenomenological
side is expressed as

$\displaystyle\Pi^{\rm{rot}}_{\text{phen}}(Q^{2})=\frac{2g_{\Omega}f_{0}}{(Q^{2}+m_{0}^{2})^{2}}+\frac{1}{\pi}\int^{\infty}_{s_{0}}ds\frac{\mathrm{Im}\Pi^{\text{rot}}_{\text{I}}(s)}{(s+Q^{2})},$
(27)

where the continuum part is re-expressed in terms of
$\Pi^{\text{rot}}_{\text{I}}$ using Eq.(13). By connecting this expression to
the OPE side and performing Borel transformation, we can obtain a new sum
rules for $\Pi^{\text{rot}}$ as

$\displaystyle\bar{\mathcal{M}}^{\text{rot}}$
$\displaystyle\equiv\mathcal{B}\big{[}\Pi^{\text{rot}}_{\text{OPE}}(Q^{2})\big{]}-\int^{\infty}_{s_{0}}dse^{-s/M^{2}}\mathrm{Im}\Pi^{\text{rot}}_{\text{I}}(s)$
$\displaystyle=\frac{2\pi
g_{\Omega}f_{0}m_{0}^{2}}{M^{2}}e^{-m_{0}^{2}/M^{2}}.$ (28)

Here, Borel transformation is defined by

$\displaystyle\mathcal{B}\equiv\lim_{\begin{subarray}{c}Q^{2}/n\rightarrow
M^{2},\\\
n,Q^{2}\rightarrow\infty\end{subarray}}\frac{\pi(Q^{2})^{n+1}}{n!}\left(-\frac{d}{dQ^{2}}\right)^{n},$
(29)

where $M$ is the Borel mass parameter. Similarly, the vacuum sum rules is
expressed as

$\displaystyle\bar{\mathcal{M}}^{\rm{vac}}$
$\displaystyle\equiv\mathcal{B}\big{[}\Pi^{\text{vac}}_{\text{OPE}}(Q^{2})\big{]}-\int^{\infty}_{s_{0}}dse^{-s/M^{2}}\mathrm{Im}\Pi^{\rm{vac}}_{\text{I}}(s)$
$\displaystyle=\pi f_{0}e^{-m_{0}^{2}/M^{2}}.$ (30)

Then $g_{\Omega}$ for the ground state can be extracted from the following
equation,

$\displaystyle
g_{\Omega}(M,s_{0})=-\frac{M^{2}}{2}\frac{\bar{\mathcal{M}}^{\rm{rot}}}{\partial\bar{\mathcal{M}}^{\rm{vac}}/\partial(1/M^{2})}.$
(31)

In an actual analysis, we need to specify appropriate threshold and a reliable
range of Borel mass, so-called Borel window which is determined by competition
between OPE convergence and pole-dominance. Here, we will just employ the
values for $\sqrt{s_{0}}$ and Borel window determined from the conventional
analysis for $m_{0}$ through the vacuum sum rules in the presence of
$\alpha_{s}$-correction[26] as the rotating part is related to the vacuum part
through the universal relation. The ground state mass can be extracted in a
similar manner from Eq.(30),

$\displaystyle
m_{0}(M,s_{0})=\sqrt{-\frac{\partial\bar{\mathcal{M}}^{\text{vac}}/\partial(1/M^{2})}{\bar{\mathcal{M}}^{\text{vac}}}}.$
(32)

Specifically, we determine $\sqrt{s_{0}}$ and Borel window by making $m_{0}$
least sensitive to $M$; See [27, 28] for more details. Here, to determine the
Borel window, we assume the sum of $\alpha_{s}$-correction and $G_{0}$
contribution is less than 30$\%$ of the total vacuum OPE and pole contribution
is dominant by above 60$\%$. In this work, we use
$m_{c}(p^{2}=-m_{c}^{2})=1.262\,\text{GeV}$, $\alpha_{s}(8m_{c}^{2})=0.21$ for
charmonia and $m_{b}(p^{2}=-m_{b}^{2})=4.12\,\text{GeV}$,
$\alpha_{s}(8m_{b}^{2})=0.158$ for bottomonia following [28]. Values for
$\sqrt{s_{0}}$ and Borel window used in this work are listed in Table 1.

Table 1: $\sqrt{s_{0}}$ and Borel window for spin-1 quarkonia222We follow the $\alpha_{s}$-correction for the axial vector given in [29]. |  $J/\psi$ |  $\chi_{c1}$ |  $\Upsilon$(1S) |  $\chi_{b1}$
---|---|---|---|---
$\sqrt{s_{0}}$ [GeV] |  3.5 |  4.0 |  10.3 |  11
($M_{\text{min}}$,$M_{\text{max}}$) [GeV] |  (1,2.3) |  (1.4,2.3) |  (3,5.5) |  (3.6,4.9)

Then we finally estimate the spin contents of spin-1 quarkonia by averaging
each component of $g_{\Omega}$ over the Borel window. We also calculate the
variance to estimate the uncertainty. The average values and their
uncertainties(subscript) are listed in Table 2.

Table 2: The spin contents of spin-1 quarkonia | Vector | Axial vector
---|---|---
|  $J/\psi$ |  $\Upsilon$(1S) |  $\chi_{c1}$ |  $\chi_{b1}$
$S_{q}$ |  $0.88_{1.8\text{e-}4}$ |  $0.92_{7.6\text{e-}5}$ |  $0.40_{8.2\text{e-}5}$ |  $0.43_{1.1\text{e-}5}$
$L_{k}$ |  $0.11_{4.9\text{e-}4}$ |  $0.076_{7.8\text{e-}5}$ |  $0.61_{5.8\text{e-}6}$ |  $0.57_{1.0\text{e-}5}$
$L_{p}$ |  $2.0\text{e-}3_{2.9\text{e-}6}$ |  $3.5\text{e-}5_{3.0\text{e-}10}$ |  $8.2\text{e-}4_{2.3\text{e-}8}$ |  $\textbf{--}1.0\text{e-}5_{3.4\text{e-}10}$
$J_{g}$ |  $8.0\text{e-}3_{5.9\text{e-}5}$ |  $1.5\text{e-}4_{7.3\text{e-}9}$ |  $\textbf{--}0.015_{5.2\text{e-}5}$ |  $\textbf{--}5.2\text{e-}5_{2.3\text{e-}8}$

In all cases the sum of the four components is exactly 1 as we expected, but
the spin contents are quite different from each other. While the quark spin
has the most dominant contribution for the vector quarkonia, the kinetic part
of the orbital angular momentum has larger contribution for the axial vector
quarkonia. The bottomonia results are somewhat comparable with the non-
relativistic quark model picture but the spin contents start to deviate from
this picture as the quark mass becomes lighter due to the significant increase
in the contribution of $\hat{L}_{k}$. For example, the $J/\psi$ is
traditionally considered as an S-wave particle but now we find that the quark
spin does not carry all of the total spin as in the case of the proton
spin[30]. Furthermore, the overall contributions of $\hat{L}_{p}$ and
$\hat{J}_{g}$ are quite small but indispensable to make $g_{\Omega}=1$
exactly. In the perturbative region, $\hat{L}_{p}$ and $\hat{J}_{g}$ start to
contribute at the next-to-leading order diagrams, so more accurate analysis
requires the computation of $\alpha_{s}$-correction for $\Pi^{\text{rot}}$.

In conclusion, we have proven that $g_{\Omega}=1$ for a heavy composite
particle in QCD and simultaneously have identified how the angular momenta of
quarks and gluons add up to its total spin. The most crucial finding in this
work is the discovery of the universal relation between the rotating part and
the corresponding vacuum part in the OPE. Because this relation indicates that
the total spin of the system is equal to the total angular momentum of its
constituents, the same formalism can be applied to any other systems to
elucidate their spin contents. Therefore, it will definitely play a crucial
role when investigating the spin contents of more complicated hadrons
especially for nucleons.

## Acknowledgements

This work was supported by Samsung Science and Technology Foundation under
Project Number SSTF-BA1901-04, and by the Korea National Research Foundation
under the grant number No.2020R1F1A1075963 and No.2019R1A2C1087107.

## References

  * Adamczyk  _et al._ [2017] L. Adamczyk _et al._ (STAR), Nature 548, 62 (2017), arXiv:1701.06657 [nucl-ex] .
  * Acharya  _et al._ [2020] S. Acharya _et al._ (ALICE), Phys. Rev. Lett. 125, 012301 (2020), arXiv:1910.14408 [nucl-ex] .
  * Abdallah  _et al._ [2022] M. Abdallah _et al._ (STAR), (2022), arXiv:2204.02302 [hep-ph] .
  * Acharya  _et al._ [2022] S. Acharya _et al._ (ALICE), (2022), arXiv:2204.10171 [nucl-ex] .
  * Liang and Wang [2005] Z.-T. Liang and X.-N. Wang, Phys. Rev. Lett. 94, 102301 (2005), [Erratum: Phys.Rev.Lett. 96, 039901 (2006)], arXiv:nucl-th/0410079 .
  * Betz  _et al._ [2007] B. Betz, M. Gyulassy, and G. Torrieri, Phys. Rev. C 76, 044901 (2007), arXiv:0708.0035 [nucl-th] .
  * Becattini  _et al._ [2008] F. Becattini, F. Piccinini, and J. Rizzo, Phys. Rev. C 77, 024906 (2008), arXiv:0711.1253 [nucl-th] .
  * Becattini and Piccinini [2008] F. Becattini and F. Piccinini, Annals Phys. 323, 2452 (2008), arXiv:0710.5694 [nucl-th] .
  * Abelev  _et al._ [2007] B. I. Abelev _et al._ (STAR), Phys. Rev. C 76, 024915 (2007), [Erratum: Phys.Rev.C 95, 039906 (2017)], arXiv:0705.1691 [nucl-ex] .
  * Adam  _et al._ [2018] J. Adam _et al._ (STAR), Phys. Rev. C 98, 014910 (2018), arXiv:1805.04400 [nucl-ex] .
  * Accardi  _et al._ [2016] A. Accardi _et al._ , Eur. Phys. J. A 52, 268 (2016), arXiv:1212.1701 [nucl-ex] .
  * Mashhoon [1988] B. Mashhoon, Phys. Rev. Lett. 61, 2639 (1988).
  * de Oliveira and Tiomno [1962] C. G. de Oliveira and J. Tiomno, Nuovo Cim. 24, 672 (1962).
  * Hehl and Ni [1990] F. W. Hehl and W.-T. Ni, Phys. Rev. D 42, 2045 (1990).
  * De Oliveira and Tiomno [1962] C. De Oliveira and J. Tiomno, Il Nuovo Cimento (1955-1965) 24, 672 (1962).
  * Huang [1994] J. C. Huang, Annalen der Physik 506, 53 (1994).
  * Papini [2002] G. Papini, Phys. Rev. D 65, 077901 (2002), arXiv:gr-qc/0201098 .
  * Cai and Papini [1991] Y. Q. Cai and G. Papini, Phys. Rev. Lett. 66, 1259 (1991).
  * Ryder [1998] L. Ryder, Journal of Physics A: Mathematical and General 31, 2465 (1998).
  * Ryder [2008] L. Ryder, General Relativity and Gravitation 40, 1111 (2008).
  * Mashhoon [1989] B. Mashhoon, Physics Letters A 139, 103 (1989).
  * Kapusta  _et al._ [2020] J. I. Kapusta, E. Rrapaj, and S. Rudaz, Phys. Rev. D 102, 125028 (2020), arXiv:2009.12010 [hep-th] .
  * Buzzegoli and Kharzeev [2021] M. Buzzegoli and D. E. Kharzeev, Phys. Rev. D 103, 116005 (2021), arXiv:2102.01676 [hep-th] .
  * Balitsky and Ji [1997] I. Balitsky and X.-D. Ji, Phys. Rev. Lett. 79, 1225 (1997), arXiv:hep-ph/9702277 .
  * Ji [1997] X.-D. Ji, Phys. Rev. Lett. 78, 610 (1997), arXiv:hep-ph/9603249 .
  * Reinders  _et al._ [1985] L. J. Reinders, H. Rubinstein, and S. Yazaki, Phys. Rept. 127, 1 (1985).
  * Gubler and Satow [2019] P. Gubler and D. Satow, Prog. Part. Nucl. Phys. 106, 1 (2019), arXiv:1812.00385 [hep-ph] .
  * Morita and Lee [2010] K. Morita and S. H. Lee, Phys. Rev. D 82, 054008 (2010), arXiv:0908.2856 [hep-ph] .
  * Reinders  _et al._ [1981] L. J. Reinders, H. R. Rubinstein, and S. Yazaki, Nucl. Phys. B 186, 109 (1981).
  * Ashman  _et al._ [1988] J. Ashman _et al._ (European Muon), Phys. Lett. B 206, 364 (1988).