# Algebraic approach to quarkyoniclike configuration and stable diquarks in
dense matter

Aaron Park<EMAIL_ADDRESS>Department of Physics and Institute of
Physics and Applied Physics, Yonsei University, Seoul 03722, Korea  Su Houng
Lee<EMAIL_ADDRESS>Department of Physics and Institute of Physics and
Applied Physics, Yonsei University, Seoul 03722, Korea

###### Abstract

We study the color-spin interaction energy of a quark, a diquark and a baryon
with their surrounding baryons and/or quark matter. This is accomplished by
classifying all possible flavor and spin states of the resulting multiquark
configuration in both the flavor SU(2) and SU(3) symmetric cases. We find that
while the baryon has the lowest interaction energy when there is only a single
surrounding baryon, the quark has the lowest interaction energy when the
surrounding has more than three baryons or becomes a quark gas. As the short
range nucleon-nucleon interactions are dominated by the color-spin
interactions, our finding suggests that the baryon modes near other baryons
are suppressed due to larger repulsive energy compared to that of a quark and
thus provides a quark model basis for the quarkyoniclike phase in dense
matter. At the same time, when the internal interactions are taken into
account, and the matter density is high so that the color-spin interaction
becomes the dominant interaction, the diquark becomes the lowest energy
configuration and will thus appear in both the dense baryonic and/or quark
matter.

Introduction: Recent inputs from multi-messenger astrophysics on neutron star
equation of state suggest that there is a sudden increase in pressure as one
approaches the core of the neutron starKojo:2021wax . Such behaviour prompted
the appearance of the quarkyonic matter type of equation of state(EOS) for the
nuclear matter when the baryon density increases to several times that of the
nuclear saturation densityMcLerran:2007qj ; McLerran:2018hbz ; Jeong:2019lhv ;
Sen:2020peq ; Duarte:2020xsp ; Duarte:2020kvi ; Zhao:2020dvu ; Sen:2020qcd ;
Margueron:2021dtx .

Recently, using the quark model with color-spin interaction, we have shown
that the short distance repulsion between the quark and baryon is smaller than
that between two baryons in the lowest energy channel, and that such
ingredients naturally lead to quarkyoniclike picture for the EOS when the
baryon density reaches four to five times nuclear matter densityPark:2021hqb .
The relevance of color-spin interaction to describe the dynamics at short
distance can be verified in several contexts. In a previous
publicationPark:2019bsz , we showed that the short distance part of the
baryon-baryon interactions for various quantum numbers from the recent lattice
calculationIshii:2006ec ; Inoue:2010hs  can be well reproduced using a
constituent quark model with color-spin interaction. The importance of color-
spin interaction in nucleon-nucleon repulsion was observed earlier within the
quark cluster modelOka:2000wj . The mass splitting between hadrons with
different spin orientations, such as the delta and the nucleon or the
psuedoscalar and vector meson, are due to the color-spin interactionLee:2007tn
. Also studies on possible multiquark configurations are based on the relative
strength of the color-spin interaction of the multiquark system compared to
that of its lowest hadron thresholdPark:2017mdp ; Noh:2021lqs .

Therefore, when the quark density becomes large so that the baryons start to
overlap, the color-spin interaction together with proper application of Pauli
principle should be the important dynamics that determines the properties of
the dense matter. In this work, we will use the quark model with color-spin
interaction to study the energy of a quark, diquarks and a baryon in the dense
matter composed of baryons and in the quark matter. This is accomplished by
classifying all possible flavor and spin states of the resulting multiquark
configurations in both the flavor SU(2) and SU(3) symmetric cases, and
studying the color-spin interactions of these configurations.

We find that while the baryon has the lowest interaction energy when there is
only a single surrounding baryon, the quark has the lowest interaction energy
when the surrounding is composed of more than three baryons or becomes a quark
gas. This is an improvement over our precious calculationPark:2021hqb  as all
possible configurations are considered and the surroundings are generalized to
study cases with more than one baryons and free quarks. Our finding implies
that the baryon modes near other baryons are suppressed due to larger
repulsive energy compared to that of a quark, which is a property when
implemented into phenomenoligical EOS’s composed of quarks and baryons leads
to the appearance of quarkyoniclike phase in dense matterJeong:2019lhv ;
Sen:2020peq ; Park:2021hqb . At the same time, when the internal interactions
are taken into account, and the matter density is high, similar to the quark
density inside a baryon so that the color-spin interaction becomes the
dominant interaction, the diquark becomes the lowest energy configuration and
will thus appear in both the dense baryonic and/or quark matter.

Color-spin interaction: In the flavor SU(3) symmetric case, the color-spin
interaction factor is determined by the following form.

$\displaystyle H_{CS}$
$\displaystyle=-\sum_{i<j}^{n}\lambda^{c}_{i}\lambda^{c}_{j}\sigma_{i}\cdot\sigma_{j}$
(1) $\displaystyle=n(n-10)+\frac{4}{3}S(S+1)+4C_{F}+2C_{C},$ (2)
$\displaystyle 4C_{F}$
$\displaystyle=\frac{4}{3}(p_{1}^{2}+p_{2}^{2}+3p_{1}+3p_{2}+p_{1}p_{2}),$ (3)

where $\lambda^{c}_{i}$ is the color SU(3) Gell-Mann matrices, $C_{F}$ is the
first kind of the Casimir operator of the flavor SU(3) and $p_{i}$ is the
number of columns containing $i$ boxes in a column in Young diagram. If there
is no strange quarks, Eq. (2) reduces to the following formula.

$\displaystyle H_{CS}=\frac{4}{3}n(n-6)+\frac{4}{3}S(S+1)+4I(I+1)+2C_{C},$ (4)

where $I$ is the total isospin.

In the quark model for a hadron, there will be an extra spatial potential that
depends on the distance between two particles multiplying the individual
color-spin factor in Eq (1) so that when calculating the total contribution of
the relevant interactions to the hadron mass, the corresponding spatial
expectation value has to be taken into account for each pair. On the other
hand, assuming that all the quarks occupy similar spatial dimension and/or the
quarks are uniformly distributed, the spatial parts will be universal for all
quark pairs. Taking into account the additional factor that is inversely
proportional to the quark masses of the pair one finds that, to fit the delta-
nucleon mass difference, the overall constant should be around
$C_{B}/m_{u}^{2}\sim$ 18 MeVLee:2007tn , which should be multiplied to the
factors in this work to estimates their approximate magnitude when the quark
density is similar to that inside a nucleon. Therefore, to study the relative
stability of different configurations, assuming that the density of quarks are
uniform in the flavor symmetric limit, one can just compare the color-spin
matrix elements given in Eq. (2).

Here, we calculate the color-spin factor experienced by a probe when the
surroundings are $n$ baryons, with $n=1,2,3$, or a quark gas, all at a
constant baryon density. The probe will be a quark, a diquark, or a baryon. As
for the diquark, we only consider the most attractive diquark where the flavor
and spin are both antisymmetric. The color of the quark and diquarks are in
the triplet or antitriplet states, respectively, for which the free energy
will be infinity in the confining phase where the surroundings are baryons.
Hence, we also consider the color singlet three diquark state in flavor SU(3),
whose quantum number is the same as the H dibaryonJaffe:1976yi .

Flavor, color and spin states of a multiquark system: Let us consider the
multiquark system consisting of $n$ baryons with a probe. In this work, we
assume that the orbital part of the wave function is totally symmetric and the
surrounding baryons are in flavor octet states. Since a baryon is a color
singlet, the color state of $3n+n_{p}$ quarks, where $n_{p}$ is the number of
quark in the probe, should be in the color state of the probe. For example, if
the probe is a baryon and $n=1$, the flavor-spin coupling state of six quarks
should be the conjugate of the color singlet state to satisfy the Pauli
exclusion principle.

Color :  , Flavor $\otimes$ Spin :
Then, we can decompose it into the flavor and the spin state as follows.

$\displaystyle[3,3]_{FS}=$
$\displaystyle[6]_{F}\otimes[3,3]_{S}+[5,1]_{F}\otimes[4,2]_{S}+[4,2]_{F}\otimes[5,1]_{S}$
$\displaystyle+[4,2]_{F}\otimes[3,3]_{S}+[4,1,1]_{F}\otimes[4,2]_{S}$
$\displaystyle+[3,3]_{F}\otimes[6]_{S}+[3,3]_{F}\otimes[4,2]_{S}$
$\displaystyle+[3,2,1]_{F}\otimes[5,1]_{F}+[3,2,1]_{F}\otimes[4,2]_{S}$
$\displaystyle+[2,2,2]_{F}\otimes[3,3]_{S}.$ (5)

Independently, we can determine the possible flavor states of multiquark
system using the outer product.

Flavor states of 2-baryon :

$\displaystyle\mathbf{8}\times\mathbf{8}=\mathbf{1}+\mathbf{8}_{(m=2)}+\mathbf{10}+\overline{\mathbf{10}}+\mathbf{27}.$
(6)

Combining Eq. (5) and Eq. (6), we can determine the possible flavor and spin
states of six quarks. The allowed flavor and spin states for all cases that we
consider are summarized in Table 1.

1b+1b | $\mathbf{1}(S=0),\mathbf{8}(S=1),\mathbf{10}(S=1),\overline{\mathbf{10}}(S=1),\mathbf{27}(S=0)$
---|---
2b+1b | $\mathbf{1}(S=\frac{3}{2}),\mathbf{8}(S=\frac{1}{2},\frac{3}{2}),\mathbf{10}(S=\frac{3}{2}),\overline{\mathbf{10}}(S=\frac{3}{2}),$
$\mathbf{27}(S=\frac{1}{2},\frac{3}{2}),\mathbf{35}(S=\frac{1}{2}),\overline{\mathbf{35}}(S=\frac{1}{2}),\mathbf{64}(S=\frac{3}{2})$
3b+1b | $\mathbf{1}(S=0),\mathbf{8}(S=1,2),\mathbf{10}(S=1),\overline{\mathbf{10}}(S=1),$
$\mathbf{27}(S=0,2),\overline{\mathbf{35}}(S=1),\overline{\mathbf{28}}(S=0)$
1b+1q | $\mathbf{3}(S=0,1),\overline{\mathbf{6}}(S=0,1),\mathbf{15}(S=0,1)$
2b+1q | $\mathbf{3}(S=\frac{1}{2},\frac{3}{2}),\overline{\mathbf{6}}(S=\frac{1}{2},\frac{3}{2}),\mathbf{15}(S=\frac{1}{2},\frac{3}{2}),$
$\mathbf{15^{\prime}}(S=\frac{1}{2},\frac{3}{2}),\mathbf{24}(S=\frac{1}{2},\frac{3}{2}),\mathbf{42}(S=\frac{1}{2},\frac{3}{2})$
3b+1q | $\mathbf{3}(S=0,1,2),\overline{\mathbf{6}}(S=0,1,2),\mathbf{15}(S=0,1,2),$
$\mathbf{15^{\prime}}(S=1),\overline{\mathbf{21}}(S=0),\overline{\mathbf{24}}(S=0,1,2),$
$\mathbf{42}(S=0,1,2),\overline{\mathbf{60}}(S=1)$
1b+1d | $\overline{\mathbf{3}}(S=\frac{1}{2}),\mathbf{6}(S=\frac{1}{2}),\overline{\mathbf{15}}(S=\frac{1}{2})$
2b+1d | $\overline{\mathbf{3}}(S=0,1),\mathbf{6}(S=1),\overline{\mathbf{15}}(S=0,1)$
$\overline{\mathbf{15^{\prime}}}(S=1),\overline{\mathbf{24}}(S=0,1),\overline{\mathbf{42}}(S=0,1)$
3b+1d | $\overline{\mathbf{3}}(S=\frac{1}{2},\frac{3}{2}),\mathbf{6}(S=\frac{1}{2},\frac{3}{2}),\overline{\mathbf{15}}(S=\frac{1}{2},\frac{3}{2}),$
$\overline{\mathbf{15^{\prime}}}(S=\frac{1}{2},\frac{3}{2}),\overline{\mathbf{24}}(S=\frac{1}{2},\frac{3}{2}),\overline{\mathbf{42}}(S=0,1),$
$\overline{\mathbf{48}}(S=\frac{1}{2})$
1b+3d | $\mathbf{8}(S=\frac{1}{2})$
2b+3d | $\mathbf{1}(S=0),\mathbf{8}(S=1),\overline{\mathbf{10}}(S=1),\mathbf{10}(S=1),\mathbf{27}(S=0)$
3b+3d | $\mathbf{8}(S=\frac{1}{2}),\overline{\mathbf{10}}(S=\frac{3}{2})$
Table 1: List of possible flavor and spin states in flavor SU(3) symmetry. b,
q and d represent a baryon, a quark and a diquark, respectively.

Now, we investigate the relative magnitude of the interaction which a quark
inside the probe sees from the surrounding $n$ baryons using the following
formula. A similar formula was used previously by us to compare the nucleon
repulsion with that of the lattice calculationPark:2019bsz .

$\displaystyle\Delta H_{CS}^{n\mathrm{b}+p}$
$\displaystyle=H_{CS}^{n\mathrm{b}+p}-H_{CS}^{n\mathrm{b}}-H_{CS}^{p},$ (7)
$\displaystyle\Delta H_{CS}^{\mathrm{avg}}$
$\displaystyle=\frac{1}{n_{p}n\sum_{C,F,S}d_{CFS}}\sum_{C,F,S}d_{CFS}\Delta
H_{CS}^{n\mathrm{b}+p}.$ (8) $\displaystyle d_{CFS}$
$\displaystyle=d_{C}d_{F}d_{S}m_{FS}.$ (9)

Here, $n\mathrm{b}$ and $p$ in the superscripts represent $n$ external baryons
and the probe, respectively. The probe will be a baryon, a quark, a diquark or
three correlated diquarks. $n_{p}$ is the number of quarks in the probe. We
will investigate cases with $n=1,2,3$, and also consider the case where
$n\mathrm{b}$ is replaced by a single quark so as to study the deconfined
phase. $d_{C},d_{F}$ and $d_{S}$ are the dimensions of the color, flavor and
spin states of $3n+n_{p}$ quarks, respectively, $m_{FS}$ is the multiplicity
of the flavor and spin states, and the summation is taken for all possible
states. Here, we divide it by $n_{p}$ to normalize the result with respect to
the single quark case. We also divide by $n$ to keep the surrounding baryon at
constant density for comparison at the same density.

Free quark gas: Finally, we consider the case where the surrounding is a free
quark gas. In such a case, we assume that the surrounding free quarks are not
correlated with each other, but are correlated with the interacting object to
satisfy the Pauli principle. Therefore, we only need to consider the average
value of the color-spin interactions for all possible diquark configurations.
There are four diquark states satisfying the Pauli principle. We represent it
for the color SU($N_{C}$) in the Table 2. If we compare it with the results
for baryons and a quark, then we should multiply it by 3 to ensure comparison
at the same density.

| $q_{i}q_{j}$
---|---
Flavor | $A$ | $S$ | $A$ | $S$
Color | $A(\overline{3})$ | $A(\overline{3})$ | $S(6)$ | $S(6)$
Spin | $A(1)$ | $S(3)$ | $S(3)$ | $A(1)$
$-\lambda_{i}\lambda_{j}\sigma_{i}\cdot\sigma_{j}$ | $-6-\frac{6}{N_{C}}$ | $2+\frac{2}{N_{C}}$ | $-2+\frac{2}{N_{C}}$ | $6-\frac{6}{N_{C}}$
$\lambda_{i}\lambda_{j}$ | $-2-\frac{2}{N_{C}}$ | $-2-\frac{2}{N_{C}}$ | $2-\frac{2}{N_{C}}$ | $2-\frac{2}{N_{C}}$
Table 2: Classification of two quark interaction due to the Pauli exclusion
principle. We denote the antisymmetric and symmetric state as $A$ and $S$,
respectively. The symbols inside the parenthesis represent the multiplet
state.

For a diquark and a free quark, since
$\mathbf{3}\times\overline{\mathbf{3}}=\mathbf{1}+\mathbf{8}$, there are two
possible color states of three quarks, which will come with the flavor-spin
configuration as below.

1\. Color :  , Flavor $\otimes$ Spin :

2\. Color :  , Flavor $\otimes$ Spin :
Considering both cases, we can calculate the average value of the color-spin
interaction.

For the color-spin interaction between three correlated diquarks and a free
quark, the color and flavor $\otimes$ spin coupling state is as follows.

Color :  , Flavor $\otimes$ Spin :

Flavor states of 7 quarks : $\mathbf{1}\times\mathbf{3}=\mathbf{3}$

Flavor and spin :  $\mathbf{3}(S=\frac{1}{2})$ .

The $H_{CS}$ factors for this 7 quarks state and the three correlated diquarks
are -12 and -24, respectively. In order to compare the result with the
interaction factor when one baryon looks at one quark, we need to divide by 2,
which corresponds to taking $n_{p}=6,n=1/3$ in Eq. (2).

For the color-spin interaction between a baryon and a quark, we can use the
result in the previous discussions.

Color-color interaction: The spin independent color-color type of interaction
is typically responsible for the confining and Coulomb type of interactions.
While such interactions are important at large separation between color states
within a color singlet configuration, the interaction of a colored object with
a color singlet configuration is small. In fact, it is zero if the color
singlet configuration has no color polarizations: that is all quarks have the
same spatial distribution. This is so because the color of the immersed object
will be the same as the color of the multiquark configuration it makes with
the surrounding color singlet baryons. Consider the following color-color
interaction factor composed of $n_{p}$ quarks from the probe and $3n$ quarks
from the $n$ baryons.

$\displaystyle\sum_{i<j}^{n_{p}+3n}\lambda_{i}^{c}\lambda_{j}^{c}$
$\displaystyle=\frac{1}{2}\bigg{[}\big{(}\sum_{i}^{n_{p}+3n}\lambda_{i}^{c}\big{)}^{2}-\sum_{i}^{n_{p}}(\lambda_{i}^{c})^{2}-\sum_{i}^{3n}(\lambda_{i}^{c})^{2}\bigg{]}$
$\displaystyle=\frac{1}{2}\bigg{[}\big{(}\sum_{i}^{n_{p}}\lambda_{i}^{c}\big{)}^{2}-\sum_{i}^{n_{p}}(\lambda_{i}^{c})^{2}\bigg{]}-\frac{1}{2}\bigg{[}\sum_{i}^{3n}(\lambda_{i}^{c})^{2}\bigg{]}$
$\displaystyle=\sum_{i<j}^{n_{p}}\lambda_{i}^{c}\lambda_{j}^{c}+\sum_{i<j}^{3n}\lambda_{i}^{c}\lambda_{j}^{c}.$
(10)

In the second line, the first square bracket is the interaction within the
probe while the second that between the color singlet baryons. That is, the
contributions between the quarks in the probe and those in the color singlet
configuration cancel out. The result is valid even for color SU($N_{C}$) with
any $N_{c}$.

Results:

SU(2)F | 1b | 2b | 3b | Free quarks
---|---|---|---|---
quark | 8 | 8.533 | 6.133 | 4.364
diquark | 8 | 8 | 8 | 8
baryon(octet) | 7.111 | 7.111 | 7.111 | 8

SU(3)F | 1b | 2b | 3b | Free quarks
---|---|---|---|---
quark | 6 | 6.446 | 4.644 | 2.823
diquark | 6 | 6.176 | 5.551 | 6.3
three correlated diquarks | 6 | 6 | 6 | 6
baryon(octet) | 5.714 | 5.78 | 4.944 | 6

Table 3: $\Delta H_{CS}^{\mathrm{avg}}$ for different probes (column) in
various surroundings (row) The upper and lower tables are for flavor SU(2) and
SU(3), respectively.

As the color-color interactions cancel each other out, to understand the most
stable structure at high density where the baryons start to overlap, we
consider the color-spin interaction energy of a quark, a diquark or a baryon
with the surrounding baryons. It should be noted that if the immersed object
is an isolated colored object, it will require large energy to bring the
isolated color charge to the position: this is why the thermal Wilson line is
infinity in the confining phase. Therefore, we also consider the color singlet
correlated three diquark state in flavor SU(3) case to search for the minimum
energy configuration in the confining phase. However, when the quark density
is large and color confining effects disappear or one is in the quark gas
phase, the color of the immersed object can be neglected.

First, we represent $\Delta H_{CS}^{\mathrm{avg}}$ values between the probe
and the surrounding for all cases in the Table 3. As can be seen in the table,
when the surrounding is a single baryon, the baryon has the smallest
interaction energy with the surrounding for both the flavor SU(2) and SU(3)
cases. From the values in the Table 3, it can be seen that in the case of
flavor SU(3), the magnitude of the color-spin interaction of baryon is smaller
than that of others by approximately 5 MeV per single quark if the density is
similar to that inside a nucleon so that the previously estimated overall
factor is used. However, when the surrounding becomes three baryon states or
the free quark gas, a quark has the lowest interaction energy. In a series of
workJeong:2019lhv ; Sen:2020peq ; Duarte:2020xsp ; Duarte:2020kvi , K.S.J. et
al. has shown that a phenomenological model involving quarks and baryons that
leads to quarkyonic EOS can be constructed when the short range baryon-baryon
repulsion, through the excluded volume effect, are introduced. More recently,
we have shown that in a realistic quark model, the short range quark-baryon
repulsion is smaller than the baryon-baryon repulsion in the lowest energy
channel, and that such difference leads to quarkyoniclike EOS. Here, we have
shown that the quark-baryon repulsion is in general smaller than the baryon-
baryon repulsion at high baryon density or free quark gas. Our results show
that the probe that feels the lowest color-spin interaction that dominates the
repulsion at short distance is quark rather than baryon, which provides a
theoretical justification for quarkyoniclike phase.

Second, to determine which probe has lower energy in dense matter, one has to
consider the internal color-spin factor of the immersed object. The octet
baryon and the diquark both have color-spin factors of -8. The additional
kinetic energy within the probe typically cancels the color-color interaction
strength rendering the color-spin interaction to be the only relevant
interaction strength within the probePark:2018wjk ; Karliner:2017qjm . Table 4
shows the results after considering these internal color-spin factors. It can
be seen that the lowest energy configuration in both the baryon matter and
quark matter is the diquark for flavor SU(2) and the color singlet three
correlated diquarks for flavor SU(3). Therefore, considering confinement, the
color singlet three correlated diquarks is the most stable configuration in
dense baryonic matter, when the density becomes similar to the quark density
within the baryon so that only the color-spin interaction becomes relevant.
Such configurations will appear in dense baryonic matter. Hence, while H
dibaryon do not seem to exist in the vacuumPark:2016cmg , similar
configurations will appear in dense baryonic matter. When deconfinement takes
place, considering the additional strange quark mass in the color singlet
three diquarks, flavor SU(2) diquarks will also appear contributing as
important configuration in dense matter: hence, a diquark matter.
Additionally, if one can show that the interaction between correlated diquarks
is relatively more attractive than the corresponding values involving baryons
or quarks, then the so-called diquarkyonic matter consisting of diquarks and
baryons or diquark condensation may emerge.

SU(2)F | 1b | 2b | 3b | Free quarks
---|---|---|---|---
quark | 8 | 8.533 | 6.133 | 4.364
diquark | 4 | 4 | 4 | 4
baryon(octet) | 4.444 | 4.444 | 4.444 | 5.333

SU(3)F | 1b | 2b | 3b | Free quarks
---|---|---|---|---
quark | 6 | 6.446 | 4.644 | 2.823
diquark | 2 | 2.176 | 1.551 | 2.3
three correlated diquarks | 2 | 2 | 2 | 2
baryon(octet) | 3.048 | 3.113 | 2.277 | 3.333

Table 4: Same as Table III after internal interactions within the probes are
considered.

## Acknowledgments

This work was supported by Samsung Science and Technology Foundation under
Project Number SSTF-BA1901-04. The work of A.P. was supported by the Korea
National Research Foundation under the grant number 2021R1I1A1A01043019. The
authors would like to thank Kiesang Jeong and Hyungjoo Kim for useful
discussions.

## References

  * (1) T. Kojo, G. Baym and T. Hatsuda, [arXiv:2111.11919 [astro-ph.HE]].
  * (2) L. McLerran and R. D. Pisarski, Nucl. Phys. A 796, 83 (2007)
  * (3) L. McLerran and S. Reddy, Phys. Rev. Lett. 122, no. 12, 122701 (2019)
  * (4) K. S. Jeong, L. McLerran and S. Sen, Phys. Rev. C 101, no.3, 035201 (2020)
  * (5) S. Sen and N. C. Warrington, Nucl. Phys. A 1006, 122059 (2021)
  * (6) D. C. Duarte, S. Hernandez-Ortiz and K. S. Jeong, Phys. Rev. C 102, no.2, 025203 (2020)
  * (7) D. C. Duarte, S. Hernandez-Ortiz and K. S. Jeong, Phys. Rev. C 102, no.6, 065202 (2020)
  * (8) S. Sen and L. Sivertsen, [arXiv:2011.04681 [astro-ph.HE]].
  * (9) T. Zhao and J. M. Lattimer, Phys. Rev. D 102, no.2, 023021 (2020)
  * (10) J. Margueron, H. Hansen, P. Proust and G. Chanfray, [arXiv:2103.10209 [nucl-th]].
  * (11) A. Park, K. S. Jeong and S. H. Lee, Phys. Rev. D 104, no.9, 094024 (2021) doi:10.1103/PhysRevD.104.094024 [arXiv:2105.12948 [hep-ph]].
  * (12) A. Park, S. H. Lee, T. Inoue and T. Hatsuda, Eur. Phys. J. A 56, no.3, 93 (2020)
  * (13) N. Ishii, S. Aoki and T. Hatsuda, Phys. Rev. Lett. 99, 022001 (2007)
  * (14) T. Inoue et al. [HAL QCD Collaboration], Prog. Theor. Phys. 124, 591 (2010)
  * (15) M. Oka, K. Shimizu and K. Yazaki, Prog. Theor. Phys. Suppl. 137, 1-20 (2000) doi:10.1143/PTPS.137.1
  * (16) S. H. Lee, S. Yasui, W. Liu and C. M. Ko, Eur. Phys. J. C 54, 259-265 (2008) doi:10.1140/epjc/s10052-007-0516-z [arXiv:0707.1747 [hep-ph]].
  * (17) A. Park, W. Park and S. H. Lee, Phys. Rev. D 96, no.3, 034029 (2017) doi:10.1103/PhysRevD.96.034029 [arXiv:1706.10025 [hep-ph]].
  * (18) S. Noh, W. Park and S. H. Lee, Phys. Rev. D 103, 114009 (2021) doi:10.1103/PhysRevD.103.114009 [arXiv:2102.09614 [hep-ph]].
  * (19) R. L. Jaffe, Phys. Rev. Lett. 38, 195 (1977).
  * (20) W. Park, S. Noh and S. H. Lee, Nucl. Phys. A 983, 1-19 (2019) doi:10.1016/j.nuclphysa.2018.12.019 [arXiv:1809.05257 [nucl-th]].
  * (21) M. Karliner and J. L. Rosner, Phys. Rev. Lett. 119, no.20, 202001 (2017) doi:10.1103/PhysRevLett.119.202001 [arXiv:1707.07666 [hep-ph]].
  * (22) W. Park, A. Park and S. H. Lee, Phys. Rev. D 93, no.7, 074007 (2016) doi:10.1103/PhysRevD.93.074007 [arXiv:1602.05017 [hep-ph]].


# Cusp singularities in the distribution of orientations of asymmetrically-
pivoted hard discs on a lattice

Sushant Saryal<EMAIL_ADDRESS>(corresponding author) Department of
Physics, Indian Institute of Science Education and Research, Pune 411008,
India  Department of Chemistry, Princeton University, Princeton, New Jersey
08544, USA  Deepak Dhar<EMAIL_ADDRESS>Department of Physics, Indian
Institute of Science Education and Research, Pune 411008, India

###### Abstract

We study a system of equal-sized circular discs each with an asymmetrically
placed pivot at a fixed distance from the center. The pivots are fixed at the
vertices of a regular triangular lattice. The discs can rotate freely about
the pivots, with the constraint that no discs can overlap with each other. Our
Monte Carlo simulations show that the one-point probability distribution of
orientations have multiple cusp-like singularities. We determine the exact
positions and qualitative behavior of these singularities. In addition to
these geometrical singularities, we also find that the system shows order-
disorder transitions, with a disordered phase at large lattice spacings, a
phase with spontaneously broken orientational lattice symmetry at small
lattice spacings, and an intervening Berezinskii-Kosterlitz-Thouless phase in
between.

##  I Introduction

In many molecular solids, the melting transition from the low-temperature
crystalline solid phase to the high-temperature liquid phase does not occur in
a single step, but one finds a multiplicity of mesophases. These are called
‘liquid crystals’, or ‘plastic solids’. In the former, the periodic three-
dimensional crystal structure is absent, but varying amount of orientational
order may be present. In the latter, the average positions of the centers of
masses of molecules do lie on a three-dimensional crystalline lattice, but
there is none or only partial orientational order. These were originally
called plastic solids, as they can be deformed easily using much less force,
compared to ‘hard’ crystals. Some examples of common materials that show
plastic solid phases are nitrogennitrogen , carbon tetrabromide powell ,
formylferrocene kaneko . The currently favored nomenclature for these is
orientationally disordered crystals. In recent years, these have attracted a
lot of interest, because of their promising applications in diverse areas such
as solid electrolytesbatteries , drug deliverydrug-delivery ,
optoelectronicsopto , barocaloricsrefrigeration , piezoelectricspiez , etc.
For a recent review of the applications, see Das el alDas .

Pauling in 1930 derived a rough criteria for strongly hindered rotational
motion of molecles in crystalline solids pauling . But it was Timmermans who
systematized the phenomenological study of plastic crystals starting from the
1930stimmermans . On the theoretical front, Pople and Karasz pople-karasz
extended the two-lattice model of Lennard-Jones and DevonshireLennard  to
account for the order-disorder transition in the orientation of molecular
crystals. A minimal model for these would be to assume the constituents as
rigid objects each identically pivoted on a lattice and free to rotate
provided no objects overlap with each other. Casey and Runnelsrunnels1  and
Freasiers and Runnels runnels2  examined a system of hard squares with centers
ﬁxed on the 1d lattice. We have recently discussed this model, and called it
rigid hard rotors on a lattice as model of multiple phases shown by plastic
crystals to describe the transitions between them dd1 ; dd2 ; dd3 . Note that
since the lattice is always present, the model does not have a ‘liquid’ phase
with no crystalline order.

In an earlier paper dd3 , we determined the exact functional form of the one-
point probability distribution function of orientations at a site for a range
of lattice spacings when a particular condition, called the at-most one
overlap (AOO) condition, holds. In this paper we particularly examine a system
of hard discs asymmetrically pivoted on a triangular lattice and study the
one-point probability distribution function (PDF) $P(\theta)$ of orientations
$\theta$ beyond the AOO condition. We find that the distribution function
$P(\theta)$ shows cusp singularities. We determine the position and
qualitative behaviour beyond the AOO condition exactly. Singularities in the
pair distribution function were studied earlier by Stillingerstill1  and
numerically observed in the the probability distribution of bond-pair angles
in a system of hard spheres still2 , but non-trivial singularities in the one-
point function have not been discussed before. We also numerically verify our
findings with the help of Monte Carlo simulations.

This paper is organized as follows. In section II we define our model. In
section III we show that there exist multiple cusp singularities in the one-
point probability distribution of orientations and exactly determine their
nature and positions. In section IV we verify our findings using Monte Carlo
simulations. Section V contains some concluding remarks.

##  II Model

Figure 1: System of hard plates pivoted asymmetrically on the triangular
lattice.

We consider a system of identical unit-radius circular discs each with an
asymmetrically placed pivot at a distance $\epsilon$ from the center, as shown
in Fig 1. The pivots form a regular triangular lattice with lattice spacing
$a$, and the discs can rotate freely about the pivots, with the constraint
that no discs can overlap with each other. Orientation of a disc pivoted at
lattice site $\vec{r}$ is specified by an angle $\theta(\vec{r})$, measured
between long-axis (passing through the pivot and center) and $x$-axis. From
elementary geometry it is clear that if $a>2(1+\epsilon)$, the discs rotate
freely, and the closest-packing limit is reached when $a=2$. We write lattice
spacing $a$ as,

$\displaystyle a=2+\epsilon x,$ (1)

with $0\leq x\leq 2$.
The partition function of the system having $N$ discs is

$\displaystyle\mathcal{Z}_{N}=\left[\prod_{\bf{r}}\int_{-\pi}^{\pi}\frac{d\theta(\bf{r})}{2\pi}\right]\prod_{\mathclap{\begin{subarray}{c}\mathbf{x},\mathbf{y}\\\
\mathbf{x}\neq\mathbf{y}\end{subarray}}}\left[1-\eta\big{(}\theta(\bf{x}),\theta(\bf{y})\big{)}\right],$
(2)

where $\eta\big{(}\theta(\mathbf{x}),\theta(\mathbf{y})\big{)}$ is an
indicator function which is one when discs at sites $\mathbf{x}$ and
$\mathbf{y}$ with orientations $\theta(\bf{x})$ and $\theta(\bf{y})$
respectively, overlap and zero otherwise. We define the entropy per site
$s(x,\epsilon)$ by

$\displaystyle
s(x,\epsilon)=\lim_{N\rightarrow\infty}\frac{\log\mathcal{Z}_{N}}{N}$ (3)

It is expected that the qualitative behaviour of $s(x,\epsilon)$ is same for
all values of $\epsilon$. In fact, as $\epsilon$ tends to zero, the function
$s(x,\epsilon)$ has a well-defined non-trivial limit

$\displaystyle s(x)=\lim_{\epsilon\rightarrow 0}s(x,\epsilon).$ (4)

In this limit, the no-overlap condition between two neighboring sites
simplifies. For two neighboring rotors, if the line joining the pivots makes
an angle $\phi$ with the $x$-axis, the no-overlap condition to first order in
$\epsilon$ becomes

$\displaystyle x-\cos(\theta-\phi)+\cos(\theta^{{}^{\prime}}-\phi)\geq 0.$ (5)

The limit $\epsilon\rightarrow 0$ has the advantage that the number of
parameters specifying the model is reduced to $1$, In the following, for the
sake of simplicity, we restrict our discussion to this case. The case of more
general $\epsilon$ presents no additional special features.

This model can also be thought of as a system of planar spins
$\\{\theta(\bf{r})\\}$ on the vertices of a triangular lattice, with nearest
neighbor interaction Hamiltonian $\mathcal{H}$ given by,

$\displaystyle\mathcal{H}=J\sum_{\bf{r}}\sum_{j=0}^{2}\Theta\left[\cos(\theta(\mathbf{r})-j\pi/3)-\cos(\theta_{j}^{{}^{\prime}}(\mathbf{r})-j\pi/3)-x\right]$

where $\theta_{j}^{{}^{\prime}}(\bf{r})$ is the neighboring spins of the spin
$\theta(\bf{r})$ in the lattice direction $j\pi/3$, and $\Theta(x)$ is the
Heavyside step function of $x$. The hard-core limit corresponds to setting $J$
to $+\infty$.

Figure 2: $\eta(\theta,\theta^{{}^{\prime}})$ for nearest-neighbor discs along
x-axis for different values of $x$ in
($\dfrac{\theta}{2\pi}$,$\dfrac{\theta^{{}^{\prime}}}{2\pi}$) plane. Shaded
area (red) corresponds to overlap region where $\eta$ is $1$. In the unshaded
area $\eta=0$. For other nearest-neighbor pairs,
$\eta(\theta,\theta^{{}^{\prime}})$ can be similarly obtained using the
symmetries of the triangular lattice.

This model is of the same form as the model of hard-core spins studied earlier
by Sommers et al sommers . These authors studied the case where our condition
Eq.(5) is replaced by $|\theta-\theta^{\prime}|\leq\alpha$. The qualitative
behavior of the models is similar. The main difference where our model differs
from theirs is the explicit breaking of isotropy in the spin space by the
lattice-direction dependent interaction. In particular, the cusp singularities
we discuss here are not present in the hard-core spin model.

##  III Probability distribution of orientations

Let $P(\theta)d\theta$ denote the probability that the disc pivoted to a
randomly picked site in the equilibrium ensemble described by the partition
function in Eq. (2) is at an orientation between $\theta$ and
$\theta+d\theta$.

The simplest case is when $x$ lies in the range $[1+\sqrt{3}/2,2]$. In this
case, it is easily seen that in any configuration of orientations of discs,
any disc can overlap with at-most one other disc. In was called the AOO
condition in our earlier paperdd3 . In this case, the calculation of partition
function simplifies to the calculation of the partition function of dimers and
vacancies on the same lattice. In this case, the activity of dimers is given
by

$\displaystyle
z=-\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\int_{-\pi}^{\pi}\frac{d\theta^{{}^{\prime}}}{2\pi}\eta(\theta,\theta^{{}^{\prime}}).$
(7)

Under the AOO condition, $P(\theta)$ is given by the expressiondd3

$\displaystyle
P(\theta)=(1-2\bar{n}(z))[1+\frac{\bar{n}(z)}{3z(1-2\bar{n}(z))}\sum_{i=1}^{6}f_{i}(\theta)]$
(8)

where $f_{i}(\theta)$ is given by

$\displaystyle
f_{i}(\theta)=\int_{-\pi}^{\pi}\frac{d\theta_{i}}{2\pi}\eta(\theta,\theta_{i})$
(9)

and $\bar{n}(z)$ is the dimer number density at activity $z$, which is given
by the low density expansion

$\displaystyle\bar{n}(z)=\sum_{n=1}^{\infty}(-1)^{n-1}a_{n}z^{n},$ (10)

where $a_{n}$ is the number of heaps made of dimers vien . For the triangular
lattice, we have

$\displaystyle\bar{n}(z)=3[z-11z+144z^{2}-2047z^{3}+30526z^{4}....\,]$ (11)

In our problem, the explicit expression for the function
$\eta(\theta,\theta_{1})$ is

$\displaystyle\eta(\theta,\theta_{1})=\begin{cases}1,&\text{if
}|\theta|\leq\arccos(x-1)\,\,\,\text{and}\,\,\,|\theta_{1}|\geq\arccos[\cos(\theta)-x].\\\
0,&\text{otherwise}.\end{cases}$ (12)

Then, it is easily seen that $f_{1}(\theta)$ is given by

$\displaystyle
f_{1}(\theta)=\begin{cases}2(\pi-\arccos[\cos(\theta)-x]),\,\,|\theta|\leq\arccos(x-1).\\\
0,\,\,\,\,\,\,\,\,\,\,\,\text{otherwise}.\end{cases}$ (13)

Other $f_{i}$’s can be easily found using the symmetries of the underlying
lattice. From this, it is easily seen that $P(\theta)$ has square-root cusp
singularities (see Appendix A for more details) at

$\displaystyle\theta_{\text{cusp}}=\frac{j\pi}{3}\,\pm\,\arccos(x-1)$ (14)

where $j=0,1,2,....5$. When $x=1+\sqrt{3}/2-\delta$, with $\delta$ positive,
but small, the AOO condition is no longer satisfied. In this case, one can
still define the graphical expansion of Eq.(2) in terms of configurations of
dimers, but now, the configurations where two or more dimers are incident on a
vertex have non-zero weight. However, if $\delta$ is small, the weights of
such vertices is small. This suggests that we can organize the terms of this
series in the following form

$\displaystyle\mathcal{Z}=\mathcal{Z}_{0}+\mathcal{Z}_{1}+\mathcal{Z}_{2}+...$
(15)

In this expansion, $\mathcal{Z}_{r}$ is the sum over terms having $r$ dimer
pairs such that the dimers in each pair have a common vertex. We may associate
an extra factor $y$ with each such overlapping dimer pair, and consider the
partition sum

$\displaystyle\mathcal{Z}(y)=\sum_{r=0}^{\infty}y^{r}\mathcal{Z}_{r}.$ (16)

We treat $y$ as a small parameter, and assuming the sum converges for small
enough $y$, treat it is as a perturbation series in $y$. If we put $y=0$ in
this series, we have a sum over the all dimer configurations. In this
ensemble, one can define the one-point function $P(\theta)$ as

$\displaystyle
P(\theta)=\sum_{\bf{x}}\langle\delta(\theta(\bf{x})-\theta)\rangle$ (17)

where the angular brackets denote the average over the ensemble. It is easy to
see that in the unperturbed ensemble $\mathcal{Z}_{0}$, Eq. (8) continues to
remain valid, but now in this ensemble there are ranges of $\theta$ where more
than one of the $f$-terms contributes in the equation. And the positions of
the cusp singularities in $P(\theta)$ are still given by Eq.(14).

Consider that the terms in $\mathcal{Z}_{1}$, involving two specified dimers
meeting at a specific site $\bf{x}$. Say the dimers are covering the bonds
$(\bf{x},\bf{y})$ and $(\bf{x},\bf{z})$. This weight can be written as a
product of two terms $T_{1}$ and $T_{2}$, with

$\displaystyle T_{1}=\int d\theta(\mathbf{x})\int d\theta(\mathbf{y})\int
d\theta(\mathbf{z})~{}~{}~{}\eta(\bf{x},\bf{y})\eta(\bf{x,}\bf{z})$ (18)

It can be shown that for small positive $\delta$, $T_{1}$ varies asdd3
$\delta^{3/2}$. $T_{2}$ is a polynomial in $z$, the sum over all possible
partial dimer coverings of the lattice, not involving sites
$\bf{x},\bf{y},\bf{z}$. Similar statement is valid for higher $r$.

If we expand the function $P(\theta)$ in powers of $y$, each term in the
perturbation sum is well-behaved, and only singularities in $\theta$ come from
integration of the functions $\eta$, and hence the positions are same as in
Eq.(14) and are robust. By analyic continuation on $y$, we expect these
results to hold for all $y$, and hence at $y=1$. Thus, we conjecture that the
cusp singularities $P(\theta)$ are given by Eq. (14) for all $x$, in the range
$0<x<2$. This is shown in Fig. 3.

Figure 3: The change in the position of cusp singularities in $P(\theta)$ as
we vary $x$. The positions vary linearly with the variable $\arccos(x-1).$

##  IV Monte Carlo Results

Now we present our findings of Monte Carlo simulations. Our simulations were
done on lattices of size varying from $50\times 50$ to $100\times 100$. We use
a single site update scheme: we pick a site $\vec{r}$ at random, and try to
change the value of $\theta_{\vec{r}}$ to $\theta_{\vec{r}}+\Delta\theta$,
where $\Delta\theta$ is a random variable with a uniform distribution from
$-\Delta_{0}$ to $+\Delta_{0}$. We accept the move if the new value does not
result in any overlap. And repeat. We average over times of order $7.2$
million MCS, after rejecting the first $8\times 10^{5}$ steps.

Figure 4: One-point PDF under AOO regime. (a) $x=1.92$ and (b)
$x=1+\frac{\sqrt{3}}{2}$. Dashed vertical lines represent the positions of
cusp singularities given in Eq.(14). Theoretical predictions, Eq.(8) and
Eq.(14), and Monte Carlo results are in very good agreement. Figure 5: One-
point PDF beyond AOO regime. (a) $x=1.7$, (b) $x=1.5$, (c) $x=1.25$, and (d)
$x=1$. Dashed vertical represent the positions of cusp singularities given in
Eq.(14).  Figure 6: Comparison of one-point PDF obtained using Monte Carlo
simulations and our ansatz given in Eq. (19). (a) $x=1.7$, $m=1.405$ and (b)
$x=1.5$, $m=1.32$.

In Fig. 4 one-point PDF is plotted using Monte Carlo simulations along with
our theoretical prediction, Eq. (8), in the AOO regime $1+\sqrt{3}/2\leq x\leq
2$. One can see that Monte Carlo findings are in very good agreement with our
theoretical prediction Eq. (8). Dashed vertical lines corresponds to the
position of the cusp singularities in Eq. (14) and it is also in excellent
agreement with the Monte Carlo simulations. Note that for $x=1+\sqrt{3}/2$,
cusp singularities merge in pairs, producing only 6 singular points.

In general, there are twelve singularities for each value of $x$, except at
special points where the singularities merge in pairs. We have determined the
cusp positions from the Monte Carlo data for several values of $x$, shown in
Fig. 5. These are in very good agreement with the predicted values.

We parameterize the distribution function for values of $x$ outside the AOO
regime using a fitting form with only one fitting parameter $m$, by

$\displaystyle
P_{\text{approx}}(\theta)=\mathcal{N}\,\prod_{i=0}^{5}[1+m~{}f(\theta-i\pi/3)],$
(19)

where $\mathcal{N}$ is the normalization constant. Note that in the AOO
regime, this expression is exact, and the parameter $m$ can be written in
terms of the density of dimers. Outside the AOO regime, the expression is only
approximate, but incorporates the known exact position of the cusp
singularities, and its analytical structure is suggested by the solution of
the model of interacting rods on the Bethe lattice dd2 . The plot shown in
FIG. 6 compares the Monte Carlo data beyond the AOO regime, for $x=0.7$ and
$0.5$, with value of $m$ chosen to provide the best fit. We see that the
ansatz provides an good qualitative description for the one-point PDF. The
deviations from this form occur only in the intervals of $\theta$ for which
the failure of the AOO condition is possible.

Our Monte Carlo simulations also revealed the presence of both the
Berezinskii–Kosterlitz–Thouless (BKT) phase and the orientational-symmetry-
broken phase in addition to the disordered phase. Clock-model also exhibits
similar kind of behaviour where there is a intervening BKT phase between the
phases with broken symmetry and the disordered phase nelson77 . A phase where
the symmetry of orientational distribution function is broken under lattice
rotations was observed for values of $x$ less than 0.5, as illustrated in FIG.
7, where the one-point probability density function is peaked at $\pi/6$ and
decays rapidly to zero as we move away from it.

Figure 7: One-point PDF obtained using Monte Carlo simulations for $x=0.45$
and (b) $x=0.3$. Figure 8: (a)-(c) panels show spatial heat map of
orientations of discs in a 600$\times$600 triangular lattice for different
$x$. Color coding is from -1 to 1 as orientations are divided by $\pi$. (d)
log-log plot of mean-square magnetization vs $L$ for different $x$.

In the intermediate region, roughly in the range $0.5<x<1$, we observe BKT
phase which is usually characterized by a power law decay of correlation
function with the exponent that depends on $x$ and no long-range order. In
FIGS. 8 (a)-(c) spatial heat maps of orientations of discs are plotted for
different values of $x$ chosen from three different phases.

BKT phase can also identified by looking into the system size dependence of
mean-square magnetization. On a $L\times L$ lattice, leading order dependence
of mean-square magnetization on $L$ for different phases varies as,

$\displaystyle\langle M_{x}^{2}\rangle+\langle
M_{y}^{2}\rangle\propto\begin{cases}L^{2}&\text{in the disordered phase }\\\
L^{4}&\text{in the ordered phase }\\\ L^{4-\eta(x)}&\text{in the BKT
phase}.\end{cases}$ (20)

where $M_{x}=\sum_{\mathbf{r}}\cos\theta(\mathbf{r}),$
$M_{y}=\sum_{\mathbf{r}}\sin\theta(\mathbf{r}),$ and $0\leq\eta(x)\leq 2$.
Angular bracket represents ensemble average. In FIG. 8(d) mean-square
magnetization vs $L$ is plotted for various values of $x$. One can see that
the system size scaling of mean-square magnetization shows that there is a
intervening BKT phase between ordered and disordered phases.

##  V Concluding remarks

The model discussed in this paper is of interest for several reasons. Firstly,
it is a model simple to define, which shows a large number of phases, and
phase transitions, just by varying one parameter. It thus provides a minimal
model for describing the multitude of phases seen in plastic solids.

Secondly, we have argued that to characterize the many phases, it is
convenient to use not a single order parameter, but the whole distribution
function of orientations. The full distribution as an order parameter has been
discussed in the context of spin-glasses parisi-mezard . Also, we have shown
that this distribution function has robust geometrical singularities, whose
qualitative behavior is easy to determine theoretically. To the best of our
knowledge, this is first system with continuous degrees of freedom that shows
non-trivial singularities in the one-point function, whose position changes
when the coupling constant is varied.

Thirdly, we could determine the angular dependence of the distribution
function $P(\theta)$ for a range of values on $x$, when the AOO condition is
satisfied. Outside this range, we showed that the distribution function can be
expanded in a perturbation series in a variable $y$. In general, problems
where the position of singularity varies with the perturbation strength are
difficult to construct. For example, the function $P(\theta)$ is expected to
have cusp singularity of the form,

$\displaystyle P(\theta)=A(\theta,x)+B(\theta,x)|\theta-\theta_{c}(x)|^{1/2},$
(21)

near a cusp singularity $\theta=\theta_{c}(x)$, where $A(\theta,x)$ and
$B(\theta,x)$ are smooth functions of $x$ (different on different sides of the
cusp). A naive perturbation series in $\delta=x-x^{*}$ for $P(\theta)$ about a
point $x^{*}$ would generate spurious singularities of the type
$|(\theta-\theta_{c}(x^{*})|^{-1/2}\delta$. Our perturbation parameter $y$
avoids this problem, as the $x$ dependent cusp singularity structure is built
in the perturbation series.

Fourthly, the connection of this problem to the hard-disc problem is also of
interest. Typical configurations generated in this model are visually not
easily distinguished from the configurations of the hard-disc model at same
density. In the limit $x$ tending to zero, we get the close-packed crystalline
solid. In our model, the centers of the discs cannot move freely, but each
center is restricted to a circle of radius $\epsilon$. The restricted model
seems to have qualitative behavior similar to the original model. For small
$x$, we do not have crystalline order, and if we construct a local coarse-
grained variable giving the average orientation of the lattice locally, this
will change slowly in space, as in the hard-disc model. In particular, for
$0<x<0.5$, one can show that there are no vortices possible. For intermediate
range of $x$, we have seen that there is Kosterlitz-Thouless phase, with
power-law decay of angular correlations.

In addition to these, the system shows a series of ordering transitions. In
fact, for low $x$, our Monte Carlo simulations also show a transition to a
glassy phase with very large correlation times. We have not discussed these
here. This seems to be a good direction for further studies.

## ACKNOWLEDGMENTS

S.S. acknowledges financial support from the Department of Physics, IISER
Pune. DD’s work is supported by the Indian government under the Senior
Scientist Scheme of the National Academy of Sciences of India. We thank
National Supercomputing Mission (NSM) for providing computing resources of
‘PARAM Brahma at IISER Pune’, which is implemented by C-DAC and supported by
the Ministry of Electronics and Information Technology (MeitY) and Department
of Science and Technology (DST), Government of India.

## Appendix A: Singularities in $P(\theta)$

We will show that the $f_{1}(\theta)$ has square-root cusp singularities at
$\theta^{\pm}_{\text{cusp}}=\pm\cos^{-1}(x-1)$. If
$\theta=\theta^{+}_{\text{cusp}}-\xi$ where $0<\xi\ll 1$, then we have,

$\displaystyle\cos\theta-x=-1+x\sqrt{2-x}\,\,\xi+O(\xi^{2}).$ (A1)

Let us define,

$\displaystyle y=\cos^{-1}[\cos{\theta}-x].$ (A2)

This can also be written as $\cos
y=\cos{\theta}-x=-1+x\sqrt{2-x}\,\,\xi+O(\xi^{2})$. If $\xi\to 0$ then
$y\to\pi$. Therefore we get,

$\displaystyle y=\pi-\sqrt{2x}\,\,(2-x)^{1/4}\,\sqrt{\xi}+\text{higher order
terms in}\,\,\,\xi$

Putting this in Eq. (13) we finally get,

$\displaystyle
f_{1}(\theta)=\sqrt{2x}\,\,(2-x)^{1/4}\sqrt{\theta^{+}_{\text{cusp}}-\theta}\,\,\,\,\,\,\,\theta\to\theta^{+}_{\text{cusp}}$
(A4)

Following same lines of argument one can show that cusp singularity also
exists at $\theta=\theta^{-}_{\text{cusp}}=-\cos^{-1}(x-1)$. Other
$f_{i}(\theta)$ can also be found from $f_{1}(\theta)$ using the symmetries of
triangular lattice, thus $f_{i}(\theta)=f_{1}\left(\theta-(i-1)\pi/3\right)$.
Similarly, two cusp singularities exist for each $f_{i}(\theta)$ whose
positions can be found using the symmetries of the triangular lattice. As a
consequence, there exist twelve cusp singularities in $P(\theta)$ whose
positions are given by,

$\displaystyle\theta_{\text{cusp}}=\frac{j\pi}{3}\,\pm\,\arccos(x-1)\,\,\,\text{where}\,\,\,j=0,1,2,....,5.$

## References

  * (1) L. Tassini, F. Gorelli, and L. Ulivi, High temperature structures and orientational disorder in compressed solid nitrogen, J. Chem. Phys. 122, 074701 (2005).
  * (2) B. M. Powell; V. F. Sears; G. Dolling, Neutron scattering from orientationally disordered solids, AIP Conference Proceedings 89, 221–229 (1982).
  * (3) Y. Kaneko and M. Sorai, Heat capacity and phase transitions of the plastic crystal formylferrocene, Phase Transitions, Vol. 80, Nos 6–7, June–July 2007, 517–528.
  * (4) J. M. Pringle, P. C. Howlett , D. R. MacFarlane, and M. Forsyth, Organic ionic plastic crystals: recent advances, J. Mater. Chem., 2010, 20, 2056-2062.
  * (5) E. Shalaev, K. Wu, S. Shamblin, J.F. Krzyzaniak, M. Descamps Crystalline mesophases: structure, mobility, and pharmaceutical properties Adv. Drug Deliv. Rev., 100 (2016), pp. 194-211.
  * (6) Z. Sun, T. Chen, X. Liu, M. Hong and J. Luo, Plastic Transition to Switch Nonlinear Optical Properties Showing the Record High Contrast in a Single-Component Molecular Crystal, J. Am. Chem. Soc., 2015, 137, 15660–15663.
  * (7) Li, B., Kawakita, Y., Ohira-Kawamura, S. et al. Colossal barocaloric effects in plastic crystals. Nature 567, 506–510 (2019).
  * (8) J. Harada, Y. Kawamura, Y. Takahashi, Y. Uemura, T. Hasegawa, H. Taniguchi and K. Maruyama, Plastic/Ferroelectric Crystals with Easily Switchable Polarization: Low-Voltage Operation, Unprecedentedly High Pyroelectric Performance, and Large Piezoelectric Effect in Polycrystalline Forms, J. Am. Chem. Soc., 2019, 141, 9349–9357.
  * (9) S. Das, A. Mondal, and C.L. Reddy, Chem. Soc. Rev., $\mathbf{49}$, 8878 (2020).
  * (10) L. Pauling, The Rotational Motion of Molecules in Crystals, Phys. Rev. $\mathbf{36}$, 430 (1930).
  * (11) J. Timmermans, Plastic crystals: a historical review, J. Phys. Chem. Solids, 18 (1961), pp. 1-8.
  * (12) J. A. Pople and F. E. Karasz, A Theory of fusion of molecular crystals: I. The effect of hindered rotation, J. Phys. Chem. Solids, 18, (1961) pp 28-39.
  * (13) J.E Lennard-Jones and A.F. Devonshire (1939), Critical and co-operative phenomena IV. A theory of disorder in solids and liquids and the process of melting. Proc. R. Soc. Lond. A170464–484.
  * (14) L. M. Casey and L. K. Runnels, Model for Correlated Molecular Rotation, J. Chem. Phys. $\mathbf{51}$, 5070 (1969).
  * (15) B. C. Freasier and L. K. Runnels , Classical rotators on a linear lattice, J. Chem. Phys. $\mathbf{58}$, 2963 (1973).
  * (16) S. Saryal, J. U. Klamser, T. Sadhu and D. Dhar, Multiple Singularities of the Equilibrium Free Energy in a One-Dimensional Model of Soft Rods, Phys. Rev. Lett. $\mathbf{121}$, 240601(2018).
  * (17) J. U. Klamser, T. Sadhu, and D. Dhar, Sequence of phase transitions in a model of interacting rods, Phys. Rev. E $\mathbf{106}$, L052101.
  * (18) S. Saryal and D. Dhar, Exact results for interacting hard rigid rotors on a d-dimensional lattice, J. Stat. Mech. (2022) 043204.
  * (19) F. H. Stillinger, Pair distribution in the classical rigid disk and sphere systems, Journal of Computational Physics, $\mathbf{7}$, 367 (1971).
  * (20) A. Donev, S. Torquato, and F. H. Stillinger, Pair correlation function characteristics of nearly jammed disordered and ordered hard-sphere packings, Phys. Rev. E $\mathbf{71}$, 011105 (2005).
  * (21) G. M. Sommers, B. Placke , R. Moessner, and S. L. Sondhi, From hard spheres to hard-core spins, Phys. Rev. B 103 104407 (2021).
  * (22) G. X. Viennot, in Proc. of the Colloque de Combinatoire Enumerative, Lecture Notes in Mathematics Vol. 1234, pp 321-350 (Springer-Verlag, Berlin, 1986).
  * (23) J. V. José , L. P. Kadanoff, S. Kirkpatrick, and D.R. Nelson, Renormalization, vortices, and symmetry-breaking perturbations in the two-dimensional planar mode, Phys. Rev. B $\mathbf{16}$, 1217(1977).
  * (24) M. Mezard, G. Parisi and M. Virasoro, Spin Glass Theory and Beyond, World Scientific (1986).