IEEEexample:BSTcontrol

# Low-Complexity Models for Acoustic Scene Classification Based on Receptive
Field Regularization and Frequency Damping

###### Abstract

Deep Neural Networks are known to be very demanding in terms of computing and
memory requirements. Due to the ever increasing use of embedded systems and
mobile devices with a limited resource budget, designing low-complexity models
without sacrificing too much of their predictive performance gained great
importance. In this work, we investigate and compare several well-known
methods to reduce the number of parameters in neural networks. We further put
these into the context of a recent study on the effect of the Receptive Field
(RF) on a model’s performance, and empirically show that we can achieve high-
performing low-complexity models by applying specific restrictions on the RFs,
in combination with parameter reduction methods. Additionally, we propose a
filter-damping technique for regularizing the RF of models, without altering
their architecture and changing their parameter counts. We will show that
incorporating this technique improves the performance in various low-
complexity settings such as pruning and decomposed convolution. Using our
proposed filter damping, we achieved the 1st rank at the DCASE-2020 Challenge
in the task of Low-Complexity Acoustic Scene Classification.111Code available
at: https://github.com/kkoutini/cpjku_dcase20

Index Terms—  low-complexity, acoustic scene classification, receptive-field
regularization, pruning, network decomposition

##  1 Introduction

The recent advances in machine learning have been mainly due to the
unprecedented successes of deep neural networks with millions or even billions
of trainable parameters that can learn from a large amount of data and solve
complex problems. Although in deep learning, the main attention has been on
reaching the highest performance, as these models continue to expand their
applications from research into industry, the memory efficiency, energy
consumption, and latency of these models become more and more important.

To address these problems, a new line of research has been established to
design low-complexity neural networks that are capable of reaching the
performance of the large models, while having many orders of magnitude fewer
parameters. In this area, three main approaches have been followed in the
literature. A first approach is knowledge distillation [1], which trains a
smaller network known as the student, by using the output of a bigger network
(the teacher) as the training signal. Second are methods that focus on
proposing efficient neural architectures that by design have fewer parameters
[2, 3]. For example, strategies such as depth-wise separable convolutions [4]
or careful tuning of the width and depth of the networks [3] aim at producing
efficient networks with lower complexity. Another approach to designing low-
complexity architectures is to use decomposition layers; given the fact that
neural network parameters are mostly represented with high-dimensional weight
tensors, several approaches have been developed to decompose these into
smaller tensors to reduce the computational load, without sacrificing the
model’s performance. In [5, 6, 7], the authors propose different methods to
decompose a single convolutional layer into multiple smaller ones, resulting
in more compact models with computational speedup. Similarly, [8, 9, 10]
investigate the use of $1\times 1$ convolutions to reduce the number of
channels before applying more expensive operations with larger filter sizes. A
third class of approaches aims at starting with a large high-performance
model, to then remove a large part of the parameters by network pruning [11,
12, 13, 14]. Parameter Pruning is the process of compressing a neural network
by zeroing out some of its parameters. Some methods focus on removing the
weights of pre-trained models [15, 16, 12], while more recent approaches
incorporate iterative training, by pruning and resetting only the non-pruned
weights [13], or pruning without the use of any training data [17]. We will
focus our analysis on the width and depth restriction methods (Section 3.1),
on decomposed CNNs (Section 3.2), and on parameter pruning (Section 3.3),
which have gained more popularity and interest among the scientific community,
and leave aside the Knowledge Distillation approaches due to their high
complexity and slow nature.

Recently CNNs have been successfully used for end-to-end Acoustic Scene
Classification (ASC) [18, 19, 20, 21], outperforming previous approaches and
setting new state-of-the-art. Recent studies on ASC with CNNs have revealed
that regularizing the RF of CNNs can significantly improve their
generalisation [22, 20, 23, 21, 24]. Further, authors in [18] provide a
systematic way for controlling the receptive-field of CNNs by adapting the
architectural design of the networks. Although such regularizations improve
generalisation, they affect the architectural design and the number of
parameters used in a model. Hence, finding methods that can achieve high
performance by following the RF-regularization principles, while having
minimal complexity becomes a challenging and important task.

In this paper, we aim at connecting the concept of RF-regularization, with the
low-complexity CNNs, and investigate the relationship between the RF size,
architecture complexity, and generalisation performance of CNNs for the task
of ASC. To this end, we analyse the performance of different low-complexity
ASC methods, under various RFs, and show how RF affects the performance in
low-complexity settings. We empirically evaluate different approaches to low-
complexity ASC, and analyse their performance under various maximum receptive
fields to reveal the connection between generalisation in a low-complexity
setup, and the RF of the models. Further, we propose a novel RF-regularization
technique called “Damping” which regularizes the RF of any CNN, without a need
to alter the topological design of the network. We show that our Damping RF
regularization achieves the best performance with both pruning and
decomposition architectures, and hence is a suitable approach for improving
generalisation of models in low-complexity settings.

##  2 Architectures

Previous work has shown the success of RF-regularized CNNs in various acoustic
tasks [22, 23]. Therefore, we base our work on the RF-regularized ResNet
architecture introduced in [18]. Furthermore, we introduce a new technique for
further restricting the effective receptive field of the network and provide
empirical evidence on its success in ASC.

###  2.1 Baseline ResNet Architecture: _CP-ResNet_

The details of the RF-regularized ResNet architecture are explained in [20],
where the authors introduce the hyperparameter $\rho$ in the architecture
design, such that the RF of the architecture can be changed by varying $\rho$.
Since previous work [18, 22, 24] showed the optimal range of the RF for
different ASC datasets to be approximately between $75$ and $150$ (for the
input spectograms explained in Section 4), we restrict our experiments to
$\rho$ values in range $3$-$12$. Furthermore, we remove the tailing $1\times
1$ convolutional layers by removing the last 5 residual blocks,222The residual
blocks from 8 to 12 as explained in Table 1 of [20] in order to reduce the
number of parameters in the baseline. As shown in architecture RN1 in [18],
these tailing layers have a minor effect on the performance of the model. We
refer to this RF-regularized ResNet architecture as _CP-ResNet_ throughout
this paper.

###  2.2 Frequency Damping: _Damped CP-ResNet_

Previous work [18] has shown that restricting the RF of deep CNNs, especially
over the frequency dimension, results in better generalization on different
ASC datasets. While [18] introduces a method for systematically tuning the RF
of CNNs, the proposed approach requires changes to the architecture, which as
a result changes the number of parameters. To address these drawbacks, we
propose a novel method to restrict the Effective Receptive Field (ERF) of CNNs
– the part of the RF that has the most influence on the output activation, as
detailed in [25, 18] – by damping the convolutional filter weights over the
frequency dimensions. Each convolutional neuron has a limited receptive field
of its layer input. Damping works in such a way that the further the input is
from the center of a neuron’s RF, the less influence it will have on that
neuron.

In practice, we damp the filters of a convolutional layer by applying an
element-wise multiplication between the convolution filter weights and a non-
trainable constant matrix $C\in\mathbb{R}^{T\times F}$ (damping matrix). The
damping matrix matches the spatial shape of the filters. It decays linearly
away from the center, so that the outermost elements of the filter over the
frequency dimension have a smaller influence on the activation.

The resulting network is called _damped CNN_ , where every convolution
operation $O_{n}=W_{n}*Z_{n-1}+B_{n}$ is replaced by $O_{n}=(W_{n}\odot
C_{n})*Z_{n-1}+B_{n}$, $*$ is the convolution operator, $\odot$ is the
element-wise multiplication operator, $Z_{n-1}$ is the output of the previous
layer, $W_{n}$ is the filter trainable weight, and $B_{n}$ is the bias. The
matrix has a value of 1 in the center, and decays linearly to reach a value
$\lambda$; we used $\lambda=0.1$ in our experiments. This approach has shown
empirical improvement over _CP-ResNet_ in different ASC datasets. We refer to
this architecture as _Damp_ throughout this paper.

##  3 Model Complexity reduction approaches

In this section, we investigate 3 different approaches to reduce the number of
parameters of CNNs. We follow the principle of reducing the number of
parameters, while keeping the final receptive field of the network constant,
which allows the comparison with the baseline models in each receptive field
setting.

###  3.1 Width and Depth restriction

Reducing the width (number of channels) and depth (number of layers) of CNNs
is a simple technique to reduce the network size. Tan and Le [3] showed that
after an optimal width of the network is reached, increasing the width further
will result in only a minimal performance gain at best. The width of the
network has quadratic influence on the number of weights, while the depth has
a lower effect as the number of parameters grows linearly with the number of
layers. Based on this fact, and the goal of our experiments that is comparing
different RF setups, we focus our efforts on different network widths.

As explained in Section 2.1, we remove the tailing residual blocks with
$1\times 1$ convolutions since we are only studying networks with $\rho$
values in range $3-12$. Therefore, we decrease the depth and the number of
parameters from the baseline network in [18, 20], with only a minimal
performance impact. However, removing these tailing layers with $1\times 1$
convolutions result in a significant reduction of the number of parameters in
the baseline (reported in [24] with $\rho=7$) from $3956K$ to $1715K$.

We change the width of our CNNs (Section 2) by changing the number of channels
in initial layers from 128 (in the baseline architecture) to $64$ and $32$.
This results in reducing the baseline parameter count from $1715K$ to $431K$
and $109K$ respectively.

###  3.2 Decomposed Convolutions

Inspired by the use of singular-value-decomposition (SVD) for convolutional
neural networks [7], we propose to directly train decomposed convolutional
layers instead of decomposing and then fine-tuning a pretrained model, as done
in [7] . Given a regular convolutional layer with dimensionality

$C_{in}\times C_{out}\times k\times k,$ (1)

with $C_{in}$ and $C_{out}$ being the number of input and output filters
respectively, and $k$ being the kernel size. Such a layer can be decomposed
into three convolutional layers using a compression factor $Z$:

$\begin{split}C_{in}\times(C_{out}/Z)\times 1\times 1\\\
(C_{out}/Z)\times(C_{out}/Z)\times k\times k\\\ (C_{out}/Z)\times
C_{out}\times 1\times 1\end{split}$ (2)

For example, a $128\times 128\times 3\times 3$ convolution has $147456$
parameters (neglecting the bias). Using a compression factor $Z=4$, we get
three convolutions $128\times 32\times 1\times 1$, $32\times 32\times 3\times
3$ and $32\times 128\times 1\times 1$, resulting in 17408 parameters. Similar
approaches to decompose convolutions for the purpose of parameter reduction
are explored in [8, 9, 10]. Our proposed decomposition shares the same
structure as the ”bottleneck” building block described in [10], but is
different from them as we do not use non-linearities and batch-normalization
within the decomposed block, which proved to yield better performance.

###  3.3 Parameter Pruning

Pruning is well studied in neural network literature and several methods and
adaptation have been proposed [13, 14, 15, 16, 12, 17]. We use magnitude
pruning [12, 14] with iteratively increasing the pruning ratio until reaching
the desired number of parameters.

Authors in [17] show that pruning approaches perform very similarly, if the
compression-ratio is up to 1%. Since in our low-complexity setting, we are not
targeting compression-rates lower than 1%, we choose to use the more simple
magnitude-pruning approach. We ramp up the number of pruned parameters with an
exponential decay from 0 to the final desired number in 100 epochs, to get
models with 250K, 300k, 400k, and 500k parameters. This allows us to remove
more weights at the beginning of the training, and fewer weights in later
stages as the model converges.

Figure 1: Testing accuracy of CP-ResNet on DCASE’18 with/without filter
damping on the frequency dimension using different CNN width setups as
explained in Section 3.1.

##  4 Experimental Setup

We evaluate our proposed approaches introduced in Section 3 on two acoustic
scene classification (ASC) datasets.

DCASE’20 Low Complexity ASC [26]: The dataset contains 40 hours of recordings
from 12 European cities in 10 different acoustic scenes that are summarized
into three categories, _indoor_ , _outdoor_ and _transportation_. We follow
the training/test split provided by the task organizers where recordings from
one city are not seen during training.

DCASE’18 ASC [27]: The dataset comprises recordings from six European cities
with similar acoustic scenes as the DCASE’20 dataset. We choose this dataset
because it has 10 classes which makes it a more challenging task compared to
DCASE’20, therefore we can better differentiate the performance between
different approaches with the 500KB limit (and more) that was used in the
DCASE’20 challenge. This dataset has a relatively smaller size resulting in
shorter training time. Therefore, we use DCASE’18 for our detailed analysis.

For feature extraction, we compute perceptually weighted Mel-spectrograms from
the provided 10 second audio snippets which we down-sample to 22.05kHz similar
to [18]. Each input channel of the stereo audio is processed independently,
normalized using the training set statistics and provided as a two-channel-
spectrogram input to the CNN. We use the same experiments and training setup
as explained in [20].

##  5 Results

Figure 2: Comparison of different complexity reduction methods in different
receptive field setups on DCASE’18. All methods apply frequency-damping. The
line color represents the accuracy. Figure 3: Comparing different complexity
reduction methods in different receptive field setups with and without damping
on DCASE’18. All methods use 128 base channels. Decomp uses a compression
factor $Z=4$. ND refers to no damping method being applied to the original
Network. The size of the dots relates to the number of model parameters.
Figure 4: Comparing different complexity reduction approaches on the DCASE’18
dataset. All networks are damped. The size of the dots relates to the number
of model parameters.

Table 1 summarizes the results of our best performing models on the DCASE’20
and DCASE’18 datasets. On both datasets we observe an improvement for the
damped ResNet architecture, compared to the baseline. Furthermore, we observe
that a simple width reduction, i. e., reducing the number of base channels and
thus decreasing the overall number of parameters, is not as effective as
pruning the weights of a bigger network. Following the pruning strategy
explained in Section 3.3, we can achieve models with a significantly lower
number of non-zero parameters, with the same performance on DCASE’20 or only
minor performance degradation on DCASE’18.

| DCASE’20 | DCASE’18
---|---|---
Model | Non-Zero | Acc. | Non-Zero | Acc.
CP-ResNet | 1061.6K | .9726 | 1715.6K | .7852
Damp. | 1061.6K | .9745 | 1715.6K | .7915
Damp.-Restricted | 268.6K | .9722 | 431.8K | .7701
Damp.-Pruned | 250K | .9745 | 400K | .7854
Damp.-Decomp. | 361K | .9700 | 417K | .7696
Table 1: Comparison of the performance of the best models of the discussed
approaches on the DCASE‘20 and DCASE’18 datasets. We report the accuracy
(Acc.) and the number of non-zero parameters. _Damp._ refers to a damped
ResNet as introduced in Section 2.2, _Damp.-Restriced_ is a width restricted
ResNet with 64 base channels (cf. Section 3.1), _Damp.-Pruned_ is a pruned
ResNet (cf. Section 3.3) and _Damp.-Decomp._ is a ResNet where all
convolutions are replaced with decomposed convolution and a compression factor
$Z=4$ (cf. Section 3.2). The table shows the mean of the accuracy over the
last 10 epochs for each model. All DCASE’20 and DCASE’18 models use $\rho=4$
and $\rho=7$, respectively.

###  5.1 Frequency Damping Comparison

To further investigate the effect of _damping_ , we compare the performance of
several baseline CP-ResNets with and without _damping_ , using different
maximum RF setups on DCASE’18. Fig. 1 shows that _damping_ not only improves
the overall performance, but is also more robust to the choice of the maximum
RF. We also observe that the damped ResNet with 128 base channels (Damped)
yields the best overall performance using $\rho=7$ to regularize the RF. The
model with 64 base channels (Damped_64) decreases the performance, while
significantly reducing the number of parameters (shown in Figure 2), and 32
base channels seems to result in an insufficient amount of parameters for this
task. As can be seen in Fig. 3, Damped CP-ResNets outperforms non-damped CP-
ResNets with the same number of parameters, both in pruning and decomposition
cases.

###  5.2 Parameter Pruning vs. Network Decomposition

Comparing all parameter reduction methods in Figs. 2 and 4, we see that
pruning is the closest model to the frequency-damped baseline (Damp) in terms
of performance, and achieves the highest accuracy with limited parameters.
Fig. 4 further shows that the pruned model achieves an accuracy close to the
baseline, while significantly reducing the number of non-zero parameters.
Using decomposed convolutions, the performance degradation is significantly
higher than for pruning. Similarly, simply reducing the width of the network
results in a large loss in accuracy.

##  6 Conclusion

In this paper, we analysed various low-complexity CNN-based approaches for
ASC, and studied the relationship between the RF and the performance in each
approach. We showed that pruning achieves better performance compared to
decomposition and width reduction methods. We proposed a filter damping
technique that can be used in low-complexity settings to regularize the RF of
models, without altering the architectures. We showed that using filter
damping improves the performance on all architectures and datasets we
evaluated, hence is a simple and effective technique for improving
generalisation of models. Our results on two datasets for ASC suggests that
the proposed filter damping can achieve state-of-the-art performance in low-
complexity ASC. Using this approach, we achieved the 1st rank at the
DCASE-2020 Challenge Low-Complexity Acoustic Scene Classification (Task1.b).

##  7 ACKNOWLEDGMENT

This work has been supported by (1) the COMET-K2 Center of the Linz Center of
Mechatronics (LCM) funded by the Austrian federal government and the federal
state of Upper Austria, and (2) the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program (grant agreement
number 670035, project ”Con Espressione”).

## References

  * [1] G. Hinton, O. Vinyals, and J. Dean, “Distilling the knowledge in a neural network,” _arXiv preprint arXiv:1503.02531_ , 2015.
  * [2] A. G. Howard, M. Zhu, B. Chen, D. Kalenichenko, W. Wang, T. Weyand, M. Andreetto, and H. Adam, “Mobilenets: Efficient convolutional neural networks for mobile vision applications,” _arXiv preprint arXiv:1704.04861_ , 2017.
  * [3] M. Tan and Q. V. Le, “Efficientnet: Rethinking model scaling for convolutional neural networks,” in _Proc. of the 36th Int. Conf. on Machine Learning (ICML)_ , 2019, pp. 6105–6114.
  * [4] A. G. Howard, M. Zhu, B. Chen, D. Kalenichenko, W. Wang, T. Weyand, M. Andreetto, and H. Adam, “Mobilenets: Efficient convolutional neural networks for mobile vision applications,” _CoRR_ , vol. abs/1704.04861, 2017\.
  * [5] V. Lebedev, Y. Ganin, M. Rakhuba, I. Oseledets, and V. Lempitsky, “Speeding-up Convolutional Neural Networks using fine-tuned CP-Decomposition,” in _Proc. of the Int. Conf. on Learning Representations (ICLR)_ , 2015.
  * [6] M. Jaderberg, A. Vedaldi, and A. Zisserman, “Speeding up convolutional neural networks with low rank expansions,” in _British Machine Vision Conf. (BMVC)_ , 2014.
  * [7] Y.-D. Kim, E. Park, S. Yoo, T. Choi, L. Yang, and D. Shin, “Compression of deep convolutional neural networks for fast and low power mobile applications,” _arXiv preprint arXiv:1511.06530_ , 2016.
  * [8] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich, “Going Deeper with Convolutions,” in _Proc. of the IEEE Conf. on Computer Vision and Pattern Recognition (CVPR)_ , 2015, pp. 1–9.
  * [9] C. Szegedy, V. Vanhoucke, S. Ioffe, J. Shlens, and Z. Wojna, “Rethinking the Inception Architecture for Computer Vision,” in _Proc. of the IEEE Conf. on Computer Vision and Pattern Recognition (CVPR)_ , 2016, pp. 2818–2826.
  * [10] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proc. of the IEEE Conf. on Computer Vision and Pattern Recognition (CVPR)_ , 2016, pp. 770–778.
  * [11] Y. LeCun, J. S. Denker, and S. A. Solla, “Optimal brain damage,” in _Advances in Neural Information Processing Systems 2_ , 1990, pp. 598–605.
  * [12] S. Janowsky, “Pruning versus clipping in neural networks.” _Physical review. A, Atomic, molecular, and optical physics_ , vol. 39, no. 12, pp. 6600–6603, 1989.
  * [13] J. Frankle and M. Carbin, “The lottery ticket hypothesis: Finding sparse, trainable neural networks,” in _7th Int. Conf. on Learning Representations (ICLR)_ , 2019.
  * [14] S. Han, J. Pool, J. Tran, and W. Dally, “Learning both weights and connections for efficient neural network,” in _Advances in Neural Information Processing Systems 28_ , C. Cortes, N. D. Lawrence, D. D. Lee, M. Sugiyama, and R. Garnett, Eds. Curran Associates, Inc., 2015, pp. 1135–1143.
  * [15] Y. LeCun, J. S. Denker, and S. A. Solla, “Optimal brain damage,” in _Advances in Neural Information Processing Systems 2_ , 1989, pp. 598–605.
  * [16] M. Mozer and P. Smolensky, “Skeletonization: A technique for trimming the fat from a network via relevance assessment,” in _Advances in Neural Information Processing Systems 1_ , 1988, pp. 107–115.
  * [17] H. Tanaka, D. Kunin, D. L. K. Yamins, and S. Ganguli, “Pruning neural networks without any data by iteratively conserving synaptic flow,” _CoRR_ , vol. abs/2006.05467, 2020.
  * [18] K. Koutini, H. Eghbal-zadeh, M. Dorfer, and G. Widmer, “The Receptive Field as a Regularizer in Deep Convolutional Neural Networks for Acoustic Scene Classification,” in _Proc. of the European Signal Processing Conf. (EUSIPCO)_ , 2019.
  * [19] H. Chen, Z. Liu, Z. Liu, P. Zhang, and Y. Yan, “Integrating the data augmentation scheme with various classifiers for acoustic scene modeling,” DCASE2019 Challenge, Tech. Rep., 2019.
  * [20] K. Koutini, H. Eghbal-zadeh, and G. Widmer, “Receptive-field-regularized CNN variants for acoustic scene classification,” in _Proc. of the Detection and Classification of Acoustic Scenes and Events 2019 Workshop (DCASE2019)_ , 2019.
  * [21] S. Suh, S. Park, Y. Jeong, and T. Lee, “Designing Acoustic Scene Classification Models with CNN Variants,” DCASE2020 Challenge, Tech. Rep., 2020.
  * [22] K. Koutini, H. Eghbal-zadeh, and G. Widmer, “CP-JKU submissions to DCASE’19: Acoustic scene classification and audio tagging with receptive-field-regularized CNNs,” DCASE2019 Challenge, Tech. Rep., 2019.
  * [23] K. Koutini, S. Chowdhury, V. Haunschmid, H. Eghbal-Zadeh, and G. Widmer, “Emotion and theme recognition in music with Frequency-Aware RF-Regularized CNNs,” in _MediaEval Benchmark Workshop_ , 2019.
  * [24] K. Koutini, F. Henkel, H. Eghbal-zadeh, and G. Widmer, “CP-JKU Submissions to Dcase’20: Low-Complexity Cross-Device Acoustic Scene Classification with RF-Regularized CNNs,” DCASE2020 Challenge, Tech. Rep., 2020.
  * [25] W. Luo, Y. Li, R. Urtasun, and R. Zemel, “Understanding the Effective Receptive Field in Deep Convolutional Neural Networks,” in _Advances in Neural Information Processing Systems 29_ , 2016, pp. 4898–4906.
  * [26] T. Heittola, A. Mesaros, and T. Virtanen, “Acoustic scene classification in dcase 2020 challenge: generalization across devices and low complexity solutions,” in _Proc. of the Detection and Classification of Acoustic Scenes and Events 2020 Workshop (DCASE2020)_ , 2020.
  * [27] A. Mesaros, T. Heittola, and T. Virtanen, “A multi-device dataset for urban acoustic scene classification,” in _Proc. of the Detection and Classification of Acoustic Scenes and Events 2018 Workshop (DCASE2018)_ , 2018, pp. 9–13.


# Two-stage melting of an inter-component Potts long-range order in two
dimensions

Feng-Feng Song1 and Guang-Ming Zhang1,2,3<EMAIL_ADDRESS>1State Key
Laboratory of Low-Dimensional Quantum Physics and Department of Physics,
Tsinghua University, Beijing 100084, China.
2Collaborative Innovation Center of Quantum Matter, Beijing, China.
3Frontier Science Center for Quantum Information, Beijing 100084, China.

###### Abstract

Interplay of topology and competing interactions can induce new phases and
phase transitions at finite temperatures. We consider a weakly coupled two-
dimensional hexatic-nematic XY model with a relative $Z_{3}$ Potts degrees of
freedom, and apply the matrix product state method to solve this model
rigorously. Since the partition function is expressed as a product of two-
legged one-dimensional transfer matrix operator, an entanglement entropy of
the eigenstate corresponding to the maximal eigenvalue of this transfer
operator can be used as a stringent criterion to determine various phase
transitions precisely. At low temperatures, the inter-component $Z_{3}$ Potts
long-range order (LRO) exists, indicating that the hexatic and nematic fields
are locked together and their respective vortices exhibit quasi-LRO. In the
hexatic regime, below the BKT transition of the hexatic vortices, the inter-
component $Z_{3}$ Potts LRO appears, accompanying with the binding of nematic
vortices. In the nematic regime, however, the inter-component $Z_{3}$ Potts
LRO undergoes a two-stage melting process. An intermediate Potts liquid phase
emerges between the Potts ordered and disordered phases, characterized by an
algebraic correlation with formation of charge-neutral pairs of both hexatic
and nematic vortices. These two-stage phase transitions are associated with
the proliferation of the domain walls and vortices of the relative $Z_{3}$
Potts variable, respectively. Our results thus provide a prototype example of
two-stage melting of a two-dimensional long-range order, driven by multiple
topological defects.

Introduction. -The concept of topological defects has been one of the
cornerstones in the study of phase transitions in two-dimensional (2D)
systems. The known example is the Berezinskii-Kosterlitz-Thouless (BKT) phase
transitionBerezinsky (1971); Kosterlitz and Thouless (1973); Kosterlitz (1974)
where unbinding of integer vortices occurs in the absence of spontaneous
breaking of continuous symmetry. The followed renormalization group
Kosterlitz-Thouless-Halperin-Nelson-Young theoryHalperin and Nelson (1978);
Nelson and Halperin (1979); Young (1979) describes a two-stage melting process
of 2D crystals: dissociation of dislocation pairs induces a BKT-type
transition from a crystal phase to a hexatic phase with a quasi-long-range
(quasi-LRO) orientational order and a further transition is resulted from
unbinding of disclinations, as a dislocation can be represented as a coupled
pair of disclinations. Despite the four decade research, it is still
challenging to find a simple microscopic model that exhibit these phenomena
theoreticallyKosterlitz (2016); Ryzhov et al. (2017).

Recently it has known that the interaction of multiple topological excitations
also plays a central role in many physical systems ranging from
superfluids/superconductors to condensed atom-molecular mixturesStrandburg
(1988); Bernard and Krauth (2011); Chae et al. (2012); Radzihovsky et al.
(2008); Shahbazi and Ghanbari (2006); de Parny et al. (2016); Serna et al.
(2017); Kobayashi et al. (2019). Among these studies, a prototype model is a
2D coupled hexatic-nematic XY spin model proposed for an unusual melting
processBruinsma and Aeppli (1982); Aeppli and Bruinsma (1984) and a hidden-
order phase transition of isotropic liquid-crystal thin filmsJiang et al.
(1993, 1996); Drouin-Touchette et al. (2022). As shown in Fig.1(a), the
hexatic degrees of freedom describe a six-fold bond-orientational field
represented by $\Theta=|\Theta|\exp(i6\vartheta)$ with $\vartheta$ as the bond
angle linking the centers of mass of neighboring moleculesHalperin and Nelson
(1978), while the nematic degrees of freedom arising from the herringbone
order in the crystalline phaseBruinsma and Aeppli (1982) are described by
$\Phi=|\Phi|\exp(i2\varphi)$. Three inequivalent herringbone patterns are
displayed in Fig.1(b), from which a different herringbone pattern with the
same orientation can be obtained by translating over a lattice vector. The
Hamiltonian is described by a coupled bilayer system

$\displaystyle H$ $\displaystyle=$ $\displaystyle-J_{2}\sum_{\langle
i,j\rangle}\cos(2\varphi_{i}-2\varphi_{j})-J_{6}\sum_{\langle
i,j\rangle}\cos(6\vartheta_{i}-6\vartheta_{j})$ (1)
$\displaystyle-K\sum_{i}\cos(6\vartheta_{i}-6\varphi_{i}),$

where $\varphi_{i}$ and $\vartheta_{i}\in[0,2\pi]$ are two $U(1)$ phase
fields, $J_{2}$ and $J_{6}$ are their respective nearest-neighbour intra-field
couplings, and the inter-component coupling $K$ denotes the minimal hexatic-
nematic coupling allowed by the relative symmetryBruinsma and Aeppli (1982).

Figure 1: (a) The bond-orientational angle $\vartheta$ joining the centers of
neighboring molecules with respect to a laboratory axis. (b) Three possible
herringbone patterns. (c) The double-layer tensor network of the partition
function and the local tensor. (d) Eigen-equation of the 1D transfer operator.

The study of this 2D coupled XY spin model in the context of nematic/hexatic
liquid crystalline systems has a very long history, and is plagued with
complexities, confusions and unexplained puzzling results, which have been
reviewed in Ref.Drouin-Touchette et al. (2022). In the strongly coupled model,
large-scale Monte Carlo simulations had determined the possible $Z_{3}$ Potts
transition relative to the BKT transitionJiang et al. (1993, 1996). In the
hexatic regime, the inter-component $Z_{3}$ Potts LRO is developed below the
BKT transition of the hexatic vortices, however, in the nematic regime the
formation of the $Z_{3}$ Potts LRO coincides with the BKT vortex binding
transitionJiang et al. (1993, 1996). A recent numerical work suggests that the
$Z_{3}$ Potts LRO in the nematic regime may form above the BKT
transitionDrouin-Touchette et al. (2022). Due to the lack of sharp
thermodynamic signatures associated with the binding of topological defects,
the scenario of a single transition cannot be fully excluded. The phase
structure in the nematic regime is expected to be revealed only in the weakly
coupled model, which reminds largely unexplored because of the requirements of
significantly larger system sizes with enlarged vortex cores.

In order to resolve this long-standing issue, we apply a state-of-art tensor
network methodVerstraete et al. (2008); Orús (2014); Haegeman and Verstraete
(2017) to surmount those difficulties. As the partition function of a 2D
statistical model can be represented as a product of two-legged 1D transfer
operator (Fig.1(c)), its eigen-equation (Fig.1(d)) can be solved by the
algorithm of variational uniform matrix product states (MPS) in the
thermodynamic limitZauner-Stauber et al. (2018); Fishman et al. (2018);
Vanderstraeten et al. (2019a, b). According to the singularity displayed by
the entanglement entropy for the 1D quantum analogue, the various phase
transitions in the global phase diagram can be precisely determinedLi et al.
(2020); Song and Zhang (2022a). In the hexatic regime, we obtain the similar
results to the previous onesJiang et al. (1993, 1996); Drouin-Touchette et al.
(2022). To our surprise, in the nematic regime, we discover that the inter-
component Potts LRO itself undergoes a two-stage melting process. An
intermediate Potts liquid phase emerges simultaneously with the formation of
charge-neutral pair of hexatic and nematic vortices. Such two-stage
transitions are associated with the separated proliferations of the domain
walls and vortices of the Potts variable, respectively. Furthermore, we show
that, when increasing the hexatic-nematic coupling strength, the two-stage
transitions gradually merge into a single transition, the result for the
strongly coupled modelJiang et al. (1993, 1996).

Tensor Network Method. -We first simplify the model Hamiltonian Eq.(1) by
introducing new variables $\phi=2\varphi$ and $\theta=6\vartheta$,

$\displaystyle H/J$ $\displaystyle=$ $\displaystyle-\Delta\sum_{\langle
i,j\rangle}\cos(\phi_{i}-\phi_{j})-(2-\Delta)\sum_{\langle
i,j\rangle}\cos(\theta_{i}-\theta_{j})$ (2)
$\displaystyle-\lambda\sum_{i}\cos(\theta_{i}-3\phi_{i}),$

where $\lambda=K/J$, $J=(J_{2}+J_{6})/2$, and $\Delta=J_{2}/J$ as the relative
interacting strength of the hexatic and nematic fields. In the absence of the
hexatic-nematic coupling, the model has $U(1)\times U(1)$ symmetry and
exhibits two independent BKT phase transitions. For $\lambda\neq 0$, the model
is invariant under a $U(1)\times Z_{3}$ transformation,
$\phi_{i}\to\phi_{i}+\alpha/3+{2\pi}k_{i}/3$ and
$\theta_{i}\to\theta_{i}+\alpha$, where $k_{i}=0,1,2$ correspond to a $Z_{3}$
degrees of freedom. Specially, a sufficient large $\lambda$ tends to lock the
hexatic and nematic fields, $\theta_{i}=3\phi_{i}$, and the model is reduced
to a generalized XY modelGranato and Kosterlitz (1986); Poderoso et al.
(2011); Canova et al. (2014, 2016).

In the tensor network framework, after a duality transformation, the partition
function is expressed as a tensor contraction over all auxiliary links,

$Z=\mathrm{tTr}\prod_{i}O_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}(i),$
(3)

where “tTr” denotes the tensor contraction, which forms a double-layer tensor
network shown in Fig. 1(c). Each local tensor $O$ is given by

$\displaystyle O_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}$
$\displaystyle=$
$\displaystyle\sum_{k}\left(\prod_{l=1}^{4}I_{n_{l}}(\beta\Delta)I_{m_{l}}(\beta(2-\Delta))\right)^{1/2}$
(4) $\displaystyle\times
I_{k}(\beta\lambda)\delta_{n_{1}+n_{2}}^{n_{3}+n_{4}+3k}\delta_{m_{1}+m_{2}+k}^{m_{3}+m_{4}},$

where $I_{n}(x)$ is the modified Bessel function of the first kind, $n$ and
$m$ are integers, and $\beta=1/(k_{B}T)$. The global $U(1)\times Z_{3}$
invariance is encoded in the local tensor. In the thermodynamic limit, the
partition function is determined by the dominant eigenvalues of the 1D quantum
transfer operator $\hat{T}$, whose eigen-equation (Fig. 1(d))
$\hat{T}|\Psi(A)\rangle=\Lambda_{\max}|\Psi(A)\rangle$ can be accurately
solved by the algorithm of variational uniform MPSZauner-Stauber et al.
(2018); Fishman et al. (2018); Vanderstraeten et al. (2019a, b). The
corresponding 1D quantum Hamiltonian $\hat{H}_{1D}=-(1/\beta)\ln\hat{T}$, and
the leading eigenvector $|\Psi(A)\rangle$ is represented by a MPS, whose
precision is controlled by the auxiliary bond dimension $D$ of the local
tensors.

From the maximal eigenvalue and its eigenvector of the 1D quantum transfer
operator, various physical quantities can be estimated accuratelyLi et al.
(2020); Song and Zhang (2021, 2022b, 2022a). As the phase transitions are
concerned, the quantum entanglement entropy is the most efficient measureVidal
et al. (2003); Pollmann et al. (2009), which can be directly determined via
the Schmidt decomposition: $S_{E}=-\sum_{\alpha=1}^{D}s_{\alpha}^{2}\ln
s_{\alpha}^{2}$ with $s_{\alpha}$ as the singular values from the bipartition
of the state $|\Psi(A)\rangle$. Various two-point correlation functions can be
evaluated by the trace of an infinite sequence of channel operators containing
two local impurity tensorsSM (See Supplementary Materials for the tensor
network method and detailed numerical calculations.).

The superfluid response is described by the spin stiffness as the second
derivative of the free energy density $f=-(1/N\beta)\ln Z$ with respect to a
twist $v$ along a reference direction:
$\rho_{s}=\frac{\partial^{2}f}{\partial^{2}v}|_{v=0}$. The twist needs to be
imposed in a way that respects the joint $U(1)$ invariance of the hexatic and
nematic fields as

$(\phi_{i},\theta_{i})\rightarrow(\phi_{i}+\vec{v}\cdot\vec{r}_{i},\theta_{i}+3\vec{v}\cdot\vec{r}_{i}),$
(5)

where $\vec{r}_{i}$ is the position vector for the lattice site $i$ and
$\vec{v}$ is a constant vector to increase the phase difference across each
neighboring hexatic spins $3$ times larger than nematic spins. So the jump of
spin stiffness is altered from the BKT predictions when the unbinding of the
$\theta$ and $\phi$ vortices happens separatelyHübscher and Wessel (2013).

The hexatic-nematic XY model has a rich physics, and the most intriguing
results are expected in the weakly coupled case. For a typical value of
$\lambda=0.1$, a global phase diagram is derived in Fig. 2. All phase
boundaries are determined by the singularity displayed in the entanglement
entropy $S_{E}$ for the 1D quantum transfer operator. Since the succession of
phases crucially depends on the intra-component coupling ratio, our results
are discussed in hexatic $\Delta<1$ and nematic $\Delta>1$ regimes,
separately.

Figure 2: The global phase diagram of the weakly coupled hexatic-nematic XY
model with a typical coupling $\lambda=0.1$. (a) Schematic pictures of
different topological excitations. (b) In the Potts liquid phase, the inter-
component Potts variable has quasi-LRO, the vortex-antivortex pairs in $\phi$
fields are formed, and the dominant topological defects of the $\theta$ fields
are composite vortex pairs with charge $q_{\theta}=\pm 3$. (c) In the
fractional vortex paired phase, the $\theta$ vortices have quasi-LRO, and the
$\phi$ vortices are fractionalized as paired vortices with charge
$q_{\phi}=1/3$. (d) In the Potts ordered phase, vortices in both types are
bound in pairs, accompanying the inter-component Potts LRO.

Phase transitions in the hexatic regime. -For a typical value $\Delta=0.8$,
our numerical results show two singular peaks in the entanglement entropy at
$T_{c1}\simeq 0.805J$ and $T_{c2}\simeq 1.145J$, respectively, as shown in
Fig. 3(a). The peak positions are nearly unchanged with increasing bond
dimensions $D=100,200,300$ so the transition points are determined with high
accuracy. The corresponding specific heat is shown in Fig. 3(b), exhibiting a
sharp divergence at $T_{c1}$ and a small rounded bump around $T_{c2}$. The
maximum of the hump is above $T_{c2}$, a typical feature of the BKT
transitionKosterlitz and Thouless (1973); Kosterlitz (1974), while the low-$T$
singular behavior is fitted by $C_{V}\propto|T-T_{c1}|^{\alpha}$ with
$\alpha\simeq 1/3$. Further evidence is provided by a $Z_{3}$ order parameter
$M_{\sigma}=\langle\cos(\sigma_{i})\rangle$ with
$\sigma_{i}=\theta_{i}/3-\phi_{i}$. As shown in Fig. 3(c), $M_{\sigma}$
becomes finite at $T_{c1}$, suggesting that a true LRO is established and the
relative phase between the $\phi$ and $\theta$ fields is fully locked. The
magnetization satisfies the scaling form $M_{\sigma}\propto(T_{c1}-T)^{\beta}$
with the critical exponent $\beta\simeq 1/9$. These critical exponents are
perfectly in agreement with the $Z_{3}$ Potts transitionWu (1982).

Figure 3:  The numerical results for $\Delta=0.8$ and $\lambda=0.1$. (a) The
entanglement entropy. (b) The specific heat. (c) The magnetization of the
Potts variable. In the inset $t$ is the reduced temperature
$t=|T_{c1}-T|/T_{c1}$. (d) The superfluid stiffness. (e) and (f) At $T=1.0J$,
the correlation function $G_{\phi}(r)$ for $q_{\phi}=1$ vortices decays
exponentially, but the correlation function $G_{3\phi}(r)$ for $q_{\phi}=1/3$
vortices decays in power law.

As shown in Fig. 3(d), the spin stiffness starts to dramatically increase from
zero at the BKT transition $T_{c2}$. When the temperature further decreases, a
small jump appears to enhance the spin stiffness at the Potts transition
$T_{c1}$ precisely. Such a small increase of spin stiffness is resulted from
bindings of integer $\phi$ vorticesCanova et al. (2014); Drouin-Touchette et
al. (2022). When cooling the system from the disordered phase, the correlation
function $G_{\theta}(r)=\langle\cos(\theta_{i}-\theta_{i+r})\rangle$ exhibits
a power law decay at $T_{c2}$, indicating the binding of hexatic vortices. In
the intermediate temperature regime ($T_{c1}<T<T_{c2}$), we found that
$G_{\phi}(r)=\langle\cos(\phi_{i}-\phi_{i+r})\rangle$ decays exponentially,
but $G_{3\phi}(r)$ exhibits an algebraic correlation. In Fig. 3(e) and (f), a
direct comparison at $T=1.0J$ suggests that the vortices are fractionalized
into the $q_{\phi}=1/3$ vortices. As the correlation length $\xi_{\phi}$
extracted from $G_{\phi}(r)$ follows an exponentially divergence above
$T_{c1}$, we thus conclude that the low-temperature transition is a hybrid BKT
and Potts transitionSM (See Supplementary Materials for the tensor network
method and detailed numerical calculations.).

Phase transitions in the nematic regime. -For a typical value $\Delta=1.2$,
our numerical calculations show two peaks in the entanglement entropy at
$T_{c1}\simeq 0.915J$ and $T_{c2}\simeq 1.15J$, respectively, as seen in Fig.
4(a). Unlike the strongly coupled modelDrouin-Touchette et al. (2022), the
specific heat $C_{V}$ exhibits a pronounced peak at $T_{c1}$ and a bump at
$T_{c2}$, as displayed in Fig. 4(b). The bump appears lightly above $T_{c2}$
as the usual BKT transition, but the peak locates below $T_{c1}$.
Surprisingly, the magnetization of the relative Potts variable displays a two-
step feature below $T_{c2}$, as shown in Fig. 4(c). The artificial finite
Potts magnetization in the intermediate phase is due to the finite-bond
dimension effect, which should decrease slowly to zero with increasing bond
dimensions in MPS approximations for $U(1)$ phasesVanderstraeten et al.
(2019b); Li et al. (2020). To figure out the transition at $T_{c1}$, we
calculate the superfluid stiffness, which also has a two-step jump, as
displayed in Fig. 4(d). In contrast to the hexatic regime, the stiffness drops
more dramatically at $T_{c1}$ than at $T_{c2}$. Since the phase twist
$\vec{v}$ applied to the hexatic $\theta$ field is three times larger than
that to the nematic $\phi$ field, the first drop of the stiffness is nearly
six times larger than the second drop.

Figure 4:  The numerical results for $\Delta=1.2$ with $\lambda=0.1$. (a) The
entanglement entropy. (b) The specific heat. (c) The Potts magnetization. (d)
The superfluid stiffness.  Figure 5:  (a) and (c) The correlation functions
$G_{\theta}(r)$ in both the Potts ordered phase ($T=0.8J$) and the Potts
liquid phase ($T=1.0J$). (e) The correlation $G_{\theta}(r)$ at $r=100$ and
the exponents in the inset. (b), (d) and (f) The correlation functions
$G_{\sigma}(r)$ in the Potts ordered, Potts liquid, and disordered phases.

The abrupt drop of spin stiffness at $T_{c1}$ is rather rare, and the binding
of charge-neutral vortex pairs of $\theta$ field can be ruled out from the
correlation function
$G_{\theta}(r)=\langle\cos(\theta_{i}-\theta_{i+r})\rangle$. On two sides of
$T_{c1}$, the $G_{\theta}(r)$ correlations have the power law behavior as
shown in Fig. 5(a) and (c), indicating that the quasi-LRO persists through the
transition. As displayed in Fig. 5(e), a direct comparison between the
amplitudes of $G_{\theta}(r)\sim r^{-\eta_{\theta}}$ at $r=100$ indicates that
the correlation above $T_{c1}$ is greatly suppressed by three orders of
magnitude, and the exponent $\eta_{\theta}$ varies with temperature depicted
in the inset. The extremely weak correlation above $T_{c1}$ may account for
the dramatic jump in total stiffness at $T_{c1}$. Below $T_{c2}$, the onset of
algebraic correlations of the $\phi$ vortices also induces the quasi-LRO of
the $\theta$ vortices. So the dominant topological excitations between
$T_{c1}$ and $T_{c2}$ are composite vortex pairs of $q_{\theta}=3$ with an
internal structure of three bound vortices.

By calculating the correlations of the relative $Z_{3}$ Potts variable, we
find that the inter-component Potts order-disorder phase transition splits
into two different transitions separated by an intermediate liquid phase with
quasi-LRO correlation. As shown in Fig. 5(b), (d) and (f), the correlation
function $G_{\sigma}(r)$ exhibits a LRO in the Potts ordered phase, an
algebraic decay in the intermediate phase, and an exponential decay in the
Potts disordered phase, respectively. Actually the Potts quasi-LRO occurs
simultaneously with the formation of charge-neutral pairs of hexatic and
nematic vortices. Moreover, the correlation length for Potts variables
displays a similar structure as the $p$-state clock model with $p>4$ (Ref.SM
(See Supplementary Materials for the tensor network method and detailed
numerical calculations.)).

Figure 6: The entanglement entropy for different $\lambda$ in the hexatic
regime (a) $\Delta=0.8$ and the nematic regime (b) $\Delta=1.2$. For a good
comparison, each curve is relatively shifted by $\Delta{S_{E}}=0.2$. The
displayed results are obtained under MPS bond dimension $D=100$.

Unlike the 2D Ising model, the excitations in the $Z_{3}$ Potts model include
both the loop-like domain walls and vortices, and the usual Potts transition
is resulted from a coupled proliferation of two topological defectsEinhorn et
al. (1980). However, a quasi-LRO phase appears in a generalized $Z_{3}$ Potts
model by artificially raising the vortex core energyBhattacharya and Ray
(2016). We notice that the effective core energy of $Z_{3}$ vortices in the
coupled hexatic-nematic model is related to the hexatic-nematic coupling. A
larger $\lambda$ makes the excitations of free vortices in the $\theta$ field
more costly, and the domain wall energy increases dramatically to forbid the
pre-excitations of domain walls. However, a small $\lambda$ is required to
achieve a relatively high ratio of the core energy of the $Z_{3}$ vortices. To
illustrate this physics, we calculated the entanglement entropy as a function
of temperature for different values of $\lambda$. As shown in Fig. 6, in the
hexatic regime ($\Delta=0.8$), the entanglement entropy results indicate that
the the model always has two phase transitions similar to those of strongly
coupled model, while in the nematic regime ($\Delta=1.2$), two-stage phase
transitions are allowed only for a small $\lambda$ value. As $\lambda$
increases, the intermediate Potts liquid phase obscures because the lower
transition gets vague, and two transitions are merged together.

Conclusion and Outlook. -We have studied a weakly coupled hexatic-nematic XY
model by using the tensor-network method. It is found that, the hexatic regime
shares the similar phase structure to that of the strongly coupled modelJiang
et al. (1993, 1996); Drouin-Touchette et al. (2022). In the nematic regime,
however, an inter-component $Z_{3}$ Potts liquid phase emerges as the
intermediate phase of the two-stage melting of the Potts LRO. Actually this
emergent Potts liquid phase opens up a promising route towards better
understanding the hidden structure of the phase transitions driven by
different topological defects. Such an intriguing splitting of phase
transitions should exist in other systems with multiple topological
excitations, such as multi-stage melting process of 2D liquid crystal
filmsPindak et al. (1981); Huang et al. (1981); Chou et al. (1998). Above the
melting of positional order, the orientational order itself possesses a
further structure. Future study can be extended to the coupled $U(1)\times
Z_{n}$ systems with larger $n$ where a richer phase diagram is expectedNelson
and Halperin (1980). Moreover, in various XY systems like condensed atom-
molecular mixtures and multi-component superfluids/superconductorsDonley et
al. (2002); Chin et al. (2010); Fernandes et al. (2010, 2019), such
unconventional phenomena should also be realizable.

Acknowledgments. The research is supported by the National Key Research and
Development Program of MOST of China (2017YFA0302902).

##  I Appendix A: Tensor network representation of our model

A good wealth of fascinating physics have been found in many-body systems but
the numerical treatment still remains challenging due to the nature of strong
correlations. Tensor networks (TN) have proven to be a very efficient tool to
overcome these challenges based on the variational ansatz wave functions. Due
to the unique advantages of faithfully capturing the entanglement structure of
many-body states, they are increasingly becoming a standard tool for not only
strongly correlated quantum systems but also many-body classical problems.

The central but often nontrivial task in condensed matter physics is to
evaluate the partition function of a system. To implement the TN method in the
coupled hexatic-nematic XY model, the first step is to convert the classical
lattice model with nearest-neighbor local interactions into a TN
representation. The coupled hexatic-nematic XY model on a 2D square lattice is
defined by the Hamiltonian

$H=-J_{2}\sum_{\langle
i,j\rangle}\cos(2\varphi_{i}-2\varphi_{j})-J_{6}\sum_{\langle
i,j\rangle}\cos(6\vartheta_{i}-6\vartheta_{j})-K\sum_{i}\cos(6\vartheta_{i}-6\varphi_{i}),$
(6)

where $\langle i,j\rangle$ refers to nearest neighbors, $\varphi_{i}$ and
$\vartheta_{i}\in[0,2\pi]$ are phase angles associated to the lattice site
$i$, $J_{2}$ and $J_{6}$ are their respective intra-component coupling
strength, and $K$ denotes the inter-component coupling. To lowest order in
$\varphi$ and $\vartheta$ by introducing $\phi=2\varphi$ and
$\theta=6\vartheta$, the simplified Hamiltonian in the same universality class
can be written as

$H/J=-\sum_{\langle
i,j\rangle}\Delta\cos(\phi_{i}-\phi_{j})-(2-\Delta)\sum_{\langle
i,j\rangle}\cos(\theta_{i}-\theta_{j})-\lambda\sum_{i}\cos(\theta_{i}-3\phi_{i}),$
(7)

where $J=\frac{1}{2}(J_{2}+J_{6})$, $\lambda=K/J$ and $\Delta=J_{2}/J\in[0,2]$
tuning the relative strength of hexatic and nematic interactions.

To derive the TN representation, we express the partition function on its
original lattice and the Boltzmann weights are represented as a tensor product
of local bonds as shown in Fig. 7(a)

$Z=\iint\frac{\mathrm{d}\phi_{i}\mathrm{d}\theta_{i}}{(2\pi)^{2}}\,\prod_{\langle
i,j\rangle}W_{U}(\phi_{i},\phi_{j})W_{I}(\phi_{i},\theta_{i})W_{L}(\theta_{i},\theta_{j}),$
(8)

where

$\displaystyle W_{U}(\phi_{i},\phi_{j})$
$\displaystyle=\mathrm{e}^{\beta\Delta\cos(\phi_{i}-\phi_{j})},$ (9)
$\displaystyle W_{I}(\phi_{i},\theta_{i})$
$\displaystyle=\mathrm{e}^{\beta\lambda\cos(\theta_{i}-3\phi_{i})},$ (10)
$\displaystyle W_{L}(\theta_{i},\theta_{j})$
$\displaystyle=\mathrm{e}^{\beta(2-\Delta)\cos(\theta_{i}-\theta_{j})}$ (11)

can be viewed as infinite matrices whose indices are the continuous $\phi$ and
$\theta$ variables. The partition function is now cast into a double layer TN
representation, where the integrations of $\int d\phi/2\pi$ and $\int
d\theta_{i}/2\pi$ are denoted as blue and green dots and the matrix indices
take the same values at the joint points.

The duality transformation changes the local tensors from continuous $U(1)$
variables onto discrete basis for further numerical analysis. As shown in Fig.
7(d), we perform the character expansion for the symmetric matrix $W$

$\displaystyle W_{U}(\phi_{i},\phi_{j})$
$\displaystyle=\sum_{n}U_{\phi_{i},n}\,I_{n}(\beta\Delta)\,U_{\phi_{j},n}^{*},$
(12) $\displaystyle W_{I}(\phi_{i},\theta_{i})$
$\displaystyle=\sum_{k}U_{\phi_{i},-3k}\,I_{k}(\beta\lambda)\,V_{\theta_{i},k}^{*},$
(13) $\displaystyle W_{L}(\theta_{i},\theta_{j})$
$\displaystyle=\sum_{m}V_{\theta_{i},m}\,I_{m}(\beta(2-\Delta))\,V_{\theta_{j},m}^{*}$
(14)

where $U_{\phi_{i},n}=\mathrm{e}^{in\phi_{i}}$,
$V_{\theta_{i},m}=\mathrm{e}^{im\theta_{i}}$, and the diagonal
$I_{n}(x),I_{m}(x),I_{k}(x)$ are the modified Bessel functions of the first
kind.

Then, as displayed in Fig. 7(d), we can simply integrate out the $U(1)$ phase
variables at each site by a Fourier transformation

$\displaystyle\int\frac{d\phi_{i}}{2\pi}U_{\phi_{i},n_{1}}U_{\phi_{i},n_{2}}U_{\phi_{i},n_{3}}^{*}U_{\phi_{i},n_{4}}^{*}U_{\phi_{i},-3k}$
$\displaystyle=\delta_{n_{1}+n_{2}}^{n_{3}+n_{4}+3k}\equiv\delta_{U},$ (15)
$\displaystyle\int\frac{d\theta_{i}}{2\pi}V_{\theta_{i},m_{1}}V_{\theta_{i},m_{2}}V_{\theta_{i},m_{3}}^{*}V_{\theta_{i},m_{4}}^{*}V_{\theta_{i},k}$
$\displaystyle=\delta_{m_{1}+m_{2}+k}^{m_{3}+m_{4}}\equiv\delta_{L}.$ (16)

In this way, the continuous variables $\phi$ and $\theta$ are now transformed
into the discrete bond indices $n,m$ and $k$. Moreover, we evenly divide the
diagonal $I$ tensors and take a contraction of the $\sqrt{I}$ tensors
connecting to the respective $\delta$ tensors. As a result, the double layer
TN representation of the partition function in the main text is obtained as
shown in Fig. 7(b). The local tensor $O$ can then be obtained by grouping the
corresponding tensor indices in the upper and lower layers and summing out the
inter-layer $k$ indices

$O_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}=\sum_{k}\left(\prod_{l=1}^{4}I_{n_{l}}(\beta\Delta)I_{m_{l}}(\beta(2-\Delta))\right)^{1/2}I_{k}(\beta\lambda)\,\delta_{n_{1}+n_{2}}^{n_{3}+n_{4}+3k}\delta_{m_{1}+m_{2}+k}^{m_{3}+m_{4}}.$
(17)

And finally the partition function is now successfully converted into a
uniform TN on the 2D square lattice as displayed in Fig. 7(c)

$Z=\mathrm{tTr}\prod_{i}O_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}(i)$
(18)

where “tTr” denotes the sum over all auxiliary bond indices. The symmetries of
the original model is well preserved in the TN representation. It is evident
that the global $U(1)\times Z_{3}$ invariance of the bilayer model is encoded
in each local tensor since
$O_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}\neq 0$ only if
$n_{1}+3m_{1}+n_{2}+3m_{2}=n_{3}+3m_{3}+n_{4}+3m_{4}$.

As a matter of fact, the equivalence of the Hamiltonian (6) and (7) can be
easily verified from the perspective of local tensors. The duality
transformation can be performed on $\varphi$ and $\vartheta$ in the same way
for Hamiltonian (6), and the integration will generate the same local tensor
as

$\displaystyle\int\frac{d\varphi_{i}}{2\pi}\mathrm{e}^{i2n_{1}\varphi_{i}}\mathrm{e}^{i2n_{2}\varphi_{i}}\mathrm{e}^{-i2n_{3}\varphi_{i}}\mathrm{e}^{-i2n_{4}\varphi_{i}}\mathrm{e}^{-i6k\varphi_{i}}$
$\displaystyle=\delta_{2n_{1}+2n_{2}}^{2n_{3}+2n_{4}+6k}=\delta_{n_{1}+n_{2}}^{n_{3}+n_{4}+3k}\equiv\delta_{U},$
(19)
$\displaystyle\int\frac{d\vartheta_{i}}{2\pi}\mathrm{e}^{i6m_{1}\vartheta_{i}}\mathrm{e}^{i6m_{2}\vartheta_{i}}\mathrm{e}^{-i6m_{3}\vartheta_{i}}\mathrm{e}^{-i6m_{4}\vartheta_{i}}\mathrm{e}^{-i6k\vartheta_{i}}$
$\displaystyle=\delta_{6m_{1}+6m_{2}+6k}^{6m_{3}+6m_{4}}=\delta_{m_{1}+m_{2}+k}^{m_{3}+m_{4}}\equiv\delta_{L}.$
(20)

Therefore all the physical quantities deducted from the partition function
should be the same for both (6) and (7).

Figure 7:  (a) Tensor network representation of the partition function with
Boltzmann weights represented as a tensor product of local bonds. (b) The
double layer TN derived in the main text. (c) The 2D uniform tensor network
representation of the partition function. (d) The eigenvalue decompositions of
the symmetric matrix $W$ and the transformation to discrete degrees of
freedom.

##  II Appendix B: Algorithm of uniform matrix product states

Although it is often straightforward to write down the TN representation for
the many-body problem, the real challenge for numerical simulations comes from
the contraction of the whole tensor network which is proved to be NP hard when
the network is complicated. Fortunately, a lot of algorithms have been
proposed to contract an infinite translation-invariant TN. One of the best
practices to address the challenge is the algorithm of uniform matrix product
statesZauner-Stauber et al. (2018); Fishman et al. (2018); Vanderstraeten et
al. (2019a), which greatly speeds up the process of finding the leading
eigenvector of the transfer matrix based on a variational ansatz.

Figure 8:  (a) Eigen-equation for the fixed-point uMPS of the transfer
operator $T$. (b) The uniform representation and two equivalent mixed
canonical forms of the fixed-point uMPS. (c) The canonical conditions of the
fixed-point local tensors and the isometric gauge transformation between them.
(d) Eigen-equations to update the left and right environmental eigenvectors of
the channel operators. (e) Eigen-equations to obtain the central tensors based
on the new environment. (f) Polar decompositions to get the left and right
canonical local tensors from the central tensors. (g) The reduced density
matrix obtained by tracing out the right part of the uMPS. (h) Expectation
value of a local observable by contracting the leading eigenvectors of the
channel operators. (i) Two-point correlation function expressed by contracting
a sequence of channel operators.

Here we provide a brief review of this efficient algorithm based on a set of
optimized eigen-solvers. In the thermodynamic limit, the partition function is
determined by the dominant eigenvalues of the 1D quantum transfer operator
$\hat{T}$, whose eigen-equation is shown in Fig. 8(a)

$\hat{T}|\Psi(A)\rangle=\Lambda_{\max}|\Psi(A)\rangle.$ (21)

For a translation-invariant system, the leading eigenvector of the transfer
matrix can be expressed as the uniform matrix product states (uMPS):

$|\Psi(A)\rangle=\sum_{\\{(n_{i}m_{i})\\}}\mathrm{Tr}\left(\cdots
A^{(n_{1}m_{1})}A^{(n_{2}m_{2})}A^{(n_{3}m_{3})}\cdots\right),$ (22)

which is constructed by infinitely many repetitions of rank-3
$A_{\alpha\beta}^{(n_{i}m_{i})}$ with auxiliary bond dimension
$\alpha,\beta=1,2,\ldots,D$. Here, $D$ is the upper bond dimension which
controls the accuracy of the approximations. As shown in Fig. 8(b), to fix the
gauge of the uMPS, we should bring the uMPS into the equivalent mixed
canonical forms as

$|\Psi\rangle=\sum_{\alpha,\beta=1}^{D}C_{\alpha,\beta}|\Psi_{\alpha}^{[-\infty,(n_{j}m_{j})]}\rangle\otimes|\Psi_{\beta}^{[(n_{j+1}m_{j+1}),+\infty]}\rangle,$
(23)

and

$|\Psi\rangle=\sum_{\alpha,\beta=1}^{D}\left(A_{C}^{(n_{j}m_{j})}\right)_{\alpha,\beta}|\Psi_{\alpha}^{[-\infty,(n_{j-1}m_{j-1})]}\rangle\otimes|\Psi_{\beta}^{[(n_{j+1}m_{j+1}),+\infty]}\rangle,$
(24)

where $|\Psi_{\alpha}^{[-\infty,(n_{j}m_{j})]}\rangle$ and
$|\Psi_{\beta}^{[(n_{j+1}m_{j+1}),+\infty]}\rangle$ are the left and right
orthonormal basis comprised of the left and right canonical local tensors

$\displaystyle|\Psi_{\alpha}^{[-\infty,(n_{j}m_{j})]}\rangle$
$\displaystyle=\sum_{\\{(n_{i}m_{i})\\}}\mathrm{Tr}\left(\cdots
A_{L}^{(n_{j-1}m_{j-1})}A_{L}^{(n_{j}m_{j})}\right),$ (25)
$\displaystyle|\Psi_{\beta}^{[(n_{j+1}m_{j+1}),+\infty]}\rangle$
$\displaystyle=\sum_{\\{(n_{i}m_{i})\\}}\mathrm{Tr}\left(A_{R}^{(n_{j+1}m_{j+1})}A_{R}^{(n_{j+2}m_{j+2})}\cdots\right).$
(26)

The left and right orthonormal tensors $A_{L}$ and $A_{R}$ satisfy the
isometric constraints as displayed in Fig. 8(c)

$\displaystyle\sum_{(nm)}\left(A_{L}^{(nm)}\right)^{\dagger}A_{L}^{(nm)}=\mathcal{I},$
(27)
$\displaystyle\sum_{(nm)}A_{R}^{(nm)}\left(A_{R}^{(nm)}\right)^{\dagger}=\mathcal{I}.$
(28)

Another fixed point equation is obtained by pulling the central tensor from
right to left

$A_{C}^{(nm)}=A_{L}^{(nm)}C=CA_{R}^{(nm)}.$ (29)

The VUMPS algorithm contains three key steps summarized in Fig. 8(d)-(f) to
find the fixed-point tensors $A_{L}$, $A_{R}$ and $C$ of the transfer operator
$\hat{T}$. These three steps (d)-(f) are repeated sequentially until the error
is less than a given convergence threshold.

(1). As shown in Fig. 8(d), we calculate the left and right environments
$F_{L}$ and $F_{R}$ by solving the eigenvalue equation with the Arnoldi
method:

$\displaystyle\mathbb{T}_{L}F_{L}$ $\displaystyle=\lambda F_{L},$ (30)
$\displaystyle\mathbb{T}_{R}F_{R}$ $\displaystyle=\lambda F_{R},$ (31)

where $\mathbb{T}_{L}$ and $\mathbb{T}_{R}$ are the channel operators of
$A_{L}$ and $A_{R}$.

(2). As shown in Fig. 8(e), we calculate the central tensors $A_{C}$ and $C$
based on the new environment by the Arnoldi method:

$\displaystyle H_{A_{C}}A_{C}$ $\displaystyle=\lambda A_{C},$ (32)
$\displaystyle H_{C}C$ $\displaystyle=C,$ (33)

where $H_{A_{C}}$ and $H_{C}$ are the effective environment comprised of the
updated $F_{L}$ and $F_{R}$.

(3). As shown in Fig. 8(f), we use the left and right polar decompositions of
new $A_{C}$ and $C$ to update $A_{L}$ and $A_{R}$ satisfying the isometric
conditions

$\displaystyle A_{C}$
$\displaystyle=U_{A_{C}}^{[l]}P_{A_{C}}^{[l]}=P_{A_{C}}^{[r]}U_{A_{C}}^{[r]},$
(34) $\displaystyle C$
$\displaystyle=U_{C}^{[l]}P_{C}^{[l]}=P_{C}^{[r]}U_{C}^{[r]},$ (35)

where $U_{A_{C}}^{[l]}$, $U_{A_{C}}^{[r]}$, $U_{C}^{[l]}$ and $U_{C}^{[r]}$
are unitary matrices, and the matrices $P_{A_{C}}^{[l]}$, $P_{A_{C}}^{[r]}$,
$P_{C}^{[l]}$ and $P_{C}^{[r]}$ are hermitian and positive. Finally, we obtain

$\displaystyle A_{L}$ $\displaystyle=U_{A_{C}}^{[l]}U_{C}^{[l]\dagger},$ (36)
$\displaystyle A_{R}$ $\displaystyle=U_{C}^{[r]\dagger}U_{A_{C}}^{[r]}.$ (37)

Actually, the row-to-row transfer matrix should play the same role as the
matrix product operator for the 1D quantum spin chains whose logarithmic form
can be mapped to a 1D quantum system with complicated interactions. With such
a correspondence, the finite-temperature properties of the 2D statistical
problem are related to a 1D quantum model at zero temperature and all the low-
temperature physics can be achieved from the fixed-point uMPS.

##  III Appendix C: Calculations of physical quantities

Once the fixed-point uMPS is achieved, various physical quantities can be
precisely calculated because the 2D network can be easily squeezed into a 1D
chain of channel operators based on the fixed uMPS.

By mapping the transfer matrix to a 1D quantum transfer operator, we bring to
it modern concepts of quantum entanglement, while also taking advantage of the
sharp criterion for detecting various phase transitions. As shown in Fig.
8(g), we perform a bipartition on the MPS and trace out the right part using
the right canonical condition. Since the
$|\Psi^{[-\infty,(n_{j}m_{j})]}\rangle$ also forms a unitary basis, the
reduced density matrix $\rho_{L}$ of the left part is expressed as

$\rho_{L}=\mathrm{Tr}_{R}|\Psi\rangle\langle\Psi|=|\Psi^{[-\infty,(n_{j}m_{j})]}\rangle
CC^{\dagger}\langle\Psi^{[-\infty,(n_{j}m_{j})]}|\sim CC^{\dagger}.$ (38)

And the entanglement entropyVidal et al. (2003) is readily obtained by

$S_{E}=-\mathrm{Tr}\,\left(\rho_{L}\ln\rho_{L}\right)=-\sum_{\alpha=1}^{D}s_{\alpha}^{2}\ln
s_{\alpha}^{2},$ (39)

where $s_{\alpha}$ are the singular values of the $C$ matrix.

With the fixed-point uMPS, the contraction of an infinite 2D tensor network
can be reduced to the trace of an infinite 1D chain of channel operators. The
expectation value of of a single-site observable $m(\phi_{j},\theta_{j})$

$\langle
m(\phi_{j},\theta_{j})\rangle=\frac{1}{Z}\prod_{i}\iint\frac{\mathrm{d}\phi_{i}\mathrm{d}\theta_{i}}{\left(2\pi\right)^{2}}\prod_{\langle
i,j\rangle}\mathrm{e}^{\beta\Delta\cos(\phi_{i}-\phi_{j})}\mathrm{e}^{\beta(2-\Delta)\cos(\theta_{i}-\theta_{j})}\mathrm{e}^{\beta\lambda\cos(\theta_{i}-3\phi_{i})}m(\phi_{j},\theta_{j})$
(40)

can be obtained by the contraction of a reduced network containing an impurity
tensor thanks to the fixed point environment tensors $F_{L}$ and $F_{R}$ as
displayed in Fig. 8(h)

$\langle
m\rangle=\mathrm{Tr}\left(\cdots\mathbb{T}_{L}\mathbb{T}_{L}\mathbb{T}_{M}\mathbb{T}_{R}\mathbb{T}_{R}\cdots\right)=\langle
F_{L}|\mathbb{T}_{M}|F_{R}\rangle$ (41)

where $\mathbb{T}_{M}=\mathrm{tTr}(A_{C}\otimes M\otimes\bar{A}_{C})$. And the
corresponding impurity tensor $M$ can be achieved by

$M_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}=\sum_{k}\left(\prod_{l=1}^{4}I_{n_{l}}(\beta\Delta)I_{m_{l}}(\beta(2-\Delta))\right)^{1/2}I_{k}(\beta\lambda)\iint\frac{\mathrm{d}\phi\mathrm{d}\theta}{(2\pi)^{2}}\,\mathrm{e}^{i\phi(n_{1}+n_{2}-n_{3}-n_{4}-3k)}\mathrm{e}^{i\theta(m_{1}+m_{2}+k-m_{3}-m_{4})}m(\phi,\theta).$
(42)

For instance, the corresponding impurity tensors for
$m_{\phi}=\mathrm{e}^{i\phi}$, $m_{\theta}=\mathrm{e}^{i\theta}$ and
$m_{\sigma}=\mathrm{e}^{i(\theta/3-\phi)}$ are

$\displaystyle\left(M_{\phi}\right)_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}$
$\displaystyle=\sum_{k}\left(\prod_{l=1}^{4}I_{n_{l}}(\beta\Delta)I_{m_{l}}(\beta(2-\Delta))\right)^{1/2}I_{k}(\beta\lambda)\,\delta_{n_{1}+n_{2}+1}^{n_{3}+n_{4}+3k}\delta_{m_{1}+m_{2}+k}^{m_{3}+m_{4}},$
(43)
$\displaystyle\left(M_{\theta}\right)_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}$
$\displaystyle=\sum_{k}\left(\prod_{l=1}^{4}I_{n_{l}}(\beta\Delta)I_{m_{l}}(\beta(2-\Delta))\right)^{1/2}I_{k}(\beta\lambda)\,\delta_{n_{1}+n_{2}}^{n_{3}+n_{4}+3k}\delta_{m_{1}+m_{2}+k+1}^{m_{3}+m_{4}},$
(44)
$\displaystyle\left(M_{\sigma}\right)_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}$
$\displaystyle=\sum_{k}\left(\prod_{l=1}^{4}I_{n_{l}}(\beta\Delta)I_{m_{l}}(\beta(2-\Delta))\right)^{1/2}I_{k}(\beta\lambda)\,\delta_{n_{1}+n_{2}}^{n_{3}+n_{4}+3k+1}\frac{\sqrt{3}}{2\pi}\frac{(-1)^{l}}{l+1/3},$
(45)

where $l=m_{1}+m_{2}+k-m_{3}-m_{4}$.

Moreover, the two-point correlation function $G(r)=\langle
m_{j}m_{j+r}\rangle$ between local observables $m_{j}(\phi_{j},\theta_{j})$
and $m_{j+r}(\phi_{j+r},\theta_{j+r})$

$G(r)=\frac{1}{Z}\prod_{i}\iint\frac{\mathrm{d}\phi_{i}\mathrm{d}\theta_{i}}{\left(2\pi\right)^{2}}\prod_{\langle
i,j\rangle}\mathrm{e}^{\beta\Delta\cos(\phi_{i}-\phi_{j})}\mathrm{e}^{\beta(2-\Delta)\cos(\theta_{i}-\theta_{j})}\mathrm{e}^{\beta\lambda\cos(\theta_{i}-3\phi_{i})}m_{j}(\phi_{j},\theta_{j})m_{j+r}(\phi_{j+r},\theta_{j+r}).$
(46)

can be calculated in the same way as a contraction of a 1D chain with two
local impurity tensors $M_{j}$ and $M_{j+r}$ as shown in Fig. 8(i)

$G(r)=\langle
F_{L}|\mathbb{T}_{M_{j}}\underbrace{\mathbb{T}_{R}\cdots\mathbb{T}_{R}}_{r-1}\mathbb{T}_{M_{j+r}}|F_{R}\rangle$
(47)

where $\mathbb{T}_{M_{j}}=\mathrm{tTr}(A_{C}\otimes M_{j}\otimes\bar{A}_{C})$
and $\mathbb{T}_{M_{j+r}}=\mathrm{tTr}(A_{R}\otimes
M_{j+r}\otimes\bar{A}_{R})$.

##  IV Appendix D: Particular calculation of spin stiffness

Although there is no local order parameter due to the absence of the $U(1)$
symmetry breaking, the onset of the superfluidity of the hexatic-nematic XY
model can be described by the spin stiffness $\rho_{s}$Fisher et al. (1973);
Nelson and Kosterlitz (1977), which characterizes the change of free energy in
response to a small uniform phase twist $v$ along a reference direction

$\rho_{s}=\left.\frac{\partial^{2}f_{v}}{\partial
v^{2}}\right|_{v=0}=-\frac{1}{N\beta}\left[\frac{1}{Z_{v}}\frac{\partial^{2}Z_{v}}{\partial
v^{2}}-\left(\frac{1}{Z_{v}}\frac{\partial Z_{v}}{\partial
v}\right)^{2}\right]_{v=0}$ (48)

where $f_{v}=-\frac{1}{N\beta}\ln Z_{v}$ is the free energy per site.

Since the two components of $\phi$ and $\theta$ are coupled together via a
relevant inter-component term, the contribution of the total spin stiffness of
both fields should be taken into account. In order to preserve the $U(1)\times
Z_{3}$ invariance of the model, the phase twist should be applied as

$(\phi_{i},\theta_{i})\to(\phi_{i}+\vec{v}\cdot\vec{r}_{i},\theta_{i}+3\vec{v}\cdot\vec{r}_{i})$
(49)

where $\vec{r}_{i}$ is the position vector for the lattice site $i$ and
$\vec{v}$ is a constant vector to increase the phase difference across each
neighboring hexatic spins $3$ times larger than nematic spins. When a global
twist is applied along $y$-axis, the phase differences between two nearest-
neighbor sites in $y$-direction are increased by $v$ for $\phi$ field and $3v$
for $\theta$ field, and the corresponding Hamiltonian is changed as

$\displaystyle\begin{split}H_{v}=&-\Delta\sum_{\langle
i,j\rangle_{x}}\cos(\phi_{i}-\phi_{j})-(2-\Delta)\sum_{\langle
i,j\rangle_{x}}\cos(\theta_{i}-\theta_{j})-\lambda\sum_{i}\cos(\theta_{i}-3\varphi_{i})\\\
&-\Delta\sum_{\langle
i,j\rangle_{y}}\cos(\phi_{i}-\phi_{j}+v)-(2-\Delta)\sum_{\langle
i,j\rangle_{x}}\cos(\theta_{i}-\theta_{j}+3v).\end{split}$ (50)

In the same time, the TN representation of the partition function is modified
as

$Z_{v}=\mathrm{tTr}\prod_{i}(O_{v})_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}(i),$
(51)

where an additional phase factor is appended on the vertical leg of the local
tensor

$(O_{v})_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}=\mathrm{e}^{i(n_{1}+3m_{1})v}O_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}.$
(52)

The precision of the direct second differentiation on the partition function
often suffers from a finite differentiation step $v$. Alternatively, the
calculations of the spin stiffness can be improved by directly contracting the
TN containing the differentiated terms independent of $v$. The first
derivative of the partition function with respect to $v$ introduces an
impurity tensor $R$ in the original TN as displayed in Fig. 9(a), where

$R_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}=\left.\frac{\partial
O_{v}}{\partial
v}\right|_{v=0}=i(n_{1}+3m_{1})O_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}.$
(53)

Because the partition function is expressed as a tensor product of $O_{v}$
tensors, the second derivative is composed by two parts: two tensors
separately differentiated once and one tensor differentiated twice. The
configurations for the second derivative are summed in Fig. 9(b), where

$S_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}=\left.\frac{\partial^{2}O_{v}}{\partial
v^{2}}\right|_{v=0}=-(n_{1}+3m_{1})^{2}O_{n_{1}m_{1},n_{2}m_{2}}^{n_{3}m_{3},n_{4}m_{4}}.$
(54)

In the previous section, we evaluate the expectation value of two-point
observables when two impurity tensors locate along a horizontal or vertical
line using the fixed point MPS. However, for the spin stiffness, because two
$R$ tensors are located at different rows and columns, we cannot simply
contract the whole network by squeezing the 2D TN into a 1D chain.
Fortunately, another contraction strategyVanderstraeten et al. (2015, 2016)
has been proposed which is based on the channel environment where the linear
fixed points can be bent into corner-shaped transfer matrices as the effective
“corner environment” in approximation to infinite quarter-planes as depicted
in Fig. 9(c). The fixed point equation for the top left corner is shown in
Fig. 9(d), where the corner-shaped environment is approximated by a bent MPS
comprised of two half infinite boundary MPS connected by an corner matrix $K$.
Once the corner-shaped fixed points in different quarter-planes are obtained,
arbitrary two-point expectation values can be calculated by contracting the
network containing $R_{i}$ and $R_{j}$ from four corners. Finally the network
is reduced into a shape of cross as shown in Fig. 9(e) with the help of the
fixed-point eigenvectors $F_{L}$ and $F_{R}$.

Figure 9:  (a) The tensor network representation of the first derivative of
the partition function. (b) The tensor network representation of the second
derivative of the partition function. (c) The contraction of an infinite
tensor network using corner environments from four corners. (d) The eigen-
equation of the top left corner-shaped fixed point. (e) The computation of a
two-point function expressed by contracting the channel operators and the
corner tensors.

All contributions of different configurations of the $R$ tensors should be
taken into consideration which is achieved by moving the $R$ tensors
independently in horizontal and vertical directions. The sum of the infinite
geometric series of two observables is a bit more involved. For exponentially
decaying observables, long-range observables can be evaluated by inverting the
corresponding channel operator as a sum of channel operators

$\sum_{n=0}^{\infty}\mathbb{T}_{L}^{n}\simeq\left[\mathcal{I}-(\mathbb{T}_{L}-|F_{r}\rangle\langle
F_{L}|)\right]^{-1},$ (55)

where $F_{r}$ is the right dominant eigenvector of $\mathbb{T}_{L}$ whose
dominant eigenvalue is renormalized to $1$.

The discontinuous jump of the spin stiffness is a characteristic feature for
unbinding of topological defects above the BKT transition. In the main text,
we show the two-step jump of the superfluid stiffness for both hexatic
($\Delta<1$) and nematic ($\Delta>1$) regime. In contrast to the case for
small $\Delta$ where the stiffness develops a cusp point with a small jump at
the lower transition temperature due to the unbinding of vortices with charge
$q_{\phi}=1$, in the large $\Delta$ regime the stiffness drops dramatically at
the lower transition temperature $T_{c1}$. The different behavior of the
stiffness jump at lower transition temperatures between small and large
$\Delta$ cases should be related to the dissociations of different kinds of
vortices in the $\phi$ and $\theta$ fields. Since the phase twist $v$ is
applied to the hexatic $\theta$ field $3$ times larger than to the nematic
$\phi$ field, the ratio of the contribution to the total stiffness is roughly
$\Delta:9(2-\Delta)$ between $\phi$ and $\theta$ fields for a weak coupling
$\lambda$Drouin-Touchette et al. (2022). That is why the first drop of the
stiffness is almost $6$ times larger than the second drop upon increasing
temperature along $\Delta=1.2$. Such a great drop is different from the
strongly coupled case with large $\lambda$ where the topological defects do
not result in a jump of the total stiffness but a gradual decrease at the
lower crossover. Further investigations into the correlation properties tells
us that the extremely weak coherence between vortices in $\theta$ field above
$T_{c1}$ should account for the dramatic jump in total stiffness even though
the $\theta$ field still remains quasi-long-range ordered (quasi-LRO).

##  V Appendix E: Correlation functions and their correlation lengths

To reveal the nature of the topological excitations in different phases, we
calculate various two-point correlation functions and the corresponding
correlation lengths. Here, we present more details for the correlation
properties for intra- and inter-component observables. Our results are
discussed in the hexatic ($\Delta<1$) and nematic ($\Delta>1$) regimes,
separately.

The correlation functions for $\Delta=0.8$ are summarized in Fig. 10 and the
corresponding correlation lengths are displayed in the first column of Fig.
12.

Figure 10:  Various two-point correlation functions are calculated under the
intra-component coupling ratio $\Delta=0.8$ and the inter-component coupling
$\lambda=0.1$. Different columns denote different values of $T=0.6J$, $1.0J$
and $1.4J$ as indicated.

In the $T=0.6J$ column below $T_{c1}$, the intra-component correlation
functions $G_{\phi}(r)$ in Fig. 10 (a) and $G_{\theta}(r)$ in Fig. 10 (g)
display the algebraic behavior, indicating the bindings of $q_{\phi}=1$ and
$q_{\theta}=1$ vortices in the $\phi$ and $\theta$ fields, respectively.
Moreover, the power law behavior of $G_{\phi\theta}(r)$ in Fig. 10 (j) implies
that the vortices and anti-vortices also bind in the inter-field. Therefore a
fully phase-coherent state of the coupled system is established in the low
temperature phase. Besides, the inter-component Potts variable develops a true
LRO and the correlation function $G_{\sigma}(r)$ in Fig. 10 (m) becomes a
constant for large $r$.

In the $T=1.0J$ column of the intermediate temperatures, the most interesting
physics happens in the $\phi$ field. When cooling down the system from high
temperature disordered phase, the buildup of correlations in the $\theta$
field first occurs at $T_{c2}$ with an algebraically decaying correlation
function $G_{\theta}(r)$ in Fig. 10 (h), corresponding to the binding of
hexatic vortices with charge $q_{\theta}=1$. For the fractional vortex paired
phase between $T_{c1}$ and $T_{c2}$, the correlation function $G_{\phi}(r)$
decays exponentially, while the correlation function $G_{3\phi}(r)$ of the
fractionalized vortices exhibit an algebraic behavior.

As displayed in Fig. 10(b) and (e), the direct comparison between difference
correlation behavior demonstrates that the integer nematic vortices
$q_{\phi}=1$ in the $\phi$ field are fractionalized into fractional nematic
vortices with $q_{\phi}=1/3$ due to the presence of the inter-component
$Z_{3}$ degrees of freedom. Since the fractional nematic vortices are
pointlike topological defects around which the phase angles of spins wind by
$2\pi/3$, each fractional nematic vortices should be at the end of a domain-
wall string across which the $\phi$ field twists by $2\pi/3$. The existence of
multiple Potts domains destroys the phase coherence in $\phi$ field and leads
to disorder in inter-component $\sigma$ variables. As we can see, both
$G_{\phi\theta}(r)$ in Fig. 10(k) and $G_{\sigma}(r)$ in Fig. 10 (n) display
the exponential decay. We should point out that the correlation length
$\xi_{\phi}$ extracted from $G_{\phi}(r)$ could be fitted by an exponentially
divergent form $\xi(T)\propto\exp(b/\sqrt{T-T_{C}})$ above $T_{c1}$ as shown
in Fig. 12 (a) with $b>0$ as a key feature of the BKT transition, which
implies that the transition at $T_{c1}$ is a hybrid BKT and Potts transition.
Furthermore, as shown in Fig. 12 (e), the correlation length $\xi_{\sigma}$
for the inter-component Potts variable extracted from

$\mathrm{e}^{-r/\xi_{\sigma}}\propto\langle\cos(\sigma_{i}-\sigma_{i+r})\rangle-\langle\cos(\sigma_{i})\rangle\langle\cos(\sigma_{i+r})\rangle$
(56)

is well fitted by $\xi_{\sigma}\propto 1/|T-T_{c1}|^{5/6}$ in agreement with
2D Potts universality class.

In the $T=1.4J$ column above $T_{c2}$, all the correlation functions decay
exponentially and the system is disordered. When approaching $T_{c2}$ from
high temperatures, the relevant correlation length $\xi_{\theta}$ is extracted
and displayed in Fig. 12(c), which can be well fitted by an exponentially
divergent form of the BKT transition.

It is interesting to see that the correlation functions of $G_{\theta}(r)$ and
$G_{3\phi}(r)$ share the same exponents while $G_{\phi}(r)$ and
$G_{\sigma}(r)$ share the same exponents. Such behavior can be interpreted as
another proof of the relevance of the inter-component couplings. Due to the
relative locking, the $\theta$ and $3\phi$ fields have the same long-
wavelength behavior and the correspondence $q_{\theta}=3q_{\phi}$ between
vortex charges is established.

However, the correlation functions for $\Delta=1.2$ are summarized in Fig. 11
and the corresponding correlation lengths are displayed in the second column
of Fig. 12.

Figure 11:  Various two-point correlation functions are calculated under the
intra-component coupling ratio $\Delta=1.2$ and the inter-component coupling
$\lambda=0.1$. Different columns denote different values of $T=0.6J$, $1.0J$
and $1.4J$ as indicated.

In the $T=0.6J$ column below $T_{c1}$, all intra-component correlation
functions show a power law behavior and the correlation function for inter-
component Potts variable becomes a constant for long distance as displayed in
Fig. 11(m). In fact, the low temperature phases in the hexatic and nematic
regime share the same physics.

In the $T=1.0J$ column of the intermediate temperatures, the power law
behavior of $G_{\theta}(r)$ in Fig. 11(h) demonstrates that the quasi-LRO in
the $\theta$ field survives through the transition at $T_{c1}$. A direct
comparison between the amplitude of $G_{\theta}(r)$ at two sides of $T_{c1}$
in the main text shows that the coherence between spins in the $\theta$ field
is greatly suppressed by $3$ orders of magnitude across the transition at
$T_{c1}$. Although the $\theta$ field is still quasi-LRO, the extremely weak
coherence above $T_{c1}$ could account for the dramatic jump in total
stiffness at $T_{c1}$. The reason for the absence of unbinding transition in
the $\theta$ field is that the onset of algebraic correlations in the $\phi$
field in Fig. 11(b) leads to the quasi-LRO in $\theta$ field because they are
coupled by the relevant inter-component term
$\lambda\cos(\theta_{i}-3\phi_{i})$. Since the vortices with charge
$q_{\phi}=1$ in the $\phi$ field is bound below $T_{c2}$, the dominant
topological excitations in the $\theta$ field between $T_{c1}$ and $T_{c2}$
should be composite vortex pairs with $q_{\theta}=3$. Moreover, the composite
$q_{\theta}=3$ vortices should further fractionalize into three $q_{\theta}=1$
vortices bound together as a larger extended vortex core to lower the gradient
energy cost, especially for weakly coupled cases of smaller $\lambda$. In this
way, the correlation $G_{\theta}(r)$ should be viewed as a higher order
expansion of the correlation function among composite vortices which has the
same long-wavelength behavior as $G_{\phi}(r)$.

In the $T=1.4J$ column above $T_{c2}$, all the correlation functions decay
exponentially and the system is disordered. When approaching $T_{c2}$ from
high temperatures, the corresponding correlation length $\xi_{\phi}$ for
$G_{\phi}(r)$ as well as $\xi_{\theta}$ for $G_{\theta}(r)$ are extracted and
displayed in Fig. 12(b) and (d), respectively. Both $\xi_{\phi}$ and
$\xi_{\theta}$ above $T_{c2}$ can be well fitted by an exponentially divergent
form of the BKT transition.

The most interesting physics is revealed in further investigations into the
inter-component Potts variable $\sigma_{i}$. In the strongly coupled case of
large $\lambda$, the LRO of $\sigma$ is found to be lost at $T_{c2}$ together
with the destruction of the quasi-LRO in both the hexatic and nematic fields
through a BKT transition. However, for the weekly coupled case, the order-
disorder transition in $\sigma$ variables splits into two transitions
separated by an intermediate phase with Potts quasi-LRO which is regarded as
the Potts liquid phase. As shown in Fig. 11(m), (n) and (o), the correlation
function $G_{\sigma}(r)$ exhibits a constant value, an algebraic decay and an
exponential decay at $T=0.6J$, $T=1.0J$ and $T=1.4J$, respectively. Such
unconventional two-stage phase transitions should be associated with the
separated proliferations of the $Z_{3}$ domain walls and $Z_{3}$ vortices,
similar to two continuous phase transitions in the $Z_{p}$ clock models with
$p>4$. The correlation length $\xi_{\sigma}$ is also extracted from
$\mathrm{e}^{-r/\xi_{\sigma}}\propto\langle\cos(\sigma_{i}-\sigma_{i+r})\rangle-\langle\cos(\sigma_{i})\rangle\langle\cos(\sigma_{i+r})\rangle$
as displayed in Fig. 12(f). The similar behavior of the correlation length was
also observed in $p$-state clock models ($p>4$)Li et al. (2020). For the clock
model, the correlation lengths scale exponentially with $T$ when the quasi-LRO
phase is approached from both the ordered and disordered sides, which means
that both the transitions at lower and higher temperatures belong to the BKT
class. However, as displayed in the inset of Fig. 12(f), the correlation
length for $T>T_{c2}$ (right inset) agrees well with the BKT behavior while
for $T<T_{c1}$ (left inset) the correlation length cannot fit well in the
exponential form. Therefore, it is reasonable to infer that the phase
transition at $T_{c1}$ slightly deviates from the standard BKT class.

Figure 12:  Different correlation lengths $\xi_{\phi}$, $\xi_{\theta}$ and
$\xi_{\sigma}$ are calculated for nematic $\phi$, hexatic $\theta$ and intra-
component $\sigma$ under $\lambda=0.1$, respectively. Different columns denote
different values of $\Delta=0.8$ and $\Delta=1.2$ as indicated.

##  VI Appendix F: Topological defects in $Z_{3}$ variables

To illustrate the splitting of the order-disorder transition in the $Z_{3}$
Potts variables, the corresponding topological defects are schematic depicted
in Fig. 13. Unlike the 2D Ising model where a phase transition is driven by
loop-like domain walls, the $Z_{3}$ Potts model allows the loop-like domain
walls as well as $Z_{3}$ vortices. The $Z_{3}$ Potts variables
$\sigma_{i}=0,1,2$ are defined at each vertex on a square lattice (black
lines). The topological defects, $Z_{3}$ domain walls and vortices, are
defined on the dual lattice (orange dotted line). Each segment of domain wall
corresponds to a mismatch in the Potts variables across the edge. The $Z_{3}$
vortices are determined by the vorticity around the vertex of the dual lattice

$\omega_{p}=\frac{1}{3}\sum_{\langle i,j\rangle\in\Box_{p}}\Delta_{ij}$ (57)

where $\Delta_{ij}=\sigma_{i}-\sigma_{j}$ are wrapped within $[-1,+1]$
surrounding the plaquette anticlockwise.

For the conventional Potts model, only a single Potts transition is observed
because the proliferation of two kinds of topological defects are coupled
together. However, we find that a quasi-LRO intermediate phase can appear in
the $Z_{3}$ Potts variables due to the separated excitations of different
topological defects in the weakly coupled model. The low temperature
transition is driven by the proliferation of domain walls while the high
temperature transition is resulted from excitations of free $Z_{3}$ vortices.
The effective core energy of $Z_{3}$ vortices in the $\sigma$ field of the
coupled hexatic-nematic model could be increased by a finite inter-component
coupling. Actually, a finite $\lambda$ makes it more costly for the
excitations of free vortices in the $\theta$ field which would proliferate
through a BKT transition at $T_{c1}$ if $\lambda=0$. On the other hand, for
large $\lambda$ limit, the energy of domain walls increases dramatically which
forbids the excitations of domain walls before the proliferation of free
$Z_{3}$ vortices and only one transition is left. Since the $\sigma$ variable
is determined by $\theta$ and $\phi$ field simultaneously, the relation
between the inter-component coupling $\lambda$ and core energy of $Z_{3}$
vortices is more complicated. An estimation made from an effective long-
wavelength modelDrouin-Touchette et al. (2022) gives that the energies of
domain walls are proportional to $\sqrt{\lambda\Delta}$ while the core
energies of the vortices in the phase field are proportional to $\Delta$.
Therefore a small $\lambda$ is required in the coupled hexatic-nematic model
to achieve a relatively high ratio of the core energy of the $Z_{3}$ vortices
compared to domain walls. An investigation of the changes of the phase
structure on different $\lambda$ is presented in the main text.

Figure 13:  Schematic representations of the topological defects of the
$Z_{3}$ variables in (a) the Potts ordered phase, (b) the Potts liquid phase,
and (c) the Potts disordered phase, respectively. The Potts variables
($\sigma_{i}=0,1,2$) are define on the original lattice while domain walls
(orange thick lines), vortices (red dots), and antivortices (blue dots) are
defined on the dual lattice.

## References

  * Berezinsky (1971) V. Berezinsky, Sov. Phys. JETP 32, 493 (1971).
  * Kosterlitz and Thouless (1973) J. M. Kosterlitz and D. J. Thouless, Journal of Physics C: Solid State Physics 6, 1181 (1973), URL https://doi.org/10.1088.
  * Kosterlitz (1974) J. M. Kosterlitz, Journal of Physics C: Solid State Physics 7, 1046 (1974), URL https://doi.org/10.1088.
  * Halperin and Nelson (1978) B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. 41, 121 (1978), URL https://link.aps.org/doi/10.1103/PhysRevLett.41.121.
  * Nelson and Halperin (1979) D. R. Nelson and B. I. Halperin, Phys. Rev. B 19, 2457 (1979), URL https://link.aps.org/doi/10.1103/PhysRevB.19.2457.
  * Young (1979) A. P. Young, Phys. Rev. B 19, 1855 (1979), URL https://link.aps.org/doi/10.1103/PhysRevB.19.1855.
  * Kosterlitz (2016) J. M. Kosterlitz, Reports on Progress in Physics 79, 026001 (2016), URL https://dx.doi.org/10.1088/0034-4885/79/2/026001.
  * Ryzhov et al. (2017) V. N. Ryzhov, E. E. Tareyeva, Y. D. Fomin, and E. N. Tsiok, Physics-Uspekhi 60, 857 (2017), URL https://dx.doi.org/10.3367/UFNe.2017.06.038161.
  * Strandburg (1988) K. J. Strandburg, Rev. Mod. Phys. 60, 161 (1988), URL https://link.aps.org/doi/10.1103/RevModPhys.60.161.
  * Bernard and Krauth (2011) E. P. Bernard and W. Krauth, Phys. Rev. Lett. 107, 155704 (2011), URL https://link.aps.org/doi/10.1103/PhysRevLett.107.155704.
  * Chae et al. (2012) S. C. Chae, N. Lee, Y. Horibe, M. Tanimura, S. Mori, B. Gao, S. Carr, and S.-W. Cheong, Phys. Rev. Lett. 108, 167603 (2012), URL https://link.aps.org/doi/10.1103/PhysRevLett.108.167603.
  * Radzihovsky et al. (2008) L. Radzihovsky, P. B. Weichman, and J. I. Park, Annals of Physics 323, 2376 (2008), ISSN 0003-4916, URL https://www.sciencedirect.com/science/article/pii/S0003491608000742.
  * Shahbazi and Ghanbari (2006) F. Shahbazi and R. Ghanbari, Phys. Rev. E 74, 021705 (2006), URL https://link.aps.org/doi/10.1103/PhysRevE.74.021705.
  * de Parny et al. (2016) L. d. F. de Parny, A. Rançon, and T. Roscilde, Phys. Rev. A 93, 023639 (2016), URL https://link.aps.org/doi/10.1103/PhysRevA.93.023639.
  * Serna et al. (2017) P. Serna, J. T. Chalker, and P. Fendley, Journal of Physics A: Mathematical and Theoretical 50, 424003 (2017), URL https://dx.doi.org/10.1088/1751-8121/aa89a1.
  * Kobayashi et al. (2019) M. Kobayashi, M. Eto, and M. Nitta, Phys. Rev. Lett. 123, 075303 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.075303.
  * Bruinsma and Aeppli (1982) R. Bruinsma and G. Aeppli, Phys. Rev. Lett. 48, 1625 (1982), URL https://link.aps.org/doi/10.1103/PhysRevLett.48.1625.
  * Aeppli and Bruinsma (1984) G. Aeppli and R. Bruinsma, Phys. Rev. Lett. 53, 2133 (1984), URL https://link.aps.org/doi/10.1103/PhysRevLett.53.2133.
  * Jiang et al. (1993) I. M. Jiang, S. N. Huang, J. Y. Ko, T. Stoebe, A. J. Jin, and C. C. Huang, Phys. Rev. E 48, R3240 (1993), URL https://link.aps.org/doi/10.1103/PhysRevE.48.R3240.
  * Jiang et al. (1996) I. M. Jiang, T. Stoebe, and C. C. Huang, Phys. Rev. Lett. 76, 2910 (1996), URL https://link.aps.org/doi/10.1103/PhysRevLett.76.2910.
  * Drouin-Touchette et al. (2022) V. Drouin-Touchette, P. P. Orth, P. Coleman, P. Chandra, and T. C. Lubensky, Phys. Rev. X 12, 011043 (2022), URL https://link.aps.org/doi/10.1103/PhysRevX.12.011043.
  * Verstraete et al. (2008) F. Verstraete, V. Murg, and J. Cirac, Advances in Physics 57, 143 (2008), eprint https://doi.org/10.1080/14789940801912366, URL https://doi.org/10.1080/14789940801912366.
  * Orús (2014) R. Orús, Annals of Physics 349, 117 (2014), ISSN 0003-4916, URL http://www.sciencedirect.com/science/article/pii/S0003491614001596.
  * Haegeman and Verstraete (2017) J. Haegeman and F. Verstraete, Annual Review of Condensed Matter Physics 8, 355 (2017), eprint https://doi.org/10.1146/annurev-conmatphys-031016-025507, URL https://doi.org/10.1146/annurev-conmatphys-031016-025507.
  * Zauner-Stauber et al. (2018) V. Zauner-Stauber, L. Vanderstraeten, M. T. Fishman, F. Verstraete, and J. Haegeman, Phys. Rev. B 97, 045145 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.045145.
  * Fishman et al. (2018) M. T. Fishman, L. Vanderstraeten, V. Zauner-Stauber, J. Haegeman, and F. Verstraete, Phys. Rev. B 98, 235148 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.98.235148.
  * Vanderstraeten et al. (2019a) L. Vanderstraeten, J. Haegeman, and F. Verstraete, SciPost Phys. Lect. Notes p. 7 (2019a), URL https://scipost.org/10.21468/SciPostPhysLectNotes.7.
  * Vanderstraeten et al. (2019b) L. Vanderstraeten, B. Vanhecke, A. M. Läuchli, and F. Verstraete, Phys. Rev. E 100, 062136 (2019b), URL https://link.aps.org/doi/10.1103/PhysRevE.100.062136.
  * Li et al. (2020) Z.-Q. Li, L.-P. Yang, Z. Y. Xie, H.-H. Tu, H.-J. Liao, and T. Xiang, Phys. Rev. E 101, 060105 (2020), URL https://link.aps.org/doi/10.1103/PhysRevE.101.060105.
  * Song and Zhang (2022a) F.-F. Song and G.-M. Zhang, Phys. Rev. Lett. 128, 195301 (2022a), URL https://link.aps.org/doi/10.1103/PhysRevLett.128.195301.
  * Granato and Kosterlitz (1986) E. Granato and J. M. Kosterlitz, Phys. Rev. B 33, 4767 (1986), URL https://link.aps.org/doi/10.1103/PhysRevB.33.4767.
  * Poderoso et al. (2011) F. C. Poderoso, J. J. Arenzon, and Y. Levin, Phys. Rev. Lett. 106, 067202 (2011), URL https://link.aps.org/doi/10.1103/PhysRevLett.106.067202.
  * Canova et al. (2014) G. A. Canova, Y. Levin, and J. J. Arenzon, Phys. Rev. E 89, 012126 (2014), URL https://link.aps.org/doi/10.1103/PhysRevE.89.012126.
  * Canova et al. (2016) G. A. Canova, Y. Levin, and J. J. Arenzon, Phys. Rev. E 94, 032140 (2016), URL https://link.aps.org/doi/10.1103/PhysRevE.94.032140.
  * Song and Zhang (2021) F.-F. Song and G.-M. Zhang, Phys. Rev. B 103, 024518 (2021), URL https://link.aps.org/doi/10.1103/PhysRevB.103.024518.
  * Song and Zhang (2022b) F.-F. Song and G.-M. Zhang, Phys. Rev. B 105, 134516 (2022b), URL https://link.aps.org/doi/10.1103/PhysRevB.105.134516.
  * Vidal et al. (2003) G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Phys. Rev. Lett. 90, 227902 (2003), URL https://link.aps.org/doi/10.1103/PhysRevLett.90.227902.
  * Pollmann et al. (2009) F. Pollmann, S. Mukerjee, A. M. Turner, and J. E. Moore, Phys. Rev. Lett. 102, 255701 (2009), URL https://link.aps.org/doi/10.1103/PhysRevLett.102.255701.
  * SM (See Supplementary Materials for the tensor network method and detailed numerical calculations.) (See Supplementary Materials for the tensor network method and detailed numerical calculations.).
  * Hübscher and Wessel (2013) D. M. Hübscher and S. Wessel, Phys. Rev. E 87, 062112 (2013), URL https://link.aps.org/doi/10.1103/PhysRevE.87.062112.
  * Wu (1982) F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982), URL https://link.aps.org/doi/10.1103/RevModPhys.54.235.
  * Einhorn et al. (1980) M. B. Einhorn, R. Savit, and E. Rabinovici, Nuclear Physics B 170, 16 (1980), ISSN 0550-3213, URL https://www.sciencedirect.com/science/article/pii/0550321380904733.
  * Bhattacharya and Ray (2016) S. Bhattacharya and P. Ray, Phys. Rev. Lett. 116, 097206 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.097206.
  * Pindak et al. (1981) R. Pindak, D. E. Moncton, S. C. Davey, and J. W. Goodby, Phys. Rev. Lett. 46, 1135 (1981), URL https://link.aps.org/doi/10.1103/PhysRevLett.46.1135.
  * Huang et al. (1981) C. C. Huang, J. M. Viner, R. Pindak, and J. W. Goodby, Phys. Rev. Lett. 46, 1289 (1981), URL https://link.aps.org/doi/10.1103/PhysRevLett.46.1289.
  * Chou et al. (1998) C.-F. Chou, A. J. Jin, S. W. Hui, C. C. Huang, and J. T. Ho, Science 280, 1424 (1998), URL https://www.science.org/doi/abs/10.1126/science.280.5368.1424.
  * Nelson and Halperin (1980) D. R. Nelson and B. I. Halperin, Phys. Rev. B 21, 5312 (1980), URL https://link.aps.org/doi/10.1103/PhysRevB.21.5312.
  * Donley et al. (2002) E. A. Donley, N. R. Claussen, S. T. Thompson, and C. E. Wieman, Nature 417, 529 (2002), URL https://doi.org/10.1038/417529a.
  * Chin et al. (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010), URL https://link.aps.org/doi/10.1103/RevModPhys.82.1225.
  * Fernandes et al. (2010) R. M. Fernandes, L. H. VanBebber, S. Bhattacharya, P. Chandra, V. Keppens, D. Mandrus, M. A. McGuire, B. C. Sales, A. S. Sefat, and J. Schmalian, Phys. Rev. Lett. 105, 157003 (2010), URL https://link.aps.org/doi/10.1103/PhysRevLett.105.157003.
  * Fernandes et al. (2019) R. M. Fernandes, P. P. Orth, and J. Schmalian, Annual Review of Condensed Matter Physics 10, 133 (2019), eprint https://doi.org/10.1146/annurev-conmatphys-031218-013200, URL https://doi.org/10.1146/annurev-conmatphys-031218-013200.
  * Fisher et al. (1973) M. E. Fisher, M. N. Barber, and D. Jasnow, Phys. Rev. A 8, 1111 (1973), URL https://link.aps.org/doi/10.1103/PhysRevA.8.1111.
  * Nelson and Kosterlitz (1977) D. R. Nelson and J. M. Kosterlitz, Phys. Rev. Lett. 39, 1201 (1977), URL https://link.aps.org/doi/10.1103/PhysRevLett.39.1201.
  * Vanderstraeten et al. (2015) L. Vanderstraeten, M. Mariën, F. Verstraete, and J. Haegeman, Phys. Rev. B 92, 201111 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.92.201111.
  * Vanderstraeten et al. (2016) L. Vanderstraeten, J. Haegeman, P. Corboz, and F. Verstraete, Phys. Rev. B 94, 155123 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.94.155123.