# Singing Voice Synthesis with Vibrato Modeling and Latent Energy
Representation

Yingjie Song1∗, Wei Song2, Wei Zhang2, Zhengchen Zhang2, Dan Zeng1, Zhi Liu1,
Yang Yu3 ∗ Work was done during internship at JD AI Research. 1 School of
Communication and Information Engineering, Shanghai University, Shanghai,
China
Email: {jie0222<EMAIL_ADDRESS><EMAIL_ADDRESS>2JD AI Research,
Beijing, China
Email: {songwei11<EMAIL_ADDRESS><EMAIL_ADDRESS>3Shanghai
Conservatory of Music, Shanghai, China
Email<EMAIL_ADDRESS>
###### Abstract

This paper proposes an expressive singing voice synthesis system by
introducing explicit vibrato modeling and latent energy representation.
Vibrato is essential to the naturalness of synthesized sound, due to the
inherent characteristics of human singing. Hence, a deep learning-based
vibrato model is introduced in this paper to control the vibrato’s likeliness,
rate, depth and phase in singing, where the vibrato likeliness represents the
existence probability of vibrato and it would help improve the singing voice’s
naturalness. Actually, there is no annotated label about vibrato likeliness in
existing singing corpus. We adopt a novel vibrato likeliness labeling method
to label the vibrato likeliness automatically. Meanwhile, the power
spectrogram of audio contains rich information that can improve the
expressiveness of singing. An autoencoder-based latent energy bottleneck
feature is proposed for expressive singing voice synthesis. Experimental
results on the open dataset NUS48E show that both the vibrato modeling and the
latent energy representation could significantly improve the expressiveness of
singing voice. The audio samples are shown in the demo
website111https://mango321321.github.io/ExpressiveSing/.

###### Index Terms:

Singing Voice Synthesis, Pitch Model, Vibrato, Energy Representation

##  I Introduction

Singing voice synthesis (SVS) systems [1, 2, 3, 4, 5, 6, 7] generate singing
voices from musical scores which contain music information such as lyrics,
tempo, pitch, etc. SVS is similar to the text-to-speech (TTS) task [8, 9, 10,
11, 12, 13] in terms of generating speech from text. The difference is that
SVS should reproduce emotional ebbs and flows beyond generating a natural-
sounding voice [14, 15, 16]. The state-of-the-art singing voice synthesis
systems already provide multiple solutions for the application, while
synthesized voices lack expressiveness compared to real voices.

Over the past few decades, researchers have made many efforts to enhance
expressiveness. In [17], Umbert et al. presented the expressiveness is related
to a set of voice parameters and classifies the expression control parameters
into four categories including melody, dynamics, rhythm, and timbre. Rooted in
[17], many subsequent works for expressive singing voice synthesis focus on
improvement about these control parameters, especially melody, dynamics a.k.a
pitch and energy. For example, to predict pitch feature better, Yi et al. [18]
utilized deep autoregressive network to capture the dependencies among the
consecutive acoustic features. Zhuang et al. [1] separated the pitch feature
from the acoustic feature to avoid the interdependence between these pitch
features and the timbre features. Ren et al. [9] introduced the pitch and
energy information into the speech generation task and presented variance
adaptors to make the generated audio expressive. In this article, we enhance
expressiveness from these two aspects as well.

Pitch conveys melody according to score and contains score-independent dynamic
movements in a singing voice, such as vibrato. Because vibrato is not
available in music score and it is unknown where the vibrato will appear, we
introduced vibrato likeliness into pitch model to predict where the vibrato
will occur. Vibrato likeliness represents the existence probability of vibrato
but there is no annotated data in the existing dataset. Additional labeling is
feasible but expensive. To address this problem, we propose a novel automatic
labeling model called Vibrato Likeliness Labeling Network (VLLN), which is
trained using the simulated data.

Many studies show that energy is related to volume and prosody of the singing
voice which affect expressiveness crucially. The energy prediction can be
recast as the task of recovering the feature of power spectrogram. In previous
work [1, 9, 19], the L2-Norm of the amplitude spectrogram is taken as energy
information to promote expressiveness. However, this $1$-dimensional energy
representation loses a lot of spectrogram information. In this work, an
autoencoder-based subnetwork is used to learn a compact and representative
latent energy representation from the input musical score information.

The main contributions of this paper are summarized as follows:

  1. 1.

An SVS system named ExpressiveSing is proposed which is adapted from an End-
to-End TTS framework, and contains an accurate pitch model and a latent energy
representation prediction module.

  2. 2.

The latent energy representation is proposed to utilize energy feature
generating expressive singing voices.

  3. 3.

A pitch model with accurate vibrato modeling is proposed by taking the vibrato
likeliness into account. The vibrato likeliness is annotated by vibrato
likeliness labeling method automatically.

The remainder of this paper is organized as follows: Section II describes the
related work of pitch model and energy representation in SVS systems. Section
III introduces the proposed singing voice synthesis system in detail. The
experimental setup is given in Section IV and results are presented in Section
V. Conclusion follows in Section VI.

##  II Related Work

###  II-A Pitch Model in SVS

The traditional pitch models include rule-based method, template-based method,
unit concatenation method and HMM-based method. The rule-based method [20]
designed rules for each singing style based on professional musical knowledge.
It mainly depended on designer’s experience, while template-based method
proposed a simpler solution. Yukara et al. [21] chose typical examples from
commercial recordings to form a library of parametrized templates. The unit
concatenation method [22] used real contours extracted from singing data and
concatenated them to generate pitch contours. Compared with the above
mentioned methods, the HMM-based method [23] is more flexible. It generated
pitch contours using hidden Markov models (HMMs). The disadvantage of this
approach was that it was poor to retain the detailed information, resulting in
the over-smoothing effect.

Recently, data-driven approaches have been introduced into pitch modeling in
SVS. They emphasized the vibrato during the pitch modeling process to further
improve the expressiveness. Yi et al. [18] utilized the dependencies among the
acoustic features of consecutive frames to help the pitch model understand
vibrato. However, it still suffered the over-smoothing effect to some extent.
Bonada et al. [24] proposed a neural-parametric pitch model which solved the
problem of over-smoothing effectively. It parametered the vibrato with phase,
depth, and rate. Actually, there was noise mixed in vibrato component. Bonada
et al. directly constrained vibrato depth to be zero at note transitions,
which was inelegant and could not eliminate noise completely.

To solve the noise problem, we propose an additional parameter, namely vibrato
likeliness, and automatically annotate the label of vibrato likeliness by the
proposed Vibrato Likeliness Labeling Network.

###  II-B Energy Representation in SVS

Typically, the previous SVS systems implicitly modeled energy during
synthesizing acoustic features such as [25, 5, 26]. The TTS system,
Fastspeech2 [9] firstly proposed variance adaptor to model energy explicitly.
It was proved to be beneficial for expressiveness. Concretely, the L2-norm of
the amplitude of each short-time Fourier transform (STFT) frame was regarded
as the energy and then the quantized energy is mapped to the embedding. In
other words, the energy was compressed to 1-dimension vector firstly. Because
the size of amplitude spectrogram is always thousands of dimensions, the
operation of compression to 1-dimension vector would cause massive information
loss. Following Fastspeech2, LiteSing [1] makes the according modification in
SVS system. The problem of information loss still existed. We propose an
autoencoder-based latent spectrogram bottleneck feature to ease the
information loss of compression. The experiment in Section V shows that our
proposed energy representation is superior to the one used in LiteSing and
Fastspeech2.

Figure 1: The overall architecture of ExpressiveSing, where ‘c’ denotes
concatenation, and ‘+’ denotes addition. ‘f’ refers to formula (3). ‘FC’ is
short for fully connected layer and ‘LR’ denotes length regulator. ‘Rep.’
stands for representation. The dashed arrow in subfigure (a) indicates that
the pitch model is trained separately.

##  III Proposed Framework

###  III-A Network Architecture

The overall architecture of proposed ExpressiveSing is shown in Fig. 1(a). The
network architecture of ExpressiveSing is a DIAN [27]-based network, which
contains an encoder for extracting phoneme representation, and a decoder for
the acoustic feature prediction. The encoder takes singer ID and phonemes at
phoneme level as input and the following length regulator expands the
generated phoneme representation to frame level. The decoder predicts acoustic
features at frame level, i.e. MGC, BAP, and pitch, which can be converted to
audio by WORLD vocoder [28]. Compared to DIAN, ExpressiveSing designs pitch
model, energy representation predictor and power spectrogram decoder
additionally. The pitch model predicts human-like pitch contour from music
notes. As for the input of pitch model, we use the note, phonemes at frame
level and singer ID, so pitch model outputs pitch sequence at frame level. The
energy representation predictor generates latent energy representation and the
power spectrogram decoder reproduces the power spectrogram from the latent
energy representation.

###  III-B Latent Energy Representation

An autoencoder-based network architecture is used in ExpressiveSing to predict
the latent energy representation as shown in Fig. 1(a). The pitch embedding
and phoneme representation are concatenated as the input of the energy
representation predictor $P_{e}$. The energy representation predictor learns a
compact and representative bottleneck feature $E\in R^{T\times N}$ of power
spectrogram, where $T$ is the frame count and $N$ denotes the dimension of the
energy representation. The power spectrogram decoder $D_{e}$ composed of $1$D
convolutions is trained to recover power spectrogram $S$ from the latent
representation $E$. The power spectrogram decoder is jointly trained with
acoustic feature prediction. The loss function for the power spectrogram
decoder is shown as follows:

$\mathcal{L}(\textbf{S},\hat{\textbf{S}})=||\textbf{S}-D_{e}(\textbf{E})||_{2}$
(1)

where $\hat{S}$ denotes the predicted power spectrogram.

###  III-C Pitch Model

Pitch contour is decomposed into pitch intonation and vibrato [29, 30].
Intonation is responsible for singing in tune according to music note. Vibrato
shows the skill of singers. We parameter the vibrato with phase, depth, rate,
and the specially-designed likeliness. Vibrato depth and rate are the width
and speed of the pitch fluctuation respectively [31, 32]. Phase refers to the
initial phase of vibrato, and likeliness stands for the possibility of the
appearance of vibrato. The training target of The vibrato depth, rate, and
phase are extracted by Hilbert transformation. The vibrato likeliness is
annotated by the proposed Vibrato Likeliness Labeling Network, which is
described in Section III-D.

Figure 1(b) illustrates the network architecture of the pitch model. The
phoneme, music notes and singer ID are embedded as the input.

A Transformer [33]-based network is used to predict the vibrato phase, depth,
rate, likeliness and intonation at the frame level. As for the loss function,
MSE loss is used in intonation, rate, depth, and phase prediction while the
cross entropy loss is adopted for vibrato likeliness prediction.

Following [18], to avoid out of tune, the predicted intonation i is post-
processed by

${\textbf{{i}}}_{post}={\textbf{{i}}}-\tilde{\textbf{{i}}}+{\textbf{{i}}}_{n}$
(2)

where ${\textbf{{i}}}_{n}$ means the music notes at frame level, and
$\tilde{\textbf{{i}}}$ is the smoothed curve of i. A triangular window is used
to convolve with i to get the smoothed result $\tilde{\textbf{{i}}}$.

For the predicted vibrato likeliness in pitch model, the average value
$l_{mean}$ is calculated on each note. If $l_{mean}$ is greater than the
threshold $\epsilon$ ($\epsilon=0.5$ in our experiments), the note is
considered to be a vibrato section, otherwise not. Finally, the vibrato
${v}_{t}~{}(t\in\\{1,2,\cdots,T\\})$ is synthesized by function:

${v}_{t}=\left\\{\begin{aligned}
&{l}_{t}\cdot{e}_{t}\cdot\cos(\frac{2\pi{r}_{t}\cdot
t}{f_{s}}+{\phi}),&&{l_{mean}\geq\epsilon}\\\
&0,&&{l_{mean}<\epsilon}\end{aligned}\right.$ (3)

where ${e}_{t},~{}{r}_{t},~{}{l}_{t}$ denote the $t$-th frame of predicted
depth, rate, and likeliness of vibrato, respectively. $~{}{\phi}$ represents
the predicted initial phase. ${f_{s}}$ stands for the sample rate. Then the
final output of pitch model is the sum of intonation ${\textbf{{i}}}_{post}$
and predicted vibrato ${v}_{t}$.

###  III-D Vibrato Likeliness Labeling Network

Training pitch model requires vibrato likeliness label, i.e., the ground truth
of vibrato existence probability. However, no such label exists in existing
singing corpus. To solve this problem, we propose a novel labeling method to
get the ground truth of vibrato likeliness automatically. We randomly
integrate vibrato to smoothed pitch contour to simulate the real pitch contour
and then VLLN is trained on the simulated data.

Concretely, the real pitch is extracted from the singing voice firstly and
then the pitch contour is convolved with a triangular window to remove
vibrato, resulting in a smoothed pitch contour. According to our observation,
vibrato commonly appears in notes with a longer duration. Hence, we select
notes with duration longer than one second to add simulated vibrato. In order
to simulate vibrato on the selected music notes, we exploit a sinusoid curve
in 6 Hz and vary its depth in range 0 to 2 in midi randomly and gradually.
Then the synthesized vibrato is added to the smoothed pitch contour to
generate the simulated data. Then frame-level labels of vibrato likeliness (0
or 1) can be obtained because we know which segment has vibrato.

Figure 2: The Vibrato Likeliness Labeling Network.

The Vibrato Likeliness Labeling Network is shown in Fig. 2, which is a
Transformer-based network. Sigmoid function is used as activation function in
the last layer. The cross entropy loss is adopted as the training objective
for vibrato likeliness prediction. The details of training and inference
procedures are shown in Fig. 3. During train stage, we employ the simulated
data with known vibrato likeliness label to train VLLN. Afterwards, in
inference stage, VLLN is used to label the vibrato likeliness for the real
pitch of singing corpus. The color on the pitch contour represents vibrato
likeliness. The darker the color, the higher the vibrato likeliness.

Figure 3: The training and inference procedures of Vibrato Likeliness Labeling
Network.

##  IV Experimental Setup

###  IV-A Dataset

NUS48E corpus [34] is used to evaluate the proposed system, which consists of
48 popular English songs, and only 20 songs are unique. The lyrics are sung
and read out by 12 singers (4 songs each singer). Only the singing data is
used for the proposed model training, which contains 115 minutes recordings.

The training data are converted to 16kHz in our experiments. All the audios
are split into segments shorter than 10 seconds and redundant silence of each
utterance is removed. We annotated 36 songs with music scores manually because
the NUS48E corpus does not contain music score labeling. The pitch model is
trained by 36 annotated songs while the acoustic model uses all 48 songs. We
keep $10\%$ as test set and use others for training. Pitch is extracted using
the WORLD vocoder [28] and the 3$\sim$8Hz part extracted by the band-pass
filter is considered to be vibrato, the hop size for acoustic feature
extraction is set to 10 milliseconds in our experiments.

###  IV-B Network Configuration and Implementation

We follow the network architecture configuration in DIAN [27], and the energy
representation predictor follows the variance adaptor setup in FastSpeech2 [9]
except that the output layer predicts a 256-dimension vector. The power
spectrogram decoder contains a stack of 4 convolutional layers and a linear
projection layer, the kernel size is 5 and channel size is 512 for each
convolution layer, batch normalization [35] and ReLU [36] activation is used
for each convolution layer. The dimension of target power spectrogram is 1025.
In training process of ExpressiveSing, we use real pitch as input. In
inference stage, the predicted pitch from pitch model is taken as input. Table
I shows the dimension of features used in ExpressiveSing.

The proposed pitch model follows the duration model network configuration in
[27] and has the following modifications: (1) music notes are also taken as
input as well as the phoneme input; (2) the pitch model contains 5 output
layers, which are used to predict intonation, phase, depth, rate and
likeliness. The loss weight of intonation prediction is 1, and the loss weight
of vibrato rate, depth, phase, and likeliness are all 1e-2.

The Vibrato Likeliness Labeling Network’s configuration is the same as pitch
model, except that there is no singer embedding in Vibrato Likeliness Labeling
Network and only contains 1 output layer.

Adam optimizer [37] with batch size 32 is used for optimization with
$\beta_{1}$ = 0.9 and $\beta_{2}$ = 0.999. ExpressiveSing is trained with
learning rate that starts from 1e-3. The learning rate exponentially decays to
1e-5 in the first 50K steps, then the training continued with a constant
learning rate of 1e-5 until 60K steps. The pitch model and Vibrato Likeliness
Labeling Network are trained with Noam[33] learning rate schedule with a peak
learning rate of 1e-4, with a total of 40k steps for each model.

TABLE I: The dimension of features used in ExpressiveSing. Feature | Dimension
---|---
Phone Representation |  512
Pitch Embedding |  32
Latent Energy Representation |  256
Mel-Generalized Cepstral |  60
Band Aperiodicity |  4
Pitch |  1

##  V Experimental Results

TABLE II: Mean opinion scores (MOS) of different systems. Pronunciation accuracy (Pronun acc.), sound quality (Quality), and naturalness are evaluated separately. System | Pronun acc. | Quality | Naturalness
---|---|---|---
Baseline |  3.120 |  2.373 |  2.460
USVC[38] |  3.160 |  2.993 |  2.707
WGANSing[39] |  3.367 |  2.960 |  3.133
UCSVC[40] |  3.900 |  3.547 |  3.780
ExpressiveSing |  3.780 |  3.413 |  3.467

###  V-A Overall Performance

The system performance is evaluated using subjective mean opinion scores
(MOS). A total of $15$ subjects participated in the evaluation. The listeners
were asked to mark the singing samples from three aspects: pronunciation
accuracy, sound quality and naturalness. The scope of the score is from 1
(extremely bad) to 5 (extremely good). Five systems are compared, and ten
samples generated by each system are prepared for the evaluation. The baseline
in our experiment is built using the network structure proposed in DIAN [27].
The difference between baseline and ExpressiveSing is that the baseline system
has no pitch model nor latent energy representation module. The proposed
ExpressiveSing is also compared with the state-of-the-art systems including
two singing voice conversion systems USVC [38] and UCSVC [40], as well as a
singing synthesis system WGANSing [39]. These systems are trained on the same
NUS48E dataset. It is notable that the singing voice conversion systems have
real human singing pitch as input, which would make the generated voices more
natural. The samples of
USVC222https://enk100.github.io/Unsupervised_Singing_Voice_Conversion/ and
UCSVC333https://singing-conversion.github.io/ are downloaded from their demo
web pages. We implemented the WGANSing method, and trained it on the same
training set as ExpressiveSing. Because WGANSing uses the real pitch
originally, to ensure fairness, WGANSing is conditioned with the predicted
pitch from our pitch model in this experiment. The evaluation results are
shown in Table II. We can see that ExpressiveSing (ES) outperforms WGANSing,
USVC and the baseline system. And ES obtains comparable scores with UCSVC,
even though UCSVC uses real pitch as input and exploits additional speech
datasets for training.

###  V-B Ablation Study

We conduct ablation studies to evaluate the effects of the proposed pitch
model, and the latent energy representation. Both subjective and objective
evaluations are performed. For subjective evaluation, the listeners are
instructed to do the A/B preference tests to compare the naturalness of the
singing samples generated by different systems. Two pairs of systems are
compared. One is the comparison between ExpressiveSing (ES) and ExpressiveSing
without vibrato likeliness (ES-V), where ES-V uses the pitch model without the
parameter of vibrato likeliness. The other is the comparison between
ExpressiveSing (ES) and ExpressiveSing without latent energy representation
(ES-E1). ES-E1 compressed energy to 1-dimension vector firstly and then a
lookup table is used to convert this single valued energy to embedding
representation, which is the same as LiteSing and Fastspeech2. The evaluation
results are shown in Fig. 4. One can find that the listeners strongly prefer
the ES model to either ES-V or ES-E1, which proves the effectiveness of the
proposed pitch model and latent energy representation.

Figure 4: A/B preference evaluation results of ablation study. ‘No pref.’
denotes no preference for the systems.

We further evaluate the pitch model with objective metrics including F0 Root
Mean Square Error (RMSE) and Correlation (CORR). The results are shown in
Table III. The synthesized F0 sequences are compared with the human singers’
F0 sequences and the score’s pitch notations, which are noted as ‘Natural’ and
’Score’ respectively. The more natural the generated singing is, the lower F0
RMSE and higher CORR are achieved. We can find that ES achieved lower F0 RMSE
and higher CORR than ES-V. To compare the proposed pitch model with the state-
of-the-art methods, we further implemented the pitch model proposed in [24],
which is named Parametric F0. The F0 RMSE and CORR of Parametric F0 are also
presented in Table III. It can be seen that ES achieved better performance
than Parametric F0.

TABLE III: Objective evaluation of systems with different pitch models. ExpressiveSing (ES), ExpressiveSing without vibrato likeliness (ES-V), and the model proposed in [24] are compared.  System | F0 RMSE (Hz) | F0 CORR
---|---|---
Natural | Score | Natural | Score
Parametric F0 [24] |  22.759 |  23.764 |  0.973 |  0.980
ES-V |  21.837 |  14.459 |  0.977 |  0.993
ES |  20.900 |  12.633 |  0.978 |  0.996

An example is illustrated in Fig. 5 to demonstrate the effect of the pitch
model. Compared to ES-V, ES generates fewer fluctuations at the positions of
note changing with the help of vibrato likeliness prediction. Hence, the noise
in the non-vibrato section is avoided, and the generated sound becomes more
steady.

Figure 5: An example of pitch contour generated by ES and ES-V, compared with
human singing voices (Natural) and the note of score (Music Note).

Mel-cepstral distortion (MCD) metric is used to verify if the proposed latent
energy representation is helpful. Three systems, ES, ES-E1, and ES-w/o, are
evaluated, where ES-w/o is the EpressiveSing model without energy information.
The MCD scores are 6.600, 6.787, and 6.800 respectively. ES obtained the
lowest MCD, which means the generated singing voice is more similar to human
singing. The result is consistent with the A/B test results. ES-w/o obtained
the highest MCD, which proved that the energy feature is helpful to expressive
singing voice synthesis.

##  VI Conclusion

According to the inherent characteristics of human singing, accurate vibrato
modeling and latent energy representation are introduced into the singing
voice synthesis system in this paper. The MOS results show that the proposed
system can synthesize expressive singing voice. The A/B test results show that
both the proposed vibrato modeling and latent energy representation could
improve the naturalness. The objective metrics of MCD, F0 RMSE and CORR are
consistent with the A/B test results. Both the subjective and objective
results verify the effectiveness of vibrato modeling of pitch and latent
energy representation.

## References

  * [1] X. Zhuang, T. Jiang, S.-Y. Chou, B. Wu, P. Hu, and S. Lui, “Litesing: Towards Fast, Lightweight and Expressive Singing Voice Synthesis,” in _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021, pp. 7078–7082.
  * [2] Y. Hono, K. Hashimoto, K. Oura, Y. Nankaku, and K. Tokuda, “Singing Voice Synthesis Based on Generative Adversarial Networks,” in _ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2019, pp. 6955–6959.
  * [3] M. Blaauw and J. Bonada, “Sequence-to-Sequence Singing Synthesis Using the Feed-Forward Transformer,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 7229–7233.
  * [4] J. Shi, S. Guo, N. Huo, Y. Zhang, and Q. Jin, “Sequence-To-Sequence Singing Voice Synthesis With Perceptual Entropy Loss,” in _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021, pp. 76–80.
  * [5] Y. Ren, X. Tan, T. Qin, J. Luan, Z. Zhao, and T.-Y. Liu, “DeepSinger: Singing Voice Synthesis with Data Mined From the Web,” in _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 2020, pp. 1979–1989.
  * [6] X. Gao, X. Tian, Y. Zhou, R. K. Das, and H. Li, “Personalized Singing Voice Generation Using WaveRNN,” in _Proc. Odyssey 2020 The Speaker and Language Recognition Workshop_ , 2020, pp. 252–258.
  * [7] K. Oura, A. Mase, T. Yamada, S. Muto, Y. Nankaku, and K. Tokuda, “Recent development of the hmm-based singing voice synthesis system—sinsy,” in _Seventh ISCA Workshop on Speech Synthesis_ , 2010.
  * [8] Y. Wang, R. Skerry-Ryan, D. Stanton, Y. Wu, R. J. Weiss, N. Jaitly, Z. Yang, Y. Xiao, Z. Chen, S. Bengio _et al._ , “Tacotron: Towards End-to-End Speech Synthesis,” _arXiv preprint arXiv:1703.10135_ , 2017.
  * [9] Y. Ren, C. Hu, X. Tan, T. Qin, S. Zhao, Z. Zhao, and T.-Y. Liu, “FastSpeech 2: Fast and High-Quality End-to-End Text to Speech,” _arXiv preprint arXiv:2006.04558_ , 2020.
  * [10] E. Song, J.-S. Kim, K. Byun, and H.-G. Kang, “Speaker-adaptive neural vocoders for parametric speech synthesis systems,” in _2020 IEEE 22nd International Workshop on Multimedia Signal Processing (MMSP)_. IEEE, 2020, pp. 1–5.
  * [11] R. Gupta and T. H. Falk, “Physiological quality-of-experience assessment of text-to-speech systems,” in _2016 IEEE 18th International Workshop on Multimedia Signal Processing (MMSP)_. IEEE, 2016, pp. 1–2.
  * [12] D. Jiang, I. Ravyse, H. Sahli, and Y. Zhang, “Accurate visual speech synthesis based on diviseme unit selection and concatenation,” in _2008 IEEE 10th Workshop on Multimedia Signal Processing_. IEEE, 2008, pp. 906–909.
  * [13] M. Schnell, M. Kustner, O. Jokisch, and R. Hoffmann, “Text-to-speech for low-resource systems,” in _2002 IEEE Workshop on Multimedia Signal Processing._ IEEE, 2002, pp. 259–262.
  * [14] H.-S. Pang, “On the use of the maximum likelihood estimation for analysis of vibrato tones,” _Applied Acoustics_ , vol. 65, no. 1, pp. 101–107, 2004\.
  * [15] T. L. Nwe and H. Li, “Exploring Vibrato-Motivated Acoustic Features for Singer Identification,” _IEEE Transactions on Audio, Speech, and Language Processing_ , vol. 15, no. 2, pp. 519–530, 2007.
  * [16] Y.-J. Luo, C.-C. Hsu, K. Agres, and D. Herremans, “Singing voice conversion with disentangled representations of singer and vocal technique using variational autoencoders,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 3277–3281.
  * [17] M. Umbert, J. Bonada, M. Goto, T. Nakano, and J. Sundberg, “Expression Control in Singing Voice Synthesis: Features, approaches, evaluation, and challenges,” _IEEE Signal Processing Magazine_ , vol. 32, no. 6, pp. 55–73, 2015.
  * [18] Y.-H. Yi, Y. Ai, Z.-H. Ling, and L.-R. Dai, “Singing Voice Synthesis Using Deep Autoregressive Neural Networks for Acoustic Modeling,” _arXiv preprint arXiv:1906.08977_ , 2019.
  * [19] L. Hantrakul, J. H. Engel, A. Roberts, and C. Gu, “Fast and Flexible Neural Audio Synthesis,” in _ISMIR_ , 2019.
  * [20] G. Berndtsson, “The kth rule system for singing synthesis,” _Computer Music Journal_ , vol. 20, no. 1, pp. 76–91, 1996.
  * [21] Y. Ikemiya, K. Itoyama, and H. G. Okuno, “Transferring vocal expression of f0 contour using singing voice synthesizer,” in _International Conference on Industrial, Engineering and Other Applications of Applied Intelligent Systems_. Springer, 2014, pp. 250–259.
  * [22] M. Umbert, J. Bonada, and M. Blaauw, “Generating singing voice expression contours based on unit selection,” in _Proc. SMAC_ , 2013.
  * [23] T. Nose, M. Kanemoto, T. Koriyama, and T. Kobayashi, “Hmm-based expressive singing voice synthesis with singing style control and robust pitch modeling,” _Computer Speech & Language_, vol. 34, no. 1, pp. 308–322, 2015\.
  * [24] J. Bonada and M. Blaauw, “Hybrid Neural-Parametric F0 Model for Singing Synthesis,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 7244–7248.
  * [25] Y. Cao, S. Liu, S. Kang, N. Hu, P. Liu, X. Liu, D. Su, D. Yu, and H. Meng, “Exploring Cross-lingual Singing Voice Synthesis Using Speech Data,” in _2021 12th International Symposium on Chinese Spoken Language Processing (ISCSLP)_. IEEE, 2021, pp. 1–5.
  * [26] Y. Gu, X. Yin, Y. Rao, Y. Wan, B. Tang, Y. Zhang, J. Chen, Y. Wang, and Z. Ma, “ByteSing: A Chinese Singing Voice Synthesis System Using Duration Allocated Encoder-Decoder Acoustic Models and WaveRNN Vocoders,” in _2021 12th International Symposium on Chinese Spoken Language Processing (ISCSLP)_. IEEE, 2021, pp. 1–5.
  * [27] W. Song, X. Yuan, Z. Zhang, C. Zhang, Y. Wu, X. He, and B. Zhou, “Dian: Duration Informed Auto-Regressive Network for Voice Cloning,” in _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021, pp. 8598–8602.
  * [28] M. Morise, F. Yokomori, and K. Ozawa, “World: A Vocoder-Based High-Quality Speech Synthesis System for Real-Time Applications,” _IEICE TRANSACTIONS on Information and Systems_ , vol. 99, no. 7, pp. 1877–1884, 2016.
  * [29] J. Sundberg, E. Prame, and J. Iwarsson, “Replicability and accuracy of pitch patterns in professional singers,” _Vocal fold physiology, controlling complexity and chaos_ , pp. 291–306, 1996.
  * [30] J. I. Shonle and K. E. Horan, “The pitch of vibrato tones,” _The Journal of the Acoustical Society of America_ , vol. 67, no. 1, pp. 246–252, 1980.
  * [31] N. Migita, M. Morise, and T. Nishiura, “A Study of Vibrato Features to Control Singing Voices,” _Proc. ICA_ , pp. 23–27, 2010.
  * [32] H.-Y. Gu and Z.-F. Lin, “Singing-voice Synthesis Using ANN Vibrato-parameter Models.” _J. Inf. Sci. Eng._ , vol. 30, no. 2, pp. 425–442, 2014.
  * [33] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention Is All You Need,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [34] Z. Duan, H. Fang, B. Li, K. C. Sim, and Y. Wang, “The nus sung and spoken lyrics corpus: A quantitative comparison of singing and speech,” in _2013 Asia-Pacific Signal and Information Processing Association Annual Summit and Conference_. IEEE, 2013, pp. 1–9.
  * [35] S. Ioffe and C. Szegedy, “Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift,” in _International conference on machine learning_. PMLR, 2015, pp. 448–456.
  * [36] A. F. Agarap, “Deep Learning using Rectified Linear Units (ReLu),” _arXiv preprint arXiv:1803.08375_ , 2018.
  * [37] D. P. Kingma and J. Ba, “Adam: A Method for Stochastic Optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [38] E. Nachmani and L. Wolf, “Unsupervised Singing Voice Conversion,” in _Interspeech 2019, 20th Annual Conference of the International Speech_. ISCA, 2019, pp. 2583–2587.
  * [39] P. Chandna, M. Blaauw, J. Bonada, and E. Gómez, “Wgansing: A Multi-Voice Singing Voice Synthesizer Based on the Wasserstein-GAN,” in _2019 27th European Signal Processing Conference (EUSIPCO)_. IEEE, 2019, pp. 1–5.
  * [40] A. Polyak, L. Wolf, Y. Adi, and Y. Taigman, “Unsupervised Cross-Domain Singing Voice Conversion,” in _Interspeech 2020, 21st Annual Conference of the International Speech_. ISCA, 2020, pp. 801–805.


#  $Z_{N}$ symmetry in $SU(N)$ gauge theories

Minati Biswal<EMAIL_ADDRESS>Indian Institute of Science Education and
Research, Mohali 140306, India  Sanatan Digal<EMAIL_ADDRESS>The Institute
of Mathematical Sciences, Chennai 600113, India  Homi Bhabha National
Institute, Training School Complex, Anushakti Nagar, Mumbai 400085, India
Vinod Mamale<EMAIL_ADDRESS>The Institute of Mathematical Sciences,
Chennai 600113, India  Homi Bhabha National Institute, Training School
Complex, Anushakti Nagar, Mumbai 400085, India  Sabiar Shaikh
<EMAIL_ADDRESS>The Institute of Mathematical Sciences, Chennai
600113, India  Homi Bhabha National Institute, Training School Complex,
Anushakti Nagar, Mumbai 400085, India

###### Abstract

We study $Z_{N}$ symmetry in $SU(N)$ gauge theories in the presence of matter
fields in the fundamental representation, by restricting the lattice partition
function integration to matter fields which are uniform in spatial directions
and gauge fields with vanishing spatial components. In this approximation the
gauge matter field interaction effectively reduces to a 1-dimensional gauged
chain. This makes analytical calculations of the matter field contribution to
the Polyakov loop free energy possible. We show that in the limit of large
number of temporal sites the explicit breaking of $Z_{N}$ symmetry in this
free energy vanishes, driven by dominance of the density of states. We argue
that the spatial links as well as the spatial modes of the matter fields
determine the boundaries separating regions where $Z_{N}$ symmetry is realised
from rest of the phase diagram.

††preprint: IMSc/2021/02/02

##  I Introduction

Gauge theories such as quantum chromodynamics (QCD), the standard model (SM)
etc. are crucial for understanding evolutions of the early universe as well as
the quark-gluon plasma (QGP) formed in relativistic heavy-ion collisions.
Studies of phases and phase transitions in these theories are important as the
main aim of experimental observations is to look for their signatures. One of
the phase transitions common to all these theories is the confinement
deconfinement transition(s) (CD) at finite temperatures. In the pure gauge
limit of these theories the CD transition is described by the center $Z_{N}\in
SU(N)$ symmetry and the Polyakov loop which plays the role of an order
parameter Kuti:1980gh ; Creutz:1980zw ; McLerran:1980pk ; Weiss:1980rj ;
Polonyi:1982wz ; Svetitsky:1982gs ; Yaffe:1982qf ; Green:1983sd ; Celik:1983wz
. Above the critical temperature, the Polyakov loop acquires a non-zero
thermal average value which leads to the spontaneous symmetry breaking (SSB)
of the $Z_{N}$ symmetry Weiss:1980rj ; Yaffe:1982qf ; Celik:1983wz ;
Svetitsky:1985ye . The SSB leads to $N$ degenerate states in the deconfined
phase and global topological defects such as strings and domain walls in
physical space Gross:1980br ; Weiss:1981ev ; Balachandran:2001qn ;
Gupta:2010pp ; Ignatius:1991nk .

The $Z_{N}$ symmetry arises from the fact that the allowed gauge
transformations in Euclidean space are periodic up to a factor $z\in Z_{N}$.
These gauge transformations are responsible for the $Z_{N}$ symmetry in pure
gauge theory. In turn the $Z_{N}$ symmetry plays an important role on the
nature of the CD transition Kogut:1979wt . The presence of matter fields in
the above gauge theories spoil the $Z_{N}$ symmetry. The requirement that the
matter fields satisfy either periodic or anti-periodic boundary conditions in
the temporal direction forces the gauge transformations to be periodic
Weiss:1981ev ; Belyaev:1991np ; Green:1983sd ; Karsch:2000zv ; Biswal:2015rul
; Biswal:2016xyq . However gauge field configurations related by $Z_{N}$ gauge
transformations can contribute to the partition function, though not equally.
Since the matter fields cannot be subjected to these gauge transformations the
different configurations will have different actions. This situation appears
similar to Ising model in the presence of external field Banks:1983me . In the
present context, the gauge-matter field interaction plays the role of explicit
$Z_{N}$ symmetry breaking term Banks:1983me ; Green:1983sd .

There are several studies of explicit breaking of $Z_{N}$ symmetry due to
matter fields over the years Weiss:1981ev ; Green:1983sd ; Belyaev:1991np ;
Ignatius:1991nk ; Dixit:1991et ; Deka:2010bc ; Biswal:2019xju . Perturbative
one loop calculations show that the $Z_{N}$ symmetry explicit breaking
increases with decrease in mass of the matter fields Weiss:1981ev . The
explicit breaking is found to increase with temperature. Recently there are
extensions of loop calculations to higher order with similar trend Guo:2018scp
. The non-perturbative studies, which are mostly around the CD transition
regions, show decrease in explicit breaking with the number of temporal
lattice points ($N_{\tau}$) Satz:1985js ; Biswal:2015rul ; Biswal:2016xyq .
These studies show a trend that in the continuum limit likely there will be
reemergence of $Z_{N}$ symmetry.
The exact calculation of the partition function to validate the reemergence of
$Z_{N}$ symmetry is almost impossible. An explanation of this analytically, is
highly desirable even with some simplifications. In this paper, we attempt an
exact calculation of the lattice partition function after restricting the
spatial gauge fields to zero and matter fields uniform in the spatial
directions. With these restrictions only the terms of the action which
explicitly break the $Z_{N}$ gauge transformations remain and the problem
effectively reduces to a 1-dimensional model. The matter fields can then be
integrated out exactly after a suitable gauge choice for the temporal links.
The integrations can be carried out for any arbitrary value of $N_{\tau}$.
From these calculations the free energy for the given gauge field background
or Polyakov loop is obtained. The results show that the explicit symmetry
breaking vanishes in the limit of large $N_{\tau}$ even when the relevant
couplings are finite and the action breaks the $Z_{N}$ symmetry. These results
suggest that for the parameters for which the $Z_{N}$ symmetry realised, the
free energy is dominated by the density of states. Recent studies of
$Z_{2}+$Higgs have shown that the histogram of the explicit breaking term
exhibits $Z_{2}$ symmetry Biswal:2021fde . The modes/aspects of the fields
neglected here will affect the $Z_{N}$ symmetry in parts of the phase diagram
by driving the system away from the point where the density of states dominate
the thermodynamics. This is observed in Monte Carlo simulations of the
partition function with the full action Biswal:2021fde .
The paper is organised as follows. In section II, we discuss the $Z_{N}$
symmetry in the presence of matter fields in the fundamental representation.
In section III, we calculate the partition function for $SU(N)+$Higgs which is
followed by the calculations for $SU(N)+$fermions in section IV. The
discussions and conclusions are presented in section V.

##  II $Z_{N}$ symmetry in $SU(N)$ gauge theories

The action for a minimally coupled $SU(N)$ gauge theory of fermions and bosons
in $3+1$ Euclidean space is given by

$\displaystyle S=\int_{V}d^{3}x\int_{0}^{\beta}d\tau\left[{1\over
2}\left\\{Tr\left(F^{\mu\nu}F_{\mu\nu}\right)+|D_{\mu}\Phi|^{2}+{m_{b}^{2}}\Phi^{\dagger}\Phi\right\\}+{\bar{\Psi}}(\not{D}+m_{f})\Psi\right]$
(1) $\displaystyle
F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}+g[A_{\mu},A_{\nu}],~{}D_{\mu}\Phi=(\partial_{\mu}+igA_{\mu})\Phi,~{}\not{D}\Psi=(\not{\partial}+ig\not{A})\Psi.$

Here $A_{\mu}$, $\Phi$ and $\Psi$ are the gauge, Higgs and the fermion fields
respectively. $\Phi$ and $\Psi$ are in the fundamental representation. $g$ is
the gauge coupling strength, $m_{b}(m_{f})$ is the mass of $\Phi(\Psi)$ fields
and $\beta$ is the inverse of temperature, i.e $\beta=1/T$. The corresponding
partition function takes the form

${\cal
Z}=\int[DA][D\Phi][D\Phi^{\dagger}][D\Psi][D\bar{\Psi}]~{}\rm{Exp}[-S].$ (2)

The fields contributing to the partition function satisfy the following
temporal boundary conditions,

$A_{\mu}(\tau=0)=A_{\mu}(\tau=\beta),~{}\Phi(\tau=0)=\Phi(\tau=\beta),\Psi(\tau=0)=-\Psi(\tau=\beta).$
(3)

These fields transform under a gauge transformation $V\in SU(N)$ as

$A_{\mu}\longrightarrow VA_{\mu}V^{-1}+{1\over
g}(\partial_{\mu}V)V^{-1},~{}\Phi\longrightarrow V\Phi,~{}\Psi\longrightarrow
V\Psi.$ (4)

In absence of the matter fields one can consider $V$ which is not necessarily
periodic, i.e $V(\tau=0)=zV(\tau=\beta)$, where $z$ is an element of the
center $Z_{N}$ of $SU(N)$. While the pure gauge action is invariant under this
transformation, the Polyakov loop,

$L(\vec{x})={1\over N}{\rm{Tr}}\left[{\rm{P}}\left\\{{\rm
Exp}{\left(-ig\int_{0}^{\beta}A_{0}d\tau\right)}\right\\}\right],$ (5)

transforms as $L\longrightarrow zL$. This transformation property of the
Polyakov loop is crucial for it playing the role of an order parameter for the
CD phase transition and also SSB of the $Z_{N}$ symmetry in the deconfined
phase. In the presence of the matter fields ($\Phi,\Psi$), the boundary
conditions in Eq.3, restrict the gauge transformations to be periodic in
$\tau$. Since the non-periodic gauge transformations can not act on the matter
fields, two configurations with Polyakov loops $L(\vec{x})$ and $zL(\vec{x})$
do not necessarily contribute equally to the partition function. Hence the
$Z_{N}$ symmetry is explicitly broken. If $L(\vec{x})$ belongs to the identity
sector of $Z_{N}$ then the configuration corresponding to $zL(\vec{x})$ will
have higher action.

Symmetry in the action automatically leads to symmetry in the free energy.
However explicit breaking at the level of action does not necessarily mean the
same is true at the level of free energy. This is because the free energy
difference between $L(\vec{x})$ and $zL(\vec{x})$ can only be decided after
the matter fields are integrated out. Integrating the matter field is very
difficult task in four dimensions. Therefore in the following we consider
$SU(N)$ gauge theory of Higgs and fermions separately by restricting to the
spatial gauge fields and spatial variations of the matter fields to zero. As
mentioned previously in this approximation the part of the action which breaks
the $Z_{N}$ symmetry remains, while the system effectively reduces to
$1-$dimensional model of gauged chain, making analytical calculations
possible.

##  III Gauged $1-d$ chain of $SU(N)+$Higgs

The $SU(N)+$Higgs part of Eq.1 on the Euclidean lattice can be written as
Biswal:2016xyq ,

$S=\beta_{g}\sum_{p}\left[1-{1\over
2}(U_{p}+U^{\dagger}_{p})\right]-b\sum_{n,\mu}(\Phi_{n}^{\dagger}U_{n,\mu}\Phi_{n+\hat{\mu}}+h.c)+a\sum_{n}\Phi^{\dagger}_{n}\Phi_{n}.$
(6)

$\beta_{g}$ is the gauge coupling constant, $a={1\over 2}$ and the coupling
$b=(m_{b}^{2}+8)^{-1}$, the Higgs mass $m_{b}$ is expressed in lattice units
Biswal:2016xyq . For unit spatial links and $\Phi$ uniform in the spatial
directions, the action reduces to,

$S=a\sum_{i=1}^{N_{\tau}}\Phi_{i}^{\dagger}\Phi_{i}-b\sum_{i=1}^{N_{\tau}}(\Phi_{i}^{\dagger}U_{i}\Phi_{i+1}+h.c.),$
(7)

apart from an overall spatial volume factor. The pure gauge part is dropped as
the effect of matter fields on the $Z_{N}$ symmetry in the back ground of
temporal gauge links is being considered. For convenience the subscripts of
the field variables have been replaced by $i$. $N_{\tau}$ denotes the number
of temporal sites. $\Phi$ satisfies periodic boundary condition, i.e
$\Phi_{N_{\tau}+1}=\Phi_{1}$. We consider a gauge choice in which
$U_{i}=\mathbb{I}$ for $i=1,2,...,N_{\tau}-1$ and $U_{N_{\tau}}=U$. The
Polyakov loop is $L=Tr(U)/N$. In order to derive the free energy $V(L)$, only
the $\Phi_{i}$ fields in the partition function ${\cal{Z}}_{L}$ are to be
integrated out,

${\cal{Z}}_{L}=\int\prod_{i=1}^{N_{\tau}}d\Phi_{i}^{\dagger}d\Phi_{i}~{}{\rm
Exp}[-S].$ (8)

For convenience the action is written as $S=S_{1}+S_{2}$ as in the following,

$S_{1}=a\Phi_{1}^{\dagger}\Phi_{1}-b\left(\Phi_{N_{\tau}}^{\dagger}U\Phi_{1}+h.c.\right),~{}S_{2}=a\sum_{i=2}^{N_{\tau}}\Phi_{i}^{\dagger}\Phi_{i}-b\sum_{i=1}^{N_{\tau}-1}(\Phi_{i}^{\dagger}\Phi_{i+1}+h.c.).$
(9)

At first, the fields $\Phi_{2}$ to $\Phi_{N_{\tau}-1}$ are integrated out
sequentially, i.e,

${\cal{Z}}=\int\prod_{i=2}^{N_{\tau}-1}d\Phi_{i}^{\dagger}d\Phi_{i}~{}{\rm
Exp}[-S_{2}]$ (10)

Afterwards the remaining integration of $\Phi_{1}$ and $\Phi_{N_{\tau}}$ can
be carried out to obtain the partition function,

${\cal{Z}}_{L}=\int
d\Phi_{1}^{\dagger}d\Phi_{1}d\Phi^{\dagger}_{N_{\tau}}d\Phi_{N_{\tau}}\left({\cal{Z}}\times{\rm
Exp}[-S_{1}]\right).$ (11)

The integration of $\Phi_{1}$ and $\Phi_{N_{\tau}}$ requires evaluation of
determinant of a matrix of size $4N\times 4N$. The integrations of $\Phi_{2}$
to $\Phi_{N_{\tau}-1}$ greatly simplify the problem, otherwise one would have
to deal with evaluating of matrix whose size depends on $N_{\tau}$.

In the integration of ${\cal{Z}}$ in Eq.10, due to the gauge choice mentioned
above the different components as well as the real and imaginary parts of
$\Phi_{i}$’s do not mix. Therefore it can be written as,

${\cal{Z}}=\prod_{r=1}^{2N}{\cal{I}}(\Phi_{N_{\tau},r}),$ (12)

where $\Phi_{{}_{N_{\tau},r}}$ is the r-th component of $\Phi_{N_{\tau}}$ and
${\cal{I}}(\Phi_{N_{\tau},r})$ is obtained by integrating out r-th component
of $\Phi_{2}$ to $\Phi_{N_{\tau}-1}$. Denoting the r-th component by $\phi$ we
can write

${\cal{I}}(\phi_{N_{\tau}})=\int\prod_{i=2}^{N_{\tau}-1}d\phi_{i}~{}{\rm
Exp}[-S_{2}^{\prime}],$ (13)

where

$S_{2}^{\prime}=a\sum_{i=2}^{N_{\tau}-1}\phi_{i}^{2}-2b\sum_{i=1}^{N_{\tau}-1}\phi_{i}\phi_{i+1}.$
(14)

The integration ${\cal{I}}$, in Eq.13, can be also be written as,

${\cal{I}}(\phi_{N_{\tau}})=\int\prod_{i=3}^{N_{\tau}-1}d\phi_{i}e^{-S_{3}^{\prime}}\int
d\phi_{2}~{}{\rm
Exp}\left[-a\phi_{2}^{2}+2\phi_{2}\left(b\phi_{1}+b\phi_{3}\right)\right]$
(15)

$S_{3}^{\prime}$ is obtained by taking out terms from $S_{2}^{\prime}$ which
are dependent on $\phi_{2}$. After $\phi_{2}$ is integrated out,

${\cal{I}}(\phi_{N_{\tau}})=\int\prod_{i=3}^{N_{\tau}-1}d\phi_{i}e^{-S^{\prime}_{3}}\sqrt{\pi\over
a}~{}{\rm Exp}\left[{1\over a}\left(b\phi_{1}+b\phi_{3}\right)^{2}\right]$
(16)

which can also be written as

$\displaystyle{\cal{I}}(\phi_{N_{\tau}})$ $\displaystyle=$
$\displaystyle\int\prod_{i=4}^{N_{\tau}-1}d\phi_{i}e^{-S^{\prime}_{4}}$ (17)
$\displaystyle\times$ $\displaystyle\sqrt{\pi\over a}\int d\phi_{3}~{}{\rm
Exp}\left[-\left(a-{b^{2}\over
a}\right)\phi_{3}^{2}+2\phi_{3}\left({b^{2}\over
a}\phi_{1}+b\phi_{4}\right)+{b^{2}\over a}\phi_{1}^{2}\right].$

Here again $S_{4}^{\prime}$ is $S_{2}^{\prime}$ without terms containing
$\phi_{2}$ and $\phi_{3}$. Given the forms of ${\cal{I}}(\phi_{N_{\tau}})$ in
Eq.15 and Eq.17 one easily write down the would be form of
${\cal{I}}(\phi_{N_{\tau}})$ after integration of $\phi_{k-1}$ as,

$\displaystyle{\cal{I}}(\phi_{N_{\tau}})$ $\displaystyle=$ $\displaystyle\int
d\phi_{k+1}....d\phi_{N_{\tau}-1}e^{-S^{\prime}_{k+1}}$ (18)
$\displaystyle\times$ $\displaystyle I_{k}\int d\phi_{k}~{}{\rm
Exp}\left[-A_{k}\phi_{k}^{2}+2\phi_{k}\left(B_{k}\phi_{1}+b\phi_{k+1}\right)+E_{k}\phi_{1}^{2}\right].$

Carrying out the $\phi_{k}$ integration will result in,

$\displaystyle{\cal{I}}(\phi_{N_{\tau}})$ $\displaystyle=$ $\displaystyle\int
d\phi_{k+2}....d\phi_{N_{\tau}-1}e^{-S^{\prime}_{k+2}}$ (19)
$\displaystyle\times$ $\displaystyle I_{k+1}\int d\phi_{k+1}~{}{\rm
Exp}\left[-A_{k+1}\phi_{k+1}^{2}+2\phi_{k+1}\left(B_{k+1}\phi_{1}+b\phi_{k+2}\right)+E_{k+1}\phi_{1}^{2}\right].$

From equations 18 and 19, one can read off the following recursion relations,

$I_{k+1}=\sqrt{\pi\over A_{k}},~{}A_{k+1}=a-{b^{2}\over
A_{k}},~{}B_{k+1}={bB_{k}\over A_{k}},~{}E_{k+1}=E_{k}+{B_{k}^{2}\over
A_{k}},$ (20)

with $I_{2}=1$, $A_{2}=a$, $B_{2}=b$ and $E_{2}=0$. Using these recursion
relations we can workout the integration, ${\cal{I}}(\phi_{N_{\tau}})$
completely. Using ${\cal{I}}(\phi_{N_{\tau}})$’s one can write the partition
function as,

$\displaystyle{\cal{Z}}_{L}$ $\displaystyle=$ $\displaystyle Q\int
d\Phi_{1}^{\dagger}d\Phi_{1}d\Phi_{N_{\tau}}^{\dagger}d\Phi_{N_{\tau}}$ (21)
$\displaystyle~{}\rm
Exp\left[-A_{N_{\tau}}\Phi_{N_{\tau}}^{\dagger}\Phi_{N_{\tau}}-C_{N_{\tau}}\Phi_{1}^{\dagger}\Phi_{1}+\left(\Phi_{N_{\tau}}^{\dagger}(B_{N_{\tau}}\mathbb{I}+bU)\Phi_{1}+H.C.\right)\right],$

where

$Q=\prod_{k=2}^{N_{\tau}}I_{k}^{{}^{n}},~{}n=2N.$ (22)

$n$ corresponds to the number of components of the $\Phi$ field. The
coefficient $C_{N_{\tau}}=a-E_{N_{\tau}}$. After the integration of the
remaining fields $\Phi_{1}$ and $\Phi_{N_{\tau}}$ the partition function takes
the form,

${\cal{Z}}_{L}=Q\sqrt{\pi^{8}\over Det(M)}.$ (23)

$M$ is $(4N\times 4N)$ given by,

$\begin{pmatrix}A_{N_{\tau}}&B_{N_{\tau}}+bU\\\
B_{N_{\tau}}+bU^{\dagger}&C_{N_{\tau}}\end{pmatrix}$

The exact form of $Det(M)$ for arbitrary $N$ is difficult to find. However the
sequential integration has greatly simplified the problem. For arbitrary
$N_{\tau}$ we need to deal with a matrix of finite size. In the following we
consider $N=2$ and evaluate ${\cal{Z}}_{L}$ explicitly for an arbitrary $U$,
which can be parametrised as,

$U=\alpha_{0}+i{\bf\alpha.\sigma},~{}\alpha=(\alpha_{1},\alpha_{2},\alpha_{3}),$
(24)

where $\sigma_{i}$’s are the Pauli matrices. The corresponding matrix $M$ is
given by,

$\begin{pmatrix}A_{N_{\tau}}&0&0&0&B_{1}&b\alpha_{3}&-b\alpha_{2}&b\alpha_{1}\\\
0&A_{N_{\tau}}&0&0&-b\alpha_{3}&B_{1}&-b\alpha_{1}&-b\alpha_{2}\\\
0&0&A_{N_{\tau}}&0&b\alpha_{2}&b\alpha_{1}&B_{1}&-b\alpha_{3}\\\
0&0&0&A_{N_{\tau}}&-b\alpha_{1}&b\alpha_{2}&b\alpha_{3}&B_{1}\\\
B_{1}&-b\alpha_{3}&b\alpha_{2}&-b\alpha_{1}&C_{N_{\tau}}&0&0&0\\\
b\alpha_{3}&B_{1}&b\alpha_{1}&b\alpha_{2}&0&C_{N_{\tau}}&0&0\\\
-b\alpha_{2}&-b\alpha_{1}&B_{1}&b\alpha_{3}&0&0&C_{N_{\tau}}&0\\\
b\alpha_{1}&-b\alpha_{2}&-b\alpha_{3}&B_{1}&0&0&0&C_{N_{\tau}}\end{pmatrix}$,

where $B_{1}=-(b\alpha_{0}+B_{N_{\tau}})$. The determinant of $M$ is given by,

$Det{M}=\left(B^{2}_{N_{\tau}}-A_{N_{\tau}}C_{N_{\tau}}+2bB_{N_{\tau}}\alpha_{0}+b^{2}\right)^{4}$
(25)

$Z_{2}$ rotation of the Polyakov loop changes $\alpha_{0}\to-\alpha_{0}$. So
in the determinant the explicit symmetry breaking of $Z_{2}$ is
$2bB_{N_{\tau}}\alpha_{0}$. It is observed that $B_{N_{\tau}}$ rapidly
decreases, vanishing in the larger $N_{\tau}$ limit restoring the $Z_{2}$
symmetry. Even for higher $N$ one can see the realisation of $Z_{N}$ symmetry
as the off diagonal elements of the matrix $M$ turn out to be just $U$ and
$U^{\dagger}$ due to vanishing of $B_{N_{\tau}}$. Effecting a $Z_{N}$
transformation, ie $U\to zU$, the factor $z$ in $U$ will cancel with $z^{*}$
in $U^{\dagger}$ leaving the determinant unchanged.

From the above results, it can be shown that there is realisation of $Z_{N}$
symmetry for fixed temperature $T$ and the Higgs mass. The number of temporal
sites $N_{\tau}$ and $T$ are related by $T=1/(a_{l}N_{\tau})$, where $a_{l}$
is the lattice spacing . Fixing $T$ amounts to $a_{l}\propto 1/N_{\tau}$.
$a_{l}$ enters into the calculations through the parameter
$b=(m_{b}^{2}+8)^{-1}$ which depends on the Higgs mass, $m_{b}$, expressed in
lattice units. For a fixed Higgs physical mass, the corresponding mass in
lattice units must decrease, i.e $m_{b}\propto a_{l}$. Consequently, in the
$N_{\tau}\to\infty$ limit the dimensionless parameter $b$ increases to
$b=1/8$. This increase, however, does not affect the realisation of $Z_{N}$
symmetry, since $bB_{N_{\tau}}$(Eq.25) vanishes in the $N_{\tau}\to\infty$
limit.

In the following we consider the effects of staggered fermion fields on the
$Z_{N}$ symmetry.

##  IV  Gauged $1-d$ chain of $SU(N)+$fermions

The lattice action for $SU(N)$ staggered fermions is given by Kilcup:1986dg

$S=\beta_{g}\sum_{p}\left[1-{1\over
2}(U_{p}+U^{\dagger}_{p})\right]+2m_{f}\sum_{n}\bar{\Psi}_{n}\Psi_{n}+\sum_{n,\mu}\eta_{n,\mu}\left[\bar{\Psi}_{n}U_{n,\mu}\Psi_{n+\mu}-\bar{\Psi}_{n}U^{\dagger}_{n-\mu,\mu}\Psi_{n-\mu}\right]$
(26)

Here the fermion mass as well as the fields are expressed in lattice unit. The
analog of Eq.7 in this case turns out to be,

$S=2m_{f}\sum_{i=1}^{N_{\tau}}\bar{\Psi}_{i}\Psi_{i}+\sum_{i=1}^{N_{\tau}-1}\left(\bar{\Psi}_{i}\Psi_{i+1}-\bar{\Psi}_{i+1}\Psi_{i}\right)-{\bar{\Psi}}_{N_{\tau}}U\Psi_{1}+\bar{\Psi}_{1}U^{\dagger}\Psi_{N_{\tau}},$
(27)

The change in the sign of the last two terms is due to the anti-periodicity of
$\Psi$. Here we have considered the $KS$ phase $\eta_{0}$ to be $+1$
Susskind:1976jm ; KlubergStern:1983dg , however the results/conclusions do not
depend on $\eta_{0}$. As in the previous section we work in the gauge in which
all temporal links except the last one are set to identity. The last link is
denoted by $U$. The corresponding Polyakov loop is $L=Tr(U)/N$. To find out
the free energy $V(L)$ we need to integrate out only the fermion fields. For
convenience we write $S=S_{1}+S_{2}$ where

$\displaystyle
S_{1}=2m_{f}\bar{\Psi}_{1}\Psi_{1}-{\bar{\Psi}}_{N_{\tau}}U\Psi_{1}+\bar{\Psi}_{1}U^{\dagger}\Psi_{N_{\tau}},$
(28) $\displaystyle
S_{2}=2m_{f}\sum_{i=2}^{N_{\tau}}\bar{\Psi}_{i}\Psi_{i}+\sum_{i=1}^{N_{\tau}-1}\left(\bar{\Psi}_{i}\Psi_{i+1}-\bar{\Psi}_{i+1}\Psi_{i}\right).$
(29)

Initially we integrate the fields $\Psi_{2}$, $\bar{\Psi}_{2}$ to
$\Psi_{N_{\tau}-1}$, $\bar{\Psi}_{N_{\tau}-1}$ sequentially just as in the
previous section. Afterwards $\Psi_{1}$, $\bar{\Psi}_{1}$ and
$\Psi_{N_{\tau}}$, $\bar{\Psi}_{N_{\tau}}$ are integrated out to obtain the
partition function,

${\cal{Z}}_{L}=\int
d\bar{\Psi}_{1}d\bar{\Psi}_{N_{\tau}}d{\Psi}_{1}d\Psi_{N_{\tau}}~{}{\rm
Exp}[-S_{1}]{\cal{Z}},$ (30)

where ${\cal{Z}}$ is given by

${\cal{Z}}=\int\prod_{i=2}^{N_{\tau}-1}d\bar{\Psi}_{i}d{\Psi}_{i}~{}{\rm
Exp}[-S_{2}].$ (31)

Since $S_{2}$ is diagonal in colour space we consider a particular colour of
$\Psi_{i}$ and denote it by $\psi_{i}$. For this choice the relevant integral
is,

${\cal{I}}_{\psi}=\int\prod_{i=2}^{N_{\tau}-1}d\bar{\psi}_{i}d{\psi}_{i}~{}{\rm
Exp}[-S_{2,\psi}].$ (32)

After integrating $\psi_{2}$, $\bar{\psi}_{2}$ and $\psi_{3}$,
$\bar{\psi}_{3}$ the integral takes the form,

$\displaystyle{\cal{I}}_{\psi}$
$\displaystyle=\int\prod_{i=4}^{N_{\tau}-1}d\bar{\psi}_{i}d{\psi}_{i}~{}{\rm
Exp}[-S^{4}_{2,\psi}]\times$ (33)
$\displaystyle\left[1+4m_{f}^{2}-2m_{f}\bar{\psi}_{1}\psi_{1}-2m_{f}\bar{\psi}_{4}\psi_{4}+\bar{\psi}_{4}\psi_{1}-\bar{\psi}_{1}\psi_{4}+\bar{\psi}_{4}\psi_{4}\bar{\psi}_{1}\psi_{1}\right]$

$S^{4}_{2,\psi}$ is obtained by dropping terms which depend on $\psi_{2}$,
$\bar{\psi}_{2}$ and $\psi_{3}$, $\bar{\psi}_{3}$. The sequential integration
$\psi_{4}$ up to $\psi_{N_{\tau}-1}$ and their conjugates leads to

$\displaystyle{\cal{I}}_{\psi}=A_{N_{\tau}}-B_{N_{\tau}}\bar{\psi}_{1}\psi_{1}-C_{N_{\tau}}\bar{\psi}_{N_{\tau}}\psi_{N_{\tau}}+\bar{\psi}_{N_{\tau}}\psi_{1}+D_{N_{\tau}}\bar{\psi_{1}}\psi_{N_{\tau}}+E_{N_{\tau}}\bar{\psi}_{N_{\tau}}\psi_{N_{\tau}}\bar{\psi}_{1}\psi_{1}$
(34)

where the coefficients $A_{N_{\tau}}$ to $E_{N_{\tau}}$ can be obtained by
recursion as

$A_{k+1}=2m_{f}A_{k}+C_{k},~{}B_{k+1}=2m_{f}B_{k}+E_{k},~{}C_{k+1}=A_{k},~{}D_{k+1}=(-1)^{k},~{}E_{k+1}=B_{k},$
(35)

with

$A_{4}=(1+4m_{f}^{2}),~{}B_{4}=2m_{f},~{}C_{4}=2m_{f},~{}E_{4}=1,$ (36)

Taking the ${\cal{I}}$ integrals into account we can write the partition
function as

$\displaystyle{\cal{Z}}_{L}=\int
d\bar{\psi}_{1}d{\psi}_{1}d\bar{\psi}_{N_{\tau}}d{\psi}_{N_{\tau}}~{}\rm{Exp}\left[\bar{\psi}_{N_{\tau}}U\psi_{1}-\bar{\psi}_{1}U^{\dagger}\psi_{N_{\tau}}\right]\times$
$\displaystyle\prod_{r}\left(1-2m_{f}\bar{\psi}^{r}_{1}\psi^{r}_{1}-2m_{f}\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{N_{\tau}}+4m_{f}^{2}\bar{\psi}^{r}_{1}\psi_{1}\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{N_{\tau}}\right)\times$
$\displaystyle\left(A_{N_{\tau}}-B_{N_{\tau}}\bar{\psi}^{r}_{1}\psi^{r}_{1}-C_{N_{\tau}}\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{N_{\tau}}+\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{1}+D_{N_{\tau}}\bar{\psi}^{r}_{1}\psi^{r}_{N_{\tau}}+E_{N_{\tau}}\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{N_{\tau}}\bar{\psi}^{r}_{1}\psi^{r}_{1}\right).$
(37)

Note that $\psi_{i}^{r}$ denotes the colour $r$ of the field $\Psi_{i}$ at the
temporal site $i$. This expression can be simplified as,

$\displaystyle{\cal{Z}}_{L}=\int
d\bar{\psi}_{1}d{\psi}_{1}d\bar{\psi}_{N_{\tau}}d{\psi}_{N_{\tau}}~{}\rm{Exp}\left[\bar{\psi}_{N_{\tau}}U\psi_{1}-\bar{\psi}_{1}U^{\dagger}\psi_{N_{\tau}}\right]\times$
$\displaystyle\prod_{r}\left(\tilde{A}-\tilde{B}\bar{\psi}^{r}_{1}\psi^{r}_{1}-\tilde{C}\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{N_{\tau}}+\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{1}+\tilde{D}\bar{\psi}^{r}_{1}\psi^{r}_{N_{\tau}}+\tilde{E}\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{N_{\tau}}\bar{\psi}^{r}_{1}\psi^{r}_{1}\right),$
(38)

where $\tilde{A}=A_{N_{\tau}}$, $\tilde{B}=(2m_{f}A_{N_{\tau}}+B_{N_{\tau}})$,
$\tilde{C}=(2m_{f}A_{N_{\tau}}+C_{N_{\tau}})$, $\tilde{D}=D_{N_{\tau}}$ and
$\tilde{E}=E_{N_{\tau}}+2m_{f}C_{N_{\tau}}+2m_{f}B_{N_{\tau}}+4m^{2}_{f}A_{N_{\tau}}$.
The superscript $r$ denotes the colour of the fermion field. For large
$m_{f}$, in ${\cal Z}_{L}$ the ratio of the leading and sub-leading term
scales as $\sim m_{f}$. The sub-leading terms contain matrix elements of $U$.
In the limit of $m_{f}\rightarrow\infty$, therefore the fermions decouple from
the gluons and exact $Z_{N}$ symmetry is recovered for any $N_{\tau}$. For
$N=2$, integration of the rest of the fields in Eq.38 leads to,

$\displaystyle{\cal{Z}}_{L}$
$\displaystyle=\tilde{E}^{2}+2\tilde{E}\tilde{A}|U_{11}|^{2}+\tilde{A}^{2}+2\tilde{B}\tilde{C}|U_{12}|^{2}+2(1-\tilde{D}Re(U^{2}_{11}))$
(39) $\displaystyle+(\tilde{E}+\tilde{A})(1-\tilde{D})tr(U).$

As one can see the $Z_{2}$ explicit breaking term is linear in
$\tilde{E}+\tilde{A}$. For non zero $m_{f}$, in the free energy $V(L)$ the
first four terms of ${\cal{Z}}_{L}$ dominate over $\tilde{E}+\tilde{A}$. The
dominance only grows with $N_{\tau}$, hence in the limit of large $N_{\tau}$
the $Z_{2}$ symmetry is recovered. For higher $N$ it is difficult to evaluate
${\cal{Z}}_{L}$ for a general $U$. To proceed further we assume the $U$ to be
$U_{rs}=\lambda_{r}\delta_{rs}$. After the exponential in Eq.38 is written as
a polynomial,

$\displaystyle{\cal{Z}}_{L}$ $\displaystyle=\int
d\bar{\psi}_{1}d{\psi}_{1}d\bar{\psi}_{N_{\tau}}d{\psi}_{N_{\tau}}\prod_{r}\left(1+\lambda_{r}\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{1}-\lambda^{*}_{r}\bar{\psi}^{r}_{1}\psi^{r}_{N_{\tau}}+\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{N_{\tau}}\bar{\psi}^{r}_{1}\psi^{r}_{1}\right)\times$
(41)
$\displaystyle\left(\tilde{A}-\tilde{B}\bar{\psi}^{r}_{1}\psi^{r}_{1}-\tilde{C}\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{N_{\tau}}+\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{1}+\tilde{D}\bar{\psi}^{r}_{1}\psi^{r}_{N_{\tau}}+\tilde{E}\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{N_{\tau}}\bar{\psi}^{r}_{1}\psi^{r}_{1}\right)$
$\displaystyle=\int
d\bar{\psi}_{1}d{\psi}_{1}d\bar{\psi}_{N_{\tau}}d{\psi}_{N_{\tau}}\times$
$\displaystyle\prod_{r}(A-B\bar{\psi}^{r}_{1}\psi^{r}_{1}-C\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{N_{\tau}}+F_{r}\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{1}+D_{r}\bar{\psi}^{r}_{1}\psi^{r}_{N_{\tau}}+E_{r}\bar{\psi}^{r}_{N_{\tau}}\psi^{r}_{N_{\tau}}\bar{\psi}^{r}_{1}\psi^{r}_{1}),$

where $A=\tilde{A}$, $B=\tilde{B}$, $C=\tilde{C}$,
$D_{r}=\tilde{D}-\lambda_{r}^{*}\tilde{A}$,
$E_{r}=\tilde{E}-\lambda_{r}\tilde{D}+\lambda^{*}_{r}+\tilde{A}$ and
$F_{r}=(1+\lambda_{r}\tilde{A})$. After the fields are integrated out we get
the following result for the partition function,

${\cal{Z}}_{L}=\prod_{r}E_{r}$ (42)

The corresponding free energy is

$V(L)=-T\sum_{r}\left\\{\rm{log}\left(\tilde{E}+\tilde{A}\right)+\rm{log}\left(1-{\lambda_{r}\tilde{D}-\lambda^{*}_{r}\over\tilde{E}+\tilde{A}}\right)\right\\}.$
(43)

For non-zero $m_{f}$, the second term vanishes in the limit of large
$N_{\tau}$. Hence, in this limit, the free energy $V(L)$ is independent of $L$
leading to realisation of the $Z_{N}$ symmetry. The form of $U$ considered
above includes $\lambda_{r}=\lambda$ for all $r$, with
$\lambda=~{}\rm{Exp}(i2\pi n/N)$ with $n=0,1,2,...N-1$. We mention here that
for this case one would have expected the explicit breaking of $Z_{N}$ to be
maximal.

To see the realisation of $Z_{N}$ symmetry for a fixed temperature and
physical fermion mass, the behaviour of $\tilde{E}$ and $\tilde{A}$ must be
studied in the limit $N_{\tau}\to\infty$ while scaling the fermion mass in
lattice units as $m_{f}\propto 1/N_{\tau}$. Unlike in the case of bosons, it
is not possible to carry this out analytically as the polynomial coefficients
in $\tilde{E}$ and $\tilde{A}$ change with $N_{\tau}$. We have numerically
checked that $\tilde{E}$ and $\tilde{A}$ monotonically increase with
$N_{\tau}$ even when the lattice fermion mass scales as $m_{f}\propto
1/N_{\tau}$. The increase though is slower compared to the case when $m_{f}$
is held fixed. This suggests that even for a fixed temperature and non-zero
physical fermion mass there will be realisation of $Z_{N}$ symmetry.

##  V Discussions and conclusions

In this paper we report on the explicit breaking of $Z_{N}$ symmetry in the
presence of bosons and fermions. We show that analytical treatment of the
problem is possible by simplifying the functional integral where the spatial
links are set to unity and matter fields uniform in spatial direction. In this
simplification most of the terms of the original action drop out except for
the ones which break the $Z_{N}$ symmetry. Also the problem reduces to
1-dimensional chain of gauged bosons/fermions making analytical calculations
possible. To derive the free energy $V(L)$ for the Polyakov loop $L$, the
partition function is evaluated for a given background of temporal gauge
links. The calculations become simple in the gauge where we can set all the
links except the last one to unity. Then the matter fields are integrated out
sequentially except for the two fields connected to the last gauge link. The
integration of the last two fields result in determinant of a finite sized
matrix for arbitrary $N_{\tau}$. In the case of Higgs, the $Z_{N}$ symmetry is
realised in the partition function for $N_{\tau}\to\infty$. In the same limit,
for fermions the explicit breaking terms drop out when the free energy $V(L)$
is calculated. The vanishingly small explicit breaking of $Z_{N}$ for
$N_{\tau}\to\infty$ can be attributed to the dominance of the density of the
states over the action. In this limit, the density of states is found to have
the $Z_{N}$ symmetry Biswal:2021fde .

The present calculations leave out the effect of the spatial links and non-
zero spatial modes of the matter fields. Consideration of these modes will
require that spatial interaction terms be included in the action, though these
terms are not directly responsible for the explicit breaking $Z_{N}$ symmetry.
It appears that these terms decide the onset of the explicit breaking via the
Higgs and the chiral transitions, which are entropy driven. It is expected
that in the phase diagram where the action dominates over the entropy or the
density of states, the $Z_{N}$ symmetry will be explicitly broken. Recent
Monte Carlo simulations in the presence of Higgs show that this is indeed the
case Biswal:2015rul ; Biswal:2021fde . The $Z_{N}$ symmetry is explicitly
broken in the Higgs broken phase even for large $N_{\tau}$.

###### Acknowledgements.

We thank A. P. Balachandran, S. Datta and S. Sharma for valuable discussions
and suggestions.

REFERENCES

## References

  * (1) J. Kuti, J. Polonyi and K. Szlachanyi, doi:10.1016/0370-2693(81)90987-4
  * (2) M. Creutz, Phys. Rev. D 21, 2308-2315 (1980) doi:10.1103/PhysRevD.21.2308
  * (3) L. D. McLerran and B. Svetitsky, doi:10.1016/0370-2693(81)90986-2
  * (4) N. Weiss, Phys. Rev. D 24, 475 (1981) doi:10.1103/PhysRevD.24.475
  * (5) J. Polonyi and K. Szlachanyi, Phys. Lett. B 110, 395-398 (1982) doi:10.1016/0370-2693(82)91280-1
  * (6) B. Svetitsky and L. G. Yaffe, Nucl. Phys. B 210, 423-447 (1982) doi:10.1016/0550-3213(82)90172-9
  * (7) L. G. Yaffe and B. Svetitsky, Phys. Rev. D 26, 963 (1982) doi:10.1103/PhysRevD.26.963
  * (8) F. Green and F. Karsch, Nucl. Phys. B 238, 297-306 (1984) doi:10.1016/0550-3213(84)90452-8
  * (9) T. Celik, J. Engels and H. Satz, Phys. Lett. B 125, 411-414 (1983) doi:10.1016/0370-2693(83)91314-X
  * (10) B. Svetitsky, Phys. Rept. 132, 1-53 (1986) doi:10.1016/0370-1573(86)90014-1
  * (11) D. J. Gross, R. D. Pisarski and L. G. Yaffe, Rev. Mod. Phys. 53, 43 (1981) doi:10.1103/RevModPhys.53.43
  * (12) N. Weiss, Phys. Rev. D 25, 2667 (1982) doi:10.1103/PhysRevD.25.2667
  * (13) A. P. Balachandran and S. Digal, Int. J. Mod. Phys. A 17, 1149-1158 (2002) doi:10.1142/S0217751X02005864 [arXiv:hep-ph/0108086 [hep-ph]].
  * (14) U. S. Gupta, R. K. Mohapatra, A. M. Srivastava and V. K. Tiwari, Phys. Rev. D 82, 074020 (2010) doi:10.1103/PhysRevD.82.074020 [arXiv:1007.5001 [hep-ph]].
  * (15) J. Ignatius, K. Kajantie and K. Rummukainen, Phys. Rev. Lett. 68, 737-740 (1992) doi:10.1103/PhysRevLett.68.737
  * (16) J. B. Kogut, Rev. Mod. Phys. 51, 659 (1979) doi:10.1103/RevModPhys.51.659
  * (17) V. M. Belyaev, I. I. Kogan, G. W. Semenoff and N. Weiss, Phys. Lett. B 277, 331-336 (1992) doi:10.1016/0370-2693(92)90754-R
  * (18) F. Karsch, E. Laermann, A. Peikert, C. Schmidt and S. Stickan, Nucl. Phys. B Proc. Suppl. 94, 411-414 (2001) doi:10.1016/S0920-5632(01)00988-4 [arXiv:hep-lat/0010040 [hep-lat]].
  * (19) M. Biswal, S. Digal and P. S. Saumia, Nucl. Phys. B 910, 30-39 (2016) doi:10.1016/j.nuclphysb.2016.06.025 [arXiv:1511.08295 [hep-lat]].
  * (20) M. Biswal, M. Deka, S. Digal and P. S. Saumia, Phys. Rev. D 96, no.1, 014503 (2017) doi:10.1103/PhysRevD.96.014503 [arXiv:1610.08265 [hep-lat]].
  * (21) T. Banks and A. Ukawa, Nucl. Phys. B 225, 145-155 (1983) doi:10.1016/0550-3213(83)90016-0
  * (22) V. Dixit and M. C. Ogilvie, Phys. Lett. B 269, 353-356 (1991) doi:10.1016/0370-2693(91)90183-Q
  * (23) M. Deka, S. Digal and A. P. Mishra, Phys. Rev. D 85, 114505 (2012) doi:10.1103/PhysRevD.85.114505 [arXiv:1009.0739 [hep-lat]].
  * (24) M. Biswal, S. Digal and P. S. Saumia, Phys. Rev. D 102, no.7, 074020 (2020) doi:10.1103/PhysRevD.102.074020 [arXiv:1907.07981 [hep-ph]].
  * (25) Y. Guo and Q. Du, JHEP 05, 042 (2019) doi:10.1007/JHEP05(2019)042 [arXiv:1810.13090 [hep-ph]].
  * (26) H. Satz, Phys. Lett. B 157, 65-69 (1985) doi:10.1016/0370-2693(85)91213-4
  * (27) M. Biswal, S. Digal, V. Mamale and S. Shaikh, [arXiv:2102.11091 [hep-lat]].
  * (28) G. W. Kilcup and S. R. Sharpe, Nucl. Phys. B 283, 493-550 (1987) doi:10.1016/0550-3213(87)90285-9
  * (29) H. Kluberg-Stern, A. Morel, O. Napoly and B. Petersson, Nucl. Phys. B 220, 447-470 (1983) doi:10.1016/0550-3213(83)90501-1
  * (30) L. Susskind, Phys. Rev. D 16, 3031-3039 (1977) doi:10.1103/PhysRevD.16.3031