action one could without problems use the property of integration by parts.
Moreover, for the standard Riemann tensor with covariant derivatives we could
use also the differential Bianchi identities (possibly also in the contracted
forms) to simplify and reduce some scalar invariant terms in the action. These
were the advantages for GR-adapted basis. We remark that these nice properties
are not anymore present in the formalism with all conformal objects. We cannot
so simply integrate by parts conformal covariant derivatives under the
spacetime volume integrals, there are correction terms with the Schouten
tensor which on this side take us out of the conformal formalism. Similarly,
we cannot use any form of the conformal Bianchi identities (with conformal
derivatives on the conformal Weyl tensor) to reduce terms further. Still we
can use first Bianchi identities (cyclicity) because they work equally well
for the Weyl tensors and also originally for Riemann tensors as in standard
differential geometry.

Therefore, roughly understanding the construction of invariant scalar
Lagrangian densities with conformal objects is hampered and one does not have
many options for reduction of terms, and one would conclude that the number of
possible different scalar invariants is bigger than in the corresponding GR
case. However, this is only a naive thinking since all the examples show that
the usage of conformal objects instead of just GR-adapted ones is beneficial
for construction of conformal operators, tensors, expressions and conformally
invariant actions as well. Of course, the last sentence is a very natural
conclusion: for general conformal objects use basic conformal objects as their
building elements. One could compare the lost freedom of integration by parts
and differential Bianchi identity in the conformal case with the additional
complications that one inevitably meets in the GR basis. In the last case, one
also has to consider terms with the contractions of Riemann tensor (Ricci
tensor and Ricci scalar) since the trace is not removed there from the basic
curvature Riemann tensor. The advantage here of the Weyl tensor in the
conformal case is that it is completely trace-free so we do not have to worry
about these additional contractions. Barring these differences the
construction of actions in conformal or GR case are exactly the same and they
lead to roughly the same structures (when in the conformal case we use
conformal derivatives instead of GR ones and conformal Weyl curvatures instead
of the Riemann ones). The estimate is that however, in the conformal case we
can create less number of terms and then this conformal basis for the
invariant actions is smaller and hence more effective. An example here is the
basis of all scalar invariant terms with six derivatives and constructed
exactly with three powers of curvatures. In general dimension, we have $8$
terms in the GR basis, while in the conformal one we have only $2$ terms. The
effectiveness of using only the trace-free objects is clearly seen here.

As a kind of point in the middle we can understand our new method when we
construct all possible terms in $D_{(2N)}$ as different orders in a sequence
of precisely $2N$ conformal covariant derivatives. In this way we do not
explicitly use conformal curvature tensors and we only consider combinations
of various orders of derivatives. As explained above this method has some
advantages and disadvantages. From the point of view of counting number of
unknown coefficients and of linear algebraic equations to solve, again this
method stays in the middle between the case of very efficient conformal
objects (derivatives and curvatures) and the GR case. As seen in the explicit
example for the generalization of Hamada operator here we had $5$ (after
normalization) unknown coefficients and effectively only two equations for
them giving us a remaining three-parameter freedom which was totally
consistent with the previous findings. We compare it with the case of Wunsch
minimal operator $D_{(6)}^{\min}$, where only one relative coefficient needed
to be determined from effectively one linear equation. This minimal form of
the operator could be seen as a result of commutation of conformal derivatives
and of various ordering of terms there; this is possible in every dimension
different than the critical one for a given $2N$ number of derivatives.
However, by the analysis of solutions to the system in (158), we see that it
is impossible to write such a minimal term as a combination with less than
three different orders of conformal derivatives. Again, the expression
provided by Wunsch as minimal in (131) is the shortest in the number of
conformal terms and then in this case the algebraic system to solve is the
smallest one with smallest number of unknown coefficients to determine. That
is why in that case the method is the most effective.

Here we also discuss other generalizations of the conformally invariant
actions with scalar fields. When the scalar theory with the action functional

$S_{(2N)}=\\!\int\\!d^{d}x\sqrt{|g|}\phi D_{(2N)}\phi$ (182)

is considered on flat spacetime, then the theory is free (but certainly with
higher derivatives – $2N$ derivatives in the scalar propagator from the
leading term when all curvatures vanish) and hence there are no interactions
and no divergences generated at the quantum level. This is because in this
case we have that $D_{(2N)}=\square^{N}=(\partial^{2}){}^{N}$ in Cartesian
coordinates. One could here analyze a slightly more general case, when we add
conformal interaction terms to the action with a form

$S_{{\rm int}}=\\!\int\\!d^{d}x\sqrt{|g|}\lambda\phi^{\alpha},$ (183)

where the exponent $\alpha=\frac{2d}{d-2N}=2+\frac{4N}{d-2N}$ can be defined
in any dimension different than a critical one $d_{{\rm crit}}=2N$. In (183)
$\lambda$ is a dimensionless parameter of self-interaction of the scalar
field. The action (183) or the total action (182)+(183) are both classically
conformally invariant. This way one adds non-linearities and classical
interactions to the conformal action of the theory. Above we have added non-
linearities without derivatives, where only generalized powers of the scalar
field are present, however the more complicated other ones are also still
possible here and they could potentially include derivatives and higher than
quadratic monomials. For example in $d=4$ and when $N=1$ this is a standard
two-derivative theory giving rise to classically scale-invariant Klein-Gordon
model on flat spacetime with $\phi^{4}$-type of scalar quartic self-
interactions – therefore it is a conformal Higgs model. Of course, this has
vast applications to the conformal Standard Model of particle physics, where
one knows that only the classical phenomenon of spontaneous symmetry breaking
destroys the classical conformal symmetry due to the mass of the Higgs
particle. But as the zero-order approximation (or for very high temperatures,
or very early Universe) one can take a conformal model as a very good and
symmetric classical starting point of further considerations.

##  IX Conclusions

The six-derivative conformal scalar operator was originally found by Hamada in
its critical dimension of spacetime, $d=6$. We showed that it is possible to
generalize this construction to arbitrary even dimensions of spacetime,
provided that this is bigger than the critical dimensions $d_{\rm crit}=6$. We
gave the explicit form of its coefficients of expansion in the basis of GR
scalars (invariant terms with respect to diffeomorphisms). The critical case
is the situation for the Hamada operator with six derivatives in six spacetime
dimensions. As the second main part of this work, we also proved a general
theorem that a scalar conformal operator with $n$ derivatives in $d=n-2$
dimensions is impossible to construct. One notices that the critical dimension
of the operator equals to the number of derivatives present in the leading
term $\square^{n}$ of the expansion in curvature terms. We can conjecture
extending our theorem that for any even dimension $d<n=d_{\rm crit}$ such
construction is impossible since the coefficients blow up. The operators in
question can be constructed only in dimensions $d\geqslant d_{\rm crit}$, when
also the conformal weight of the basic scalar field $w$ must be non-positive.
This is an interesting pattern for the construction of various conformal
operators that perhaps we can explore also for other vectors and tensors
(higher rank representations of the Lorentz group) [118].

In this article, we showed different ways of how to generalize the original
construction of Hamada conformal operator when acting on scalar fields. We
started in the situation with the critical dimension for this operator and
since it must contain six derivatives in the leading part, then this situation
was in $d=6$ spacetime dimensions. We extended the construction to arbitrary
dimension $d$ but also found that in dimensions $d=2$ and $d=4$ such
construction is impossible due to various singularities. Our form of the
generalized Hamada operator works naturally in spacetimes with Minkowskian
signature as well as for the completely Euclidean spaces useful for purely
geometrical applications of conformal symmetry. In the course of our article,
we first discussed the explicit and straightforward method of finding
coefficients in front of all terms containing up to six derivatives and built
with standard gravitational curvatures and standard GR-covariant derivatives
acting on them. This required solving some big system of linear algebraic
equations. Finally there we encountered the 3-parameter freedom since the
generalized form of Hamada operator is not unique and there are various
parameters or terms that could be added to its “minimal” form for free.

In the last sections of this article, we embarked on explaining and applying
the formalism of conformal covariant derivatives as constructed by Wunsch.
With this theoretical tool, we showed that the construction can be a bit
simplified for the case of six-derivative operator, while it was completely
clear for the cases of studies of Ørsted-Penrose and Paneitz operator. In the
last two simplifying cases, the usage of conformal covariant derivatives
$\overset{c}{\nabla}_{\mu}$ and the conformally covariant box operator
$\overset{c}{\square}$ give directly answers for the conformal operators
containing explicitly 2 or 4 derivatives in their leading term respectively.
However, for the case of six-derivative operator (and all higher) this
construction is not fully and completely successful since even the minimal
operators must also contain correction terms which are higher in conformal
curvatures (Weyl or Bach tensors). We also discussed the possible ambiguities
and freedom that we found for construction using these lines for
generalization of Hamada operator. By comparison with the results from the
previous sections, we found a full agreement numerically for the form of
constructed conformally covariant operators as well as for a qualitative
description of the level of freedom that we have in such constructions. Of
course, using both methods we found that the freedom for the generalized
Hamada conformal operator is contained in three arbitrary real parameters.

Finally, we showed that using a combination of different orders of six
conformal covariant derivatives, one reproduces the same results as before for
the generalization of the conformal Hamada operator. However, this last
construction is not completely successful in the critical dimension ($d_{{\rm
crit}}=6$ for Hamada operator), but works in any dimensions $d\neq 2,4,6$.
Therefore we provide some _new_ results written in a _new_ and distinct way.
We also showed the roadmap how one can think about looking for generalization
for operators containing more than six derivatives in their leading terms.

We also explained using algebraic methods and arguments why the conformal
operator with six derivatives in its action must necessarily contain curvature
correction terms (except the case of special dimension $d=3$) and why the
coefficient of these additional correcting terms arise only as a solution of
system of algebraic equations and why never they can be obtained from some
other constructions using different definition of conformal covariant
derivative or conformal box operator. For this last assertion, we used and
studied the fact that the generalization of the conformal Hamada operator
cannot be defined in $d=4$ and traced it back to the singularity of the
explicit coefficients of curvature correction terms which must contain in the
denominators factors $(d-4)$. This was of course, a special case of a general
theorem, that we have also proved here using algebraic and infinitesimal
conformal transformation methods, stating that the conformal operator with the
leading term containing precisely $2N$ derivatives cannot exist in spacetime
dimension $d=2N-2$.

Finally, we comment that the operators considered on the classical level that
we have studied in this paper, can bring a lot of applications. First with
them we can construct conformally invariant actions for scalar field theories
containing precisely $2N$ derivatives and in any dimensions [116]. These
actions are necessarily quadratic in scalar fields (so they give contributions
only to the propagator on the quantum level of quantum scalars), but we can
easily also add here some interaction terms (without derivatives, which
examples of were given above, or containing also derivatives). Such conformal
actions then define good scalar models, where the conformal symmetry is
preserved and realized explicitly on the classical level. (As it is expected
and proved by various quantum computations these pure scalars models (but
coupled to external gravitation) unfortunately cease to be conformal on the
quantum level of considerations.) Obviously, we can use them for the
definition of new conformal field theories, which are put on curved spacetime
backgrounds, where the conformal curvature does not vanish. These in turn, may
have new applications for the famous AdS/CFT correspondence in the new
framework where the boundary theory is placed, not on usual flat background,
but when it is characterized by some non-zero conformal curvature tensors.
From the mathematical point of view these new conformal operators are relevant
for the problems of coupling of scalar fields to external geometry (in
particular to conformal geometry). They are also essential for the study of
6-dimensional Yamabe problem, analysis of Huyghens principle and conformal
wave equations as well as to the problems of the analytic properties of
conformal Laplacians operators (their generalizations and their powers), their
determinants, logarithmic traces, spectra and positive eigenvalues. Some of
the mentioned topics are discussed in [66, 107, 108, 109, 110, 115, 111] and
in the further literature cited there.

On the other hand, we can also mention applications to recent studies of
conformal anomaly induced action, where these examples of generalized Paneitz
and Hamada operators play significant roles [94, 95, 96]. These actions
succeed in trying to capture all important non-local quantum effects directly
related to the conformal anomaly (and therefore neglecting other less
important effects in the realm of all full quantum corrections [120, 121, 122,
114]). In recent studies they were for example applied to study quantum
gravitational perturbative corrections in cosmology [123]. As different from
standard physical applications, the methods of conformal symmetry in Euclidean
space are also quite powerful and they are already exploited in theoretical
considerations in computer graphics (CG).

Eventually, we can discuss the possible future applications of the found form
of the Hamada operators in any dimensions $d$ of space(spacetime). Generally
such scalar conformal operators with six derivatives constitute a basis for
CFT algebras of primary scalar operators coupled to external geometry (so on
non-trivial background spacetimes) hence they will give rise to new CFT’s
coupled to gravity. Such operators may have long ranging applications from
studies of conformal anomalies (i.e. integrating various nonlocal terms like
in [94, 95, 96]), through applications of scalar fields in cosmology. There
could be applications also for CFT’s: like the description near fixed points
(FP) of RG flows of theories with six-derivative scalars coupled to background
geometry, so related to asymptotic safety program in gravity (where there is a
conformal symmetry in the vicinity of FP). Some obvious applications can be
probably found in condensed matter physics, in the situations, where there are
some fundamental scalars around like for example for an excitation of graphene
curved sheets varying in time ($d=3$ case) [91, 92, 90].

All these present (and also hopefully expanded in future) motivations ensure
us that the studies that we have undertaken in this paper about the
generalization of the Hamada conformal operator are important both from the
mathematical as well as physical points of view.

## Acknowledgements

We thank Prof. Dr. Ilya Shapiro for useful comments and initial encouragement
for studying the topic of this work. L.R. wishes to thank the Department of
Physics, Federal University of Juiz de Fora for kind hospitality. P.R.B.R.V.
is grateful to CAPES for supporting his Ph.D. projects.

## References

  * [1] H. Weyl, Annalen Phys. 54, 117-145 (1917)
  * [2] H. Weyl, Annalen Phys. 59, 101-133 (1919)
  * [3] P. Di Francesco, P. Mathieu and D. Senechal, Conformal Field Theory, Springer, New York, 2012
  * [4] A. A. Belavin, A. M. Polyakov and A. B. Zamolodchikov, Nucl. Phys. B 241, 333-380 (1984)
  * [5] J. Polchinski, Nucl. Phys. B 303, 226-236 (1988)
  * [6] L. Rachwał, Universe 4, no.11, 125 (2018) [arXiv:1808.10457 [hep-th]].
  * [7] H. Weyl, Z. Phys. 24, 230-232 (1923)
  * [8] H. Weyl, Z. Phys. 56, 330-352 (1929)
  * [9] R. Bach, Math. Z. 13, 119-133 (1922)
  * [10] R. Bach, H. Weyl, Math. Z. 13, 134-145 (1922)
  * [11] S. Deser, Annals Phys. 59, 248-253 (1970)
  * [12] F. Englert, C. Truffin and R. Gastmans, Nucl. Phys. B 117, 407-432 (1976)
  * [13] M. Kaku, P. K. Townsend and P. van Nieuwenhuizen, Phys. Lett. B 69, 304-308 (1977)
  * [14] E. S. Fradkin and A. A. Tseytlin, Nucl. Phys. B 201, 469-491 (1982)
  * [15] E. S. Fradkin and A. A. Tseytlin, Phys. Rept. 119, 233-362 (1985)
  * [16] M. Pawlowski and R. Raczka, Found. Phys. 24, 1305-1327 (1994) [arXiv:hep-th/9407137 [hep-th]].
  * [17] P. D. Mannheim and D. Kazanas, Astrophys. J. 342, 635-638 (1989)
  * [18] J. D. Barrow and S. Cotsakis, Phys. Lett. B 214, 515-518 (1988)
  * [19] G. de Berredo-Peixoto and I. L. Shapiro, Phys. Rev. D 70, 044024 (2004) [arXiv:hep-th/0307030 [hep-th]].
  * [20] P. D. Mannheim, Found. Phys. 42, 388-420 (2012) [arXiv:1101.2186 [hep-th]].
  * [21] P. Jizba, H. Kleinert and F. Scardigli, Eur. Phys. J. C 75, no.6, 245 (2015) [arXiv:1410.8062 [hep-th]].
  * [22] P. Jizba, L. Rachwał and J. Kňap, Phys. Rev. D 101, no.4, 044050 (2020) [arXiv:1912.10271 [hep-th]].
  * [23] D. M. Ghilencea, JHEP 03, 049 (2019) [arXiv:1812.08613 [hep-th]].
  * [24] T. Fulton, F. Rohrlich and L. Witten, Rev. Mod. Phys. 34, 442-457 (1962)
  * [25] R. Penrose, Gen.Rel.Grav. 43 (2011) 901-922 ,
  * [26] S. W. Hawking, A. R. King and P. J. Mccarthy, J. Math. Phys. 17, 174-181 (1976)
  * [27] R. Penrose, Cycles of Time - An Extraordinary New View of the Universe, Knopf Doubleday Publishing Group, New York, 2010
  * [28] J. D. Bekenstein, Annals Phys. 82, 535-547 (1974)
  * [29] S. Deser, Phys. Lett. B 134, 419-421 (1984)
  * [30] I. Jack and H. Osborn, Nucl. Phys. B 343, 647-688 (1990)
  * [31] A. Cappelli, D. Friedan and J. I. Latorre, Nucl. Phys. B 352, 616-670 (1991)
  * [32] J. Maldacena, [arXiv:1105.5632 [hep-th]].
  * [33] G. ’t Hooft, [arXiv:1410.6675 [gr-qc]].
  * [34] G. ’t Hooft, Int. J. Mod. Phys. D 26, no.03, 1730006 (2016)
  * [35] C. Bambi, L. Modesto, S. Porey and L. Rachwał, JCAP 09, 033 (2017) [arXiv:1611.05582 [gr-qc]].
  * [36] C. Bambi, L. Modesto, S. Porey and L. Rachwał, Eur. Phys. J. C 78, no.2, 116 (2018) [arXiv:1711.04575 [gr-qc]].
  * [37] L. Modesto and L. Rachwal, [arXiv:1605.04173 [hep-th]].
  * [38] C. Bambi, L. Modesto and L. Rachwał, JCAP 05, 003 (2017) [arXiv:1611.00865 [gr-qc]].
  * [39] S. Deser, M. J. Duff and C. J. Isham, Nucl. Phys. B 111, 45-55 (1976)
  * [40] S. Deser and A. Schwimmer, Phys. Lett. B 309, 279-284 (1993) [arXiv:hep-th/9302047 [hep-th]].
  * [41] J. Erdmenger and H. Osborn, Nucl. Phys. B 483, 431-474 (1997) [arXiv:hep-th/9605009 [hep-th]].
  * [42] I. L. Shapiro and G. Cognola, Phys. Rev. D 51, 2775-2781 (1995) [arXiv:hep-th/9406027 [hep-th]].
  * [43] J. A. de Barros and I. L. Shapiro, Phys. Lett. B 412, 242-252 (1997) [arXiv:hep-th/9706123 [hep-th]].
  * [44] M. Asorey, E. V. Gorbar and I. L. Shapiro, Class. Quant. Grav. 21, 163-178 (2003) [arXiv:hep-th/0307187 [hep-th]].
  * [45] D. F. Carneiro, E. A. Freiras, B. Goncalves, A. G. de Lima and I. L. Shapiro, Grav. Cosmol. 10, 305-312 (2004) [arXiv:gr-qc/0412113 [gr-qc]].
  * [46] I. L. Shapiro, Class. Quant. Grav. 14, 391-406 (1997) [arXiv:hep-th/9610129 [hep-th]].
  * [47] I. L. Shapiro and J. Sola, Phys. Lett. B 530, 10-19 (2002) [arXiv:hep-ph/0104182 [hep-ph]].
  * [48] R. Balbinot, A. Fabbri and I. L. Shapiro, Nucl. Phys. B 559, 301-319 (1999) [arXiv:hep-th/9904162 [hep-th]].
  * [49] M. J. Duff, Class. Quant. Grav. 11, 1387-1404 (1994) [arXiv:hep-th/9308075 [hep-th]].
  * [50] S. D. Odintsov and I. L. Shapiro, Class. Quant. Grav. 8, L57-L60 (1991)
  * [51] I. Antoniadis and E. Mottola, Phys. Rev. D 45, 2013-2025 (1992)
  * [52] S. Deser, Helv. Phys. Acta 69, no.4, 570-581 (1996) [arXiv:hep-th/9609138 [hep-th]].
  * [53] G. A. Vilkovisky, Class. Quant. Grav. 9, 895-903 (1992)
  * [54] A. O. Barvinsky, Y. V. Gusev, G. A. Vilkovisky and V. V. Zhytnikov, J. Math. Phys. 35, 3525-3542 (1994) [arXiv:gr-qc/9404061 [gr-qc]].
  * [55] A. O. Barvinsky, Y. V. Gusev, G. A. Vilkovisky and V. V. Zhytnikov, Nucl. Phys. B 439, 561-582 (1995) [arXiv:hep-th/9404187 [hep-th]].
  * [56] A. O. Barvinsky, A. G. Mirzabekian and V. V. Zhytnikov, [arXiv:gr-qc/9510037 [gr-qc]].
  * [57] R. J. Riegert, Phys. Lett. B 134, 56-60 (1984)
  * [58] E. S. Fradkin and A. A. Tseytlin, Phys. Lett. B 134, 187 (1984)
  * [59] I. L. Buchbinder, S. D. Odintsov and I. L. Shapiro, Phys. Lett. B 162, 92-96 (1985)
  * [60] R. Penrose, Proc. Roy. Soc. Lond. A 284, 159 (1965)
  * [61] R. Penrose, Annals Phys. 10, 171-201 (1960)
  * [62] B. Ørsted, J. Funct. Anal. 44, 1-23 (1981).
  * [63] E. S. Fradkin and A. A. Tseytlin, Phys. Lett. B 110, 117-122 (1982)
  * [64] S. M. Paneitz, Symmetry, Integrability and Geometry: Methods and Applications SIGMA 4, 036 (2008), preprint arxiv.org/abs/0803.4331v1 [arXiv: 0803.4331v1 [math.DG]].
  * [65] M. G. Eastwood and M. Singer, Phys. Lett. A 107, 73-74 (1985)
  * [66] T. P. Branson, B. Ørsted, in Proceedings of the American Mathematical Society, 113 (3), 669–682 (1991)
  * [67] S-Y. A. Chang, in M. Christ; C. Kenig; C. Sadorsky (eds.), Harmonic Analysis and Partial Differential Equations; Essays in honor of Alberto P. Calderon, Chicago Lectures in Mathematics, 127–150 (1999)
  * [68] H. Kengo, In Complex Geometry (Osaka 1990) Lecture Notes in Pure and Applied Math, New York, Marcel Dekker 143 67–76, (1993)
  * [69] S. Park and R. P. Woodard, Gen. Rel. Grav. 42, 2765-2783 (2010) [arXiv:0910.4756 [gr-qc]].
  * [70] I. L. Shapiro and A. G. Zheksenaev, Phys. Lett. B 324, 286-292 (1994)
  * [71] K. A. Meissner and H. Nicolai, Phys. Lett. B 660, 260-266 (2008) [arXiv:0710.2840 [hep-th]].
  * [72] M. Asorey, E. V. Gorbar and I. L. Shapiro, Braz. J. Phys. 35, 1086-1090 (2005)
  * [73] P. H. Ginsparg, [arXiv:hep-th/9108028 [hep-th]].
  * [74] A. W. W. Ludwig, In *Trieste 1992, Proceedings, Low-dimensional quantum field theories for condensed matter physicists*, 389-455, (1992)
  * [75] V. Riva and J. L. Cardy, Phys. Lett. B 622, 339-342 (2005) [arXiv:hep-th/0504197 [hep-th]].
  * [76] A. Iorio, L. O’Raifeartaigh, I. Sachs and C. Wiesendanger, Nucl. Phys. B 495, 433-450 (1997) [arXiv:hep-th/9607110 [hep-th]].
  * [77] L. Modesto and L. Rachwal, Nucl. Phys. B 889, 228-248 (2014) [arXiv:1407.8036 [hep-th]].
  * [78] L. Modesto and L. Rachwał, Nucl. Phys. B 900, 147-169 (2015) [arXiv:1503.00261 [hep-th]].
  * [79] L. Modesto, M. Piva and L. Rachwal, Phys. Rev. D 94, no.2, 025021 (2016) [arXiv:1506.06227 [hep-th]].
  * [80] A. S. Koshelev, K. Sravan Kumar, L. Modesto and L. Rachwał, Phys. Rev. D 98, no.4, 046007 (2018) [arXiv:1710.07759 [hep-th]].
  * [81] S. Weinberg, In S. W. Hawking; W. Israel (eds.). General Relativity: An Einstein centenary survey. Cambridge University Press, 790–831 (1979)
  * [82] S. Weinberg, In Zichichi A. (eds) Understanding the Fundamental Constituents of Matter. The Subnuclear Series, vol 14. Springer, Boston, MA (Erice Subnucl.Phys.1976:1)
  * [83] K. Rummukainen, M. Tsypin, K. Kajantie, M. Laine and M. E. Shaposhnikov, Nucl. Phys. B 532, 283-314 (1998) [arXiv:hep-lat/9805013 [hep-lat]].
  * [84] R. Kallosh and A. Linde, JCAP 07, 002 (2013) [arXiv:1306.5220 [hep-th]].
  * [85] M. Henkel. In: Conformal Invariance and Critical Phenomena. Texts and Monographs in Physics. Springer, Berlin, Heidelberg. 1999
  * [86] R. Utiyama and B. S. DeWitt, J. Math. Phys. 3, 608-618 (1962)
  * [87] K. j. Hamada, Prog. Theor. Phys. 105, 673-690 (2001) [arXiv:hep-th/0012053 [hep-th]].
  * [88] K. j. Hamada, Prog. Theor. Phys. 108, 399-433 (2002) [arXiv:hep-th/0203250 [hep-th]].
  * [89] K. j. Hamada and S. Horata, Prog. Theor. Phys. 110, 1169-1210 (2004) [arXiv:hep-th/0307008 [hep-th]].
  * [90] A. Iorio and G. Lambiase, Phys. Rev. D 90, no.2, 025006 (2014) [arXiv:1308.0265 [hep-th]].
  * [91] A. Iorio, Annals Phys. 326, 1334-1353 (2011) [arXiv:1007.5012 [hep-th]].
  * [92] A. Iorio, Eur. Phys. J. Plus 127, 156 (2012) [arXiv:1207.6929 [hep-th]].
  * [93] R. M. Wald, General Relativity, The University of Chicago Press, Chicago, 1984
  * [94] F. M. Ferreira, I. L. Shapiro and P. M. Teixeira, Eur. Phys. J. Plus 131, no.5, 164 (2016) [arXiv:1507.03620 [gr-qc]].
  * [95] F. M. Ferreira and I. L. Shapiro, Phys. Lett. B 772, 174-178 (2017) [arXiv:1702.06892 [hep-th]].
  * [96] F. M. Ferreira and I. L. Shapiro, Phys. Rev. D 99, no.6, 064032 (2019) [arXiv:1812.01140 [gr-qc]].
  * [97] L. Rachwał, Higher Conformal Couplings of Scalar Fields from Weyl Gauge Theory, In Proceedings of the LIII Cracow School of Theoretical Physics 2013 in Zakopane, Poland, 2013.
  * [98] S. Faci, Class. Quant. Grav. 30, 115005 (2013) [arXiv:1212.2599 [hep-th]].
  * [99] S. Faci, EPL 101, no.3, 31002 (2013) [arXiv:1206.4362 [hep-th]].
  * [100] F. Bastianelli, S. Frolov and A. A. Tseytlin, JHEP 02, 013 (2000) doi:10.1088/1126-6708/2000/02/013 [arXiv:hep-th/0001041 [hep-th]].
  * [101] M. Asorey, L. Rachwał, I. L. Shapiro and W. Cesar e Silva, Eur. Phys. J. Plus 136, no.10, 1043 (2021) [arXiv:2107.13125 [hep-th]].
  * [102] M. Safari, A. Stergiou, G. P. Vacca and O. Zanusso, [arXiv:2112.01084 [hep-th]].
  * [103] G. Anastasiou and R. Olea, Phys. Rev. D 94, no.8, 086008 (2016) [arXiv:1608.07826 [hep-th]].
  * [104] G. Anastasiou, I. J. Araya and R. Olea, JHEP 01, 134 (2021) [arXiv:2010.15146 [hep-th]].
  * [105] V. Wunsch , Some New Conformal Covariants, Zeitschrift für Analysis und ihre Anwendungen Vol. 19, No. 2, pp. 339–357 (2000)
  * [106] Volkmar Wunsch, On conformally invariant differential operators, Math. Nachr. 129 (1986), 269–281. MR 864639,
  * [107] Thomas P. Branson, Differential operators canonically associated to a conformal structure, Math. Scand. 57 (1985), no. 2, 293–345. MR 832360,
  * [108] Thomas Branson, Peter Gilkey, and Juha Pohjanpelto, Invariants of locally conformally flat manifolds, Trans. Amer. Math. Soc. 347 (1995), no. 3, 939–953. MR 1282884,
  * [109] Thomas P. Branson and Bent Orsted, Conformal indices of Riemannian manifolds, Compositio Math. 60 (1986), no. 3, 261–293. MR 869104
  * [110] Thomas P. Branson and Bent Orsted, Conformal geometry and global invariants, Differential Geom. Appl. 1 (1991), no. 3, 279–308. MR 1244447,
  * [111] H. D. Fegan, Conformally invariant first order differential operators, Quart. J. Math. Oxford Ser. (2) 27 (1976), no. 107, 371–378. MR 482879,
  * [112] C. Robin Graham, Ralph Jenne, Lionel J. Mason, and George A. J. Sparling, Conformally invariant powers of the Laplacian. I. Existence, J. London Math. Soc. (2) 46 (1992), no. 3, 557–565. MR 1190438,
  * [113] Graham, C.R. (1992), Conformally Invariant Powers of the Laplacian, II: Nonexistence. Journal of the London Mathematical Society, s2-46: 566-576.
  * [114] V. P. Gusynin and V. V. Romankov, Sov. J. Nucl. Phys. 46, 1097-1099 (1987) ITF-86-126P.
  * [115] Branson T.P., Conformally covariant equations on differential forms, Comm. Partial Differential Equations 7 (1982), 393–431
  * [116] H. Osborn and A. Stergiou, JHEP 04, 157 (2015) [arXiv:1501.01308 [hep-th]].
  * [117] T. Arakelian, D. R. Karakhanian, R. P. Manvelyan and R. L. Mkrtchian, Phys. Lett. B 353, 52-56 (1995)
  * [118] J. Erdmenger and H. Osborn, Class. Quant. Grav. 15, 273-280 (1998) [arXiv:gr-qc/9708040 [gr-qc]].
  * [119] L. Bonora, P. Pasti and M. Bregola, Class. Quant. Grav. 3, 635 (1986)
  * [120] A. O. Barvinsky, G. H. S. Camargo, A. E. Kalugin, N. Ohta and I. L. Shapiro, Phys. Rev. D 108, no.8, 086018 (2023) [arXiv:2308.05251 [hep-th]].
  * [121] G. H. S. Camargo and I. L. Shapiro, Phys. Rev. D 106, no.4, 045004 (2022) [arXiv:2206.02839 [hep-th]].
  * [122] M. Asorey, W. C. Silva, I. L. Shapiro and P. R. B. d. Vale, Eur. Phys. J. C 83, no.2, 157 (2023) [arXiv:2202.00154 [hep-th]].
  * [123] W. Cesar e Silva and I. L. Shapiro, Symmetry 13, no.1, 50 (2020) [arXiv:2012.10554 [hep-th]].
  * [124] Gover, A. Rod, Peterson, Lawrence J., Conformally Invariant Powers of the Laplacian, Q-Curvature, and Tractor Calculus, Communications in Mathematical Physics, 339-378, vol. 235, no. 2 (2003).
  * [125] S. Curry and A. R. Gover, [arXiv:1412.7559 [math.DG]].


#  Low-frequency excess vibrational modes in two-dimensional glasses

Lijin Wang1,∗ Grzegorz Szamel2 Elijah Flenner2,† 1School of Physics and
Materials Science, Anhui University, Hefei 230601, P. R. China  2Department of
Chemistry, Colorado State University, Fort Collins, Colorado 80523, USA

###### Abstract

Glasses possess more low-frequency vibrational modes than predicted by Debye
theory. These excess modes are crucial for the understanding of the low
temperature thermal and mechanical properties of glasses, which differ from
those of crystalline solids. Recent simulational studies suggest that the
density of the excess modes scales with their frequency $\omega$ as
$\omega^{4}$ in two and higher dimensions. Here, we present extensive
numerical studies of two-dimensional model glass formers over a large range of
glass stabilities. We find that the density of the excess modes follows
$D_{\text{exc}}(\omega)\sim\omega^{2}$ up to around the boson peak, regardless
of the glass stability. The stability dependence of the overall scale of
$D_{\text{exc}}(\omega)$ correlates with the stability dependence of low-
frequency sound attenuation. However, we also find that in small systems,
where the first sound mode is pushed to higher frequencies, at frequencies
below the first sound mode there are excess modes with a system size
independent density of states that scales as $\omega^{3}$.

Low-temperature glasses exhibit thermal and mechanical properties Zeller1971 ;
Anderson1972 ; Zaitlin1975 ; Phillips ; Pohl2002 ; Ozawa_PNAS2018 ;
Ketkaew_2018NC  that distinguish them from low temperature crystalline solids.
The low frequency vibrational modes in crystalline solids are plane waves.
Their density of states is well described by Debye theory and scales with
frequency $\omega$ as $\omega^{d-1}$ where $d$ is the spatial dimension. For
glasses, there are additional low-frequency modes that result in a peak in the
reduced total density of states $D(\omega)/\omega^{d-1}$ in $d$ spatial
dimensions, which is referred to as the boson peak parisi_nature ; Inoue1991 ;
Buchenau2007 ; tanaka_nm2008 . Understanding the nature of the additional
modes provides insight into the the physics behind the anomalous properties of
glasses and super-cooled liquids Widmer-Cooper_np ; rottler_prx ; chen_prl ;
manning_prl ; Zylberg_PNAS2017 ; Flenner_SM2020 ; Lerner2019JCP ;
Wang2019SMattenuation ; Xu2010EPL .

Mean field theory meanfield1 ; meanfield2  predicts that the density of the
low frequency excess modes $D_{\text{exc}}(\omega)$ grows as $\omega^{\beta}$
with ${\beta}=2$, while several phenomenological models Buchenau1991 ;
Schober1996 ; Gurevich_prb2003 ; Gurarie_prb2003 ; Kumar2021arxiv  predict
${\beta}=3$ or $4$. Fluctuating elasticity theory Schrimacher_prl2007
predicts that $D_{\text{exc}}(\omega)$ scales as $\omega^{d+1}$, which is the
same scaling that the same theory predicts for the sound attenuation
coefficient. An analysis based on a fold stability predicts
$D_{\text{exc}}(\omega)$ scales as $\omega^{3}$ in glasses approaching
marginal stability Xu2017prl . Other recent theories predict the density of
the excess modes scales as $\omega^{4}$ Ikeda2019pre ; Stanifer2018pre ;
Ji2019pre ; Bouchbinder2020 .

Simulations are a useful tool to examine low-frequency vibrational modes since
characteristics of each mode can be investigated, but studying finite systems
presents some difficulties. One is that the plane-wave-like modes occur around
discrete frequencies, which can be approximated using Debye theory, but the
density of states does not exactly follow the Debye prediction. One has to be
careful in the calculation of the density of states due to this discrete
nature of the plane-wave-like modes Bouchbinder2018 . Simulating two-
dimensional (2D) glasses adds another potential difficulty since Mermin-Wagner
Flenner2015 ; Vivek2017 ; Illing2017 ; Flenner2018  fluctuations lead to
pronounced finite size effects in some static and dynamic properties of two-
dimensional solids.

With increasing system size the glass behaves increasingly as a continuous
elastic solid, and it is expected that there are plane-wave-like modes similar
to those of Debye theory. This observation has lead researchers to distinguish
between plane-wave-like modes and additional modes. One simple way to do this
is to use the participation ratio, which is a measure of how many particles
significantly participate in the mode Mazzacurati1996 ; ikeda_pnas . A more
sophisticated approach is to introduce an order parameter that quantifies the
similarity between a low-frequency mode in an amorphous solid and a plane wave
ikeda_pnas . Although these two methods are naturally suited for large
systems, in principle they can be used for systems of any size. An alternative
approach to distinguish modes of different nature is to study small systems in
which the first plane-wave-like mode (which is found at a frequency close to
the one predicted by Debye theory) is pushed to higher frequencies
lerner_prl2016 . The low-frequency modes found in these small systems are
postulated to be the modes in excess of the Debye prediction.

Mizuno, Shiba, and Ikeda ikeda_pnas  used the participation ratio and an order
parameter to separate modes into extended and excess modes in large, over one-
million particle, two- (2D) and three-dimensional (3D) systems. In both
dimensions they found that the density of the modes with large participation
ratio obeyed Debye scaling. In 3D they found that the density of the excess
modes, which they determined are quasi-localized, scales as
$D_{loc}(\omega)\sim\omega^{4}$.

The scaling of the density of excess modes Mizuno et al. found in 3D agrees
with the scaling observed previously by studying small systems lerner_prl2016
. Subsequent work by Wang et al. Wang2019NC  confirmed that the picture
observed by Mizuno et al. in 3D is also found in glasses of various
stabilities, up to the stability of laboratory glasses. Numerical simulations
have demonstrated overwhelmingly the universality of $D(\omega)\sim\omega^{4}$
scaling in 3D model glass formers, irrespective of glass preparation
protocols, or interaction potentials spinglassPRL ; lerner_prl2016 ;
ikeda_pnas ; ikeda_pre2018 ; Shimada2018pre ; lerner2018jcp ; Wang2019NC ;
lerner2020pre ; Angelani2018 ; Lerner2020pnas ; Bonfanti2020prl ; Lopez2020prl
; Das2021prl ; Shimada2020pre ; Ji2020pre ; LernerPRE2017 . In their studies
of large 2D systems, Mizuno et al. found very few low frequency modes with
small participation ratio or with small values of plane-wave order parameter.
However, the work of Kapteijns, Bouchbinder, and Lerner Lerner2Dprl ;
Krishnan2D , who studied small systems, found modes below the first plane-
wave-like mode with density scaling as $\omega^{4}$. Notably, unlike in higher
dimensions, Kapteijns et al. found that the pre-factor for the $\omega^{4}$
scaling grew with system size as $[\log{N}]^{5/2}$ in 2D. For the much larger
system studied by Mizuno et al., it might be expected that there would be a
discernible increase in the density of states over the Debye spectrum, but the
logarithmic increase with system size would make the increase modest.

In this work, we present results for the density of excess modes,
$D_{\text{exc}}(\omega)$, in 2D model glass formers with different interaction
potentials and stability. We used two ways to calculate
$D_{\text{exc}}(\omega)$. The first method is to subtract off the infinite
size system Debye prediction. Except for very low frequencies, this should
allow one to examine how the density of modes in excess of the Debye
prediction changes with frequency, but the discrete nature of the spectrum at
low frequencies makes it hard to determine the low-frequency growth of
$D_{\text{exc}}(\omega)$. We used this concept and found that
$D_{\text{exc}}(\omega)\sim\omega^{2}$ in 2D, which differs from previous
observations. Importantly, $D_{\text{exc}}(\omega)$ is correlated with the low
frequency scaling of sound attenuation, which resembles the correlation we
found between the density of the excess modes and the sound attenuation in 3D
Wang2019SMattenuation . To make a more direct connection with previous results
we also studied small systems. Unlike previous work, we found a system size
and model independent $\omega^{3}$ scaling of modes far below the first mode
predicted by Debye theory. However, these low-frequency modes are very rare
even in poorly annealed systems and are absent in our very stable systems.

We performed extensive simulation studies of four 2D model glass formers with
spherically symmetric interactions: (I) a system of polydisperse particles
interacting via an inverse power law potential $\propto{r}^{-n}$ ($r$ is the
distance of two particles) with $n=12$ (IPL-12) berthier_2020nc ; (II) a
bidisperse system with inverse power law potential with $n=10$ (IPL-10)
lerner2018jcp ; (III) a bidisperse system with Lennard-Jones potential (LJ)
BruningLJ ; (IV) a bidisperse system with a harmonic potential (HARM)
OHern2003 . Details regarding the four models can be found in the Supplemental
Material SM .

We created zero-temperature ($T=0$) glasses by quenching instantaneously
equilibrated liquid configurations at parent temperatures $T_{p}$ to $T=0$
using the fast inertial relaxation engine fire . Equilibrated liquids at very
high parent temperatures were obtained by performing molecular dynamic
simulations using LAMMPS lammps . As discussed elsewhere, glasses obtained
using this method are not very stable. To generate stable glasses for the
IPL-12 system, we employed the swap Monte Carlo method Grigera2001 ;
berthier_prx2017 ; berthier_prl2016  to prepare equilibrated supercooled
liquids at low $T_{p}$, down to around $37\%T_{g}$, where $T_{g}\approx 0.082$
is the estimated experimental glass temperature berthier_2020nc .

The normal modes of $T=0$ glasses were obtained by diagonalizing the Hessian
matrix using ARPACK arpack  and Intel Math Kernel Library mkl . The density of
states is given by
$D(\omega)=\frac{1}{2N-2}\sum_{l=1}^{2N-2}\delta(\omega-\omega_{l})$ with
$\omega_{l}$ the frequency of mode $l$ and $N$ the number of particles. In
glasses, there are no pure plane-wave modes and the frequencies of the plane-
wave-like modes are generally clustered around the Debye predictions
Lernerarxiv2021 ; Bouchbinder2018 . Since Debye theory predicts discrete modes
in finite systems, one has to be careful when calculating the density of
states since that calculation requires a division by the bin size
$\delta\omega$. If $\delta\omega$ is not chosen correctly, the density is
inaccurate. The calculation of the cumulative density of states
$I(\omega)=\int_{0}^{\omega}D(\omega^{\prime})d\omega^{\prime}$ does not
suffer from this issue since it amounts to counting the number of states up to
$\omega$ and dividing by the total number of states. For this reason, we focus
on the cumulative density of states $I(\omega)$.

To obtain the excess modes Mizuno and coworkers ikeda_pnas  defined a
threshold of the participation ratio $P_{c}=0.01$ to divide plane-wave-like
modes and quasi-localized modes. They concluded that there are few to no low-
frequency quasi-localized modes in poorly annealed 2D glasses ikeda_pnas .
Additionally, the Debye theory accurately predicted the low-frequency density
of states, but there was still a boson peak at higher frequencies. We
attempted to use the participation ratio to separate the modes, but we found
that the scaling behavior of the excess modes in 2D stable IPL-12 model
glasses depends strongly on the choice of $P_{c}$. This strong dependence
makes it impossible to determine the scaling of $I(\omega)$ using the
participation ratio. Therefore, we utilized a different procedure by
subtracting from the cumulative density of states the Debye prediction
Schrimacher_prl2007

$I_{\text{exc}}(\omega)=I(\omega)-I_{D}(\omega),$ (1)

where $I(\omega)$ is the cumulative density of states of all modes and
$I_{D}(\omega)$ is the Debye prediction kittel ,
$I_{D}(\omega)=A_{D}\omega^{d}/d$ with Debye level $A_{D}$ determined
independently from mechanical moduli kittel . We note that this procedure does
not take into account that the mode frequencies are discrete for finite
systems.

Figure 1:  (color online) (a) Cumulative density of states of excess modes
$I_{\text{exc}}(\omega)=I(\omega)-I_{D}(\omega)$ at $T_{p}=0.400$ (circles),
$0.085$ (diamonds), and $0.030$ (squares) in $N=20000$ system in the 2D IPL-12
model. $I(\omega)$ is the total cumulative density of states while
$I_{D}(\omega)=1/2A_{D}\omega^{2}$ with $A_{D}$ Debye level. The red lines are
fits to $I_{\text{exc}}(\omega)=1/3A_{2}\omega^{3}$ while vertical lines
indicate boson peak frequencies. (b) $T_{p}$ dependence of $A_{2}$, with the
estimated glass transition temperature $T_{g}$ indicated for reference. Glass
stability increases with decreasing $T_{p}$. (Inset) $A_{2}$ against the pre-
factor $B_{3}$ in $\Gamma(\omega)=B_{3}\omega^{3}$ with $\Gamma$ the
transverse sound attenuation coefficient. Details for the calculation of
$\Gamma$ can be found in Ref.Wang2019SMattenuation . The blue line indicates a
fit to $A_{2}\sim B_{3}^{\gamma}$ with $\gamma\approx 2/3$.

Figure 1(a) shows excess cumulative density of states $I_{\text{exc}}(\omega)$
for our 2D IPL-12 model glasses for $N=20000$ at different parent temperatures
$T_{p}$. The glass stability increases with decreasing $T_{p}$ Wang2019NC ;
berthier_2020nc ; Lerner2020pnas . We use parent temperatures ranging from
$T_{p}=0.400$, which is above the onset temperature of slow dynamics
$T_{o}=0.250$, down to $T_{p}=0.030$, which is below the estimated laboratory
glass temperature $T_{g}=0.082$ berthier_2020nc . For the lowest frequencies
where we can clearly estimate a power law, we find that
$I_{\text{exc}}(\omega)$ scales as $\omega^{3}$, $I_{\text{exc}}(\omega)\simeq
A_{2}\omega^{3}/3$, which suggests that $D_{\text{exc}}(\omega)\simeq
A_{2}\omega^{2}$. We find that this scaling continues up to around the Ioffe-
Regal limit, which is around the boson peak frequency tanaka_nm2008 ,
irrespective of the glass’s stability, see Fig. 1(b). We find that the
coefficient quantifying the magnitude of the excess modes density, $A_{2}$, is
stability-dependent. $A_{2}$ is nearly constant for larger $T_{p}$, but
decreases by a factor of 13 for our lowest $T_{p}$. This indicates that there
are fewer excess modes for increasingly stable 2D glasses, which is consistent
with observations for 3D glasses Wang2019NC ; Lerner2020pnas ; Ji2020pre .

Previous work Wang2019SMattenuation ; Lerner2019JCP ; Schrimacher_prl2007
found a connection between sound attenuation and density of states of excess
low-frequency modes. Inspired by this work we examined whether
$I_{\text{exc}}(\omega)$ is related to sound attenuation in 2D. The frequency
dependence of the transverse sound attenuation coefficient $\Gamma(\omega)$ in
2D glasses follows the Rayleigh scattering scaling as
$\Gamma(\omega)=B_{3}\omega^{3}$ and $B_{3}$ decreases with increasing glass
stability Lerner2019JCP . Here, we study the relation between $B_{3}$ and
$A_{2}$, see the inset to Fig. 1b. We find $A_{2}\sim B_{3}^{\gamma}$ with
$\gamma\approx 2/3$, and thus we establish that in 2D the excess modes density
is related to sound attenuation, which is consistent with our result for 3D
glasses.

Since the method introduced in this Letter is different from methods used
before by us and others, we checked what results it produces if used for 3D
glasses where it has been firmly established that
$I_{\text{exc}}(\omega)\sim\omega^{5}$. In Fig. 2 we show
$I_{\text{exc}}(\omega)$ in 3D IPL-12 glasses for two stabilities, a very
poorly annealed glass with $T_{p}=0.200$ and a very stable glass with
$T_{p}=0.062$. These are the same glasses examined in Ref. Wang2019NC . We
find that $I_{\text{exc}}(\omega)\sim\omega^{5}$ up to a frequency close to
the boson peak for both glasses, which indicates the resulting scaling of
$I_{\text{exc}}(\omega)$ determined using Eq. 1 is consistent with that of
$I_{\text{exc}}(\omega)$ calculated with previously used procedures Wang2019NC
; ikeda_pnas ; ikeda_pre2018 . Additionally, these results suggest that in 3D
the end of the $\omega^{5}$ scaling of $I_{\text{exc}}(\omega)$, $\omega_{g}$,
is around the boson peak frequency. We note that our procedure cannot be used
for frequencies below the lowest frequency plane wave mode predicted by the
Debye theory. More importantly, it will only reveal the proper scaling if
there is a near continuum of modes Bouchbinder2018 .

Figure 2:  (color online) Cumulative density of states of excess modes
$I_{\text{exc}}(\omega)=I(\omega)-I_{D}(\omega)$ in the 3D IPL-12 model
glasses, with the boson peak frequency indicated. $I(\omega)$ is the total
cumulative density of states and $I_{D}(\omega)=A_{D}\omega^{3}/3$ with
$A_{D}$ being the Debye level. Red filled squares and circles represent
$I(\omega)$ at frequencies below around the first Debye frequency at
$T_{p}=0.200$ and $T_{p}=0.062$, respectively, in the $N=48000$ system. The
red lines correspond to $I_{\text{exc}}(\omega)\sim\omega^{5}$.

Debye theory predicts that the lowest normal mode frequency increases with
decreasing $N$. It has been argued that the scaling of the excess modes could
be obtained from the low frequency density of states for small systems since
the frequency of the plane-wave-like modes is pushed to high frequencies.
Using this reasoning, $I_{\text{exc}}(\omega)$ can be determined by
calculating the total cumulative density of states $I(\omega)$ for frequencies
below the first predicted Debye mode frequency. One may expect
$I(\omega)=I_{\text{exc}}(\omega)$ for low frequencies if the excess modes
were independent of the plane–wave like modes, which we found for 3D glasses
with different stabilities, see Fig. 2. The low-frequency tail of $I(\omega)$
is well described by a power law $I_{\text{exc}}(\omega)\sim\omega^{5}$ for
$T_{p}=0.200$ and $T_{p}=0.062$. However, in 2D glasses much below the first
Debye frequency we find $I(\omega)\sim\omega^{4}$, which suggests that
$D_{exc}(\omega)\sim\omega^{3}$. Previous studies reported
$D_{\text{exc}}(\omega)\sim\omega^{4}$ Lerner2Dprl ; Krishnan2D , which would
imply that $I_{\text{exc}}(\omega)\sim\omega^{5}$. To make sure this
observation is model independent, we calculated $I(\omega)$ for small systems
at frequencies much lower than the first Debye frequency in different models
of glass formers.

Figure 3:  (color online) The total cumulative density of states $I(\omega)$
for $N=3000$ system at $T_{p}=0.400$ with different ensemble sizes in the 2D
IPL-12 model. The solid and dashed lines represent power laws of $\omega^{4}$
and $\omega^{5}$ respectively. For reference, the lowest Debye mode frequency
is about $0.261$. (Inset) The same data plotted as $I(\omega)/\omega^{4}$ vs.
$\omega$.

In Fig. 3 we show the total cumulative density of states $I(\omega)$ for
$N=3000$ system in the 2D IPL-12 model. There is a range of frequencies below
the lowest Debye mode frequency ($\approx 0.261$) that is well described by
$I(\omega)\sim\omega^{4}$. To check the quartic scaling, we examined
$I(\omega)/\omega^{4}$, which is shown in the inset to Fig. 3, and we find
that there is a low-frequency plateau. Previous results suggest that
$I(\omega)\sim\omega^{5}$ at frequencies much below the first Debye mode in 2D
glasses Lerner2Dprl ; Krishnan2D , and that $I(\omega)$ should be system size
dependent. However, we find that the $\omega^{5}$ scaling is only valid for an
intermediate-frequency regime below the peak of $I(\omega)$, and appears to be
only a transition between the low frequency scaling and the change of the
scaling due to the emergence of plane-wave-like modes. We also find that
$I(\omega)$ is system size independent at frequencies much lower than the
first frequency predicted by Debye theory.

There are very few modes that contribute to the low frequency $\omega^{4}$
scaling of $I(\omega)$ for our least stable 2D glass. For example, on average
there is only one mode lying in the low-frequency $I(\omega)\sim\omega^{4}$
regime every one hundred configurations for the $N=3000$ system. The lower the
frequency we want to examine, the larger the ensemble size $N_{En}$ (number of
configurations) we need. However, we do not observe $N_{En}$ dependence of the
quartic scaling regime in $N=3000$ system when $N_{En}$ ranges from around 0.1
million to our maximum 0.71 million examined, see Fig. 3. The same conclusion
can also be drawn in our study of the $N=1000$ system where the maximum
$N_{En}$ is around 2.2 million SM . In addition, we checked that the
previously reported $\omega^{5}$ scaling in some systems is due to ensemble
size not being large enough, which hinders the observation of the $\omega^{4}$
scaling at much lower frequencies. We also find that the $\omega^{5}$ scaling
regime vanishes for very small systems, Fig. 4. Since the number of these low-
frequency modes decreases with increasing stability, we could not examine the
stability dependence of these modes. However, we do not exclude the
possibility that thermal relaxation can change the scaling of these low-
frequency modes LernerPRE2017 .

We find that this low-frequency quartic scaling of $I(\omega)$ in 2D is
universal, i.e. it does not depend on the model glass former, see Fig. 4. In
addition, there are common features shared by each model. First, the pre-
factor of the quartic scaling does not depend on system size. This conclusion
is different than the conclusion of Ref. Lerner2Dprl , namely that
$I_{\text{exc}}(\omega)=A_{4}\omega^{5}/5$ and the pre-factor $A_{4}$ grows as
$(\log N)^{5/2}$, i.e. there should be more excess modes in a larger system.
Second, the low-frequency scaling $I(\omega)\sim\omega^{4}$ works up to a
larger frequency with decreasing system size. It deserves further study to
check whether the upper frequency for this scaling correlates with any known
characteristic frequency. Figures 3 and 4, therefore, demonstrate that
$I(\omega)\sim\omega^{4}$ at frequencies much lower than the first frequency
predicted by Debye theory in 2D glasses and the pre-factor for this power law
shows no system size dependence.

Figure 4:  (color online) The total cumulative density of states $I(\omega)$
for different system sizes in (a) the 2D IPL-12 model, (b) 2D IPL-10 model,
(c) 2D LJ model and (d) 2D HARM model. In each panel, the solid line
corresponds to the power law of $\omega^{4}$ and represents the low frequency
data very well for all model glass formers examined. We plotted the same data
as $I(\omega)/\omega^{4}$ against $\omega$ in the Supplemental Material SM .

We utilized two methods to examine the excess density of states in 2D glasses.
In large systems, we find evidence that the excess density of states scale as
$\omega^{2}$. We find that the pre-factor $A_{2}$ of this scaling law
correlates with the sound attenuation coefficient. However, in small systems,
where the frequencies of plane-wave-like modes are pushed higher, we find that
the modes below the lowest Debye frequency have density of states scaling as
$\omega^{3}$, with a system size-independent pre-factor. This inconsistent
behavior is not found in 3D glasses using the same analysis.

Our results leave several open questions. First, why is the scaling of excess
modes different above and below the first mode predicted by Debye theory? One
possibility is that our systems are not large enough to accurately determine
$I_{\text{exc}}(\omega)$ by subtracting off the Debye contribution at low
frequencies. However, we do find an overlapping frequency range where we find
that the excess density of states scales as $\omega^{2}$ when we subtract off
the Debye contribution and where the excess density of states scales as
$\omega^{3}$ for very small systems, see Fig. 3 of the Supplemental Material
SM . Thus, it seems that the presence of plane waves influences the scaling of
the excess density of states in 2D glasses.

There may also be a gap in the excess density of states and the $\omega^{2}$
and the $\omega^{3}$ scaling does not extend to $\omega=0$, which would be
consistent with the conclusions of Ref. ikeda_pnas . It is very difficult to
numerically test these possibilities with the current computer power since the
$\omega^{2}$ and $\omega^{3}$ scaling represents very few modes at small
$\omega$. Future related theoretical work may shed some light on this issue.

Third, why is the excess density of states different from that of 3D glasses?
Fluctuation elasticity theories predict that $D_{\text{exc}}(\omega)$ depends
on spatial dimension as $\omega^{d+1}$ Schrimacher_prl2007 . Thus, the
predicted $D_{\text{exc}}(\omega)\sim\omega^{4}$ in 3D glasses is consistent
with 3D numerical observations. The predicted scaling of
$D_{\text{exc}}(\omega)\sim\omega^{3}$ of 2D glasses is consistent with what
we find in small systems. However, it has been shown that
$D_{\text{exc}}(\omega)\sim\omega^{4}$ in 2D glasses Lerner2Dprl  in contrast
to the fluctuating elasticity theory prediction. We notice that the dynamics
of 2D and 3D glass-forming liquids was reported to be fundamentally different
Flenner2015 . It may be also interesting to probe whether there is any
connection between the difference in dynamics between 2D and 3D glass-forming
liquids and the difference in the density of states of excess modes between 2D
and 3D glasses.

Finally, it is possible that the upper frequency cutoff of the low frequency
scaling, $\omega_{g}$, is below the frequency range where we found
$\omega^{2}$ scaling of $D_{\text{exc}}(\omega)$ in large 2D glasses. If this
were the case, it would differ from our finding in 3D that $\omega_{g}$ is
around the boson peak frequency. We leave it for future work to examine what
determines $\omega_{g}$ and if it depends on dimension.

We wish to thank Andrea Ninarello for generously providing equilibrated
configurations at very low parent temperatures and E. Lerner for comments on
the manuscript. L. W. acknowledges the support from National Natural Science
Foundation of China (No. 12004001), Anhui Province (No. S020218016), Hefei
City (No. Z020132009), and Anhui University (Start-up fund). E.F. and G.S.
acknowledge the support from NSF Grant CHE-1800282.

## References

  * (1)
  * (2<EMAIL_ADDRESS>  * (3<EMAIL_ADDRESS>  * (4) R. C. Zeller and R. O. Pohl, Phys. Rev. B 4, 2029 (1971).
  * (5) P. W. Anderson, B. I. Halperin and C. M. Varma, _Philos. Mag._ , 25, 1 (1972).
  * (6) W. A. Phillips, J. Low Temp. Phys. 7, 351 (1972).
  * (7) M. P. Zaitlin and A. C. Anderson, _Phys. Rev. B: Solid State_ , 12, 4475 (1975).
  * (8) R. O. Pohl, X. Liu and E. Thompson, _Rev. Mod. Phys._ , 74, 991 (2002).
  * (9) M. Ozawa, L. Berthier , G. Biroli, A. Rosso and G. Tarjus, _Proc. Natl. Acad. Sci. U.S.A._ 115, 6656 (2018).
  * (10) J. Ketkaew, W. Chen, H. Wang, A. Datye, M. Fan, G. Pereira, U. D. Schwarz, Z. Liu, R. Yamada, W. Dmowski, M. D. Shattuck, C. S. O’Hern, T. Egami, E. Bouchbinder and J. Schroers, Nat. Commun. 9, 3271 (2018).
  * (11) T. S. Grigera, V. Martin-Mayor, G. Parisi, and P. Verrocchio, Nature (London) 422, 289 (2003).
  * (12) K. Inoue, T. Kanaya, S. Ikeda, K. Kaji, K. Shibata, M. Misawa, and Y. Kiyanagi, J. Chem. Phys. 95, 195332 (1991).
  * (13) H. Shintani and H. Tanaka, Nat. Mater. 7, 870 (2008).
  * (14) U. Buchenau, A. Wischnewski, M. Ohl, and E. Fabiani, J. Phys.: Condens. Matter 19, 205106 (2007).
  * (15) K. Chen, M. L. Manning, P. J. Yunker, W. G. Ellenbroek, Z. X. Zhang, A. J. Liu, and A. G. Yodh, Phys. Rev. Lett. 107, 108301 (2011).
  * (16) A. Widmer-Cooper, H. Perry, P. Harrowell, and D. R. Reichman, Nat. Phys. 4, 711 (2008).
  * (17) S. S. Schoenholz, A. J. Liu, R. A. Riggleman, and J. Rottler, Phys. Rev. X 4, 031014 (2014).
  * (18) M. L. Manning and A. J. Liu, Phys. Rev. Lett. 107, 108302 (2011).
  * (19) J. Zylberg, E. Lerner, E. Y. Bar-Sinai and E. Bouchbinder, Proc. Natl. Acad. Sci. U.S.A. 114, 7289 (2017).
  * (20) E. Flenner, L. Wang and G. Szamel, Soft Matter, 16, 775 (2020).
  * (21) A. Moriel, G. Kapteijns, C. Rainone, J. Zylberg, E. Lerner and E. Bouchbinder, J. Chem. Phys., 151, 104503 (2019).
  * (22) L. Wang, L. Berthier, E. Flenner, P. Guan and G. Szamel, Soft Matter, 15, 7018 (2019).
  * (23) N. Xu, V. Vitelli, A. J. Liu and S. R. Nagel, Europhys. Lett., 90, 56001 (2010).
  * (24) E. DeGiuli, A. Laversanne-Finot, G. Düring, E. Lerner, and M. Wyart, Soft matter 10, 5628 (2014).
  * (25) S. Franz, G. Parisi, P. Urbani, and F. Zamponi, Proc. Nat. Acad. Sci. USA 112, 14539 (2015).
  * (26) U. Buchenau, Yu. M. Galperin, V. L. Gurevich, and H. R. Schober, Phys. Rev. B 43, 5309 (1991).
  * (27) H. R. Schober and C. Oligschleger, Phys. Rev. B 53, 11469 (1996).
  * (28) V. L. Gurevich, D. A. Parshin, and H. R. Schober, Phys. Rev. B 67, 094203 (2003).
  * (29) V. Gurarie and J. T. Chalker, Phys. Rev. B 68, 134207 (2003).
  * (30) A. Kumar, I. Procaccia, and M. Singh, arXiv:2102.12368 (2021).
  * (31) W. Schirmacher, G. Ruocco, and T. Scopigno, Phys. Rev. Lett. 98, 025501 (2007).
  * (32) N. Xu, A. Liu, and S. R. Nagel, Phys. Rev. Lett. 119, 215502 (2017).
  * (33) H. Ikeda, Phys. Rev. E 99, 050901(R) (2019).
  * (34) E. Stanifer, P. K. Morse, A. A. Middleton, and M. L. Manning, Phys. Rev. E 98, 042908 (2018).
  * (35) E. Bouchbinder, E. Lerner, C. Rainone, P. Urbani, and F. Zamponi, arXiv:2012.11558 (2020).
  * (36) W. Ji, M. Popović, T. W. J. de Geus, E. Lerner, and M. Wyart, Phys. Rev. E 99, 023003 (2019).
  * (37) E. Bouchbinder and E. Lerner, New J. Phys. 20, 073022 (2018).
  * (38) E. Flenner and G. Szamel, Nature Commun. 6, 7392 (2015).
  * (39) S. Vivek, C.P. Kelleher, P.M. Chaikin, and E.R. Weeks, Proc. Natl. Acad. Sci. USA 114, 1850 (2017).
  * (40) B. Illing, S. Fritschi, H. Kaiser, C.L. Klix, G. Maret, and P. Keim, Proc. Natl. Acad. Sci. USA 114, 1856 (2017).
  * (41) E. Flenner and G. Szamel, Proc. Natl. Acad. Sci. USA 116, 2015 (2018).
  * (42) V. Mazzacurati, G. Ruocco, M. Sampoli, Europhys. Lett. 34, 681 (1996).
  * (43) H. Mizuno, H. Shiba, and A. Ikeda, Proc. Natl. Acad. Sci. U.S.A. 114, E9767 (2017).
  * (44) E. Lerner, G. Düring, and E. Bouchbinder, Phys. Rev. Lett. 117, 035501 (2016).
  * (45) L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel and E. Flenner, Nat. Commun. 2019, 10, 26.
  * (46) M. Baity-Jesi, V. MartÅ¡Âªn-Mayor, G. Parisi, and S. Perez-Gaviro, Phys. Rev. Lett. 115, 267205 (2015).
  * (47) M. Shimada, H. Mizuno, and A. Ikeda, Phys. Rev. E 97, 022609 (2018).
  * (48) G. Kapteijns, E. Bouchbinder, and E. Lerner, J. Chem. Phys. 148, 214502 (2018).
  * (49) E. Lerner, Phys. Rev. E 101, 032120 (2020).
  * (50) L. Angelani, M. Paoluzzi, G. Parisi and G. Ruocco, Proc. Natl. Acad. Sci. USA, 2018, 115, 8700.
  * (51) C. Rainone, E. Bouchbinder, and E. Lerner, Proc. Natl. Acad. Sci. USA, 2020, 117, 5228.
  * (52) S. Bonfanti, R. Guerra, C. Mondal, I. Procaccia, and S. Zapperi, Phys. Rev. Lett. 125, 085501 (2020).
  * (53) K. G. Lopez, D. Richard, G. Kapteijns, R. Pater, T. Vaknin, E. Bouchbinder, and E. Lerner, Phys. Rev. Lett. 125, 085502 (2020).
  * (54) M. Shimada, H. Mizuno, L. berthier, and A.Ikeda, Phys. Rev. E 101, 052906 (2020).
  * (55) P. Das and I. Procaccia, Phys. Rev. Lett. 126, 085502 (2021).
  * (56) M. Shimada, H. Mizuno, M. Wyart, and A.Ikeda, Phys. Rev. E 98, 060901(R) (2018).
  * (57) E. Lerner and E. Bouchbinder, Phys. Rev. E 96, 020104(R) (2017).
  * (58) W. Ji, T. W. J. de Geus, M. Popović, E. Agoritsas, and M. Wyart, Phys. Rev. E 102, 062110 (2020).
  * (59) G. Kapteijns, E. Bouchbinder and E. Lerner, Phys. Rev. Lett., 2018, 121, 055501.
  * (60) V. V. Krishnan, K. Ramola, and S. Karmakar, arXiv:2014.0918 (2021).
  * (61) L. Berthier, P. Charbonneau, A. Ninarello, M. Ozawa, and S. Yaida, Nat. Commun. 10, 1508 (2019).
  * (62) R. Bruning, D. A. St-Onge, S. Patterson, and W. Kob, J. Phys.: Condens. Matter 21, 035117 (2009).
  * (63) C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E 68, 011306 (2003).
  * (64) See Supplemental Material at * for simualtion details of the four interaction potential models and more supporting analysis.
  * (65) E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, Phys. Rev. Lett. 97, 170201 (2006).
  * (66) S. Plimpton, J. Comput. Phys. 117, 1 (1995).
  * (67) T. S. Grigera and G. Parisi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 63, 045102(R) (2001).
  * (68) A. Ninarello, L. Berthier, and D. Coslovich, Phys. Rev. X 7, 021039 (2017).
  * (69) L. Berthier, D. Coslovich, A. Ninarello, and M. Ozawa, Phys. Rev. Lett. 116, 238002 (2016).
  * (70) http://www.caam.rice.edu/software/ARPACK/.
  * (71) https://software.intel.com/en-us/mkl/.
  * (72) G. Kapteijns, E. Bouchbinder and E. Lerner, arXiv:2106.12613.
  * (73) C. Kittel, _Introduction to Solid State Physics_ , 7th Ed. (Wiley, New York, 1996). (Wiley, New York, 1996).

See pages - of SM-V02.pdf


# Evidence for a Compact Object in the Aftermath of the Extra-Galactic
Transient AT2018cow

Dheeraj R. Pasham1∗, Wynn C. G. Ho2, William Alston3, Ronald Remillard1,
Mason Ng1, Keith Gendreau4, Brian D. Metzger5, Diego Altamirano6,
Deepto Chakrabarty1, Andrew Fabian7, Jon Miller8, Peter Bult4,9,
Zaven Arzoumanian4, James F. Steiner10 Tod Strohmayer4,
Francesco Tombesi4,9,11,12, Jeroen Homan13, Edward M. Cackett14, Alice
Harding15

1Kavli Institute for Astrophysics and Space Research, Massachusetts Institute
of Technology,
Cambridge, MA 02139, USA
2Department of Physics and Astronomy, Haverford College, 370 Lancaster Avenue,
Haverford, PA 19041, USA
3European Space Agency (ESA), European Space Astronomy Centre (ESAC),
Villanueva de la Cañada, Madrid, E-28691, Spain
4NASA Goddard Space Flight Center, Greenbelt, MD, USA
5Columbia Astrophysics Laboratory, Columbia University, New York, NY 10027,
USA
6Department of Physics and astronomy, University of Southampton,
University Road, Southampton, SO17 1BJ, UK
7Institute of Astronomy, University of Cambridge, UK
8Department of Astronomy, 311 West Hall, 1085 South University Ave.,
Ann Arbor, MI 48109-1107
9Department of astronomy, University of Maryland College Park, MD, 20742
10Center for Astrophysics, Harvard & Smithsonian, Cambridge, MA, 02138
11Department of Physics, Tor Vergata University of Rome, Via della Ricerca
Scientifica 1,
00133 Rome, Italy
12INAF – Astronomical Observatory of Rome, Via Frascati 33, 00040 Monte Porzio
Catone, Italy
13Eureka Scientific, Inc., Oakland, CA 94602, USA
14Department of Physics and Astronomy, Wayne State University, 666 W Hancock,
Detroit, MI 48201
15Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545

∗To whom correspondence should be addressed; E-mail<EMAIL_ADDRESS>
> The brightest Fast Blue Optical Transients (FBOTs) are mysterious
> extragalactic explosions that may represent a new class of astrophysical
> phenomena (?). Their fast time to maximum brightness of less than a week and
> decline over several months and atypical optical spectra and evolution are
> difficult to explain within the context of core-collapse of massive stars
> which are powered by radioactive decay of Nickel-56 and evolve more slowly
> (?, ?). AT2018cow (at redshift of 0.014) is an extreme FBOT in terms of
> rapid evolution and high luminosities (?, ?, ?, ?). Here we present evidence
> for a high-amplitude quasi-periodic oscillation (QPO) of AT2018cow’s soft
> X-rays with a frequency of 224 Hz (at 3.7$\sigma$ significance level or
> false alarm probability of 0.02%) and fractional root-mean-squared amplitude
> of $>$30%. This signal is found in the average power density spectrum taken
> over the entire 60-day outburst and suggests a highly persistent signal that
> lasts for a billion cycles. The high frequency (rapid timescale) of 224 Hz
> (4.4 ms) argues for a compact object in AT2018cow, which can be a neutron
> star or black hole with a mass less than 850 solar masses. If the QPO is the
> spin period of a neutron star, we can set limits on the star’s magnetic
> field strength. Our work highlights a new way of using high time-resolution
> X-ray observations to study FBOTs.

High-cadence sky surveys that can scan the same portions of the sky multiple
times per night have uncovered fast-evolving optical transients (?, ?, ?).
These ”fast” transients rise to their peak brightness within $\lesssim$ 10
days and fade away within a month or two (e.g., (?)). They are spatially
coincident with external galaxies but are offset from their nuclei (e.g., (?,
?)). Their optical spectra are often blue with occasional presence of Hydrogen
and Helium features. The peak luminosities of FBOTs range from the faint end
of core-collapse supernovae to the bright end of superluminous supernovae (see
Fig. 1 of (?) and Fig. 1 of (?)). A recent study (?) has shown that the
majority of FBOTs are extreme cases of core-collapse supernovae. However, a
subset of high luminosity FBOTs with peak bolometric luminosities $\sim$1044
erg s-1 cannot be explained as an extension of properties of core-collapse
supernovae (?, ?, ?). Thus, several alternate mechanisms mentioned above have
been proposed to explain the properties of luminous FBOTs. These include
emission from the interaction of the supernova shockwave with a dense
circumstellar medium (?, ?, ?), injection of energy from spin down of a young
magnetar formed either in a core-collapse supernova or a binary neutron star
merger (?, ?), accretion onto a newly formed compact object in a failed
supernova (?), mergers of binary white dwarfs (?), and intermediate-mass black
holes (IMBHs: a few$\times$104-5 $M_{\odot}$) tidally disrupting stars (?, ?).

Prior to June 2018, the majority of FBOTs were first identified in archival
images (e.g., (?, ?, ?, ?)). AT2018cow, discovered by the ATLAS sky survey (?)
in a galaxy at a distance of $\approx$ 60 Mpc (?), was discovered in real
time. Its brightness rise of more than 5.7 magnitudes in just 4 days (see Fig.
1 of (?)) was remarkable and its peak bolometric luminosity of 4$\times$1044
erg s-1 makes it the brightest FBOT known so far (?). As the discovery was
promptly reported (?), the source received a significant amount of multi-
wavelength coverage. Radio, millimeter, optical, UV, X-ray and gamma ray
properties of the source are described in various papers (?, ?, ?, ?, ?, ?,
?). However, in spite of exquisite coverage the physical origin of AT2018cow
remains elusive.

Given its high X-ray luminosity (peak value of a roughly 1043 erg/s, (?)) and
variability on timescales of a few tens of hours (?, ?, ?), compact object
(accretion) powered scenarios have been proposed for AT2018cow. These
suggestions include emission from tidal disruption of a star by an
intermediate-mass black hole with mass in the range of 104-6 $M_{\odot}$ (?),
fallback accretion in a failed supernova (?, ?), and energy injection by a
newly born neutron star in a supernova (?).

Several works over the last few decades (?) have found that when accreting
compact object (stellar-mass black hole/neutron star) X-ray binaries go into
outbursts–due to enhanced accretion–they sometimes exhibit high-frequency (a
few$\times$(10-100) Hz) quasi-periodic variability in their X-ray brightness.
There is no clear consensus on the exact mechanisms that produce these so-
called High-Frequency Quasi-Periodic Oscillations (HFQPOs) but it is generally
agreed upon that they originate from a region close to the compact object
where the dynamics of motion are dictated by the compact object’s strong
gravitational field (see (?) and references therein), and they represent a
direct evidence for the presence of a compact object. Some active galactic
nuclei have also shown evidence for QPOs which have been argued to be
analogous to HFQPOs of stellar-mass black holes (?, ?, ?). More recently,
HFQPO analogs (frequencies of a few mHz) have also been found in stellar tidal
disruption events involving $\sim$106 $M_{\odot}$ black holes (?, ?, ?)
suggesting that, perhaps, such QPOs are universal among all compact object
systems that undergo extreme changes in accretion.

To test these hypotheses for a compact object (accretion) powered scenario, we
studied AT2018cow’s X-ray (0.25-2.5 keV) variability using an average power
density spectrum derived from the entire monitoring data taken by the Neutron
Star Interior Composition Explorer (NICER) on board the International Space
Station. We find evidence for an X-ray QPO in the average PDS (left panel of
Fig. 1). The QPO signal has a centroid frequency, full width half maximum
(FWHM) and a fractional root-mean-squared amplitude of 224.4$\pm$1.0 Hz, $<$
16 Hz, and 30$\pm$3%, respectively. Using a rigorous Monte Carlo approach (see
supplementary material, SM) that takes into account the nature of the
underlying noise continuum and the search trials, we find the global false
alarm probability of this signal to be $\approx$2$\times$10-4 (or 3.7$\sigma$
equivalent for a normal distribution; see right panel of Fig. 1 and SM).

We rule out various instrumental and particle backgrounds as the origin for
this QPO signal (See SM). Based on the long term light curves derived from
NICER and Neil Gehrels Swift X-ray telescope, which has imaging capability
(see top panel of Fig. 2), and late time XMM-Newton X-ray images of
AT2018cow’s field of view (FoV), we also rule out a contaminating source as
the origin of this QPO (see Fig. 3). Furthermore, over the last three years of
NICER operations a signal of similar nature has never been found in any of the
several dozens of other targets (see, for example, Fig. LABEL:fig:agnpds).
Based on these tests, we conclude that the signal is consistent with
originating from AT2018cow.

After establishing that the signal is statistically significant and ruling out
an instrumental and a background origin, we extracted the QPO’s signal-to-
noise as a function of the accumulated exposure (See Fig. LABEL:fig:deltachi
and supplementary movie S1). It is evident that the QPO’s strength increases
gradually with increasing exposure. This suggests that the signal is
persistent and present at some level throughout the $\sim$ 2 month monitoring
period. Taken at face value, this suggests that it is stable over $\sim$60
days/4.44 ms $\gtrsim$ 109 cycles. Interestingly, the mean slope of the curve
in Fig. LABEL:fig:deltachi also steepens around the same time (near day 17)
when high-amplitude X-ray flares start to appear on days timescale (see the
blue diamonds in the upper panel of Fig. 2). Moreover, the fractional rms of
the QPO jumps around day 17 which also coincides with the time when the
optical spectrum of AT2018cow underwent dramatic changes (see (?) for
details). We separated the total exposure into two epochs, before and after
day 17. We extract an average PDS from each of these two time intervals and
the QPO is fit with a Lorentzian. The QPO’s fractional rms amplitude appears
to be higher at later times (see bottom panel of Fig. 2).

The frequency of this QPO alone can set stringent constraints on the
underlying physical mechanism producing X-rays in AT2018cow. The causality
argument suggests that the physical size of an object producing this signal
cannot be larger than the light crossing size, i.e., speed of light $\times$
(1/224 s) $\approx$ 1.3$\times$108 cm. This small size points us naturally
towards a compact object. Emission from shock interactions (?, ?) is
disfavored because CSM–shock interactions are not known to ”pulse”. If the
compact object is a black hole then assuming the emission originates from the
innermost stable circular orbit allows us to set an upper limit on the black
hole mass. These limits are 95 $M_{\odot}$ and 850 $M_{\odot}$ for a maximally
prograde spinning (spin has the same direction as the material falling in) and
a retrograde (vice versa) spinning black hole, respectively. A larger emission
radius (in units of gravitational radii) would require the object to be even
more compact. This rules out a heavy IMBH ($\gtrsim$ 850 $M_{\odot}$) in
AT2018cow (?).

HFQPOs with frequencies in the hundreds of Hz range have been seen in a
handful of stellar-mass black holes with known dynamical mass estimates (?).
They often appear in pairs with frequency ratio of 2:3 (?). In these systems
it appears that the black hole mass scales inversely with the HFQPO frequency
(see, for example, Fig. 4.17 of (?)). While the frequency and width of
AT2018cow’s QPO are similar to HFQPOs of stellar-mass black holes, the
observed fractional rms amplitude is high (25-45% compared to a few percent in
stellar-mass black holes (?)), and a harmonic is not apparent here (see sec.
LABEL:supsec:harmonics and Table LABEL:tab:rmstab for upper limits).
Furthermore, the QPO’s energy dependence is also uncertain given our narrow
X-ray band pass of 0.25-2.5 keV. Also, the absence of any red noise (even at
lower frequencies, i.e., a few$\times$mHz; see Fig. LABEL:fig:xmmpds) is
unlike accreting X-ray binaries. All the stellar-mass X-ray binaries with
HFQPOs are relatively highly inclined (?). Given AT2018cow’s high (a
few$\times$1042 erg/s) average X-ray luminosity beaming, which can boost the
apparent brightness, is likely present, and this points towards a low
inclination for this system. Thus, given all these factors, it is unclear if a
direct comparison can be made with HFQPOs of stellar-mass black hole binaries.
Nevertheless, if we assume the same scaling law applies and that 224 Hz is the
fundamental harmonic, then the implied black hole mass in AT2018cow is $\sim$
4 $M_{\odot}$.

Alternatively, the compact object could be a neutron star with the QPO
representing either its spin rate of 224 Hz or an analog of the upper kHz QPOs
which can have high fractional rms amplitudes (see Fig. 2 of (?)). We explore
multiple physical scenarios for the former case. First, we consider the case
where AT2018cow’s luminosity is driven by the spin down of an isolated
millisecond magnetar. In this case, the high and changing bolometric
luminosity would require high magnetic field strength which in turn would lead
to a rapid change in the spin frequency. This is inconsistent with the QPO’s
stability based on its width ($<16\mbox{ Hz}/60\mbox{ d}=3\times 10^{-6}\mbox{
Hz s${}^{-1}$}$). This inconsistency rules out an isolated magnetar powered
scenario. Then we consider a case in which the luminosity is powered by
accretion of matter onto the magnetic poles of a neutron star. While this
model provides a possible mechanism to induce periodicity in the light curve
and a declining luminosity, its prediction for spin period evolution is
inconsistent with the QPO’s stability (see SM for more details). Lastly, if
this QPO is due to a process similar to that operating in magnetar QPOs, the
physics and damping mechanisms would be in an entirely different regime given
that the stability of AT2018cow’s QPO is significantly ($\gtrsim$105) longer
than those of magnetar burst QPOs (?).

Within the bounds of a central engine being either a neutron star, stellar-
mass black hole or a few hundred M⊙ black hole, the QPO presented here can be
produced under different scenarios. For example, it has been proposed that a
flare like AT2018cow could arise from young star clusters containing a
stellar-mass black hole (mass $\sim$ 10 M⊙) tidally disrupting a main sequence
star within the cluster (?) or from fallback accretion of material in a failed
supernova whose core collapses into either a neutron star or a stellar-mass
black hole (?). In all these scenarios, such a QPO signal could be envisioned.
Using the Atacama Large Millimeter/submillimeter Array (ALMA) (?) spatially
resolved the molecular gas content and star formation rate within AT2018cow’s
host galaxy. They find that AT2018cow is situated between a peak in molecular
gas content and a blue star cluster. The presence of these regions indicate
active star formation around AT2018cow and favors an association with deaths
of massive stars for AT2018cow. The implications of this QPO for AT2018cow’s
existing models are summarized in Table 1.

In summary, we present evidence for a 224 Hz X-ray QPO persistent for one
billion cycles from AT2018cow statistically significant at the 3.7$\sigma$
level. Previous works (?) have suggested a compact object in AT2018cow based
on X-ray variability on a few tens of hours timescale and also its X-ray
spectral resemblance to some accreting systems. The detection of regular
variations on millisecond timescale presented here provides the most direct
way to infer the presence of a compact object, although it still remains
unclear whether that compact object is a stellar-mass black hole or a neutron
star. If AT2018cow originated from the death of a massive star, our findings
represent the birth of a compact object in a supernova. Similar signals in
future FBOTs could allow astronomers to study infant compact objects
immediately after birth. Assuming all luminous FBOTs are AT2018cow-like which
have an estimated volumetric rate of $<$10-7 yr-1 Mpc-3 (?, ?, ?, ?), NICER
can, in principle, find and study QPOs from one such system every 3 years (see
SM sec. LABEL:supsec:rates).

Table 1: Models for AT2018cow’s X-ray emission and their validity against the
224 Hz QPO reported here. This table can be considered as a modified extension
of Table 2 of (?).

Model/class of models |  References∙ |  Consistent with QPO? |  Notes
---|---|---|---
Shock interactions with circumstellar medium (CSM) |  (?, ?) |  No |  CSM interaction could explain emission at other (non X-ray) wavelengths, i.e., optical and radio but inconsistent with rapid X-ray variability
An embedded internal shock formed from interaction with dense CSM |  (?) |  No |  A compact embedded internal shock is, in principle, consistent with the size constraint provided by the 224 Hz QPO. However, all X-ray QPOs known thus far in literature are from compact objects. So an internal shock model is disfavored.
Accreting intermediate-mass black hole ($\gtrsim$103 $M_{\odot}$) |  (?, ?, ?) |  No |  Based on causality argument the compact object producing the QPO has to be less than 850 $M_{\odot}$ (see main text)
Neutron star (formed from merging white dwarfs/Supernova) |  (?, ?) |  Yes |  Constraints on magnetic field if the 224 Hz QPO represents the spin period. However, the persistence of the signal in a narrow frequency range is challenging to explain (see SM)
Stellar-mass black hole (accreting from outer layers of a failed supernova/tidally disrupting a star) |  (?, ?, ?, ?) |  Yes |  The QPO frequency is similar to those often seen in known stellar-mass black holes but the X-ray luminosity, rms, and QPO’s stability are unlike any known stellar-mass black hole systems (see main text)
|  |  |

∙Not an exhaustive list. Also see references therein.

Figure 1: (a) Average X-ray PDS of AT2018cow showing evidence for a quasi-
periodicity near 224 Hz. This PDS was computed by averaging 105 256-second
soft X-ray (0.25-2.5 keV) light curves sampled at 1/2048 s. The resulting PDS
was further re-binned by a factor of 2048 which gives a frequency resolution
of 8 Hz. The strongest excess above the Poisson noise level of 2 is around 224
Hz. The power values in the rest of the PDS continuum are consistent with
white noise (see SM sections LABEL:supsec:kstest and LABEL:supsec:fap and
Figures LABEL:fig:kstest, LABEL:fig:probplot, LABEL:fig:xmmpds). The PDS is
normalized such that the mean value surrounding 224 Hz is equal to the Poisson
value of 2. The best-fit constant + Lorentzian models are indicated by the
dashed red curve. (b) Estimates for the statistical significance of the 224 Hz
QPO. The likelihood of finding a QPO from noise fluctuations, i.e., false
alarm probability = 1-CDF, (y-axis) vs the maximum improvement in $\chi^{2}$
by fitting the simulated noise PDS with a constant+Lorentzian over modeling it
with a constant at every frequency searched during the identification of the
signal on the left (see section LABEL:supsec:mcsims and Fig. LABEL:fig:pdcdfs
for more details). The global false alarm probability of finding a QPO as
strong as the one observed is $\approx$ 2$\times$10-4 (1 in 5000;
$\approx$3.7$\sigma$).

Figure 2: (a) Comparison of NICER/XTI and Neil Gehrels Swift/XRT long-term
light curves. Both y-axes are in units of counts/sec. It is evident that
AT2018cow’s long-term soft X-ray (0.25-2.5 keV) variability as observed by the
non-imaging NICER telescope is same as that observed with Neil Gehrels
Swift/XRT (0.3-2.5 keV) which has imaging capability. This suggests that the
flux observed by NICER is dominated by AT2018cow with minimal contamination
from other nearby astrophysical sources. (b) Fractional root-mean-squared
amplitude of the 224 Hz QPO vs time. The signal appears to be stronger during
the end of the outburst. The two values are 25$\pm$4 and 45$\pm$7%
corresponding to exposures during 12.7${}^{+5.2}_{-7.3}$ and
33.2${}^{+30.5}_{-15.3}$ days, respectively. The jump in the QPO’s strength
coincides with higher levels of variability on days timescale (XRT light curve
in (a)) and also with the dramatic changes in optical spectra seen around the
same time (?).

Figure 3: XMM-Newton and Swift images of NICER’s field of view of AT2018cow
showing that there was minimal contamination from field sources. (a) XMM-
Newton/MOS1 image of AT2018cow’s field of view long after AT2018cow faded. The
position of AT2018cow is indicated by a blue circle with a radius of 33”.
There is no contaminating point source directly coincident with AT2018cow. (b)
Stacked Neil Gehrels Swift/XRT image of the field of view of AT2018cow. The
blue circle of radius 47” is centered on AT2018cow’s optical position
(16:16:00.220 +22:16:04.91; J2000.0). This particular image was extracted by
using all the archival XRT images as of July 2020. It is evident that while
there are a few point sources present in NICER’s field of view, their
contribution to overall XTI flux is negligible when compared to AT2018cow (see
also Fig. 2). In both panels the outer/red (dashed) circles show NICER/XTI’s
approximate field of view of 3.1’. The north and east arrows are each 100”
long. The colorbar show counts.

## References

  * 1. C. Inserra, Observational properties of extreme supernovae, Nature Astronomy 3, 697 (2019).
  * 2. M. R. Drout, et al., Rapidly Evolving and Luminous Transients from Pan-STARRS1, Astrophys. J. 794, 23 (2014).
  * 3. M. Pursiainen, et al., Rapidly evolving transients in the Dark Energy Survey, Mon. Not. R. Astron. Soc. 481, 894 (2018).
  * 4. L. E. Rivera Sandoval, et al., X-ray Swift observations of SN 2018cow, Mon. Not. R. Astron. Soc. 480, L146 (2018).
  * 5. R. Margutti, et al., An Embedded X-Ray Source Shines through the Aspherical AT 2018cow: Revealing the Inner Workings of the Most Luminous Fast-evolving Optical Transients, Astrophys. J. 872, 18 (2019).
  * 6. S. J. Prentice, et al., The Cow: Discovery of a Luminous, Hot, and Rapidly Evolving Transient, Astrophys. J. 865, L3 (2018).
  * 7. D. A. Perley, et al., The fast, luminous ultraviolet transient AT2018cow: extreme supernova, or disruption of a star by an intermediate-mass black hole?, Mon. Not. R. Astron. Soc. 484, 1031 (2019).
  * 8. M. Tanaka, et al., Rapidly Rising Transients from the Subaru Hyper Suprime-Cam Transient Survey, Astrophys. J. 819, 5 (2016).
  * 9. A. Rest, et al., A fast-evolving luminous transient discovered by K2/Kepler, Nature Astronomy 2, 307 (2018).
  * 10. A. Y. Q. Ho, et al., The Photometric and Spectroscopic Evolution of Rapidly Evolving Extragalactic Transients in ZTF, arXiv e-prints p. arXiv:2105.08811 (2021).
  * 11. D. A. Perley, et al., Real-time Discovery of AT2020xnd: A Fast, Luminous Ultraviolet Transient with Minimal Radioactive Ejecta, arXiv e-prints p. arXiv:2103.01968 (2021).
  * 12. D. L. Coppejans, et al., A Mildly Relativistic Outflow from the Energetic, Fast-rising Blue Optical Transient CSS161010 in a Dwarf Galaxy, Astrophys. J. 895, L23 (2020).
  * 13. A. Y. Q. Ho, et al., The Koala: A Fast Blue Optical Transient with Luminous Radio Emission from a Starburst Dwarf Galaxy at z = 0.27, Astrophys. J. 895, 49 (2020).
  * 14. S.-C. Leung, et al., A Model for the Fast Blue Optical Transient AT2018cow: Circumstellar Interaction of a Pulsational Pair-instability Supernova, Astrophys. J. 903, 66 (2020).
  * 15. O. D. Fox, N. Smith, Signatures of circumstellar interaction in the unusual transient AT 2018cow, Mon. Not. R. Astron. Soc. 488, 3772 (2019).
  * 16. D. Kasen, L. Bildsten, Supernova Light Curves Powered by Young Magnetars, Astrophys. J. 717, 245 (2010).
  * 17. B. D. Metzger, A. L. Piro, Optical and X-ray emission from stable millisecond magnetars formed from the merger of binary neutron stars, Mon. Not. R. Astron. Soc. 439, 3916 (2014).
  * 18. M. Lyutikov, S. Toonen, Fast-rising blue optical transients and AT2018cow following electron-capture collapse of merged white dwarfs, Mon. Not. R. Astron. Soc. 487, 5618 (2019).
  * 19. N. P. M. Kuin, et al., Swift spectra of AT2018cow: a white dwarf tidal disruption event?, Mon. Not. R. Astron. Soc. 487, 2505 (2019).
  * 20. M. Tanaka, et al., Rapidly Rising Transients from the Subaru Hyper Suprime-Cam Transient Survey, Astrophys. J. 819, 5 (2016).
  * 21. M. Pursiainen, et al., Rapidly evolving transients in the Dark Energy Survey, Mon. Not. R. Astron. Soc. 481, 894 (2018).
  * 22. J. L. Tonry, et al., ATLAS: A High-cadence All-sky Survey System, Publ. Astron. Soc. Pac. 130, 064505 (2018).
  * 23. S. J. Smartt, et al., ATLAS18qqn (AT2018cow) - a bright transient spatially coincident with CGCG 137-068 (60 Mpc), The Astronomer’s Telegram 11727, 1 (2018).
  * 24. A. Y. Q. Ho, et al., AT2018cow: A Luminous Millimeter Transient, Astrophys. J. 871, 73 (2019).
  * 25. K. Huang, et al., ALMA Polarimetry of AT2018cow, Astrophys. J. 878, L25 (2019).
  * 26. E. Quataert, D. Kasen, Swift 1644+57: the longest gamma-ray burst?, Mon. Not. R. Astron. Soc. 419, L1 (2012).
  * 27. P. Mohan, T. An, J. Yang, The Nearby Luminous Transient AT2018cow: A Magnetar Formed in a Subrelativistically Expanding Nonjetted Explosion, Astrophys. J. 888, L24 (2020).
  * 28. J. E. McClintock, R. A. Remillard, Black hole binaries (2006), vol. 39, pp. 157–213.
  * 29. W. N. Alston, et al., Discovery of an $\sim$2-h high-frequency X-ray QPO and iron K$\alpha$ reverberation in the active galaxy MS 2254.9-3712, Mon. Not. R. Astron. Soc. 449, 467 (2015).
  * 30. D. I. Ashton, M. J. Middleton, Searching for energy-resolved quasi-periodic oscillations in AGN, Mon. Not. R. Astron. Soc. 501, 5478 (2021).
  * 31. M. Gierliński, M. Middleton, M. Ward, C. Done, A periodicity of ~1hour in X-ray emission from the active galaxy RE J1034+396, Nature 455, 369 (2008).
  * 32. D. R. Pasham, et al., A loud quasi-periodic oscillation after a star is disrupted by a massive black hole, Science 363, 531 (2019).
  * 33. R. C. Reis, et al., A 200-Second Quasi-Periodicity After the Tidal Disruption of a Star by a Dormant Black Hole, Science 337, 949 (2012).
  * 34. D. Lin, J. A. Irwin, O. Godet, N. A. Webb, D. Barret, A ~3.8 hr Periodicity from an Ultrasoft Active Galactic Nucleus Candidate, Astrophys. J. 776, L10 (2013).
  * 35. E. M. Ribeiro, M. Méndez, M. G. B. de Avellar, G. Zhang, K. Karpouzas, The amplitude of the kilohertz quasi-periodic oscillations in 4U 1636-53 in the frequency-energy space, Mon. Not. R. Astron. Soc. 489, 4980 (2019).
  * 36. A. L. Watts, T. E. Strohmayer, Neutron star oscillations and QPOs during magnetar flares, Advances in Space Research 40, 1446 (2007).
  * 37. K. Kremer, et al., Fast Optical Transients from Stellar-Mass Black Hole Tidal Disruption Events in Young Star Clusters, arXiv e-prints p. arXiv:2012.02796 (2020).
  * 38. C. Frohmaier, et al., From core collapse to superluminous: the rates of massive stellar explosions from the Palomar Transient Factory, Mon. Not. R. Astron. Soc. 500, 5142 (2021).
  * 39. M. Taylor, et al., The Core Collapse Supernova Rate from the SDSS-II Supernova Survey, Astrophys. J. 792, 135 (2014).
  * 40. J.-H. Chen, R.-F. Shen, Tidal Disruption of a Main-sequence Star by an Intermediate-mass Black Hole: A Bright Decade, Astrophys. J. 867, 20 (2018).
  * 41. J. Dexter, D. Kasen, Supernova Light Curves Powered by Fallback Accretion, Astrophys. J. 772, 30 (2013).
  * 42. E. Quataert, D. Lecoanet, E. R. Coughlin, Black hole accretion discs and luminous transients in failed supernovae from non-rotating supergiants, Mon. Not. R. Astron. Soc. 485, L83 (2019).
  * 43. K. C. Gendreau, et al., SPIE (2016), vol. 9905 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, p. 99051H.
  * 44. N. Gehrels, et al., The Swift Gamma-Ray Burst Mission, Astrophys. J. 611, 1005 (2004).
  * 45. D. N. Burrows, et al., The Swift X-Ray Telescope, Space Sci. Rev. 120, 165 (2005).
  * 46. F. Jansen, et al., XMM-Newton observatory. I. The spacecraft and operations, Astron. Astrophys. 365, L1 (2001).
  * 47. L. Strüder, et al., The European Photon Imaging Camera on XMM-Newton: The pn-CCD camera, Astron. Astrophys. 365, L18 (2001).
  * 48. G. Prigozhin, et al., SPIE (2012), vol. 8453 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, p. 845318.
  * 49. R. A. Remillard, et al., An Empirical Background Model for the NICER X-ray Timing Instrument, arXiv e-prints p. arXiv:2105.09901 (2021).
  * 50. J. S. Kaastra, J. A. M. Bleeker, Optimal binning of X-ray spectra and response matrix design, Astron. Astrophys. 587, A151 (2016).
  * 51. M. van der Klis, Timing Neutron Stars (1988), vol. 262, pp. 27–70.
  * 52. S. Vaughan, R. Edelson, R. S. Warwick, P. Uttley, On characterizing the variability properties of X-ray light curves from active galaxies, Mon. Not. R. Astron. Soc. 345, 1271 (2003).
  * 53. S. L. Shapiro, S. A. Teukolsky, Black holes, white dwarfs, and neutron stars : the physics of compact objects (1983).
  * 54. B. D. Metzger, P. Beniamini, D. Giannios, Effects of Fallback Accretion on Protomagnetar Outflows in Gamma-Ray Bursts and Superluminous Supernovae, Astrophys. J. 857, 95 (2018).
  * 55. V. F. Shvartsman, Neutron Stars in Binary Systems Should Not Be Pulsars., Soviet Astronomy 15, 342 (1971).
  * 56. A. F. Illarionov, R. A. Sunyaev, Why the Number of Galactic X-ray Stars Is so Small?, Astron. Astrophys. 39, 185 (1975).
  * 57. V. M. Lipunov, Astrophysics of Neutron Stars (1992).
  * 58. K. Davidson, J. P. Ostriker, Neutron-Star Accretion in a Stellar Wind: Model for a Pulsed X-Ray Source, Astrophys. J. 179, 585 (1973).
  * 59. A. L. Piro, C. D. Ott, Supernova Fallback onto Magnetars and Propeller-powered Supernovae, Astrophys. J. 736, 108 (2011).
  * 60. K. Parfrey, A. Spitkovsky, A. M. Beloborodov, Torque Enhancement, Spin Equilibrium, and Jet Power from Disk-Induced Opening of Pulsar Magnetic Fields, Astrophys. J. 822, 33 (2016).
  * 61. M. A. Alpar, On Young Neutron Stars as Propellers and Accretors with Conventional Magnetic Fields, Astrophys. J. 554, 1245 (2001).
  * 62. S. M. Ransom, S. S. Eikenberry, J. Middleditch, Fourier Techniques for Very Long Astrophysical Time-Series Analysis, Astron. J. 124, 1788 (2002).
  * 63. B. C. Andersen, S. M. Ransom, A Fourier Domain “Jerk” Search for Binary Pulsars, Astrophys. J. 863, L13 (2018).
  * 64. M. Bagchi, D. R. Lorimer, S. Wolfe, On the detectability of eccentric binary pulsars, Mon. Not. R. Astron. Soc. 432, 1303 (2013).

  * Supplementary Materials.
Materials and Methods
Figs. S1 to S11
Table S1
Caption for movie S1
Movie S1
Supplementary Text
All the codes to download the data from HEASARC public archive, reduce it, and
reproduce the results in this paper are available as supplementary files