# Critical properties of various sizes of cluster
in the Ising percolation transition

Lizhu Chen  School of Physics and Optoelectronic Engineering,
Nanjing University of Information Science and Technology, Nanjing 210044,
China
<EMAIL_ADDRESS>YeYin Zhao  Key Laboratory of Quark and Lepton Physics
(MOE) and Institute of Particle Physics,
Central China Normal University, Wuhan 430079, China  Xiaobing Li  Key
Laboratory of Quark and Lepton Physics (MOE) and Institute of Particle
Physics,
Central China Normal University, Wuhan 430079, China  Zhiming Li  Key
Laboratory of Quark and Lepton Physics (MOE) and Institute of Particle
Physics,
Central China Normal University, Wuhan 430079, China  Yuanfang Wu  Key
Laboratory of Quark and Lepton Physics (MOE) and Institute of Particle
Physics,
Central China Normal University, Wuhan 430079, China

###### Abstract

It is proposed that the $O(n)$ spin and geometrical percolation models can
help to study the QCD phase diagram due to the universality properties of the
phase transition. In this paper, correlations and fluctuations of various
sizes of cluster in the Ising model are systematically studied. With a finite
size system, we demonstrate how to use the finite size scaling and fixed point
behavior to search for critical point. At critical point, the independency of
system size is found from skewness and kurtosis of the maximum, second and
third largest cluster and their correlations. It is similar to the Binder-
ratio, which has provided a remarkable identification of the critical point.
Through an explanation of the universal characteristic of skewness and
kurtosis of the order parameter, a possible application to the relativistic
heavy-ion collisions is also discussed.

###### keywords:

Percolation transition; various sizes of cluster; fixed point; QCD phase
transition.

PACS numbers: 25.75.Gz, 25.75.Nq, 24.10.Pa

##  1 Introduction

In high-energy nuclear collisions, we study properties of the excited nuclear
matter with QCD degrees of freedom and search for the signal of the QCD phase
transition. Lattice QCD has shown that it occurs a smooth crossover from
hadronic phase to the QGP phase at $\mu_{B}=0$ MeV [1]. It is expected that
the phase transition turns to be a first order at high baryon density and low
temperature [2, 3]. The endpoint of the first order phase transition to the
crossover is referred to as the critical point, which has been studied by
high-order cumulants of conserved quantities at RHIC/STAR [4, 5, 6].

Theoretically, Lattice QCD is of course the first principle to study the QCD
phase transition [1, 7, 8, 9, 10]. However, simulations for much larger
$\mu_{B}$ have sign problem [11]. Currently, the $O(1)$ (Ising model), $O(2)$
and $O(4)$ spin models are widely applied to study the QCD phase transition
due to the universality of the critical phenomena [12], although the system is
totally different. It has been pointed out that the QCD critical point belongs
to the same universality class as 3 dimensional (3D) Ising model. Using a
parametric representation or Monte Carlo simulations, the critical behavior of
the high-order cumulants [13, 14], effects of the critical slowing down [15,
16], and the finite size scaling [17, 18] are extracted and applied to study
the QCD critical point. Especially, characteristics of the non-monotonic and
sign changes in high-order cumulants are taken as the most important signals
for probing the QCD critical point in experiment [19].

Besides the $O(n)$ spin models, the 2 dimensional (2D) geometrical percolation
models can also help us to understand the QGP formation [20, 21, 22, 23, 24,
25]. With growing energy and size of the colliding nuclei, the number of color
strings grows and starts to overlap to form different sizes of the cluster on
the transverse nucleon interaction plane, very much similar to the disks in
the 2D percolation theory. Currently, the 2D percolation approach within Color
String Percolation Model (CSPM) had already given a successfully description
of most of the experimental data in the soft region. The CSPM can also
describe the QCD to hadron cross-over transition through the temperature
dependence of the shear viscosity over entropy density ratio ($\eta/s$) [24].

It should be noted, however, the critical properties in CSPM are derived from
Ref. [percolation-Stauffer], which are related to the percolation transition
with infinite size of the system. In heavy-ion collisions, the effect of the
limited size of the formed matter can not be neglected. Consequently, we
should systematically study the 2D percolation transition with finite size of
the system. By a suitable definition of the cluster in the 2D Ising model, it
has been found that the thermal magnetization transition can be directly
mapped into a geometrical percolation transition [27].

The typical quantities in geometrical percolation transition are the
percolation strength ($P$), the maximum largest size of cluster ($\rm
S_{max}$), the average cluster size ($\left<S\right>$) and so on [28]. Besides
these observables, correlation and fluctuation of the other differently ranked
sizes of cluster, such as the second largest cluster ($\rm S_{2nd}$) and the
third largest cluster ($\rm S_{3rd}$), are also meaningful. We will
systematically study correlation and fluctuations of the differently ranked
sizes of the cluster in 2D Ising percolation transition with the finite size
of the system.

Due to the finite size of the system, we should use the finite size scaling to
study the critical behavior. In particular, if a critical related observable
obey the following formula of the finite size scaling [29]:

$Q(T,L)=F_{Q}(tL^{1/\nu}).$ (1)

In this case, at $T=T_{c}$, we can get

$Q(T_{c},L)=F_{Q}(0).$ (2)

It means the temperature dependence of $Q(T,L)$ with different system sizes
will intersect to a fixed point at $T=T_{c}$. The critical point can be
directly obtained from the temperature dependence of $Q(T,L)$. In addition,
Eq. (2) doesn’t contain any critical exponent. This feature is independent of
the dimensionality $D$ of the lattice or the number of dimensions of the order
parameter. Consequently, it is universal both in 2D and 3D $O(n)$ models, i.e.
$n=1,2,4$. One most typical observable is the Binder-ratio, which has provided
a remarkable identification of the critical point in simulations. In this
paper, we will show that the fixed point behavior can also be directly
observed from skewness and kurtosis of the differently ranked sizes of the
cluster. In addition, we will prove that the fixed point behavior for skewness
and kurtosis is also universal in 2D and 3D $O(n)$ models.

This paper is organized as the following: In section 2, we will discuss the
thermal magnetization transition and geometrical percolation transition in the
Ising model. Techniques of using finite size scaling and fixed point behavior
to search for the critical point are demonstrated. Correlations of the
differently ranked sizes of cluster, the maximum and second largest cluster
($C_{12}$), the maximum and third largest cluster ($C_{13}$), the second and
third largest cluster ($C_{23}$), are studied in Section 3. In Section 4, we
will investigate critical properties of skewness and kurtosis of the maximum,
second and third largest cluster. Signatures of the fixed point, non-
monotonic, and sign changes are systematically discussed, respectively. The
universal feature of fixed point in skewness and kurtosis is discussed in
section 5. A possible application to relativistic heavy-ion collisions is also
discussed. Finally, the results are summarized in Section 6.

##  2 Percolation Transition in the Ising Model

It has been pointed out that the QCD critical point belongs to the same
universality class as 3D Ising model. In this case, magnetization $M$ is the
order parameter. Here, $M=\frac{1}{V}\left<\sum\limits_{i}s_{i}\right>$, where
$s_{i}$ is the spin at site $i$ which can take only two values $\pm 1$. Using
a parametric representation or Monte Carlo simulations, the Ising model can be
extensively applied to study the QCD critical point.

Beside the research of the above thermal magnetization transition, it can also
be applied to study the geometrical percolation transition. In this case, the
spin $s_{i}$ at site $i$ can still take only two values $\pm 1$. The thermal
equilibrium configurations by using Wolff algorithm are still produced [30].
It means that the system is already thermal when studying the geometrical
percolation transition. With a well-defined of the cluster, we can study both
the thermal magnetization transition and geometrical percolation transition in
the Ising model. In this paper, the clusters are constructed by linking the
nearest-neighbor spins with the same spin value. The size of a cluster
($n_{s}$) is the number of sites in a cluster. The maximum cluster size is
corresponding to the largest $n_{s}$ in each event. Similarly, the values of
the second and third largest cluster size are related to the second and third
largest $n_{s}$, respectively.

The simulations are performed on different square lattices of size $N=L^{2}$
with $L=50,100,150,200$ without magnetic field. The Wolff cluster algorithm is
employed to generate events [30]. With given system size and temperature, we
totally run the program for 40,010,000 Wolff Monte Carlo steps, of which 2
million states are used in the calculations. The initial 10,000 steps are used
to insure equilibrium. To make sure the generated configurations are
independent, we only use one state every 20 Monte Carlo steps. To treat the
boundaries of the lattice, various types of boundary conditions are employed
in simulations, such as the periodic boundary condition, free boundary
condition and helical boundary condition and so on. With different types of
boundary condition, the scaling formula is still valid while the scaling
function is different [31]. The periodic boundary condition is employed in our
simulations in this paper. The Bootstrap method is used to estimate the
statistical uncertainty [32, 33]. Except some results for kurtosis of the
differently ranked sizes of the cluster, the statistical errors of other
results are all smaller than the size of the symbols used in this paper.

Figure 1: (a) and (b) Distributions of the maximum and second largest cluster
size at different system sizes near the thermal critical point, respectively.

Figs. 1(a) and (b) show the distribution of the maximum and second largest
cluster sizes near the thermal magnetization transition temperature,
respectively. The size of the maximum largest cluster decreases with the
temperature near the thermal magnetization critical point. On the contrary,
the size of the second largest cluster size increases with temperature. It is
hard to extract any information about critical point from these two original
distributions. For the system with finite volume, We should use the finite
size scaling to study the critical behavior and search for the critical point.

According to the finite size scaling [31, 34], the scaling function of the
maximum cluster near the critical point is

$S_{\rm max}=L^{d-\beta_{p}/\nu_{p}}\cdot F_{s}(tL^{1/\nu_{p}}).$ (3)

Here $d=2$ is the spacial dimensionality of the system.
$t=\left(T-T_{c}\right)/T_{c}$ is the reduced temperature and $T_{c}$ is the
critical temperature. $\nu_{p}$ and $\beta_{p}$ are the percolation exponents.
Eq.(3) means that the scaled variable $S_{\rm max}/L^{d-\beta_{p}/\nu_{p}}$
with different values of $t$ and $L$ can be described by a single scaling
function $F_{s}(tL^{1/\nu_{p}})$.

Figure 2:  (Color online) (a) and (b) The finite size scaling of the maximum
and second largest cluster, respectively.

Figs. 2(a) and (b) demonstrate finite size scaling of the maximum and second
largest cluster near the critical point. The values of $\beta_{p}$ and
$\nu_{p}$ are consistent with critical exponents of the percolation transition
[35]. The exponent of the percolation correlation length $\nu_{p}$ is equal to
the exponent of the magnetization correlation length $\nu_{m}$, $\nu_{p}$ =
$\nu_{m}$ =1. And the critical percolation temperature coincides with the
critical magnetization temperature. It can be clearly observed that they both
satisfy the finite size scaling in the whole critical region. Although the
scaling function is different, the critical exponent of the second largest
cluster is the same as that of the maximum largest cluster.

Currently, the finite size scaling has been applied to study the QCD critical
point in heavy-ion collisions [17, 18]. However, as we had discussed in Ref.
[Potts-2], this method often lacks precision, since it is hard to select the
best scaling function through visual observation.

An improved method to identify the critical point is the fixed point behavior
[36]. At $T=T_{c}$, the scaling function $F_{s}$ in Eq. (3) becomes

$F_{s}(0)=S_{\rm max}/L^{d-\beta_{p}/\nu}.$ (4)

It shows $S_{\rm max}/L^{d-\beta_{p}/\nu}$ is independent of system size $L$
at $T=T_{c}$. If we plot $S_{\rm max}/L^{d-\beta_{p}/\nu}$ as a function of
$T$, all curves of different system sizes intersect exactly at a fixed point,
i.e., the critical point.

Figure 3:  (a) and (b) The fixed point extracted from the maximum largest and
second largest cluster, respectively.

Figs. 3(a) and (b) present temperature dependence of the the scaled maximum
largest cluster $S_{\rm max}/L^{d-\beta_{p}/\nu}$ and second largest cluster
$S_{\rm 2nd}/L^{d-\beta_{p}/\nu}$, respectively. With different system sizes,
the curves remarkably intersect at $T=T_{c}$. Here we want to address that
$S_{\rm max}$ and $S_{\rm 2nd}$ should be scaled by $L^{d-\beta_{p}/\nu}$ to
obtain the fixed-point. Mathematically, the goodness of the fixed point
behavior can be quantified [36]. However, it is still different with the
observable discussed in Eq. (2), which don’t need to be re-scaled by other
variable. The observables that can directly observe the fixed point will be
demonstrated in the following.

##  3 Correlations of the Differently Ranked Sizes of Cluster

As demonstrated in Fig. 1, in vicinity of the critical point, the size of the
maximum largest cluster is a decreasing function of temperature, while the
size of the second largest cluster increases as a function of temperature. It
indicates the negative correlation of the maximum and the second largest
cluster ($C_{12}$). The quantity of $C_{12}$ is defined as

$C_{12}=\frac{\langle S_{\rm{max}}S_{\rm{2nd}}\rangle}{\langle
S_{\rm{max}}\rangle\langle S_{\rm{2nd}}\rangle}-1.$ (5)

Fig. 4(a) shows $C_{12}$ as a function of $T$ in vicinity of the critical
point. With four different lattice sizes ($L=50,100,150,200$), the non-
monotonic behavior is observed and the values of $C_{12}$ are all negative. It
first drops down below $T_{c}$, and then goes up close to zero. Figs. 4(b) and
(c) show that the correlation of the maximum and the third largest cluster
($C_{13}$) and the correlation of the second and the third largest cluster
($C_{23}$) are not 0, either. When $T>T_{c}$, $C_{13}$ significantly changes
from negative to positive. In contrary, $C_{23}$ shows a sharp change from
positive to negative near the critical point. In addition, the values of
$C_{23}$ are one order of magnitude bigger than that of $C_{12}$ and $C_{13}$.

Figure 4:  (Color online) Upper panel: (a), (b) and (c) Temperature dependence
of $C_{12}$ , $C_{13}$ , and $C_{23}$, respectively. The lines are just a
guide to the eye. Lower panel: the same measurements as in corresponding upper
panel, but they are as functions of the scaling variable $tL^{1/\nu}$.

The structure of those correlations are different. However, Figs. 4(a) to (c)
show that with different system sizes, the fixed point behavior can be
directly observed at critical temperature. They do not need to be re-scaled by
other variable. Theoretically, this feature is very helpful to search for the
critical point. It is similar with the Binder-ratio, which has been
extensively used to search for the critical point [29].

Since the critical point can be identified directly, one can expect the
following finite size scaling formula:

$C(T,L)=F_{C}(tL^{1/\nu})$ (6)

Figs. 4(d) to (f) show the scaling behavior of $C_{12}$, $C_{13}$, and
$C_{23}$, respectively. With different system lattice sizes, the finite size
scaling is held for these three observables. This scaling formula is totally
the same as that in Eq. (1), though the scaling functions are different.
Consequently, correlations of the differently ranked size of cluster are all
sensitive to the critical point. They can directly identify the critical point
through the temperature dependency with various different system sizes.

##  4 Skewness and Kurtosis of the Various Sizes of Cluster Near the Critical
Point

In heavy-ion collisions, due to the finite size system and finite evolution
time, enhancement of the critical correlations and fluctuations is suppressed.
High-order cumulants are therefore proposed to study the QCD phase transition.
It has been proposed that ratio of the high-order cumulants are sensitive to
the critical point. On the other hand, skewness ($S$) and kurtosis ($\kappa$)
can also give us more information about the distribution. Skewness is a
measurement of the symmetry of a distribution whereas kurtosis measures the
peakedness of the distribution. Skewness and kurtosis with the maximum, second
and third largest sizes of cluster can be expressed as

${\rm Skewness}\rm\\_n=\frac{\langle(S_{n}-\langle
S_{n}\rangle)^{3}\rangle}{\langle(S_{n}-\langle
S_{n}\rangle)^{2}\rangle^{1.5}},$ (7)

and

${\rm Kurtosis}\rm\\_n=\frac{\langle(S_{n}-\langle
S_{n}\rangle)^{4}\rangle}{\langle(S_{n}-\langle
S_{n}\rangle)^{2}\rangle^{2}}-3.$ (8)

Figure 5:  (Color online) (a) and (b) and (c) Temperature dependence of
skewness of sizes of the maximum cluster, the second largest cluster, and the
third largest cluster, respectively. (d) and (e) and (f) Temperature
dependence of kurtosis of sizes of the maximum cluster, the second largest
cluster, and the third largest cluster, respectively. The statistical
uncertainty is obtained by Bootstrap method. The solid lines are just a guide
to the eye.

Fig. 5 demonstrates critical behavior of skewness and kurtosis of sizes of the
maximum cluster ($\rm skewness_{\\_max}$ and $\rm kurtosis_{\\_max}$), the
second largest cluster ($\rm skewness_{\\_2nd}$ and $\rm kurtosis_{\\_2nd}$),
the third largest cluster ($\rm skewness_{\\_3rd}$ and $\rm
kurtosis_{\\_3rd}$), respectively. It is found that behavior of the sign
change is not always robust, although it really exists. Fig. 5(a) shows
negative $\rm skewness{\\_max}$ below $T_{c}$. When $T>T_{c}$ it goes up and
finally becomes positive. However, the positive signal is very small comparing
to the negative phenomena. In addition, the positive value is hard to observe
if the system size is not large enough, as shown by $L=50$. When using the
high-order cumulants to study the critical point, behavior of the sign change
is not the necessary requirement.

Values of $\rm skewness_{\\_2nd}$ are systematically larger than $\rm
skewness_{\\_max}$, and $\rm kurtosis_{\\_2nd}$ are also significantly larger
than $\rm kurtosis_{\\_max}$. The previous section had shown that the values
of $C_{23}$ are one order of magnitude bigger than that of $C_{12}$ and
$C_{13}$. That is the reason why we propose studying the critical properties
of the various sizes of cluster, instead only the maximum largest cluster. On
the other hand, Fig. 5 also shows that values of the kurtosis are all one
order of magnitude bigger than the corresponding skewness with differently
ranked size of cluster. It means the higher order cumulant, the more sensitive
to the critical point. If the system size is large enough and the statistics
are sufficient, we can study as many differently ranked clusters and high-
order cumulants as possible.

It is remarkable that the critical point can be directly identified from the
temperature dependence of all these observables. Values of the skewness and
kurtosis are independent of the system size at $T=T_{c}$. They are not zero
within a broad range of temperature near the critical point.

Skewness and kurtosis of differently ranked sizes of cluster show colorful
critical phenomena near the critical point. It is interesting to study these
observables in CSPM or any other model that can give a connection between the
percolation transition and heavy-ion collisions. Here we do not present
critical behavior of ratios of the high-order cumulants of differently ranked
sizes of cluster. As we will mention in the next section, they are also
sensitive to the correlation length, but without the fixed point phenomena.

Finally, the fixed point feature of skewness and kurtosis at the critical
point is not the unique characteristic of the various sizes of cluster in 2D
geometrical percolation transition. The fixed point behavior of skewness and
kurtosis can also be observed in thermal magnetization transition in the 3D
Ising, $O(2)$, and $O(4)$ models [37, 12].

##  5 Universal Feature of Fixed Point in Skewness and Kurtosis

With a finite size system, the critical behavior of the $n-th$ order cumulant
($K_{n}$) of order parameter can be described by the following formula:

${K_{n}}=L^{nd}L^{-n\beta/\nu}F_{n}(tL^{1/\nu}).$ (9)

This formula can be implemented to 2D and 3D thermal magnetization transition
and geometrical percolation transition. Of course, the critical exponents
$\beta$ and $\nu$ are different and determined by the universality class.

Then the $n-th$ order susceptibility is measured by $K_{n}$ scaled by $L^{d}$.
It means:

$\chi_{n}=L^{(n-1)d}L^{-n\beta/\nu}F_{n}(tL^{1/\nu}).$ (10)

For example, if we want to derive the scaling formula of the susceptibility:

${\chi_{2}}=L^{d-2\beta/\nu}F_{K_{2}}(tL^{1/\nu})=L^{\gamma/\nu}F_{2}(tL^{1/\nu}).$
(11)

Here, $\gamma=d\nu-2\beta$ is derived from the following scaling relations of
the critical exponents: $2-\alpha=d\nu$ together with
$\alpha=2-2\beta-\gamma$.

Eq. (9) also tell us that the scaling forms of the skewness
($S=K_{3}/K_{2}^{3/2}$) and kurtosis ($\kappa=K_{4}/K_{2}^{2}$) are
$S=F_{S}(tL^{1/\nu})$ and $\kappa=F_{\kappa}(tL^{1/\nu})$, respectively. They
don’t diverge with the system size. Values of skewness and kurtosis are finite
and they are sensitive to the critical point with a broad range. At $T=T_{c}$,
the value of skewness or kurtosis is a system size independent constant. The
behavior of the system size independency at the critical point is independent
of the dimensionality of the lattice or the number of dimensions of the order
parameter. It is similar with the Binder-ratio, which is a powerful quantity
to identify the critical point with the finite size system. Consequently, the
key point of application of the fixed-point feature is that we need to know
the order parameter of the system. It can only be observed in skewness and
kurtosis of the order parameter.

In heavy-ion collisions, the baselines of skewness and kurtosis of conserved
quantities under many different distributions are investigated. In a normal
distribution, the values of skewness and kurtosis are both zero. But beyond
that, they are all determined by the parameters of the statistical
distributions, such as the Poisson distribution, the Skellam distribution, and
the negative binomial distribution and so on [38]. Currently, the energy
dependence of $\kappa\sigma^{2}$ of net-proton number had already shown a hint
of the non-monotonic behavior [19]. However, the measured skewness and
kurtosis still follow general trends of those statistical distributions [4, 5,
6]. It is interesting to study why $\kappa\sigma^{2}$ show the non-monotonic
signal while $\kappa$ doesn’t show the critical feature of a fixed point or a
large deviation. As we know, there are still many uncertainties in applying
the fixed point behavior to relativistic heavy-ion collisions, such as the
finite time evolution, the volume fluctuation, and effect of the conserved law
and so on. We will consider these effects in the future.

##  6 Summary

In summary, it is interesting to study critical properties of the various
sizes of cluster in geometrical percolation transition model, which can help
us to realize the formation of QGP. We suggest that we should firstly study
percolation transition under finite size of the system, before applying to
CSPM or any other model that can give a connection between the percolation
transition and heavy-ion collisions.

In this paper, we study critical properties of various sizes of cluster with
finite size of the system in 2D Ising model. We demonstrate the finite size
scaling of the maximum and second largest cluster near the critical point. We
also present how to use the fixed point behavior to search for the critical
point.

We systematically study fluctuations and correlations of the differently
ranked sizes of the cluster, especially the quantities that we can directly
identify the critical point from the temperature dependency. For the second
largest cluster, the values of skewness and kurtosis are both one order of
magnitude bigger than those extracted from the maximum largest cluster. And
the correlation $C_{23}$ is also one order of magnitude bigger than $C_{12}$
and $C_{13}$. Consequently, instead of only the maximum largest cluster, we
should study as many differently ranked sizes of cluster as possible, if the
statistics are sufficient and the system size is large enough.

Behavior of the sign change is not always significant for skewness and
kurtosis of differently ranked sizes of cluster, although it really exists.
The non-monotonic behavior is not obvious, either. It is not sufficient to
draw the conclusion only based on one typical phenomenon or one observable.
Instead, we should comprehensively study as many critical related observables
as possible.

At critical point, the fixed point is observed from skewness and kurtosis with
the maximum, second and third largest sizes of cluster, and correlations with
the differently ranked sizes of cluster ($C_{12}$, $C_{13}$, and $C_{23}$).
For skewness and kurtosis, the fixed point behavior is universal and had
already been observed in $O(1)$, $O(2)$ and $O(4)$ model with thermal dynamic
transition. It is worth pointing out that this feature is similar with that
derived from the Binder-ratio, which has provided a remarkable examination of
the critical point. Currently, there are still many uncertainties when
applying the fixed point behavior in relativistic heavy-ion collisions.
However, it is still interesting to study the reason why it is not observed in
skewness and kurtosis in heavy-ion collisions.

## Acknowledgements

This work was supported by the postdoctral science and technology project of
Hubei Province under Grant No. 2018Z27, and the Fundamental Research Funds for
the Central Universities under Grant No. CCNU19ZN019.

## References

  * [1] Y. Aoki, G. Endrodi, Z. Fodor et al. Nature, 443: 675 (2006).
  * [2] S. Ejiri, Phys. Rev. D, 78: 074507 (2008).
  * [3] E. S. Bowman and J. I. Kapusta, Phys. Rev. C, 79: 015202 (2009).
  * [4] L. Adamczyk et al. (STAR Collaboration), Phys. Rev. Lett. 112: 032302 (2014).
  * [5] L. Adamczyk et al. (STAR Collaboration), Phys. Rev. Lett. 113: 092301 (2014).
  * [6] L. Adamczyk et al. (STAR Collaboration), Phys. Lett. B 785: 551 (2018).
  * [7] A. Bazavov et al. (HotQCD Collaboration), Phys. Lett. B 795: 15 (2019).
  * [8] A. Bazavov et al. (HotQCD Collaboration), Phys. Rev. Lett. 123: 062992 (2019).
  * [9] Y. Aoki, S. Borsanyi, S. Durr et al. JHEP 06: 088 (2009).
  * [10] S. Gupta, X. Luo, B.Mohanty et al. Science 332: 1525 (2011).
  * [11] A. Bazavov et al. (HotQCD Collaboration), Phys. ReV. D 96: 074501 (2017).
  * [12] Xue Pan, Lizhu Chen, X. S. Chen et al. Nucl. Phys. A 913: 206 (2013).
  * [13] M. A. Stephanov, Phys. Rev. Lett. 102: 032301 (2009).
  * [14] M. A. Stephanov, Phys. Rev. Lett. 107: 052301 (2011).
  * [15] M. A. Stephanov, K. Rajagopal, and E. V. Shuryak, Phys. Rev. D, 60: 114028 (1999).
  * [16] B. Berdnikov and K. Rajagopal, Phys. Rev. D, 61: 105017 (2000).
  * [17] Roy A. Lacey, Phys. Rev. Lett, 114: 142301 (2015).
  * [18] Roy A. Locey, Peifeng Liu, Niseem Magdy et al. arXiv:1606.08071 [nucl-ex].
  * [19] J. Adam et al. (STAR Collaboration), Phys. Rev. Lett. 126: 092301 (2021).
  * [20] T. Celik, F. Karsch, H. Satz, Phys. Lett. B 97: 128 (1980).
  * [21] Xu Ming-Mei, Yu Mei-Ling, Liu Lian-Shou, Phys. Rev. Lett. 100: 092301 (2008).
  * [22] Yu Mei-Ling, Xu Ming-Mei, Liu Zheng-You et al. Chin. Phys. C 33: 552 (2009).
  * [23] Brijesh K Srivastava, Nucl. Phys. A 926: 142 (2014).
  * [24] M. A. Braun, J. Dias de Deus, A. S. Hircsh et al. Phys. Rep. 599: 1 (2015).
  * [25] R. P. Scharenberg, B. K. Srivastava, and C. Pajares, Phys. Rev. D 100: 114040 (2019).
  * [26] D. Stauffer, Phys. Rep. 54: 1 (1979).
  * [27] A. Coniglio, C. R. Nappi, F. Peruggi et al. J. Phys. A 10: 205 (1977).
  * [28] D. Stauffer and A.Aharony, Introduction to Percolation Theory 2nd. ed. (Taylor and Francis, London, 1992).
  * [29] K. Binder, Z. Phys. B 43, 119 (1981); K. Binder, Rep. Prog. Phys. 60: 487 (1997).
  * [30] U. Wolff, Phys. Rev. Lett. 62: 361(1989).
  * [31] V. Privman, Finite-size Scaling and Numerical Simulation of Statistical Systems (World Scientific, Singapore, 1990) p35-37.
  * [32] B. Efron, Computers and the Theory of Statistics: Thinking the Unthinkable (Society for Industrial and Applied Mathematics, 1979).
  * [33] M. S. Abdallah et al. (STAR Collaboration), arXiv:2101.12413 [nucl-ex].
  * [34] J. Cardy, Scaling and Renormalization in Statistical Physics (Cambridge University Press, Cambridge, UK, 1996).
  * [35] A. L. Stella and C. Vanderzande, Phys. Rev. Lett. 62: 1067(1989).
  * [36] Yanhua Zhang, Yeyin Zhao, Lizhu Chen et al. Phys. Rev. E 100: 052146 (2019).
  * [37] Chen Li-Zhu, Pan Xue, Chen Xiao-Song et al. Chin. Phys. C 36: 727 (2012).
  * [38] X. Luo, B. Mohanty, and N. Xu, Nucl. Phys. A 931: 808 (2014).


# Diffusive loss shaping of quantum cascade laser frequency combs

Mithun Roy<EMAIL_ADDRESS>Department of Electrical and Computer
Engineering, University of Texas at Austin, Austin, Texas, USA  Department of
Electrical Engineering, University of Notre Dame, Notre Dame, Indiana, USA
Zhenyang Xiao  Department of Electrical and Computer Engineering, University
of Texas at Austin, Austin, Texas, USA  Department of Electrical Engineering,
University of Notre Dame, Notre Dame, Indiana, USA  Sadhvikas Addamane  Center
for Integrated Nanotechnologies, Sandia National Laboratories, Albuquerque,
New Mexico 87185, USA  David Burghoff  Department of Electrical and Computer
Engineering, University of Texas at Austin, Austin, Texas, USA  Department of
Electrical Engineering, University of Notre Dame, Notre Dame, Indiana, USA

###### Abstract

Integrated optical frequency combs based on active cavities such as quantum
cascade lasers (QCLs) have emerged as promising light sources in the mid-
infrared and terahertz spectral regions. Their bandwidths are limited by two
separate, yet equally important effects: dispersion and diffusion. However,
while dispersion has been extensively engineered, diffusion—a phenomenon
originating from gain variation—has not. We show theoretically and
experimentally that the addition of carefully-engineered diffusive loss can
enhance the bandwidth of QCL combs. Adding resonant loss to the cavity of a
terahertz QCL can counteract the diffusive effect of the gain medium and
allows broader bandwidth combs to form, fully exploiting the bandwidth and
dynamic range of the gain medium. These results are well-explained by active
cavity mean-field theory. Interestingly, this strategy also permits our
structures to generate soliton-esque pulsed states of light. Our results give
a new degree of freedom for the creation of active chip-scale combs, and can
be applied to a wide array of cavity geometries.

##  I Introduction

In recent years, there has been significant interest in the generation of
equidistant frequency combs formed in cavities with distributed gain. In
particular, quantum cascade lasers (QCLs)—semiconductor lasers that can
operate both in the mid-infrared and terahertz regions and are capable of
producing watt-level output power—are of interest due to their broadband and
compact nature, especially for spectroscopic applications [1, 2, 3],
radiometry [4], metrology, and quantum information science [5]. While QCLs
have gone through tremendous improvements with respect to output power,
operating temperature [6], and frequency range, precise comb formation remains
challenging. For many years, it has been understood that intriguing phase-
locking mechanisms can occur [7], but pulse formation by conventional passive
mode-locking remained difficult [8].

One of the most salient approaches is one that was observed to occur
spontaneously inside QCL Fabry-Pérot cavities [9, 10, 11], one that was
originally believed to relate to a large third-order nonlinearity in
intersubband systems but is now recognized as a broader emergent phenomenon
that arises from the motion of the gain grating. Briefly, gain saturation
combines with the asymmetry in the field at the facets to create an effective
quasi-$\chi^{(3)}$ nonlinearity [12, 13]. This nonlinearity creates a phase
modulation dependent on the phase of the field itself, and this causes the
dynamics of the system to be governed by a phase-driven nonlinear Schrodinger
equation [13]. This is in stark contrast to Kerr combs (which use amplitude
modulation to generate phase modulation) or passive mode-locking (which uses
amplitude modulation to generate amplitude modulation). Thus, the natural
state of these combs is a frequency-modulated (FM) mode of operation, where
the frequency is strongly modulated but the amplitude is not. FM-like states
have been observed in many different semiconductor laser systems including but
not limited to QCLs [14, 15, 16, 17], and has been verified using a number of
techniques, including SWIFTS [18, 19], FACE [20], and upconversion sampling
[21, 22].

However, in contrast to pulse-forming mechanisms, the ultimate bandwidth
limits of FM modes of operation are not well-understood. It is known that FM
combs require low but nonzero dispersion in order to maintain stability, but
even if the optimal dispersion is engineered [11, 23], this is not enough. By
modifying and generalizing the Lugiato-Lefever equation to describe Fabry-
Pérot lasers [13], we previously showed that the primary limiting factor for
these combs is not dispersion, but gain curvature (variation in the gain at
its peak). The fundamental FM comb—referred to as an extendon—can have
infinite bandwidth without gain curvature. Even when the laser is far above
threshold and has ample gain across a broad bandwidth, that is not a guarantee
that the FM comb will actually be able to utilize all this gain. As the
instantaneous frequency detunes from the center of the gain peak, the
instantaneous intensity falls in concert (Fig. 1a). Other frequencies will be
able to lase in this window instead, and the competition between these two
processes leads to chaotic multimode behavior. Alternatively, if the gain is
low, the diffusive effect of gain curvature can instead cause the laser to
instead produce continuous wave (CW) light.

Figure 1: a. Effect of gain curvature on FM comb formation. For an FM comb
with gain curvature, frequency modulation gives rise to intensity modulation
that destabilizes the comb. b. By introducing diffusive loss shaping, in which
resonant structures near the cavity reduce the gain curvature and mitigate
amplitude fluctuations, unwanted amplitude fluctuations are suppressed.

In this work, we introduce the concept of diffusive loss shaping to balance
the effects of gain curvature. By adding carefully-designed resonant loss
structures that are evanescently coupled to the cavity of a terahertz QCL and
controlling their size and distance (Fig. 1b), we are able to flatten the
effective gain of the laser and broaden the comb bandwidth that can be
achieved. We show theoretically and experimentally that this strategy can
produce combs whose bandwidth is close to the absoulte gain bandwidth limit of
a medium—the ability to make a single-mode laser. We verify the coherence of
our combs and, using SWIFTS, demonstrate that our devices can produce not only
FM combs but also pulses. The ability to engineer both dispersion and
diffusion enables new exciting prospects for integrated combs, as the
diffusive loss and gain have not been well-explored in integrated active
cavities and provide new degrees of freedom.

##  II Basic concept and theory

Gain curvature is intrinsic to any laser. The peak of a gain medium is set by
the transition frequency and falls as the frequency is detuned from it, giving
rise to a negative $\frac{\partial^{2}g}{\partial\omega^{2}}$. For a
homogeneously-broadened transition, this curvature is related to the dephasing
time. FM-type combs form in semiconductor lasers (such as QCLs) due to a
combination of gain saturation and with the asymmetry in the field at the
facets. An FM comb with gain curvature can most simply be described in the
presence of gain curvature by a phase-driven nonlinear Schrodinger equation as
[13]

$\displaystyle\frac{\partial E}{\partial
t}=\frac{i}{2}\left(\beta-i\frac{D_{g}}{2}\right)\frac{\partial^{2}E}{\partial
z^{2}}+i\gamma|E|^{2}\angle EE-r(\left|E\right|^{2}-P_{0}),$

where $\beta$ represents dispersion, $D_{g}$ represents gain curvature, the
middle term represents the nonlinear phase potential, and the final term
resists amplitude modulation. Gain curvature is diffusive and therefore acts
as imaginary dispersion. For pulsed lasers, it primarily serves to broaden
pulses, but for FM combs the effect is more subtle. As the solution is already
maximally-chirped, diffusion serves primarily to modulate the output,
converting the self-phase modulation into amplitude modulation.

Figure 2: Gain curvature theory. a. Comb bandwidths that can be achieved given
a gain bandwidth, provided the optimum dispersion is chosen. The linear
relationship demonstrates the necessity of low gain curvature. b. Effect of
gain inhomogeneity on FM comb spectra. Even if a coherent comb can form, a
minor change in the gain can lead to an enormous spectral modulation. c.
Effect of gain curvature on simulated THz QCLs. Without gain curvature, a
laser biased well above threshold can still form broadband extendon combs,
with a linear chirp and an amplitude pulse. d,e. With gain curvature, either
chaos or continuous wave operation is observed (at high and low gain values,
respectively).

However, the core requirement that allows for stable extendon solutions to
form is that the amplitude variation remain small. If the gain curvature grows
too large, this will lead to unwanted amplitude variation (Fig. 1a). Even when
the dispersion is chosen optimally, gain curvature limits the comb bandwidths
that can be achieved. For example, Fig. 2a shows the maximum comb bandwidth
that can be achieved as a function of effective gain bandwidth (simulated
using mean-field theory), demonstrating a linear relationship between the two.
For typical QCL parameters this implies a maximum comb bandwidth of
approximately a third of the full-width-half-maximum (FWHM) of the gain
medium, which is usually much narrower than the range over which the gain
medium can lase. Ideally, a comb should be able to lase at any frequencies
whose gain exceeds the threshold. The effects of gain curvature can have
deleterious effects even if a comb forms. This can be dramatically observed in
the case of double-peaked gain media (Fig. 2b), which are frequently created
as an accidental byproduct of the material growth process of THz QCLs. Even
just 10% of gain modulation can lead to combs whose amplitude varies by two
orders of magnitude.

Without gain curvature, there is no theoretical limit to the bandwidth of an
FM comb. The bandwidth of an extendon is inversely proportional to dispersion
and can be made arbitrarily large simply by choosing smaller dispersion values
or adding more gain [13, 24]. For example, Fig. 2c shows the spectrum and
temporal profile of a THz QCL comb biased well above threshold and without
gain curvature. A constant-amplitude, linear chirp forms (except at the jump
point, where a pulse forms), and the spectrum is highly coherent. However,
when gain curvature is enabled (Fig. 2d), severe amplitude fluctuations
develop and the spectrum becomes chaotic. No comb forms. When the gain is
lowered but curvature is left on (Fig. 2e), the laser instead enters a CW
mode, as the diffusion counteracts the cross-steepening nonlinearity
responsible for FM comb formation. Not only does gain curvature limit the
bandwidth, but it also limits the dynamic range over which combs can form.

In order to remove the deleterious effects of gain curvature, in this work we
engineer diffusive loss structures that add positive gain curvature. We
introduce small disks on both sides of a laser cavity, covering almost the
entire length of the cavity. These disks are made of the gain medium, and they
introduce resonant loss. Therefore, through an appropriate design, it is
possible to reduce the net gain curvature of the overall gain and increase
comb bandwidth. In contrast to the heterogeneous gain medium concept, our
strategy is less susceptible to material growth uncertainties, as we are able
to tailor our cavities electromagnetically and even for a particular wafer. It
can also account for unwanted curvature in homogeneous wafers as well (such as
double-peaked media). For this work, we demonstrated our concept on terahertz
quantum cascade laser combs made in the metal-metal waveguide platform. Disks
are introduced on both sides of a Fabry-Pérot (FP) laser cavity. Each of the
disks is, in fact, the QCL gain medium in a metal-metal waveguide. Due to the
proximity of these disks to the cavity, light couples into these disks, which
results in cavity loss. The loss has resonance characteristics, which can be
controlled by varying the disk parameters. Figure 3a shows a schematic view of
a unit cell of our design, in which the whole structure is obtained by simply
replicating the unit cell periodically. The structure is symmetric with
respect to the plane bisecting the FP cavity, which mitigates scattering of
the fundamental lateral mode of the FP cavity into higher-order asymmetric
modes.

##  III Results

The active region of the QCL we chose for our work uses a regrown version of
the same design that was that used in Ref. [25]. This is a photon-phonon [26]
design that has low threshold and broadband operation, and when processed into
tunable VECSELs was able to achieve lasing over a range of 880 GHz. However,
as this operation did not all occur at the same bias, it should be considered
an upper bound for the maximum possible bandwidth of the structure. To design
our shaping structures, we assumed that the unsaturated gain was a Lorentzian
with a peak value of 40 $\mathrm{c}\mathrm{m}^{-1}$ and had a FWHM of 0.67 THz
[27].

For most of these structures, we seek to add a relatively small amount of
engineered loss, ranging from 2 to 8 $\mathrm{c}\mathrm{m}^{-1}$. Though the
overall effect on the gain profile is therefore relatively minor (Fig. 3b),
this can be sufficient to eliminate the curvature at the top of the gain peak.
Figure 3c shows the individual and total loss obtained by simulating a
structure (COMSOL) that has two disk pairs. Since the FWHM of the total loss
depends on the radii of the disks used, one might include more disk pairs of
varying radii to increase the width of the total loss. Moreover, by varying
the distance between the disks and the laser cavity (i.e., the coupling
distance), losses of different amplitude can be introduced. To verify the
efficacy of our approach, we designed diffusive loss shapers with varying
amplitudes (2, 4, and 8 $\mathrm{c}\mathrm{m}^{-1}$ of loss), and these were
interspersed with Fabry-Pérots. All of these structures were designed with
double-chirped mirrors to eliminate dispersion (designed to compensate for
dispersion of 0.1 $\mathrm{p}\mathrm{s}^{2}\mathrm{/}\mathrm{m}\mathrm{m}$)
[11]. The fabrication of these devices was done following a standard metal-
metal waveguide process.

Figure 3: Diffusive loss shaping of terahertz QCL combs. a. Basic concept of
the design, in which a periodically-arranged series of disks are evanescently
coupled to the laser. By adjusting the radius and distance of the disks, the
losses can be precisely controlled. b. Target design of gain medium for a
structure designed to add 4 $\mathrm{c}\mathrm{m}^{-1}$ of loss. c. Multiple
ring types can be combined to achieve broader loss spectra and to more finely
tune the profile. d. Structures that were designed to demonstrate efficacy of
the diffusive shaping. FP lasers were interspersed with shaping structures
designed to add 2, 4, and 8 $\mathrm{c}\mathrm{m}^{-1}$ of loss. Figure 4:
Effect of diffusive loss shaping on terahertz QCL combs, which spectra shown
on the left, corresponding beatnotes shown in the insets, and beatnote maps
shown on the right. Starting from the top, the structures have 0, 2, 4, and 8
$\mathrm{c}\mathrm{m}^{-1}$ of loss shaping. The spectrum and beatnote map
improves significantly for the device with 4 $\mathrm{c}\mathrm{m}^{-1}$ of
shaping, but the final one adds too much loss and the spectrum bifurcates, and
no stable comb regime was found.

##  IV Discussion

In Fig. 4, we show the spectra and beatnote map for a Fabry-Pérot (indicated
by 0 $\mathrm{c}\mathrm{m}^{-1}$) and three diffusive devices with 2, 4, and 8
$\mathrm{c}\mathrm{m}^{-1}$ of diffusive loss, respectively. The spectra shown
here are the broadest spectra measured with narrow beatnotes for the
respective devices, and the corresponding current biases are indicated on the
map (vertical red line). Spectra were measured with a nitrogen-purged FTIR and
pyroelectric detector, so atmospheric absorption lines are apparent and
prevent measurement above 3.6 THz, where a deep absorption line is present.
The devices tested are 3.9 mm long, 20 $\mu\mathrm{m}$ wide, and contain DCMs
with the same parameters. The beatnote map shows that the FP and devices with
2 and 4 $\mathrm{c}\mathrm{m}^{-1}$ losses operate as frequency combs for
certain bias ranges. Comparing the spectra for these devices, we find that the
comb bandwidth for the 4 $\mathrm{c}\mathrm{m}^{-1}$ loss device is the
broadest, about 200 GHz (30%) broader than that for the device without any
disks. In addition, the beatnote map of this device is by far the cleanest,
operating as a comb essentially across the full dynamic range of the laser
(without any additional tuning). A minor dip in the spectrum appears at the
designed loss frequency (3.21 THz), but the overall flatness, bandwidth, and
dynamic range make this comb much better than the reference Fabry-Pérot. The
spectrum of the 8-$\mathrm{c}\mathrm{m}^{-1}$-loss device is too strongly
affected by the addition of loss, and as a result it has a large hole in the
middle of the spectrum. Additionally, it no longer possessed any stable comb
regimes, only incoherent multimode regimes.

In order to more fully characterize these devices and to understand their
dynamics, coherence and temporal measurements were performed using SWIFTS [11,
18, 19]. A room-temperature Schottky mixer (WM-57 from Virginia Diodes, with
response beyond 4 THz) was used to detect the optical beatnote, and high-
signal-to-noise ratio (SNR) beatnotes (40 dB) could be achieved (Fig. 5a).
Despite the high SNR in the intermediate frequency (IF) chain, the DC-coupled
monitor signal is not intended to be low-noise and had limited SNR. Despite
this, in the broadest-bandwidth comb regime the device has a coherence
spectrum that matches the normal spectrum product well at the frequency ranges
where the spectrum product is above noise (Fig. 5b). To verify the ultimate
bandwidth of the comb and compare with the best theoretical performance, the
spectrum was measured with a higher-dynamic-range superconducting bolometer
(Fig. 5c). The spectrum is flat at the top and has SWIFTS temporal traces
consistent with FM modes of operation. Although the blue side of the spectrum
is difficult to precisely confirm owing to the strong atmospheric absorption
around 3.6 THz, our results establish a comb bandwidth of at least 700 GHz,
which is approximately 80% of the 880 GHz range over which a nominally
identical wafer was able to lase over [25]. While this comparison is not
ideal—even nominally identical wafers can be different, and a single-mode
laser has the benefit of bias tuning—but this illustrates that diffusive loss
shaping can achieve combs that are nearly as broad as the bandwidth of the
gain spectrum that is above threshold.

Interestingly, diffusively-shaped devices possess unexpected pulsed modes of
operation as well. While at most biases the spectra are primarily FM-like
(Fig. 5d), when the comb is biased to very near an unstable regime clear
pulses form instead (Fig. 5e). The pulses are strikingly large—approximately
3-4 times larger than the FM portion of the wave—and would generally be
considered more characteristic of solitons than of extendons. Though QCLs do
not generally form pulses readily, THz QCLs are more amenable to pulse
formation on account of their longer gain recovery times [8, 28, 29, 30, 31].
Though the limited SNR makes precise evaluation of the origin challenging, it
is likely that it originates from the boundary pulse occasionally observed in
FM combs [14, 16]. This boundary pulse appears in simulation as well (e.g.,
see the top-right panel of Fig. IId), but is usually not so pronounced. For
example, in Fig. IId the boundary pulse is twice the height of the FM portion
of the wave, not 3-4 times larger. While a full analytical theory that exactly
captures the effect of the boundary pulse has not yet been developed, it is
known from numerics that multimode behavior arises at the boundary and that
the boundary pulse is largest at the edge of stability (which is the case
here).

Figure 5: Temporal characterization of loss-shaped devices. a. Optical
beatnote of the 4-$\mathrm{c}\mathrm{m}^{-1}$-loss device as measured on a
room-temperature Schottky mixer in the broadest-bandwidth comb state. b.
SWIFTS coherence measurement using the self-referenced scheme [18] at the
broadest-bandwidth bias, demonstrating coherence across the full range over
which the spectrum product is above noise. c. High dynamic range spectrum
measured using a superconducting bolometer, demonstrating 700 GHz of
simultaneous comb bandwidth. d. SWIFTS intensity for a comb in an FM-type comb
regime, in which the intensity is only weakly modulated e. SWIFTS intensity
profile for a comb in a pulsed comb regime, with pulses 3-4 times the
continuous wave portion of the spectrum.

Going forward, our results imply new degrees of freedom in the engineering of
active cavity combs, especially QCLs. Up until now, the development of novel
comb states has primarily focused on the dispersion degree of freedom, and to
the extent that gain has been engineered it has been through the creation of
broadband heterogeneous designs. While heterogeneous designs will likely be
critical for achieving octave-spanning combs, achieving sufficiently flat gain
spectra through MBE growth alone would be a monumental task. The very act of
combining multiple stacks naturally leads to gain variations, and even minor
(few $\mathrm{c}\mathrm{m}^{-1}$) variations can lead to enormous changes in
performance. However, loss shaping allows for the possibility of trimming
these variations, not just using a fixed fabricated design but even
dynamically. Adding bias to these structures would allow the resonant loss to
become resonant gain, giving even finer control of the comb spectrum. This
strategy is also compatible with every type of cavity, including enhanced
Fabry-Pérots [32] and rings [33, 34, 35, 36].

##  V Conclusion

In conclusion, we have demonstrated theoretically and experimentally the
pivotal role of diffusion in integrated QCL combs. By introducing the concept
of diffusive loss shaping—employing resonant structures closely coupled to the
laser cavity—we eliminated gain curvature and achieved comb bandwidths near
the intrinsic limit of the medium (80% of the maximum lasing bandwidth of a
single-mode laser). We verified the coherence of the combs produced, and
surprisingly showed that these devices could yield not just frequency-
modulated comb states, but also pulsed states. This research unveils the
enormous potential harbored by the simultaneous engineering of both dispersion
and diffusion, and the introduction of diffusive loss in active cavities paves
the way for more robust, versatile, and advanced chip-scale frequency comb
systems.

##  VI Acknowledgments

D.B. acknowledges support from ONR grant N00014-21-1-2735, AFOSR grant no.
FA9550-20-1-0192, and NSF grant ECCS-2046772; this research is funded in part
by the Gordon and Betty Moore Foundation through Grant GBMF11446 to the
University of Texas at Austin to support the work of D.B.

This work was performed, in part, at the Center for Integrated
Nanotechnologies, an Office of Science User Facility operated for the US
Department of Energy (DOE) Office of Science. Sandia National Laboratories is
a multimission laboratory managed and operated by National Technology &
Engineering Solutions of Sandia, LLC, a wholly-owned subsidiary of Honeywell
International, Inc., for the US DOE’s National Nuclear Security Administration
under contract DE-NA-0003525. The views expressed in the article do not
necessarily represent the views of the US DOE or the United States Government.

### VI.1

#### VI.1.1

## References

  * Villares  _et al._ [2015] G. Villares, J. Wolf, D. Kazakov, M. J. Süess, A. Hugi, M. Beck, and J. Faist, Applied Physics Letters 107, 251104 (2015).
  * Yang  _et al._ [2016] Y. Yang, D. Burghoff, D. J. Hayton, J.-R. Gao, J. L. Reno, and Q. Hu, Optica 3, 499 (2016).
  * Burghoff  _et al._ [2019] D. Burghoff, N. Han, and J. H. Shin, Optics Letters 44, 2966 (2019).
  * Benirschke  _et al._ [2021] D. J. Benirschke, N. Han, and D. Burghoff, Nature Communications 12, 4244 (2021).
  * Gabbrielli  _et al._ [2022] T. Gabbrielli, N. Bruno, N. Corrias, S. Borri, L. Consolino, M. Bertrand, M. Shahmohammadi, M. Franckié, M. Beck, J. Faist, A. Zavatta, F. Cappelli, and P. De Natale, Advanced Photonics Research 3, 2200162 (2022).
  * Khalatpour  _et al._ [2023] A. Khalatpour, M. C. Tam, S. J. Addamane, J. Reno, Z. Wasilewski, and Q. Hu, Applied Physics Letters 122, 161101 (2023).
  * Wang  _et al._ [2007] C. Y. Wang, L. Diehl, A. Gordon, C. Jirauschek, F. X. Kärtner, A. Belyanin, D. Bour, S. Corzine, G. Höfler, M. Troccoli, J. Faist, and F. Capasso, Physical Review A 75, 031802 (2007).
  * Riccardi  _et al._ [2023] E. Riccardi, V. Pistore, S. Kang, L. Seitner, A. De Vetter, C. Jirauschek, J. Mangeney, L. Li, A. G. Davies, E. H. Linfield, A. C. Ferrari, S. S. Dhillon, and M. S. Vitiello, Nature Photonics , 1 (2023).
  * Hugi  _et al._ [2012] A. Hugi, G. Villares, S. Blaser, H. C. Liu, and J. Faist, Nature 492, 229 (2012).
  * Khurgin  _et al._ [2014] J. B. Khurgin, Y. Dikmelik, A. Hugi, and J. Faist, Applied Physics Letters 104, 081118 (2014).
  * Burghoff  _et al._ [2014] D. Burghoff, T.-Y. Kao, N. Han, C. W. I. Chan, X. Cai, Y. Yang, D. J. Hayton, J.-R. Gao, J. L. Reno, and Q. Hu, Nature Photonics 8, 462 (2014).
  * Opačak and Schwarz [2019] N. Opačak and B. Schwarz, Physical Review Letters 123, 243902 (2019).
  * Burghoff [2020] D. Burghoff, Optica 7, 1781 (2020).
  * Singleton  _et al._ [2018] M. Singleton, P. Jouy, M. Beck, and J. Faist, Optica 5, 948 (2018).
  * Hillbrand  _et al._ [2020] J. Hillbrand, D. Auth, M. Piccardo, N. Opačak, E. Gornik, G. Strasser, F. Capasso, S. Breuer, and B. Schwarz, Physical Review Letters 124, 023901 (2020).
  * Sterczewski  _et al._ [2020] L. A. Sterczewski, C. Frez, S. Forouhar, D. Burghoff, and M. Bagheri, APL Photonics 5, 076111 (2020).
  * Henry  _et al._ [2018] N. Henry, D. Burghoff, Q. Hu, and J. B. Khurgin, Optics Express 26, 14201 (2018).
  * Burghoff  _et al._ [2015] D. Burghoff, Y. Yang, D. J. Hayton, J.-R. Gao, J. L. Reno, and Q. Hu, Optics Express 23, 1190 (2015).
  * Han  _et al._ [2020] Z. Han, D. Ren, and D. Burghoff, Optics Express 28, 6002 (2020).
  * Cappelli  _et al._ [2019] F. Cappelli, L. Consolino, G. Campo, I. Galli, D. Mazzotti, A. Campa, M. S. de Cumis, P. C. Pastor, R. Eramo, M. Rösch, M. Beck, G. Scalari, J. Faist, P. D. Natale, and S. Bartalini, Nature Photonics 13, 562 (2019).
  * Täschler  _et al._ [2021] P. Täschler, M. Bertrand, B. Schneider, M. Singleton, P. Jouy, F. Kapsalidis, M. Beck, and J. Faist, Nature Photonics , 1 (2021), arxiv:2105.04870 .
  * Täschler  _et al._ [2023] P. Täschler, A. Forrer, M. Bertrand, F. Kapsalidis, M. Beck, and J. Faist, Laser & Photonics Reviews 17, 2200590 (2023).
  * Villares  _et al._ [2016] G. Villares, S. Riedi, J. Wolf, D. Kazakov, M. J. Süess, P. Jouy, M. Beck, and J. Faist, Optica 3, 252 (2016).
  * Humbard and Burghoff [2022] L. Humbard and D. Burghoff, Optics Express 30, 5376 (2022), arxiv:2110.04356 .
  * Curwen  _et al._ [2019] C. A. Curwen, J. L. Reno, and B. S. Williams, Nature Photonics , 1 (2019).
  * Amanti  _et al._ [2009] M. I. Amanti, G. Scalari, R. Terazzi, M. Fischer, M. Beck, J. Faist, A. Rudra, P. Gallo, and E. Kapon, New Journal of Physics 11, 125022 (2009).
  * Burghoff  _et al._ [2011] D. Burghoff, T.-Y. Kao, D. Ban, A. W. M. Lee, Q. Hu, and J. Reno, Applied Physics Letters 98, 061112 (2011).
  * Micheletti  _et al._ [2023] P. Micheletti, U. Senica, A. Forrer, S. Cibella, G. Torrioli, M. Frankié, M. Beck, J. Faist, and G. Scalari, Science Advances 9, eadf9426 (2023).
  * Wang  _et al._ [2020] F. Wang, V. Pistore, M. Riesch, H. Nong, P.-B. Vigneron, R. Colombelli, O. Parillaud, J. Mangeney, J. Tignon, C. Jirauschek, and S. S. Dhillon, Light: Science & Applications 9, 51 (2020).
  * Bacon  _et al._ [2016] D. R. Bacon, J. R. Freeman, R. A. Mohandas, L. Li, E. H. Linfield, A. G. Davies, and P. Dean, Applied Physics Letters 108, 081104 (2016).
  * Barbieri  _et al._ [2011] S. Barbieri, M. Ravaro, P. Gellie, G. Santarelli, C. Manquest, C. Sirtori, S. P. Khanna, E. H. Linfield, and A. G. Davies, Nature Photonics 5, 306 (2011).
  * Senica  _et al._ [2023] U. Senica, A. Dikopoltsev, A. Forrer, S. Cibella, G. Torrioli, M. Beck, J. Faist, and G. Scalari, Frequency-modulated combs via on-chip field enhancement (2023), arxiv:2305.01483 [physics] .
  * Meng  _et al._ [2022] B. Meng, M. Singleton, J. Hillbrand, M. Franckié, M. Beck, and J. Faist, Nature Photonics 16, 142 (2022).
  * Kazakov  _et al._ [2022] D. Kazakov, T. P. Letsou, M. Beiser, Y. Zhi, N. Opačak, M. Piccardo, B. Schwarz, and F. Capasso, arXiv:2206.03379 10.48550/arXiv.2206.03379 (2022), arxiv:2206.03379 [nlin, physics:physics] .
  * Opačak  _et al._ [2023] N. Opačak, D. Kazakov, L. L. Columbo, M. Beiser, T. P. Letsou, F. Pilat, M. Brambilla, F. Prati, M. Piccardo, F. Capasso, and B. Schwarz, arXiv:2304.10796 10.48550/arXiv.2304.10796 (2023), arxiv:2304.10796 [nlin, physics:physics] .
  * Khan  _et al._ [2023] M. I. Khan, Z. Xiao, S. J. Addamane, and D. Burghoff, Frequency combs in optically-injected terahertz ring quantum cascade lasers (2023), arxiv:2308.11851 [physics] .


\node (v1) at (-0.5,2) [circle, draw = black, fill = black, inner sep = 0.5mm, label=above:{$3$}] {};
\node (v2) at (-2,1.5) [circle, draw = black, fill = black, inner sep = 0.5mm, label=left:{$2$}] {};
\node (v3) at (0.5,0.5) [circle, draw = black, fill = black, inner sep = 0.5mm, label=above:{$4$}] {};
\node (v4) at (2,1.5) [circle, draw = black, fill = black, inner sep = 0.5mm, label=right:{$5$}] {};
\node (v5) at (0,-1.5) [circle, draw = black, fill = black, inner sep = 0.5mm, label=below:{$1$}] {};
\node [circle, draw = black, fill = black, inner sep = 0.5mm, label=left:{$6$}] (v6) at (-2,-1) {};
\draw (-0.5,2) -- (-2,1.5);
\draw (2,1.5) -- (-0.5,2);
\draw (0.5,0.5) -- (0,-1.5);
\draw (-2,1.5) -- (2, 1.5);
\draw[dashed] (0,-1.5) -- (-0.5,2);
\draw (-2,1.5) -- (0,-1.5) -- (2,1.5);
\draw[fill=gray, opacity=0.7] (0,-1.5) -- (0.5,0.5) -- (2,1.5) -- (0,-1.5);
\draw[fill=gray, opacity=0.5] (-0.5,2) -- (-2,1.5)-- (0.5,0.5) -- (2,1.5) -- (-0.5,2);
\draw[fill=gray, opacity=0.7] (-2,1.5) -- (0,-1.5) -- (0.5,0.5) -- (-2,1.5);
\draw (0,-1.5) -- (-2,-1) -- (-2,1.5);
\draw[fill=gray, opacity=0.5] (0,-1.5) -- (-2,-1) -- (-2,1.5);
\end{tikzpicture}$$

\node (v1) at (-1,0.5) [circle, draw = black, fill = black, inner sep = 0.5mm, label=above:{$2$}] {};
\node (v3) at (0,-1.5) [circle, draw = black, fill = black, inner sep = 0.5mm, label=below:{$1$}] {};
\node (v2) at (1.5,0.5) [circle, draw = black, fill = black, inner sep = 0.5mm, label=above:{$6$}] {};
\node (v4) at (-3.5,-0.5) [circle, draw = black, fill = black, inner sep = 0.5mm, label=left:{$5$}] {};
\node (v5) at (2,-0.5) [circle, draw = black, fill = black, inner sep = 0.5mm, label=above:{$4$}] {};
\node (v7) at (4,-0.5) [circle, draw = black, fill = black, inner sep = 0.5mm, label=above:{$3$}] {};
\draw (-1,0.5) -- (1.5,0.5) -- (0,-1.5) -- (-3.5,-0.5) -- (-1,0.5) -- (0,-1.5);
\draw (-1,0.5) -- (2,-0.5) -- (0,-1.5) -- (4,-0.5) -- (2,-0.5);
\draw (1.5,0.5) -- (4,-0.5);
\draw[dashed] (-3.5,-0.5) -- (1.5,0.5);
\draw[fill=gray, opacity=0.8] (-1,0.5) -- (-3.5,-0.5) -- (0,-1.5) -- (-1,0.5);
\draw[fill=gray, opacity=0.8] (-1,0.5) -- (0,-1.5) -- (2,-0.5) -- (-1,0.5);
\draw[fill=gray, opacity=0.8] (v5) --(0,-1.5) node (v8) {}-- (4,-0.5) -- (2,-0.5);
\draw[fill=gray, opacity=0.3] (1.5,0.5) -- (-1,0.5) node (v6) {} -- (2,-0.5) -- (4,-0.5) -- (1.5,0.5) node (v9) {};
\draw[fill=gray, opacity=0.4] (-1,0.5) -- (1,-0.18) -- (1.5,0.5) -- (-1,0.5);
\end{tikzpicture}$$

Let $\mathfrak{g}$ be a complex Lie algebra of a connected linear algebraic Lie group $G.$ If $\ad^*$ is the coadjoint representation of $\mathfrak{g},$ then we have $\mathfrak{g}(\varphi)=\{x\in\mathfrak{g}|(\ad^*x)(\varphi)=0\},$ and we have the following definition.

An element $\varphi$ of $\mathfrak{g}^*$ is said to be of reductive type if $\faktor{\mathfrak{g}(\varphi)}{\mathfrak{z}}$ is a reductive Lie algebra whose center consists of semisimple elements of $\mathfrak{g},$ where $\mathfrak{z}$ is the center of $\mathfrak{g}.$ A Lie algebra $\mathfrak{g}$ is called quasi-reductive if it has linear forms of reductive type.

Quasi-reductive Lie algebras were first introduced by Duflo (<cit.>, 1982) because of their applications in harmonic analysis. The problem of their classification has been of recent interest among those studying biparabolic (“seaweed") Lie algebras. Panyushev (<cit.>, 2005) showed that all biparabolic subalgebras of $\mathfrak{sl}(n)$ and $\mathfrak{sp}(2n)$ and some biparabolic subalgebras of $\mathfrak{so}(n)$ are quasi-reductive (cf. <cit.> and <cit.>), and the full classification of quasi-reductive biparabolics in the classical cases was completed by Duflo (<cit.>, 2012) and Panyushev and Yakimova (<cit.>, 2018). Baur and Moreau (<cit.>, 2011) identified all quasi-reductive parabolic subalgebras of the exceptional Lie algebras, and Ammari (<cit.>, 2022) then completed the full classification of biparabolics of exceptional type. However, the study of quasi-reductive Lie algebras outside the category of biparabolics has remained mostly unexplored. We now offer some insight into the classification in the Lie poset case, with the restriction to index one in mind.

The following theorem appears in (<cit.>, 2022). (Reference dissertation)

Let $\mathfrak{g}$ be a simple Lie algebra with Cartan subalgebra $\mathfrak{h}\subset\mathfrak{g}.$ If $\mathfrak{q}$ is an index-one, quasi-reductive Lie algebra, then $\mathfrak{q}$ is contact.

In fact, the main result of <cit.> states that a biparabolic Lie algebra $\mathfrak{q}$ is contact precisely when $\mathfrak{q}$ is quasi-reductive and has index one. By Definition <ref>, this result implies that, in particular, each contact biparabolic has a diagonalizable Reeb vector. Note that for each toral-pair in Theorems <ref>, <ref>, <ref>, and <ref>, the Reeb vector is diagonalizable. This observation, along with Theorem <ref>, leads us to the following conjecture.

An index-one, type-A Lie poset algebra is quasi-reductive if and only if it is contact.

\def\Node{\node [circle, fill, inner sep=0.5mm]}
\Node (64) at (-18,0){};
\Node (65) at (-18,1){};
\Node (66) at (-18,2){};
\Node (67) at (-18.5,3){};
\Node (68) at (-17.5,3){};
\draw (64)--(65)--(66)--(67);
\draw (66)--(68);

\node at (-18,-1){$\mathcal{S}_1=\mathcal{Q}_1$};

\draw[->] (-17.25,1.5)--(-16.25,1.5);
\node at (-16.75,2){$A_1$};

\Node (53) at (-14.5,0){};
\Node (54) at (-14,1){};
\Node (55) at (-15,1){};
\Node (56) at (-14,2){};
\Node (57) at (-15,2){};
\Node (58) at (-14,3){};
\Node (59) at (-15,3){};
\draw (53)--(54)--(56)--(58);
\draw (53)--(55)--(57)--(59);
\draw (54)--(57)--(58);
\draw (55)--(56)--(59);

\Node (60) at (-15.5,0){};
\Node (61) at (-15.5,1){};
\Node (62) at (-15.5,2){};
\Node (63) at (-16,3){};
\draw (60)--(61)--(62)--(63);
\draw (62)--(59);

\node at (-15,-1){$\mathcal{Q}_2$};

\draw[->] (-13.25,1.5)--(-12.25,1.5);
\node at (-12.75,2){$D_1$};

\Node (39) at (-9,1.5){};
\Node (40) at (-8.5,3){};
\Node (41) at (-10,3){};
\Node (42) at (-10.5,0){};
\draw (42)--(39)--(40);
\draw (39)--(41);

\Node (43) at (-10,1){};
\Node (44) at (-11,1){};
\Node (45) at (-10,2){};
\Node (46) at (-11,2){};
\Node (47) at (-10,3){};
\Node (48) at (-11,3){};
\draw (42)--(43)--(45)--(47);
\draw (42)--(44)--(46)--(48);
\draw (43)--(46)--(47);
\draw (44)--(45)--(48);

\Node (49) at (-11.5,0){};
\Node (50) at (-11.5,1){};
\Node (51) at (-11.5,2){};
\Node (52) at (-12,3){};
\draw (49)--(50)--(51)--(52);
\draw (51)--(48);

\node at (-10.25,-1){$\mathcal{Q}_3$};

\draw[->] (-8.25,1.5)--(-7.25,1.5);
\node at (-7.75,2){$C$};

\Node (21) at (-2.5,0){};
\Node (22) at (-2.5,1){};
\Node (23) at (-3,2){};
\Node (24) at (-3.5,3){};
\Node (25) at (-3.5,0){};
\Node (26) at (-3.5,1){};
\draw (21)--(22)--(23)--(24);
\draw (25)--(26)--(23);

\Node (27) at (-4,1.5){};
\Node (28) at (-5.5,0){};
\Node (29) at (-5,3){};
\draw (28)--(27)--(29);
\draw (27)--(24);

\Node (30) at (-6,1){};
\Node (31) at (-5,1){};
\Node (32) at (-6,2){};
\Node (33) at (-5,2){};
\Node (34) at (-6,3){};
\draw (28)--(30)--(32)--(34);
\draw (28)--(31)--(33)--(29);
\draw (30)--(33)--(34);
\draw (31)--(32)--(29);

\Node (35) at (-6.5,0){};
\Node (36) at (-6.5,1){};
\Node (37) at (-6.5,2){};
\Node (38) at (-7,3){};
\draw (35)--(36)--(37)--(38);
\draw (37)--(34);

\node at (-4.5,-1){$\mathcal{Q}_4$};

\draw[->] (-1.75,1.5)--(-0.75,1.5);
\node at (-1.25,2){$F$};

\Node (1) at (0,0){};
\Node (2) at (0,1){};
\Node (3) at (0,2){};
\Node (4) at (-0.5,3){};
\Node (5) at (0.5,3){};
\draw (1)--(3)--(4);
\draw (3)--(5);

\Node (6) at (1,0){};
\Node (7) at (0.5,1){};
\Node (8) at (1.5,1){};
\Node (9) at (0.5,2){};
\Node (10) at (1.5,2){};
\Node (11) at (1.5,3){};
\draw (6)--(7)--(9)--(5);
\draw (6)--(8)--(10)--(11);
\draw (7)--(10)--(5);
\draw (8)--(9)--(11);

\Node (12) at (2.5,1.5){};
\Node (13) at (3,3){};
\draw (6)--(12)--(11);
\draw (12)--(13);

\Node (14) at (3,0){};
\Node (15) at (4,0){};
\Node (16) at (3,1){};
\Node (17) at (4,1){};
\Node (18) at (3.5,2){};
\draw (14)--(16)--(18)--(13);
\draw (15)--(17)--(18);

\Node (19) at (4.5,1){};
\Node (20) at (4,2){};
\draw (15)--(19)--(20)--(13);
\draw (14)--(20);

\node at (2,-1){$\mathcal{Q}_5=\mathcal{P}$};
\end{tikzpicture}$$

Contact sequence of $\mathcal{P}$

\def\Node{\node [circle, fill, inner sep=0.5mm]}

\Node[label=left:{\small 1}] (1) at (0,0){};
\Node[label=left:{\small 2}] (2) at (0,1){};
\Node[label=left:{\small 3}] (3) at (0,2){};
\Node[label=above:{\small 4}] (4) at (-0.5,3){};
\Node[label=above:{\small 5}] (5) at (0.5,3){};
\draw (1)--(3)--(4);
\draw (3)--(5);

\Node[label=left:{\small 6}] (6) at (1,0){};
\Node[label=right:{\small 7}] (7) at (0.5,1){};
\Node[label=left:{\small 8}] (8) at (1.5,1){};
\Node[label=right:{\small 9}] (9) at (0.5,2){};
\Node[label=left:{\small 10}] (10) at (1.5,2){};
\Node[label=above:{\small 11}] (11) at (1.5,3){};
\draw (6)--(7)--(9)--(5);
\draw (6)--(8)--(10)--(11);
\draw (7)--(10)--(5);
\draw (8)--(9)--(11);

\Node[label=left:{\small 12}] (12) at (2.5,1.5){};
\Node[label=above:{\small 18}] (13) at (3,3){};
\draw (6)--(12)--(11);
\draw (12)--(13);

\Node[label=left:{\small 13}] (14) at (3,0){};
\Node[label=left:{\small 14}] (15) at (4,0){};
\Node[label=left:{\small 15}] (16) at (3,1){};
\Node[label=left:{\small 16}] (17) at (4,1){};
\Node[label=left:{\small 17}] (18) at (3.5,2){};
\draw (14)--(16)--(18)--(13);
\draw (15)--(17)--(18);

\Node[label=right:{\small 19}] (19) at (4.5,1){};
\Node[label=right:{\small 20}] (20) at (4,2){};
\draw (15)--(19)--(20)--(13);
\draw (14)--(20);
\end{tikzpicture}$$

Labeled Hasse diagram of $\mathcal{P}$