# A Computationally Efficient Pipeline Approach to Full Page Offline
Handwritten Text Recognition

Jonathan Chung  Amazon Inc.
Vancouver, BC
<EMAIL_ADDRESS>Thomas Delteil  Amazon Inc.
Vancouver, BC
<EMAIL_ADDRESS>
###### Abstract

Offline handwriting recognition with deep neural networks is usually limited
to words or lines due to large computational costs. In this paper, a less
computationally expensive full page offline handwritten text recognition
framework is introduced. This framework includes a pipeline that locates
handwritten text with an object detection neural network and recognises the
text within the detected regions using features extracted with a multi-scale
convolutional neural network (CNN) fed into a bidirectional long short term
memory (LSTM) network. This framework achieves comparable error rates to state
of the art frameworks while using less memory and time. The results in this
paper demonstrate the potential of this framework and future work can
investigate production ready and deployable handwritten text recognisers.

###### Index Terms:

OCR, handwritten text recognition, text analysis, recurrent neural networks,
convolutional neural networks

##  I Introduction

Due to the heterogeneous nature of handwritten text, it is much more difficult
to automatically recognise compared to printed text (two orders of magnitude
of difference in error rate, see Table 3 vs Table 5 in [1]). Recent successes
in handwriting recognition can be attributed to developments in deep neural
networks. However, due to large computational costs, the systems are usually
limited to recognising characters, words, and lines. We propose a full page
offline handwriting recognition framework that is less computationally
expensive compared to existing frameworks.

##  II Related work

###  II-A Text localisation

Text localisation is an essential component of document layout analysis and
accurate text localisation is crucial for handwriting recognition [2].
Handcrafted features that utilise blob detection, clustering, edge detection,
and histogram projections dominate in the traditional techniques. More
recently, data-driven techniques are becoming more prominent with the growth
of neural networks. Such techniques can be categorised based on the method in
which the position of the text is defined. This includes lines [3], bounding
boxes [4], or areas containing “text pixels” [2].

In this paper, we predict bounding boxes around the text using deep learning
techniques of object detection. Given an image that contains multiple objects,
object detection identifies bounding boxes that encompass the objects along
with the confidence of the class of the object. In this work, the Single Shot
MultiBox Detector (SSD) [5] framework was applied to text localisation.

###  II-B Text recognition

Handwritten digit recognition with the MNIST dataset was among the first work
in deep learning [6]. However, the learning problem was limited to images of
single digit characters. Significant advances in handwritten text recognition
were realised by the description of the multidimensional recurrent neural
networks (MD-RNN) in Graves et al. [7] and the Connectionist Temporal
Classification (CTC) loss. A number of advances based on the MD-RNN were
reported including using attribute embeddings [8], dropout [9], Tucker
decomposition [10] etc. Recent works conducted by Puigcerver [11] suggest that
the multidimensional aspects of the MD-RNN can be replaced with feeding image-
features (from a CNN) into a one-dimensional LSTM to significantly reduce the
memory requirements of the systems. The described methods are either limited
to single words or single lines of handwritten text. Bluche et al. [12, 13]
described an end-to-end system that uses an MD-RNN along with an LSTM to
encode multiple lines of text. Although the described system shows promise to
automatically recognise multiple lines, it may not be a practical solution as
it requires a large amount of computational power [11]. Winglinton et al. [14]
utilised a region proposal network to find the starting positions of text
lines and a line follower network was trained to trace the line of text. This
was followed by using a CNN-LSTM approach to recognise the characters.

Figure 1: Overview of the system described.

###  II-C Approach overview

Previous works showed that the MD-RNN requires a large amount of computation
power when it is used to recognise multiple lines of handwritten text. A less
computationally expensive alternative could be realised if multiple lines of
handwriting recognition were not directly performed. To achieve this, our
described framework is comprised of two major components: text localisation
and recognition. Text localisation identifies the positions of handwritten
text given an image of the full page. Once a passage of handwritten text was
identified, segmentation was conducted to locate each line of text. Text
recognition refers to converting an image of a line of handwritten text into a
string with the corresponding characters and denoising the string with a
language model. By limiting handwriting recognition to single lines, the
computational costs associated with this framework can be dramatically reduced
compared to previous works that utilise the MD-RNN.

##  III Methods

Rather than designing an end-to-end network, we took a modular approach
consistent with described components in the literature. This principle allows
components of the framework to be easily replaced and tested with different
ones. An overview of our system is provided in Figure 1.

###  III-A Text localisation

The purpose of text localisation is to identify bounding boxes of each line of
text given an image containing both printed and handwritten text. The text
localisation procedure consists of two stages: passage identification and line
segmentation.

####  III-A1 Passage identification

The goal of passage identification is to predict the location of the
handwritten passage (bounding boxes containing $x$, $y$ coordinates, width and
height of the bounding box in percentages of the page size). To simplify this
step, we assume that there is one passage of printed text and one passage of
handwritten text (using the IAM dataset [15] see Section IV-A for more
details). This was achieved by extracting image features from a pre-trained
truncated 34 layer residual network (ResNet34) [16] trained on ImageNet. In
the ResNet34, the weights of the first convolutional layer were averaged into
one channel to accommodate for greyscale images. The features were then fed
into three fully connected layers: two layers with 64 units and a relu
activation and one layer with 4 units and a sigmoid activation. The four units
with sigmoid activation correspond to the $x$, $y$ coordinates, width, and
height of the bounding box in percentages. The network was trained to minimise
the mean squared error.

####  III-A2 Line segmentation

Given an image containing only handwritten text, this component predicts
bounding boxes surrounding each line of text. We modelled this as an object
detection problem to detect words followed by using a clustering algorithm to
combine words into lines. A two stage approach was taken because early
experiments showed that the network was prone to missing objects when
identifying handwritten text. By detecting individual words, the chance of the
network missing an entire line of words was less likely.

In our implementation, the SSD architecture [5] was used to predict bounding
boxes relative to anchor points and predict the probability that the bounding
boxes are encompassing words. The downsampler consists of two convolutional
layers, batch normalisation layer, and a relu activation function. The class
and bounding box predictor consists of a single convolutional layer with 6 (4
positional + 2 for classes) output channels. To adapt the SSD to our
requirements, image features were extracted with a similar network described
in Section III-A1 (ResNet34). Furthermore, the anchor boxes were adapted to
resemble words (only squares and rectangles with widths $>$ height). The SSD
was trained to minimise the cross-entropy loss for the class (handwriting or
not handwriting) and the L1 loss for the bounding box. Non-maximum suppression
was performed to filter out objects overlapping bounding boxes.

After the bounding boxes of words were detected, a greedy algorithm was then
used to cluster the words into lines proposals based on the overlap in the
y-direction (see Algorithm 1).

Input: _words_ a list of bounding boxes (x, y, w, h)

Output: _lines_ a list of bounding boxes (x, y, w, h)

Sort words based on the y-axis

prev_bb = _None_ , line_bb = []

for  _bb in sorted_words_ do

if  _prev_bb is not _None__ then

if  _get_y_overlap(prev_bb, bb) > 0.4_ then

_#If the word in bb is within the same line_

line_bb.append(bb)

else

_#If the word in bb is in a different line_

output.append(get_bb(line_bb))

line_bb = []

end if

prev_bb = bb

end for

output.append(get_bb(line_bb))

Helper functions:

get_y_overlap \- calculates the percentage overlap between two bounding boxes

get_bb \- gets the encompassing x, y, w, h coordinates given a list of
bounding boxes

Algorithm 1 Word to line clustering algorithm

The following heuristics were used to evaluate the line proposals:

  * •

Lines must have a minimum area

  * •

Lines that exceed boundaries of the page are removed

  * •

Lines (excluding the last line) that are substantially shorter than the median
width of the lines are removed

  * •

Lines that are much longer than the median height are split into 2 lines
(accounts for double lines)

  * •

Lines with starting positions that significantly deviate from other lines are
removed

  * •

Remove lines that greatly overlap with other lines

Lines that are not eliminated by the heuristics algorithm are used as the
output of the text localisation stage.

###  III-B Text recognition

Text recognition takes images containing single lines of handwritten text and
recognises the corresponding characters. Our approach includes handwriting
recognition then denoising the output with a language model.

####  III-B1 Handwriting recognition

Following a similar scheme to [11], we implemented a CNN-biLSTM network. It
makes use of a multi-scale CNN for image feature extraction, then the features
are fed into a bidirectional LSTM. The network was trained to optimise the CTC
loss (shown in Figure 2). Intuitively, the CNN generates image features that
are spatially aligned to the input image. The image features are then sliced
along the direction of the text to generate a fixed number of “timesteps” and
sequentially fed into an LSTM.

The CNN used to generate image features was identical to the residual network
described in Section III-A1 (ResNet34, Figure 2-a). In order to account for
varying sizes of the input image (e.g., lines that contain only one word
compared to lines that contain seven words), multiple downsamples of the image
features are provided (Figure 2-b, identical to the downsampler in the SSD
used in Section III-A2). The image features and downsampled image features
were each fed into separate biLSTMs. The outputs of the biLSTMs were
concatenated along the time dimension and decoded into a $N\times M$ array
where $N$ is the maximum length of the sequence and $M$ is the number of
unique characters (Figure 2-c). This array is fed into the language model
denoiser.

Figure 2: Handwriting recognition CNN-biLSTM.

####  III-B2 Language model denoiser

The $N\times M$ output of the CNN-biLSTM needs to be transformed into the
output string. As the output contains $N\times M$ probabilities corresponding
to each character of the sequence, a naive solution (greedy solution) is to
take the maximum probability (argmax) of each of the $N$ slices and collapse
the characters using the CTC collapsing function. Inspired by [7], a beam
search approach can alleviate such issues by combining multiple decoding paths
to generate candidate strings. A language model can be included in the beam
search decoding to weigh the proposals based on their likeliness. The beam
search approach required substantially more computational power as our early
experiments, revealed that there was up to $11\times$ computational time
increase compared to the greedy solution).

In this paper, a language denoiser network was developed. Given a noisy input
string, the network denoises the string in a sequence-to-sequence
configuration. A previous approach [17] encodes the noisy input at the
character level and decode the clean output at the word level to ensure that
the output is only composed of in-vocabulary words. This has proved relatively
effective however falls apart for out of vocabulary words like names and
places. To circumvent this issue, a character to character encoding / decoding
scheme based on the Transformer architecture was used [18]. The denoiser was
trained on sentences from an external database of public domain novels [19].
Characters are randomly inserted and deleted from the sentences (with a
uniform distribution). Also, characters are replaced with visually similar
counterparts (e.g., ‘d’ can be replaced with ‘c’ and ‘l’) in an attempt to
model the real noisy distribution of the handwriting recognition model. The
generated noisy sentences are used to predict its original counterpart. During
inference, the output of the trained denoiser is fed into a beam search
algorithm to generate candidate strings. We make use of the following
heuristics to rank them:

  1. 1.

Pick the candidate strings with the highest proportion of in-vocabulary words.

  2. 2.

Pick the candidate strings with the lowest Levenshtein distance.

  3. 3.

Pick the candidate string with the lowest perplexity score using an off-the-
shelf pre-trained language model.

##  IV Experiment

###  IV-A Evaluation

The system was evaluated with the IAM dataset [15]. The IAM dataset contains
1539 pages of scanned documents. Each scanned document contains printed text
and 657 writers were asked to write the contents of the printed text in the
space provided. The dataset was split into train and test data, where the test
dataset includes validation 1, validation 2, and test data designated by the
authors of the dataset.

We both evaluated the system qualitatively by visually evaluating the
transcription of examples and quantitatively by computing the character error
rate (CER). Furthermore, we conducted memory and timing comparative analysis.

The CER was calculated with SCLITE [20] and the effects of the following
components were evaluated: 1) no line heuristics, 2) no language model (argmax
algorithm), 3) with beam search [7], and 4) with the denoiser described in
Section III-B2. The predicted and actual text were aligned and the average CER
was calculated line-by-line. Our method was compared to similar methods
presented in [14, 13, 12].

###  IV-B Training details

The networks were developed with Apache’s MXNet deep learning framework [21].
The networks for each component (passage identification, line segmentation,
handwriting recognition, and language model denoising) were trained separately
and the Adam optimiser was used for all the networks [22]. Data augmentation
including random translation, shearing, and occlusions were performed.
However, many typical data augmentation are not applicable to this application
(e.g., flipping and random cropping). In the word/line object recognition
component, to circumvent this issue, lines or words were randomly blanked out.
Details of the implementation can be seen here
(https://github.com/awslabs/handwritten-text-recognition-for-apache-mxnet).

##  V Results

Figure 3: Qualitative results: full page to line images. a) paragraph
segmentation, b) word segmentation, c) word to line conversion. The pipeline
of the images goes from top to bottom within each column.

Figure 3 shows actual results of paragraph segmentation and word to line
segmentation. We can observe that the paragraph segmentation algorithm mostly
predicts the bounding boxes of the handwriting component successfully,
however, the third column presents a failure case where the last line is not
encompassed by the predicted bounding box. Given an image containing only
handwritten text, the word detection algorithm can detect tight bounding boxes
for each word. However, as mentioned in Section III-A2, there are several
short words (typically with words $<$ 3 characters) that are not detected.
Despite the missing words, we can observe in Figure 3-c that all the lines
were successfully detected.

Figure 4: Handwriting recognition and language modelling. Four line images are
displayed and under each line image contains the predicted string ([AM]:
greedy argmax algorithm with no language modelling, [BS]: beam search [7],
[D]: our denoiser Section III-B2)

Figure 4 presents the selected examples to show differences between the
language model component of the described system. First, we can observe that
the greedy algorithm ([AM]) performs reasonably well and the beam search
([BS]) algorithm does not dramatically improve the results. In a), we can see
that the word “noused” was converted into “roused”, which may be based on the
preceding word “head” and the visual similarity of ‘n’ and ‘r’. In b), the
handwriting looks like “beclared” but the denoiser replaced ‘b’ with ‘d’ based
on the learnt language model. In c), the ‘t’ in “desterted” was deleted also
based on language modelling. In d), none of the algorithms were successful to
correct the sentences, and the denoiser worsened the CER.

The CER presented in Table I suggests that line heuristics dramatically
improved handwriting recognition. Qualitatively evaluating the results suggest
that the line heuristics algorithm eliminated incorrectly identified lines
that caused large disparities when aligning the predicted and correct text.
The denoiser achieved a 1.4 CER decrease compared to the greedy argmax
algorithm and beam search algorithm. When compared to previous works on
recognising cropped images (i.e., feeding a _cropped_ image containing only
the handwritten portion compared to the _full_ page with printed and
handwritten text, as indicated by Seg. in Table I), our method outperforms
Bluche [13]. However, methods described in Bluche [12] and Wigington [14] had
lower CER compared to our method.

TABLE I: CER results Method | Seg. | CER
---|---|---
Greedy argmax without line heuristics (ours) | Full | 28.3
Greedy argmax algorithm (ours) | Full | 9.90
Beam search algorithm (ours) | Full | 9.90
Denoiser (ours) | Full | 8.5
Bluche [12] | Cropped | 7.9
Bluche [13] | Cropped | 16.2
Wigington [14] | Cropped | 6.4
TABLE II: Memory and timing requirements Method | Time taken (sec) | Memory (gb)
---|---|---
Ours | 0.389 $\pm$ 0.0888 | 2.6
Bluche [12] | 21.2 $\pm$ 5.6 | N/a
Bluche [13] | 0.62 $\pm$ 0.14 | N/a
Wigington [14] | 0.546 $\pm$ 0.447 | 9.7

Table II presents the memory and timing requirements for our memory compared
to existing methods. When comparing the mean time taken to run an image, our
method requires approximately $1.5\times$ less time compared to [14] and [13].
Our method also utilises substantially less memory (approximately $3.5\times$
less memory) compared to [14] (unfortunately, the memory requirements for [12]
and [13] could not be attained). Since our memory usage is substantially
smaller, it is possible to run multiple images at the same time; effectively
reducing the time required by a third.

##  VI Conclusion

In this paper, we presented a full page offline handwritten text recognition
framework. This framework consists of a pipeline where the handwritten text is
localised (text localisation) followed by converting images of words into
strings (text recognition). Our method achieved a CER of 8.50. The main
advantage of the framework introduced is the reduced computational costs
compared to existing methods. For a tradeoff of CER$\approx$2 comparing to
[14], the throughput could be effectively $4.2\times$ when using a similar
amount of memory.

In conclusion, the framework that we presented is a computationally cheap
alternative to performing full page offline handwritten text recognition. The
results in this paper demonstrate the potential of this framework and future
work can investigate different components of the pipeline for improved
results.

## Acknowledgement

Thank you Simon Corston-Oliver, Vishaal Kapoor, Sergey Sokolov, Soji Adeshina,
Martin Klissarov, and Thom Lane for their helpful feedback for this project.

## References

  * [1] M. Yousef, K. F. Hussain, and U. S. Mohammed, “Accurate, data-efficient, unconstrained text recognition with convolutional neural networks,” 2018.
  * [2] G. Renton, Y. Soullard, C. Chatelain, S. Adam, C. Kermorvant, and T. Paquet, “Fully convolutional network with dilated convolutions for handwritten text line segmentation,” _International Journal on Document Analysis and Recognition (IJDAR)_ , pp. 1–10, 2018.
  * [3] T. Grüning, R. Labahn, M. Diem, F. Kleber, and S. Fiel, “Read-bad: A new dataset and evaluation scheme for baseline detection in archival documents,” in _2018 13th IAPR International Workshop on Document Analysis Systems (DAS)_. IEEE, 2018, pp. 351–356.
  * [4] B. Moysset, C. Kermorvant, C. Wolf, and J. Louradour, “Paragraph text segmentation into lines with recurrent neural networks,” in _Document Analysis and Recognition (ICDAR), 2015 13th International Conference on_. IEEE, 2015, pp. 456–460.
  * [5] W. Liu, D. Anguelov, D. Erhan, C. Szegedy, S. Reed, C.-Y. Fu, and A. C. Berg, “Ssd: Single shot multibox detector,” in _European conference on computer vision_. Springer, 2016, pp. 21–37.
  * [6] Y. LeCun, B. E. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. E. Hubbard, and L. D. Jackel, “Handwritten digit recognition with a back-propagation network,” in _Advances in neural information processing systems_ , 1990, pp. 396–404.
  * [7] A. Graves and J. Schmidhuber, “Offline handwriting recognition with multidimensional recurrent neural networks,” in _Advances in neural information processing systems_ , 2009, pp. 545–552.
  * [8] J. I. Toledo, S. Dey, A. Fornés, and J. Lladós, “Handwriting recognition by attribute embedding and recurrent neural networks,” in _Document Analysis and Recognition (ICDAR), 2017 14th IAPR International Conference on_ , vol. 1. IEEE, 2017, pp. 1038–1043.
  * [9] V. Pham, T. Bluche, C. Kermorvant, and J. Louradour, “Dropout improves recurrent neural networks for handwriting recognition,” in _Frontiers in Handwriting Recognition (ICFHR), 2014 14th International Conference on_. IEEE, 2014, pp. 285–290.
  * [10] H. Ding, K. Chen, Y. Yuan, M. Cai, L. Sun, S. Liang, and Q. Huo, “A compact cnn-dblstm based character model for offline handwriting recognition with tucker decomposition,” in _Document Analysis and Recognition (ICDAR), 2017 14th IAPR International Conference on_ , vol. 1. IEEE, 2017, pp. 507–512.
  * [11] J. Puigcerver, “Are multidimensional recurrent layers really necessary for handwritten text recognition?” in _Document Analysis and Recognition (ICDAR), 2017 14th IAPR International Conference on_ , vol. 1. IEEE, 2017, pp. 67–72.
  * [12] T. Bluche, “Joint line segmentation and transcription for end-to-end handwritten paragraph recognition,” in _Advances in Neural Information Processing Systems_ , 2016, pp. 838–846.
  * [13] T. Bluche, J. Louradour, and R. Messina, “Scan, attend and read: End-to-end handwritten paragraph recognition with mdlstm attention,” in _Document Analysis and Recognition (ICDAR), 2017 14th IAPR International Conference on_ , vol. 1. IEEE, 2017, pp. 1050–1055.
  * [14] C. Wigington, C. Tensmeyer, B. Davis, W. Barrett, B. Price, and S. Cohen, “Start, follow, read: End-to-end full-page handwriting recognition,” in _Proceedings of the European Conference on Computer Vision (ECCV)_ , 2018, pp. 367–383.
  * [15] U.-V. Marti and H. Bunke, “The iam-database: an english sentence database for offline handwriting recognition,” _International Journal on Document Analysis and Recognition_ , vol. 5, no. 1, pp. 39–46, 2002.
  * [16] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [17] S. Ghosh and P. O. Kristensson, “Neural networks for text correction and completion in keyboard decoding,” _arXiv preprint arXiv:1709.06429_ , 2017\.
  * [18] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” in _Advances in Neural Information Processing Systems_ , 2017, pp. 5998–6008.
  * [19] “Public domain novels http://www.textfiles.com/etext/fiction/.”
  * [20] J. Fiscus, “Sclite scoring package version 1.5,” _US National Institute of Standard Technology (NIST), URL http://www. itl. nist. gov/iaui/894.01/tools_ , 1998.
  * [21] T. Chen, M. Li, Y. Li, M. Lin, N. Wang, M. Wang, T. Xiao, B. Xu, C. Zhang, and Z. Zhang, “Mxnet: A flexible and efficient machine learning library for heterogeneous distributed systems,” _arXiv preprint arXiv:1512.01274_ , 2015\.
  * [22] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.


# Morphological evolution of disk galaxies and their concentration, asymmetry
and clumpiness (CAS) properties in simulations across Toomre’s $Q$ parameter

T. Chantavat  The Institute for Fundamental Study, Naresuan University,
Phitsanulok, 65000, Thailand  S. Yuma P. Malelohit  Department of Physics,
Faculty of Science, Mahidol University, Bangkok 10400, Thailand  T.
Worrakitpoonpon  Institute of Science, Suranaree University of Technology,
Nakhon Ratchasima 30000, Thailand T. Worrakitpoonpon<EMAIL_ADDRESS>
(Received ¡date¿ / accepted ¡date¿; )

###### Abstract

We investigate the morphological and structural evolutions of disk galaxies in
simulations for a wide range of Toomre’s $Q$ parameter. In addition to the
inspection of conventional bar modes, we compute the concentration, asymmetry
and clumpiness (CAS) parameters to enlarge the understanding of the galaxy
evolution. These parameters are widely employed to analyze the light
distribution of the observed galaxies, but the adaptation to numerical
simulations is not much considered. While the bar formation takes place in a
considerable range of $Q$ around $1$, barred galaxies originating from $Q>1$
and $Q<1$ disks yield the CAS values that differ significantly. Disks starting
with $Q<1$ develop clumps due to local gravitational instabilities along with
the bar and these clumps play a central role in enhancing the CAS values. That
process is absent in $Q>1$ counterparts in which the evolution is dominated by
linearly unstable two-armed modes that lead to lower CAS values. Likewise,
unbarred galaxies that are obtainable from disks with $Q$ far below and far
above $1$ exhibit greatly different CAS magnitudes. It turns out that the CAS
parameters can serve as indicators of the initial kinematical state and the
evolution history of a disk of any morphology. In addition, we find an
alternative mechanism of the formation of the lopsided barred galaxy when
$Q\lesssim 1$. Bars that evolve in the midst of the clumps can spontaneously
become lopsided at the end.

Galaxy dynamics(591) — Barred spiral galaxies(136) — Galaxy formation(595) —
N-body simulations(1083)

##  1 Introduction

Disk galaxies constitute an important fraction in the Hubble sequence and
their morphologies manifest great diversities such as spiral, barred-spiral,
and lenticular shapes. However, a full understanding of their origins remains
to be resolved. In the standard theoretical framework, sub-structures harbored
in the disk components are considered as the amplified perturbations of
various forms. For instance, the galactic bar is the remnant of the preceding
linearly unstable two-armed modes (Kalnajs, 1971; Contopoulos &
Papayannopoulos, 1980; Vauterin & Dejonghe, 1996; Jalali & Hunter, 2005). On
the other hand, the spiral structures can be achieved by the spiral density
wave (Lin & Shu, 1964), the sheared gravitationally unstable medium (Goldreich
& Lynden-Bell, 1965), or the sheared synchronized epicyclic motion (Julian &
Toomre, 1966; Toomre, 1981; De Rijcke et al., 2019).

Despite the well-formulated theoretical frameworks, numerical simulations
unveiled much intricacy as it was found that the final disk morphologies
depended strongly on the disk and halo intrinsic properties. A number of
simulations demonstrated that the Toomre’s $Q$ parameter (Toomre, 1964) could
determine how a circular disk morphologically evolved (Hohl, 1971;
Athanassoula & Sellwood, 1986; Curir et al., 2006; Bekki, 2023;
Worrakitpoonpon, 2023). Otherwise, the disk random kinetic energy was proved
usable as an indicator of the fate of a disk (Ostriker & Peebles, 1973;
Efstathiou et al., 1982; Athanassoula, 2008; Romeo et al., 2023). In a larger
scope, the concern of the disk stability was not only limited to the disk
initial kinematical properties, but the halo physical and kinematical
properties were also found important to the subsequent disk evolution. Many
works reported the correlation between the halo/bulge mass concentration and
the stability against bar modes (Sellwood, 1980; Sellwood & Evans, 2001;
Athanassoula & Misiriotis, 2002; Shen & Sellwood, 2004; Athanassoula et al.,
2005; Sheth et al., 2008; Sellwood, 2016; Saha & Elmegreen, 2018; Kataria &
Das, 2018, 2019; Kataria et al., 2020; Jang & Kim, 2023). The halo spin was
also found to be influential to the disk dynamical evolution (Saha & Naab,
2013; Long et al., 2014; Collier et al., 2019; Kataria & Shen, 2022; Li et
al., 2023; Joshi & Widrow, 2024; Chiba & Kataria, 2024). On the other hand,
the formation of the spiral structure rather relied on local effects such as
the local instabilities and the disk shearing (D’Onghia et al., 2013; Fujii et
al., 2018; Michikoshi & Kokubo, 2020). The cold disk environment was the
favorable condition to form spiral arms (Evans & Read, 1998; Zakharova et al.,
2023), whereas the disk thickness tended to suppress that process (Ghosh &
Jog, 2018; Bauer & Widrow, 2019). Unlike the bar component, the formation and
evolution of the spiral structures were less dependent on the choices of halo
parameters (Sellwood, 2021).

In the observational aspect, morphological studies of galaxies started from
basic visual classification to the classification via physical properties of
the galaxies. The classic Hubble sequence classified galaxies into elliptical,
spiral, and irregular galaxies based on the apparent stellar distribution of
galaxies (Hubble, 1926). The visual inspection, however, has limitations for
distant galaxies, as their surface brightness drops due to the cosmological
dimming effect. Feature like spiral arms may get too faint to be detected at
high redshifts (e.g., Papaderos et al., 2023). In fact, distant star-forming
galaxies do not share the same structure or morphology as the local spiral
galaxies (Yuma et al., 2011, 2012). Star-forming galaxies started to show a
similar structure as those in the local universe at $z<0.85$ (Takeuchi et al.,
2015).

Another approach is to fit the azimuthally averaged surface brightness profile
with a predesignate formula such as de Vaucouleurs, exponential, and Sérsic
profiles (de Vaucouleurs, 1948; Sérsic, 1963). This method is efficiently used
to distinguish an elliptical from the spiral galaxy in the local universe
(e.g., Fischer et al., 2019; Domínguez Sánchez et al., 2022) up to very high
redshifts of $z=16$ (Ono et al., 2023). However, bulge-disc decomposition by
fitting two components of Sérsic and exponential profiles to the galaxies is
sometimes physically unrealistic with the Sérsic index of the bulge as high as
$n\sim 8$ (Fischer et al., 2019).

Non-parametric measurement of light distribution in a galaxy has been
intensively introduced to diminish the effect of model or assumption (e.g.,
Bershady et al., 2000; Conselice et al., 2000; Conselice, 2003, 2014).
Concentration ($C$), asymmetry ($A$), and clumpiness ($S$) indices or so
called CAS system are among the most common non-parametric methods used to
study the galaxy structure. Besides their efficiency in classifying elliptical
and spiral galaxies, the $CAS$ parameters are related to the evolution of
galaxies (Conselice, 2003). The concentration index is correlated with the
bulge-to-total ratio and stellar mass of the galaxies, while the clumpiness
showed a strong correlation with H$\alpha$ equivalent width, which is
indicative of the star-forming activity. The CAS system has recently been used
to study the galaxy structure and its evolution at high redshifts with data
obtained with the James Webb Space Telescope (Kartaltepe et al., 2023).
Although there are various studies using the CAS parameters in numerical
simulations, they mainly focused on the merging process of the galaxies (e.g.,
Conselice, 2006; Lotz et al., 2008, 2010a, 2010b).

In this work, we proceed on the investigation of the disk morphological
evolution starting from various Toomre’s $Q$ parameter. Apart from the
inspection of the conventional bar modes, we enlarge our scope to the CAS
properties and how they can be relevant to observations. These parameters are
widely employed to examine the structural properties of observed galaxies, and
we will adopt them to analyze the numerical simulations in this study. The
article is organized as follows. First, Sec. 2 describes the numerical models
and accuracy controls. Next in Sec. 3, we introduce the parameters used in
this work, including those for probing the global non-axisymmetric features
and the CAS indices. Then, Sec. 4 reports the numerical results and
discussions are provided therein. Finally, Sec. 5 concludes this study.

##  2 Numerical simulations and accuracy controls

The self-gravitating $N$-body simulations are handled by GADGET-2 (Springel et
al., 2001; Springel, 2005). The density profile of a disk of particles follows
the exponential profile with vertical distribution that reads

$\rho_{d}(r,z)=\frac{M_{d}}{4\pi
R_{0}^{2}z_{0}}e^{-r/R_{0}}\text{sech}^{2}{\bigg{(}\frac{z}{z_{0}}\bigg{)}}$
(1)

where $M_{d}$ is the disk mass, $R_{0}$ is the disk scale radius, and $z_{0}$
is the disk scale height. We choose $M_{d}=10^{11}\ M_{\odot}$, $R_{0}=5\
\text{kpc}$, and $z_{0}=0.2\ \text{kpc}$. The disk is radially and vertically
truncated at $5R_{0}$ and $5z_{0}$, respectively. We principally investigate
the disks of $5\times 10^{6}$ particles, unless otherwise specified. To
properly imitate the disk evolution in a dark matter halo, the disk is put in
static spherical Hernquist potential given by

$\Phi_{h}(r)=-\frac{GM_{h}}{r+r_{h}}$ (2)

where $M_{h}$ and $r_{h}$ are the halo mass and the halo scale radius, which
are fixed to $2.5\times 10^{12}\ M_{\odot}$ and $75\ \text{kpc}$,
respectively. The radial $Q$ profile corresponds to the ratio of the radial
velocity dispersion $\sigma_{r}$ to the minimum required value for the local
stability according to the Toomre’s criterion (Toomre, 1964), i.e.,

$Q=\frac{\sigma_{r}\kappa}{3.36G\Sigma}$ (3)

where $\kappa$ is the epicyclic frequency calculated from the composite disk-
halo potential $\Phi_{tot}$. In the expression (3), $\Sigma$ is the disk
radial surface density which is proportional to $e^{-r/R_{0}}$. To construct a
disk of particles in dynamical equilibrium with a static halo potential, we
adopt the prescriptions of Hernquist (1993) for the disk velocity structures.
The squared radial velocity dispersion is proportional to the disk surface
density, namely

$\sigma_{r}^{2}\propto e^{-r/R_{0}}.$ (4)

The constant of proportionality in Eq. (4) is adjusted so that the $Q$ value
calculated by Eq. (3) at $2R_{0}$, i.e., the reference radius, is equal to a
chosen $Q$ value, and this $Q$ is the representative $Q$ for a case. In other
words, a specific $Q$ can be obtained by properly adjusting the constant of
proportionality in Eq. (4). The tangential velocity dispersion
$\sigma_{\theta}$ is obtained by the relation

$\sigma_{\theta}^{2}=\frac{\kappa^{2}}{4\Omega^{2}}\sigma_{r}^{2}$ (5)

where $\Omega$ is the angular frequency of circular orbit calculated from
$\Phi_{tot}$. The vertical velocity dispersion $\sigma_{z}$ relates to
$\Sigma$ as

$\sigma_{z}^{2}=\pi Gz_{0}\Sigma.$ (6)

By these choices, we have $\sigma_{z}^{2}\propto\sigma_{r}^{2}$. The mean
tangential velocity $\bar{v}_{\theta}$ as a function of radius is obtained by
the axisymmetric Jeans equation as follows

$\bar{v}_{\theta}^{2}=r\frac{d\Phi_{tot}}{dr}+\frac{r}{\Sigma}\frac{d(\sigma_{r}^{2}\Sigma)}{dr}+\sigma_{r}^{2}-\sigma_{\theta}^{2},$
(7)

whereas the mean radial and vertical velocities (or $\bar{v}_{r}$ and
$\bar{v}_{z}$, respectively) are initially zero. The three random velocity
components are drawn from the cut-off Gaussian distribution with corresponding
local velocity dispersion ellipsoid.

Calculations of spline-softened mutual forces for disk particles are
facilitated by the tree code. We adjust the opening angle to $0.7$ and the
softening length to $5\ \text{pc}$ for all particles. The integration time
step is controlled to be not greater than $0.2\ \text{Myr}$. The accuracy is
such that the deviations of both the total energy and the disk angular
momentum at the end of simulation, which is specified to be $9.6\ \text{Gyr}$,
from the initial values are not greater than $0.1\%$. We examine the cases
where $Q=0.65,0.8,0.95,1.1,1.2$ and $1.5$.

##  3 Parameters

###  3.1 Bar parameter

Non-axisymmetric features of a disk can be evaluated by the $m$-mode Fourier
amplitudes as a function of radius $\tilde{A}_{m}(r)$ defined as

$\tilde{A}_{m}(r)=\frac{\sqrt{a^{2}_{m}+b^{2}_{m}}}{A_{0}}$ (8)

where $a_{m}$ and $b_{m}$ are the $m$-mode Fourier coefficients calculated
from particles inside the annulus of radius $r$ and $A_{0}$ is the
corresponding $m=0$ amplitudes. The $m$-mode strength $A_{m}$ is designated by
the maximum $\tilde{A}_{m}$ within $r_{max}$, or

$A_{m}\equiv\max_{r<r_{max}}[\tilde{A}_{m}]$ (9)

where we fix $r_{max}$ to $10\ \text{kpc}$ for all calculations. We employ
$A_{2}$ as the bar strength.

###  3.2 CAS parameters

CAS system is a non-parametric measurement of light distribution of a galaxy.
In the simulations, we obtain the surface density instead of the surface
brightness. So we assume constant mass-to-light ratio across the entire
simulated galaxy to translate the surface density into the surface brightness.
It consists of 3 parameters: concentration ($C$), asymmetry ($A$), and
clumpiness ($S$). We adopt the concentration index originally defined by
Bershady et al. (2000) as

$C=5\log_{10}\left(\frac{r_{80}}{r_{20}}\right),$ (10)

where $r_{20}$ and $r_{80}$ are the radii governing 20% and 80% of the growth
curve within 1.5 times the Petrosian radius at $r(\eta=0.2)$, respectively.

The asymmetry index shows the fraction of galaxy component that is not
symmetric, which was originally defined as

$A=\min\left(\frac{\sum\lvert
I_{0}-I_{180}\rvert}{\sum\lvert{I_{0}}\rvert}\right)-\min\left(\frac{\sum\lvert
B_{0}-B_{180}\rvert}{\sum\lvert I_{0}\rvert}\right),$ (11)

where $I_{0}$ and $B_{0}$ are the original image and background area close to
the galaxy, respectively (Conselice et al., 2000). $I_{180}$ is the image
after rotating $180^{\circ}$ from the center of the galaxy in the line of
sight. The center is varied until we obtain the minimum value of the asymmetry
index. To reduce the effect of noise, Conselice et al. (2000) consider the
background region near the object by randomly selecting its center, rotating
the background by $180^{\circ}$, and calculating the rightmost term of
Equation 11. However, we do not have noise in the simulation. So we simply use
only the first term of Equation 11 to estimate an asymmetry index. That is

$A=\min\left(\frac{\sum\lvert
I_{0}-I_{180}\rvert}{\sum\lvert{I_{0}}\rvert}\right).$ (12)

Lastly, the clumpiness index is defined as the ratio of light in high-
frequency structure to the total light of the galaxy (Conselice, 2003). It can
be written as

$S=10\times\left[\left(\frac{\sum(I_{xy}-I^{\sigma}_{xy})}{\sum
I_{xy}}\right)-\left(\frac{\sum(B_{xy}-B^{\sigma}_{xy})}{\sum
I_{xy}}\right)\right],$ (13)

where $I_{xy}$ is is the original image, while $I^{\sigma}_{xy}$ is the
smoothed image. Likewise, $B_{xy}$ and $B^{\sigma}_{xy}$ are, respectively,
the original and blurred background with equal area to the size of the galaxy
(Conselice, 2014). The original image is smoothed with a two-dimensional
Gaussian kernel with the size of $\sigma=0.3\times r(\eta=0.2)$ (Conselice,
2003).

##  4 Numerical results

###  4.1 Evolution of bi-symmetric modes

Figure 1: Time evolution of $A_{2}$ for disks with different $Q$. Figure 2:
Disk surface density map for $Q=0.65$ (top row), $0.95$ (middle row) and $1.1$
(bottom row) at indicated times.

First of all, we inspect the time evolution of $A_{2}$ for disks with
different $Q$ in Fig. 1. It turns out that the development of the bar modes
takes place when $Q\leq 1.2$. Despite that $Q$ is less than $1$, the $m=2$
modes can still grow in the midst of the local instabilities. A disk with
lower $Q$ tends to form a bar more rapidly, in accordance with past studies
(Athanassoula & Sellwood, 1986; Hozumi, 2022). On the other hand, the value of
$Q=1.5$ is sufficient to suppress the bar instability. We further examine the
disk configuration for different $Q$ in Fig. 2 to verify the growth of the bar
modes. In the case where $Q=1.1$, we observe the buildup of weak multi-armed
modes at $5.52\ \text{Gyr}$ which are overseen by the swing amplification
(Julian & Toomre, 1966; Toomre, 1981). That mechanism engenders the spiral
arms from the synchronized epicyclic motion in combination with the disk
shearing. After $8\ \text{Gyr}$, the two-armed modes exert the dominance over
the multi-arm modes as the configuration takes the form of the two-armed
barred structure, and it remains in shape until the end. When $Q=0.95$, which
is slightly below the local stability threshold, we observe a more rapid
formation of the bar amidst the local instabilities as the bar is spotted at
$5.04\ \text{Gyr}$. Compared to the previous case, we remark clumpy sub-arms
around the central bar. As time progresses, the shape of the bar and spiral
arms deviates from bi-symmetry and they are surrounded by long-lasting clumps
which remain until the end. The clumps emerge from the gravitationally
unstable fluctuations, so they are distributed randomly in the disk
environment. We speculate that the breaking of bi-symmetry is attributed to
those local instabilities and the bar evolution in the midst of such
instabilities can be described as follows. While the two-armed modes are
growing, the unstable local over-densities progress to clumps in random places
in parallel. Afterward, the bi-symmetric component is deformed by the
interactions with clumps as observed in the last snapshot.

The case where $Q=0.65$ is fully governed by local gravitational instabilities
which lead to a different morphological evolution from the other cases. The
mass concentrations develop quickly, supposedly in the time scale of the local
free-fall time, before they are sheared by disk differential rotation. As a
consequence, the multi-armed spiral pattern is formed rapidly within $\sim
2.6\ \text{Gyr}$ which is significantly shorter than the time scale of the bar
formation. The spiral arms appear more clumpy than those arising from the
unstable two-armed modes or the swing amplification when $Q>1$. These arms are
short-lived as they are found fragmented into clumps at $4.92\ \text{Gyr}$.
The clumps are still spotted around the concentrated disk center at the end of
simulation. The fragile spiral structure is another indication of the
dominance of the local instability over the spiral mode instability, which
plays an important role in other cases. The plot of the density map informs
that, referring to Fig. 1, the growth of $A_{2}$ does not signify the
development of the bar-like structure. The developed $m=2$ mode amplitudes are
caused by the anisotropically distributed clumps. We note that the number of
the clumps in the end reduces from the number in the snapshot before which is
possibly due to either the merger between them or the falling to the disk
center. The decrease of the clump number is in accordance with the decline of
$A_{2}$ with time near the end. In addition, we note the breaking of disk
circular symmetry by clumps, which can be explained by the same arguments for
the breaking of the bi-symmetry of the barred-spiral structure in the $Q=0.95$
case.

The bar formation in a disk with $Q\geq 1$ conforms with past mainstream
studies but only a few studies visited the regime where $Q$ was below $1$
(Athanassoula, 2003; Debattista et al., 2006; Worrakitpoonpon, 2023). We have
some important remarks for the latter regime that deserve attention. Our
simulation of $Q=0.95$ disk demonstrates that, firstly, the disk is populated
by long-lasting clumps that persist until the end. Furthermore, the overall
symmetry is strongly broken compared with the $Q=1.1$ counterpart. These
anomalies are attributed to the initial local instability. About the evolution
pattern of the $Q=0.65$ disk, there was a conjecture for a disk of comparable
$Q$ to be the progenitor of lenticular (or S0) galaxies which do not host the
bar and the spiral arms (Saha & Cortesi, 2018). That hypothesis is partly
plausible for us as our simulations suggest that the S0 galaxy candidate can
also emerge from $Q=1.5$, i.e., when the disk is stable against all types of
perturbations. The differentiation between the practically similar
configurations originating from the locally stable and unstable disks needs a
suitable indicator. These related properties will be addressed in Sec. 4.2.

###  4.2 Evolutions of CAS parameters

In continuity with the clumpy and asymmetric disks for $Q<1$ observed in Sec.
4.1, we continue with a more systematic evaluation of such properties. Shown
in Fig. 3 are the time evolutions of CAS parameters for different $Q$. From
the plots, we are able to classify the evolution schemes into 3 sub-families:
the $Q<1$ disks, the $Q>1$ disks that form a bar, and the disks that are bar-
stable. Those with $Q<1$ develop highly concentrated, asymmetric, and clumpy
features as remarked by large increases of the three parameters. With a closer
look, the clumpiness $S$ increases most rapidly because of the local
gravitational instabilities residing in the initial state. It remains at peak
value for $\sim 2\ \text{Gyr}$ before it decays sharply. Afterward, it keeps
decaying slowly until the end but it is still higher than the other families.
When $Q$ is closer to $1$, the $S$ index develops more slowly and reaches a
lower peak. This is because the disk is closer to the local stability
threshold. The decline of $S$ after the peak can be explained, in parallel
with the evolution in Fig. 2, that the number of clumps reduces with time. We
recall that the case where $Q=0.65$ undergoes the formation of clumpy multi-
arm pattern in midway, in accordance with the time at which $S$ is at peak. In
the end, a high level of clumpiness still retains, although the clumpy spiral
arms dissolve. Considering the concentration $C$ and asymmetry $A$, they reach
the peaks well after the $S$ peak. This indicates that the developments of $C$
and $A$ features are the secondary processes induced by the clumpiness that
develops earlier. This order of increases is in accordance with earlier
speculation that anisotropically distributed clumps are the source of
asymmetry. Here, we validate that the enhanced concentration is also the
consequence. Considering the evolution of $A$ index, the values at peak for
$Q=0.65$ and $0.8$ are comparable and they are significantly higher than that
for $Q=0.95$. After the peaks, the $A$ index for all three subcritical cases
decays down. The rapid decay is explained that some clumps fall to the disk
center, resulting in the enhancement of $C$ and the decline of $S$
accordingly. The $C$ parameter, on the other hand, does not show any sign of
decay with time until the end of simulation, which indicates the continuous
accumulation of the density at the center by merger with clumps. In summary,
the CAS evolution of this disk family is dominated by local instabilities
marked by rapid increase of $S$, followed by the increases of $C$ and $A$. It
is worth reminded that the case where $Q$ is far below $1$ does not develop
the bar-like structure, while those starting with $Q$ slightly below $1$ can
develop a bar. These differences cannot be distinguished by using only the CAS
parameters. A direct visual inspection on the disk configuration is required.

Figure 3: Time evolution of CAS parameters from the simulation with $N=5\times
10^{6}$. Different colors in all panels represent different values of $Q$
ranging from $Q=0.65$ to $Q=1.5$. Figure 4: CAS parameters as a function of
$Q$ for $N=5\times 10^{6}$ (top panels) and $N=1\times 10^{6}$ (bottom
panels). Blue, orange, and green solid lines indicate the snapshots at $t=4.8$
Gyr, $t=7.2$ Gyr, and $t=9.6$ Gyr, respectively. Solid circles represent the
values of $Q$ we adopted in each simulation.

Disks that are bar-unstable and start with $Q>1$ exhibit moderate increases of
all CAS parameters and the order of increases differs from those for $Q<1$. It
is the $A$ parameter that increases first due to the growth of two-armed modes
that finally become the barred-spiral structure, and that established
structure induces the growths of $C$ and $S$. In this disk family, the
linearly unstable two-armed modes play a central role unlike the situation in
the previous disk family. Although the local gravitational instabilities are
suppressed, $C$ and $S$ can nevertheless be enhanced by the bar and the two-
armed modes, in line with the concentrated bar and spiral arms observed at the
end in Fig. 2, but they are not at the levels attained in a locally unstable
disk. Close to the end of simulations, the $A$ index decreases for $Q=1.1$,
which marks the onset of the decay of bar modes in coherent with the decrease
of $A_{2}$ in Fig. 1. The case where $Q=1.2$ does not yet exhibit a decay
because the bar is formed later. In the bar-stable family, all parameters
remain close to initial values. This indicates that both local and bar-mode
instabilities are effectively suppressed.

The variation of CAS parameters can otherwise be plotted as a function of $Q$
as shown in Fig. 4 at different times for $N=5\times 10^{6}$ and $10^{6}$. The
investigation of the case with a lower $N$ has a purpose to examine how the
results change by reducing $N$ since the clumps originate from the
gravitationally unstable finite-$N$ fluctuations, whose amplitudes scale with
$1/\sqrt{N}$. In both cases, we observe that $C$ and $S$ tend to decrease as
$Q$ increases at $9.6\ \text{Gyr}$, as expected from a lesser degree of local
instabilities. With a closer look, varying $N$ causes some notable differences
in evolutionary pattern and magnitude. The trivial consequence is that the $S$
parameter for $N=10^{6}$ is averagely higher because of more pronounced local
instabilities. However, the decrease with $Q$ retains for both $N$. For the
$C$ parameter, it monotonically decreases with $Q$ in all time snapshots for
$N=5\times 10^{6}$. When $N$ steps down to $10^{6}$, the value for $Q=0.65$ is
significantly lower than the value for $Q=0.8$ during $4.8-7.2\ \text{Gyr}$,
but the monotonic decrease with $Q$ is eventually achieved at $9.6\
\text{Gyr}$. This can be explained that the developed clumpy feature lasts
longer if $N$ is lower before it falls to the disk center. This speculation is
supported by a slower decay of $S$ index from $4.8-9.6\ \text{Gyr}$ than the
decay rate for $N=5\times 10^{6}$. For the $A$ parameter, the decrease with
$Q$ is more evident in the middle of evolution because of the existing clumps
and the forming barred-spiral structure whereas there is no clear tendency of
the variation with $Q$ at the end. It turns out that using the $A$ index as an
indicator of the initial $Q$ is at best when non-axisymmetric features are at
peak. The variation of $A$ with time and with $Q$ for $N=10^{6}$ is
qualitatively similar to the $N=5\times 10^{6}$ counterpart: the decrease with
$Q$ is obvious at $4.8\ \text{Gyr}$ but the curve becomes flattened in the
range of $Q\leq 1.1$ in the end. In addition, the decay of $S$ results in the
decrease of $A$ as some clumps are falling to the symmetric disk center and
become part of it.

The inspection of CAS parameters provides a new insight in the way that the
origin and evolution history of a disk of any morphology and structural detail
can be specified. In the bar-forming regime with $Q$ ranging from $0.8$ to
$1.2$, disks starting with $Q$ below and above $1$ exhibit distinct
evolutionary tracks of CAS indices and they lead to significantly different
values. More specifically, barred galaxies formed in $Q<1$ disks can be
identified by $C>4.0$ and $S>0.15$, referring to the case of $N=5\times
10^{6}$, and vice versa for bar-forming disks starting with $Q>1$. Fig. 3
suggests that the proposed $C$ border value is applicable in the time frame of
galaxy age from $7-9.6\ \text{Gyr}$. The differentiation of the disk origin by
that proposed $S$ border value is applicable to a wider time frame of galaxy
age: the value of $S=0.15$ can distinguish the initial kinematical conditions
of disk galaxies of age from $5-9.6\ \text{Gyr}$. On the other hand, the $A$
value can be useful when the asymmetric feature is at peak. The clumpy barred-
spiral disk can produce $A$ that is significantly greater than $0.6$ and it
saturates on that level for more than $4\ \text{Gyr}$. This is significantly
greater than the asymmetry caused solely by the barred-spiral pattern, which
can attain the value around $0.6$ at most for a short period of $\sim 0.5\
\text{Gyr}$ before it starts to decay. The $A$ parameter can therefore be
utilized in the following way. The observed value of $A$ greater than $0.6$
can be an evidence of the locally unstable initial disk. Otherwise, the
establishment of the barred-spiral galaxies without local instabilities cannot
yield an $A$ value greater than $0.6$ in most of their lifetime.

Considering unbarred cases, this configuration can be achieved when $Q$ is
either far above or far below $1$. Both cases yield the lenticular-like
morphology at the end but the measured CAS parameters are at remarkably
different levels. The values of $C>4.5$, $A>0.4$, and $S>0.2$ can infer a cold
progenitor of a lenticular disk, which undergoes violent evolution forming
clumpy short-lived spiral arms. Orbiting clumps that remain are another piece
of evidence of cold origin. On the contrary, those with modest CAS magnitudes
undergo smoother evolution with no perturbative structures formed because all
local and global non-axisymmetric modes are suppressed. The $C$, $A$ and $S$
values can be as low as $2.7$, $0.1$ and $0.05$, respectively.

In observations, the concentration index increases with the increasing bulge-
to-total ratio for the local galaxies (Conselice, 2003). The early-typed disk
and elliptical galaxies, which the bulge is already well established, showed
the concentration index of $C>4.0$. Although this criterion is consistent with
the simulated barred galaxies with $Q<1$ in our simulations, the concentration
index alone cannot be efficiently used to indicate the bar structure. The
asymmetry index is famously used to identify the merging system in both
observations and simulations. Lotz et al. (2010a) used the simulations to show
that the asymmetry index could reach $A>0.4$ at the stage of final merger. It
has been later used as the criterion for merging system in the observations
(Kartaltepe et al., 2023). Considering our numerical results, such high level
of the $A$ index can be achieved and maintained for long period in cases with
$Q<1$, whereas that level can be attained momentarily during the peak of the
bar modes when $Q>1$. Our important finding is that this high level of
asymmetry can be attained spontaneously in locally unstable disks. In case of
an isolated galaxy, the clumpy object can also show high asymmetry index, as
it is strongly correlated with the clumpiness (Conselice, 2003). This
correlation is also seen in our simulations (Fig. 3). Early-type spiral and
elliptical galaxies showed very low asymmetry ($A<0.1$) and clumpiness
($S<0.1$), while the late-type spiral or starburst with high star formation
rate show the larger values. The enhanced star formation in galaxies that
exhibit a higher clumpiness index is actually understandable as the formation
of the clumps also promote the local star formation. This observed correlation
emphasizes the possibility of the locally unstable disk to be a progenitor of
clumpy spiral galaxy.

In summary, we find a utility of CAS parameters apart from their correlations
with some present physical properties of galaxies such as the bulge-to-total
ratio and the star formation rate. The CAS magnitudes of an established disk
galaxy of any morphology can determine its origin and evolution history
reasonably well. Although our numerical model consists of purely self-
gravitating particles, whereas a real galaxy consists of many components, the
conclusion of the relation between the CAS values and the evolutionary history
remains applicable. This is because the key factor leading to the difference,
namely the local instability, is overseen by the gravity and the gravity is
the governing force of all components. In addition, it was documented that
observed disk galaxies with estimated $Q$ less than $1$, in which local
instabilities played a central role, tended to be gas-poor (Renaud et al.,
2021; Aditya, 2023). This ascertains the dominance of the gravitational
dynamics, which is our focus here, over the gas dynamics for this $Q<1$
regime. Another limitation is that our study considers the evolution of a disk
in isolation. A galaxy with merging history might yield the elevated $A$
without local instabilities that might not conform with our conclusion.

###  4.3 Radial heating with and without clumps

It has been justified that the evolving non-axisymmetric forces from the bar
and the spiral arms could radially heat the disk environment, i.e., the radial
heating (Jenkins & Binney, 1990; Minchev & Quillen, 2006). Shown in Fig. 5 is
the radial velocity dispersion profile for $Q=0.8$ and $1.1$ at different
times covering the process of the bar formation and its subsequent evolution.
In both cases, we capture the radial heating as the profile raises with time,
but the detailed procedures are different. For $Q=1.1$, the entire profile
heats up smoothly by global bi-symmetric forces. On the contrary, bumps in the
profile of the $Q=0.8$ case during the heating are spotted, most notably at
$3.84\ \text{Gyr}$. The bumps are lowered but still observable at $6.72\
\text{Gyr}$ before the profile is smoothed out at $9.60\ \text{Gyr}$. We note
that the final profile is flatter than the $Q=1.1$ profile at the same time.
These differences are attributed to the additional local heating by clumps
which is absent when $Q>1$, in which case the disk is heated by global bi-
symmetric forces without local effect. The bumps in the profile during the
radial heating and the flatness of the profile at the end can be indicative to
the kinematically cold disk origin, in addition to the elevated CAS values
that we addressed in the previous section.

Figure 5: Radial velocity dispersion profiles at different times for $Q=0.8$
(top panel) and $1.1$ (bottom panel).

The concept of the radial heating by clumps shares similarities with the
framework of the heating by the giant molecular cloud (GMC) proposed by, for
instance, Hänninen & Flynn (2002); Aumer et al. (2016); Fujimoto et al.
(2023). The major difference is, however, that those studies constructed the
GMC in the context of the locally stable disk, i.e., $Q\geq 1$, while we
inspect the locally unstable disk that spontaneously forms clumps. However, we
do not rule out the relevance to that hypothesis as in reality, the disk
kinematical map is not smooth and perfectly axisymmetric as imposed in the
initial conditions for simulations. The GMC can potentially be formed in a
location where local $Q$ is below $1$, although the disk has $Q\geq 1$ on
average.

###  4.4 Lopsided bar and spiral arms in subcritical $Q$ disks

Figure 6: Time evolutions of $A_{1}$ (dotted line) and $A_{2}$ (solid line)
for $Q=0.8$ (top panel), $0.95$ (middle panel) and $1.1$ (bottom panel).
Figure 7: Surface density map of disk center for $Q=0.95$ at $6.96\
\text{Gyr}$.

In continuity with the notice of an apparently lopsided bar for the case where
$Q=0.95$ in Sec. 4.1, here we carry out a more systematic inspection of that
appearance. It is a usual practice to probe the lopsidedness by $A_{1}$ and we
will inspect its evolution in time along with $A_{2}$. Shown in Fig. 6 are
time evolutions of $A_{1}$ and $A_{2}$ for various $Q$. In all cases, the
$m=2$ modes are growing faster and they lead to the barred structure when
those modes are at peak. Afterward, the $m=1$ modes are growing whereas the
$m=2$ counterparts are decaying for cases with $Q=0.8$ and $0.95$ and the
$m=1$ modes become dominant at the end. The development of $m=1$ modes
relative to the growth of $A_{2}$ for $Q=0.95$ is slower than in the $Q=0.8$
case because of less pronounced local effect. The emerging lopsidedness as
speculated by the dominance of $A_{1}$ is validated by Fig. 7 which
illustrates the lopsided barred-spiral disk around the time of $A_{1}$ peak
for $Q=0.95$. The evolution schemes of $A_{1}$ and $A_{2}$ affirm our
hypothesis in Sec. 3.1 that anisotropically distributed clumps cause the
breaking of bi-symmetry. On the other hand, the $m=1$ modes are not much
amplified after the growth of the $m=2$ modes when $Q=1.1$, and the latter
modes remain the dominant modes until the end. This is because of the lack of
the local instabilities.

That the lopsided bar can potentially be established in subcritical $Q$ regime
is an alternative explanation because past studies often related the origin of
lopsided feature to the imposed macroscopic asymmetry of the initial state.
For instance, numerical results of Levine & Sparke (1998); Łokas (2021)
demonstrated that a disk that was initially off-centered from the halo origin
could produce a lopsided bar. Other studies attributed the formation of
lopsided bar to tidal interactions (Yozin & Bekki, 2014; Łokas, 2022; Varela-
Lavin et al., 2023). There were hypotheses that other astrophysical processes
such as the gas accretion (Bournaud et al., 2005; Dupuy et al., 2019) and the
stellar feedback (Manuwal et al., 2022) could generate the lopsidedness. Our
important finding is that the established bi-symmetry can spontaneously be
broken and the disk becomes lopsided, although the initial disk is constructed
from a circularly symmetric profile, if the gravitationally unstable
fluctuations are taken into account. In our framework, the initial asymmetry
and the external forces are not required.

##  5 Conclusion

We examine the evolution of disk galaxies in simulations for various initial
$Q$ values covering those that are below the Toomre’s stability threshold and
those above it. The scope is not only limited to the evolution of bi-symmetric
features, but we also investigate the development of concentration, asymmetry
and clumpiness (CAS) properties in that range of $Q$. The CAS parameters are
widely adopted to discriminate galaxies in observations and we will apply
those parameters to our numerical results.

We find a utility of CAS indices when they are considered along with other
conventional plots such as the bar parameter or the density map as they are
able to identify the kinematical condition of the initial disk and the
evolution history. For instance, the measured CAS values of cases that form a
bar, which cover the initial $Q\in[0.8,1.2]$, are able to distinguish them
into two sub-groups. A disk starting with $Q<1$ develops a bar and clumps in
parallel, yielding high CAS indices at the end, while a $Q>1$ counterpart
develops only a bar with two-armed spiral pattern. In the latter group, the
CAS values are significantly lower. For unbarred cases which are the
candidates of the lenticular galaxy, this morphology can be obtained from an
initial disk with $Q$ either far below or far above $1$. The discrimination of
unbarred disks achieved from those initial states at the two extremities can
also be done using the CAS indices. The former case yields considerably more
elevated CAS values than those in the latter case. High CAS values in far
subcritical $Q$ disks also reflect the clumpy and asymmetric appearance which
is not complicated to be remarked visually.

In addition, the presence of long-lived clumps surrounding the barred-spiral
structure can lead to remarkable outcomes. The first consequence is the
distinct features of the radial velocity dispersion profile caused by the
local radial heating by clumps. It gives rise to the bumps in the profile
while the heating is progressing and the flatter profile at the end than that
without local heating. Secondly, the interactions between clumps and bi-
symmetric components are able to deform the bar to lopsided configuration.
This is an alternative explanation of the formation of the lopsided bar apart
from other scenarios which require an initial seed of asymmetry or external
forces.

This research has received funding support from the NSRF via the Program
Management Unit for Human Resources & Institutional Development, Research and
Innovation (grant number B16F640076). TC is supported by Naresuan University
(NU), and the National Science, Research and Innovation Fund (NSRF) Grant
number R2566B091. TW has funding support from Thailand Science Research and
Innovation (TSRI) via Suranaree University of Technology [grant number
179349]. SY is supported by Mahidol University (Fundamental Fund: fiscal year
2023 by National Science Research and Innovation Fund (NSRF)) and by office of
the Ministry of Higher Education, Science, Research and Innovation through
research grant for new scholars (RGNS63-175). Numerical simulations are
facilitated by HPC resources of Chalawan cluster of the National Astronomical
Research Institute of Thailand.

## References

  * Aditya (2023) Aditya, K. 2023, MNRAS, 522, 2543, doi: 10.1093/mnras/stad1143
  * Athanassoula (2003) Athanassoula, E. 2003, MNRAS, 341, 1179, doi: 10.1046/j.1365-8711.2003.06473.x
  * Athanassoula (2008) —. 2008, MNRAS, 390, L69, doi: 10.1111/j.1745-3933.2008.00541.x
  * Athanassoula et al. (2005) Athanassoula, E., Lambert, J. C., & Dehnen, W. 2005, MNRAS, 363, 496, doi: 10.1111/j.1365-2966.2005.09445.x
  * Athanassoula & Misiriotis (2002) Athanassoula, E., & Misiriotis, A. 2002, MNRAS, 330, 35, doi: 10.1046/j.1365-8711.2002.05028.x
  * Athanassoula & Sellwood (1986) Athanassoula, E., & Sellwood, J. A. 1986, MNRAS, 221, 213, doi: 10.1093/mnras/221.2.213
  * Aumer et al. (2016) Aumer, M., Binney, J., & Schönrich, R. 2016, MNRAS, 462, 1697, doi: 10.1093/mnras/stw1639
  * Bauer & Widrow (2019) Bauer, J. S., & Widrow, L. M. 2019, MNRAS, 486, 523, doi: 10.1093/mnras/stz478
  * Bekki (2023) Bekki, K. 2023, MNRAS, 523, 5823, doi: 10.1093/mnras/stac3097
  * Bershady et al. (2000) Bershady, M. A., Jangren, A., & Conselice, C. J. 2000, AJ, 119, 2645, doi: 10.1086/301386
  * Bournaud et al. (2005) Bournaud, F., Combes, F., Jog, C. J., & Puerari, I. 2005, A&A, 438, 507, doi: 10.1051/0004-6361:20052631
  * Chiba & Kataria (2024) Chiba, R., & Kataria, S. K. 2024, MNRAS, 528, 4115, doi: 10.1093/mnras/stae288
  * Collier et al. (2019) Collier, A., Shlosman, I., & Heller, C. 2019, MNRAS, 488, 5788, doi: 10.1093/mnras/stz2144
  * Conselice (2003) Conselice, C. J. 2003, ApJS, 147, 1, doi: 10.1086/375001
  * Conselice (2006) —. 2006, ApJ, 638, 686, doi: 10.1086/499067
  * Conselice (2014) —. 2014, ARA&A, 52, 291, doi: 10.1146/annurev-astro-081913-040037
  * Conselice et al. (2000) Conselice, C. J., Bershady, M. A., & Jangren, A. 2000, ApJ, 529, 886, doi: 10.1086/308300
  * Contopoulos & Papayannopoulos (1980) Contopoulos, G., & Papayannopoulos, T. 1980, A&A, 92, 33
  * Curir et al. (2006) Curir, A., Mazzei, P., & Murante, G. 2006, A&A, 447, 453, doi: 10.1051/0004-6361:20053418
  * De Rijcke et al. (2019) De Rijcke, S., Fouvry, J.-B., & Pichon, C. 2019, MNRAS, 484, 3198, doi: 10.1093/mnras/stz166
  * de Vaucouleurs (1948) de Vaucouleurs, G. 1948, Annales d’Astrophysique, 11, 247
  * Debattista et al. (2006) Debattista, V. P., Mayer, L., Carollo, C. M., et al. 2006, ApJ, 645, 209, doi: 10.1086/504147
  * Domínguez Sánchez et al. (2022) Domínguez Sánchez, H., Margalef, B., Bernardi, M., & Huertas-Company, M. 2022, MNRAS, 509, 4024, doi: 10.1093/mnras/stab3089
  * D’Onghia et al. (2013) D’Onghia, E., Vogelsberger, M., & Hernquist, L. 2013, ApJ, 766, 34, doi: 10.1088/0004-637X/766/1/34
  * Dupuy et al. (2019) Dupuy, J. L., Heitsch, F., & Cecil, G. 2019, MNRAS, 486, 2020, doi: 10.1093/mnras/stz975
  * Efstathiou et al. (1982) Efstathiou, G., Lake, G., & Negroponte, J. 1982, MNRAS, 199, 1069, doi: 10.1093/mnras/199.4.1069
  * Evans & Read (1998) Evans, N. W., & Read, J. C. A. 1998, MNRAS, 300, 106, doi: 10.1046/j.1365-8711.1998.01864.x
  * Fischer et al. (2019) Fischer, J. L., Domínguez Sánchez, H., & Bernardi, M. 2019, MNRAS, 483, 2057, doi: 10.1093/mnras/sty3135
  * Fujii et al. (2018) Fujii, M. S., Bédorf, J., Baba, J., & Portegies Zwart, S. 2018, MNRAS, 477, 1451, doi: 10.1093/mnras/sty711
  * Fujimoto et al. (2023) Fujimoto, Y., Inutsuka, S.-i., & Baba, J. 2023, MNRAS, 523, 3049, doi: 10.1093/mnras/stad1612
  * Ghosh & Jog (2018) Ghosh, S., & Jog, C. J. 2018, A&A, 617, A47, doi: 10.1051/0004-6361/201832988
  * Goldreich & Lynden-Bell (1965) Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125, doi: 10.1093/mnras/130.2.125
  * Hänninen & Flynn (2002) Hänninen, J., & Flynn, C. 2002, MNRAS, 337, 731, doi: 10.1046/j.1365-8711.2002.05956.x
  * Hernquist (1993) Hernquist, L. 1993, ApJS, 86, 389, doi: 10.1086/191784
  * Hohl (1971) Hohl, F. 1971, ApJ, 168, 343, doi: 10.1086/151091
  * Hozumi (2022) Hozumi, S. 2022, MNRAS, 510, 4394, doi: 10.1093/mnras/stab3704
  * Hubble (1926) Hubble, E. P. 1926, ApJ, 64, 321, doi: 10.1086/143018
  * Jalali & Hunter (2005) Jalali, M. A., & Hunter, C. 2005, ApJ, 630, 804, doi: 10.1086/432370
  * Jang & Kim (2023) Jang, D., & Kim, W.-T. 2023, ApJ, 942, 106, doi: 10.3847/1538-4357/aca7bc
  * Jenkins & Binney (1990) Jenkins, A., & Binney, J. 1990, MNRAS, 245, 305
  * Joshi & Widrow (2024) Joshi, R., & Widrow, L. M. 2024, MNRAS, 527, 7781, doi: 10.1093/mnras/stad3666
  * Julian & Toomre (1966) Julian, W. H., & Toomre, A. 1966, ApJ, 146, 810, doi: 10.1086/148957
  * Kalnajs (1971) Kalnajs, A. J. 1971, ApJ, 166, 275, doi: 10.1086/150957
  * Kartaltepe et al. (2023) Kartaltepe, J. S., Rose, C., Vanderhoof, B. N., et al. 2023, ApJ, 946, L15, doi: 10.3847/2041-8213/acad01
  * Kataria & Das (2018) Kataria, S. K., & Das, M. 2018, MNRAS, 475, 1653, doi: 10.1093/mnras/stx3279
  * Kataria & Das (2019) —. 2019, ApJ, 886, 43, doi: 10.3847/1538-4357/ab48f7
  * Kataria et al. (2020) Kataria, S. K., Das, M., & Barway, S. 2020, A&A, 640, A14, doi: 10.1051/0004-6361/202037527
  * Kataria & Shen (2022) Kataria, S. K., & Shen, J. 2022, ApJ, 940, 175, doi: 10.3847/1538-4357/ac9df1
  * Levine & Sparke (1998) Levine, S. E., & Sparke, L. S. 1998, ApJ, 496, L13, doi: 10.1086/311239
  * Li et al. (2023) Li, X., Shlosman, I., Heller, C., & Pfenniger, D. 2023, MNRAS, doi: 10.1093/mnras/stad2799
  * Lin & Shu (1964) Lin, C. C., & Shu, F. H. 1964, ApJ, 140, 646, doi: 10.1086/147955
  * Łokas (2021) Łokas, E. L. 2021, A&A, 655, A97, doi: 10.1051/0004-6361/202141604
  * Łokas (2022) —. 2022, A&A, 662, A53, doi: 10.1051/0004-6361/202142845
  * Long et al. (2014) Long, S., Shlosman, I., & Heller, C. 2014, ApJ, 783, L18, doi: 10.1088/2041-8205/783/1/L18
  * Lotz et al. (2008) Lotz, J. M., Jonsson, P., Cox, T. J., & Primack, J. R. 2008, MNRAS, 391, 1137, doi: 10.1111/j.1365-2966.2008.14004.x
  * Lotz et al. (2010a) —. 2010a, MNRAS, 404, 575, doi: 10.1111/j.1365-2966.2010.16268.x
  * Lotz et al. (2010b) —. 2010b, MNRAS, 404, 590, doi: 10.1111/j.1365-2966.2010.16269.x
  * Manuwal et al. (2022) Manuwal, A., Ludlow, A. D., Stevens, A. R. H., Wright, R. J., & Robotham, A. S. G. 2022, MNRAS, 510, 3408, doi: 10.1093/mnras/stab3534
  * Michikoshi & Kokubo (2020) Michikoshi, S., & Kokubo, E. 2020, ApJ, 897, 65, doi: 10.3847/1538-4357/ab9369
  * Minchev & Quillen (2006) Minchev, I., & Quillen, A. C. 2006, MNRAS, 368, 623, doi: 10.1111/j.1365-2966.2006.10129.x
  * Ono et al. (2023) Ono, Y., Harikane, Y., Ouchi, M., et al. 2023, ApJ, 951, 72, doi: 10.3847/1538-4357/acd44a
  * Ostriker & Peebles (1973) Ostriker, J. P., & Peebles, P. J. E. 1973, ApJ, 186, 467, doi: 10.1086/152513
  * Papaderos et al. (2023) Papaderos, P., Östlin, G., & Breda, I. 2023, A&A, 673, A30, doi: 10.1051/0004-6361/202245769
  * Renaud et al. (2021) Renaud, F., Romeo, A. B., & Agertz, O. 2021, MNRAS, 508, 352, doi: 10.1093/mnras/stab2604
  * Romeo et al. (2023) Romeo, A. B., Agertz, O., & Renaud, F. 2023, MNRAS, 518, 1002, doi: 10.1093/mnras/stac3074
  * Saha & Cortesi (2018) Saha, K., & Cortesi, A. 2018, ApJ, 862, L12, doi: 10.3847/2041-8213/aad23a
  * Saha & Elmegreen (2018) Saha, K., & Elmegreen, B. 2018, ApJ, 858, 24, doi: 10.3847/1538-4357/aabacd
  * Saha & Naab (2013) Saha, K., & Naab, T. 2013, MNRAS, 434, 1287, doi: 10.1093/mnras/stt1088
  * Sellwood (1980) Sellwood, J. A. 1980, A&A, 89, 296
  * Sellwood (2016) —. 2016, ApJ, 819, 92, doi: 10.3847/0004-637X/819/2/92
  * Sellwood (2021) —. 2021, MNRAS, 506, 3018, doi: 10.1093/mnras/stab1924
  * Sellwood & Evans (2001) Sellwood, J. A., & Evans, N. W. 2001, ApJ, 546, 176, doi: 10.1086/318228
  * Sérsic (1963) Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41
  * Shen & Sellwood (2004) Shen, J., & Sellwood, J. A. 2004, ApJ, 604, 614, doi: 10.1086/382124
  * Sheth et al. (2008) Sheth, K., Elmegreen, D. M., Elmegreen, B. G., et al. 2008, ApJ, 675, 1141, doi: 10.1086/524980
  * Springel (2005) Springel, V. 2005, MNRAS, 364, 1105, doi: 10.1111/j.1365-2966.2005.09655.x
  * Springel et al. (2001) Springel, V., Yoshida, N., & White, S. D. M. 2001, New A, 6, 79, doi: 10.1016/S1384-1076(01)00042-2
  * Takeuchi et al. (2015) Takeuchi, T. M., Ohta, K., Yuma, S., & Yabe, K. 2015, ApJ, 801, 2, doi: 10.1088/0004-637X/801/1/2
  * Toomre (1964) Toomre, A. 1964, ApJ, 139, 1217, doi: 10.1086/147861
  * Toomre (1981) Toomre, A. 1981, in Structure and Evolution of Normal Galaxies, ed. S. M. Fall & D. Lynden-Bell, 111–136
  * Varela-Lavin et al. (2023) Varela-Lavin, S., Gómez, F. A., Tissera, P. B., et al. 2023, MNRAS, 523, 5853, doi: 10.1093/mnras/stad1724
  * Vauterin & Dejonghe (1996) Vauterin, P., & Dejonghe, H. 1996, A&A, 313, 465, doi: 10.48550/arXiv.astro-ph/9603094
  * Worrakitpoonpon (2023) Worrakitpoonpon, T. 2023, ApJ, 958, 128, doi: 10.3847/1538-4357/acf657
  * Yozin & Bekki (2014) Yozin, C., & Bekki, K. 2014, MNRAS, 439, 1948, doi: 10.1093/mnras/stu075
  * Yuma et al. (2012) Yuma, S., Ohta, K., & Yabe, K. 2012, ApJ, 761, 19, doi: 10.1088/0004-637X/761/1/19
  * Yuma et al. (2011) Yuma, S., Ohta, K., Yabe, K., Kajisawa, M., & Ichikawa, T. 2011, ApJ, 736, 92, doi: 10.1088/0004-637X/736/2/92
  * Zakharova et al. (2023) Zakharova, D., Sotnikova, N. Y., Smirnov, A. A., & Savchenko, S. S. 2023, MNRAS, 525, 2069, doi: 10.1093/mnras/stad2368


# Heavy Flavoured Meson Fragmentation Functions in $e^{+}e^{-}$ annihilation
up to NNLO + NNLL 111Contribution to the proceedings of the $57^{\mathrm{th}}$
Recontres de Moriond 2023, QCD and High Energy Interactions.

Leonardo Bonino 222Speaker at the conference. ∗ Matteo Cacciari † ‡ and
Giovanni Stagnitto ∗ ∗ Physik-Institut, Universität Zürich,
Winterthurerstrasse 190, CH-8057 Zürich, Switzerland
† Sorbonne Université CNRS, Laboratoire de Physique Théorique et Hautes
Énergies, LPTHE, F-75005 Paris, France
‡ Université Paris Cité, F-75006 Paris, France

###### Abstract

In this contribution, we make use of the QCD perturbative fragmentation
function formalism to describe the one-particle inclusive fragmentation of a
heavy quark produced in $e^{+}e^{-}$ annihilation at
$\mathcal{O}(\alpha_{S}^{2})$. We perform the computation analytically in
Mellin-space. We resum soft-gluons effects in initial conditions and
coefficient functions and perform evolution up to NNLL accuracy, obtaining the
first NNLO + NNLL prediction for charm quark production. We study the impact
of different matching schemes and Landau pole prescriptions in soft-gluon
resummation, finding significant differences. We extract simple non-
perturbative fragmentation functions for $B$ and $D^{*}$ mesons by comparing
the perturbative prediction with the data from CLEO, BELLE and LEP
experiments. We find that for charm mesons the experimental results from
CLEO/BELLE and from LEP are not reconcilable with the standard DGLAP
evolution.

##  1 Introduction

The perturbative fragmentation function formalism is used in QCD to resum
observables describing the fragmentation of heavy quarks into hadrons. We
consider the one-particle inclusive production of a single heavy hadron $H$
($D$ or $B$ meson) at high energy in $e^{+}+e^{-}\to V(Q)\to H(p)+X$,
inclusive on all other final state particles $X$, and we study it to next-to-
next-to-leading order accuracy in QCD. A paper describing these findings is in
preparation [1].

###  1.1 The Fragmentation Function formalism

If the quark is heavy, the parton fragmentation function can be factored in
perturbative $D^{p}$ and non-perturbative $D^{np}$ components such that the
hadron cross section is $\sigma_{H}(x,Q^{2})=\hat{\sigma}\otimes D^{p}\otimes
D^{np}$, with $x$ energy fraction of the heavy quark, ratio of kinematic
invariants and $\hat{\sigma}$ massless partonic cross section. The process-
independent $D^{p}$ can be written as initial conditions
$D^{\mathrm{(ini)}}_{b/Q}$ evolved via DGLAP from the initial factorization
scale $\mu_{0F}$, to the final one $\mu_{F}$ [2], such that the only process-
dependent ingredients are the massless partonic coefficient functions
$C_{a}^{(proc)}$. It is customary to move to Mellin space where the cross
section reads

$\frac{1}{\sigma_{tot}}\sigma_{N}\left(Q^{2},m^{2}\right)=\frac{\sigma^{(LO)}}{\sigma_{tot}}\sum_{a,b}C_{a,N}^{(e^{+}e^{-})}\left(Q^{2},\mu^{2},\mu_{F}^{2}\right)E_{ab,N}\left(\mu_{F}^{2},\mu_{0F}^{2}\right)D_{b/\mathcal{Q},N}^{\mathrm{(ini)}}\left(\mu_{0}^{2},\mu_{0F}^{2},m^{2}\right)\,,$
(1)

with $\mu$ renormalisation scale, $\sigma_{tot}$ total heavy quark cross
section and $\sigma^{(LO)}$ massless parton leading order cross section. The
initial conditions are available analytically in Mellin space up to NNLO [3,
4, 5], as well as the coefficient functions [6, 7, 8]. We perform the
evolution at NNLL accuracy with the public library $\mathtt{MELA}$ [9].

###  1.2 Soft-Gluon Resummation

In the Sudakov region ($x\to 1$), large logarithms of soft origin appear such
that initial conditions and coefficient functions must be resummed [10, 11].
Mellin space analytical expressions for the resummation have recently become
available up to NNLL [5, 12, 13]. In our calculation, we have implemented both
additive and multiplicative matching (such as $\log$-$R$) to match to the
fixed order result. A singularity in the resummed initial conditions, giving
rise to a branch cut and related to the divergent behaviour of the running
coupling $\alpha_{S}(\mu_{0}^{2})$ near the Landau pole $N^{L}_{\mathrm{ini}}$
at $\mu_{0}\simeq\Lambda_{QCD}$, signals the onset of non-perturbative
phenomena at large $N$ 333For the charm quark at
$\sqrt{Q^{2}}=10.6\,\mathrm{GeV}$ $N^{L}_{\mathrm{ini}}\sim 7$ and for the
bottom at $\sqrt{Q^{2}}=91.2\,\mathrm{GeV}$ $N^{L}_{\mathrm{ini}}\sim 32$.. To
have phenomenologically sensible perturbative predictions a prescription for
the pole must be used. The CNO [14] prescription accounts to a shift in $N$
with a tunable parameter $f$, while the CGMP [13] prescription truncates the
exponent of the Sudakov factor 444While the CNO prescription does not
introduce power corrections larger than generally expected for the processes
in question, this is the case for the CGMP prescription.. In Section 2 we
study the effects on phenomenology of the prescription chosen.

##  2 Towards phenomenology

The non-perturbative part of the fragmentation function is inferred from a fit
to the data. NLO + NLL fits for both $B$ and $D$ hadrons are available [14]
while an effective $\alpha_{S}$ approach [15] was calling for a full NNLO
analysis. Recently the NNLO + NNLL level was reached for the bottom quark [13]
and fits were performed as well. Here [1] we present $B$-hadrons fits with our
implementation and the first fits to the $D$-mesons fragmentation functions at
NNLO + NNLL level. Many experiments throughout the years measured one-particle
inclusive observable in $e^{+}e^{-}$ collisions. For $B$-hadrons at the
$Z^{0}$-peak we used data from ALEPH [17], SLD [18], OPAL [19] and DELPHI [20]
555While the ALEPH set refers specifically to $B$-mesons, the SLD, DELPHI and
OPAL data are for all weakly decaying bottom-flavoured hadrons and include
about $10\%$ of baryons. . For $D$-mesons the data are available at
$\sqrt{Q^{2}}=91.2\,\mathrm{GeV}$ from ALEPH [21] and at
$\sqrt{Q^{2}}=10.6\,\mathrm{GeV}$ from BELLE [22] and CLEO [23, 24]. We fit
$B$ data with a two-parameters non-perturbative ansatz and $D$ data with a
four-parameters one as in [16, 14].

###  2.1 Charm ratio

A fully-perturbative observable, the ratio of charm meson moments from LEP and
from CLEO and BELLE, was studied at NLO + NLL [14], showing an inconsistency
between theory and experiments. In Fig. 1 we study this observable at NNLO +
NNLL, confirming the inconsistency. We remark that this observable and the
prediction exclusively depend upon the evolution between the $\Upsilon(4S)$
and $Z^{0}$ energies. A possible explanation for the inconsistency can be
found in the presence of large non-perturbative power suppressed effects in
the coefficient functions.

Figure 1: Charm ratio up to NNLO + NNLL. Theory band uncertainty:
$\mu_{Z/\Upsilon}=\xi M_{Z/\Upsilon}$ with $1/2<\xi<2$.

###  2.2 Effects of Landau Pole prescription

We have also studied the effects of changing the Landau pole prescriptions on
the initial conditions 666On the coefficient function the effects of the
prescription are negligible since the pole appears at larger values of $N$..
We observe a variation much larger than any scale variation and, regardless of
the prescription chosen, the result is also unstable under variation of the
matching scheme. We conclude that on the initial conditions, the
characteristics of the perturbative result are strongly influenced by the
choice of the Landau pole regularization. The impact of this effect on the
final cross section is shown in Fig. 2 left for the bottom. Due to the smaller
mass (and low experiment scale of $\sqrt{Q^{2}}=10.6\,\mathrm{GeV}$), the
effect is much more severe for the charm as in Fig. 2 right. We notice that
the CGMP prescription is not suited to describe the $D$-meson data already at
NLO level with our non-perturbative ansatz, due to the steep fall at $x\sim
0.8$. Since many data points are in that region, for all our fits the CNO
prescription was used.

Figure 2: Left: Bottom differential cross section al NLO + NLL and NNLO + NNLL
with CNO ($f=1.25$) and CGMP prescriptions for
$\mu_{F}=\sqrt{Q^{2}}=91.2\,\mathrm{GeV}$ in $\log$-$R$ matching scheme.
Right: Same for the charm for $\mu_{F}=\sqrt{Q^{2}}=10.6\,\mathrm{GeV}$.

###  2.3 Fits to the data

As expected, the impact of the effect seen in Fig. 2 on phenomenology is
significant and absorbed in a different results of the fit.

Figure 3: Left: Fits to SLD data with CNO ($f=1.25$) prescription and
$\log$-$R$ matching. Right: Same for BELLE data ($f=2.0$).

In Fig. 3 our results of the fits to the $B$ and $D$ hadrons fragmentation
functions are presented for two of the data sets. While all fits for the
$B$-hadrons are satisfactory, at NNLO + NNLL level we are not able to describe
the $D$-mesons data in the $x>0.9$ region (for all data sets) as shown in Fig.
3. We suspect the reason for this being related to the Landau pole non-
perturbative effects.

##  3 Conclusions

We have presented the first results for the NNLO + NNLL one-particle inclusive
charm fragmentation in $e^{+}e^{-}$ annihilation and relative fits to the
data. We uncover large non-perturbative effects arising from the choice of the
prescription for Landau pole. The characteristics of the perturbative result
are strongly influenced by the choice of the Landau pole regularization, which
has a direct impact on the appropriate shape of the non-perturbative component
needed to describe experimental data. Furthermore, the inclusion of one more
perturbative order confirms the charm ratio inconsistency. In conclusion, non-
perturbative effects, possibly of different origin, play a crucial role in
precision observables for charm quark fragmentation.

## Acknowledgments

LB would like to thank Thomas Gehrmann for useful discussions. This work has
received funding from the Swiss National Science Foundation (SNF) under
contract 200020-204200.

## References

## References

  * [1] L. Bonino, M. Cacciari and G. Stagnitto, In preparation.
  * [2] B. Mele and P. Nason, Phys. Rev. D 361, 626-644 (1991).
  * [3] K. Melnikov and A. Mitov, Phys. Rev. D 70, 034027 (2004).
  * [4] A. Mitov, Phys. Rev. D 71, 054021 (2005).
  * [5] F. Maltoni, G. Ridolfi, M. Ubiali and M. Zaro, JHEP 10, 027 (2022).
  * [6] P. J. Rijken and W. L. van Neerven, Nucl. Phys. B 487, 233-207 (1997).
  * [7] J. Blümlein and V. Ravindran, Nucl. Phys. B 749, 1-24 (2006).
  * [8] A. Mitov and S. Moch, Nucl. Phys. B 751, 18-52 (2006).
  * [9] V. Bertone, S. Carrazza and E. R. Nocera, JHEP 03, 46 (2015).
  * [10] S. Catani and M. Cacciari, Nucl. Phys. B 617, 253-290 (2001).
  * [11] G. Ridolfi, M. Ubiali and M. Zaro, JHEP 01, 196 (2020).
  * [12] U. Aglietti, G. Corcella and G. Ferrera, Nucl. Phys. B 775, 162-201 (2007).
  * [13] M. Czakon, T. Generet, A. Mitov and R. Poncelet, JHEP 03, 251 (2023).
  * [14] M. Cacciari, P. Nason and C. Oleari, JHEP 0604, 006 (2006).
  * [15] G. Corcella and G. Ferrera, JHEP 12, 029 (2007).
  * [16] G. Colangelo and P. Nason, Phys. Lett. B 285, 167–171 (1992).
  * [17] ALEPH Coll., Phys. Lett. B 512, 30-48 (2001).
  * [18] SLD Coll., Phys. Rev. D 65, 092006 (2002).
  * [19] OPAL Coll., Eur. Phys. J C 29, 463-478 (2003).
  * [20] DELPHI Coll., Eur. Phys. J C 71, 1557 (2011).
  * [21] ALEPH Coll., Eur. Phys. J C 16, 597-611 (2000).
  * [22] BELLE Coll., Phys. Rev. D 73, 032002 (2006).
  * [23] CLEO Coll., Phys. Rev. D 70, 112001 (2004).
  * [24] CLEO Coll., Phys. Rev. D 57, 1350-1358 (1998).