* [55] James Y. Widdicombe, Thomas Helfer and Eugene A. Lim  “Black hole formation in relativistic Oscillaton collisions”  In _JCAP_ 01, 2020, pp. 027  DOI: 10.1088/1475-7516/2020/01/027
  * [56] Zainab Nazari, Michele Cicoli, Katy Clough and Francesco Muia  “Oscillon collapse to black holes”  In _JCAP_ 05, 2021, pp. 027  DOI: 10.1088/1475-7516/2021/05/027
  * [57] Richard Easther and Eugene A. Lim  “Stochastic gravitational wave production after inflation”  In _JCAP_ 04, 2006, pp. 010  DOI: 10.1088/1475-7516/2006/04/010
  * [58] Jean Francois Dufaux et al.  “Theory and Numerics of Gravitational Waves from Preheating after Inflation”  In _Phys. Rev. D_ 76, 2007, pp. 123517  DOI: 10.1103/PhysRevD.76.123517
  * [59] Stefan Antusch, Francesco Cefala and Stefano Orani  “Gravitational waves from oscillons after inflation” [Erratum: Phys.Rev.Lett. 120, 219901 (2018)]  In _Phys. Rev. Lett._ 118.1, 2017, pp. 011303  DOI: 10.1103/PhysRevLett.118.011303
  * [60] Jing Liu, Zong-Kuan Guo, Rong-Gen Cai and Gary Shiu  “Gravitational Waves from Oscillons with Cuspy Potentials”  In _Phys. Rev. Lett._ 120.3, 2018, pp. 031301  DOI: 10.1103/PhysRevLett.120.031301
  * [61] Thomas Helfer, Eugene A. Lim, Marcos A.. Garcia and Mustafa A. Amin  “Gravitational Wave Emission from Collisions of Compact Scalar Solitons”  In _Phys. Rev. D_ 99.4, 2019, pp. 044046  DOI: 10.1103/PhysRevD.99.044046
  * [62] Takashi Hiramatsu, Evangelos I. Sfakianakis and Masahide Yamaguchi  “Gravitational wave spectra from oscillon formation after inflation”  In _JHEP_ 03, 2021, pp. 021  DOI: 10.1007/JHEP03(2021)021
  * [63] Wayne Hu, Rennan Barkana and Andrei Gruzinov  “Cold and fuzzy dark matter”  In _Phys. Rev. Lett._ 85, 2000, pp. 1158–1161  DOI: 10.1103/PhysRevLett.85.1158
  * [64] Jan Ollé, Oriol Pujolàs and Fabrizio Rompineve  “Oscillons and Dark Matter”  In _JCAP_ 02, 2020, pp. 006  DOI: 10.1088/1475-7516/2020/02/006
  * [65] Masahiro Kawasaki, Wakutaka Nakano and Eisuke Sonomoto  “Oscillon of Ultra-Light Axion-like Particle”  In _JCAP_ 01, 2020, pp. 047  DOI: 10.1088/1475-7516/2020/01/047
  * [66] Elisa G.. Ferreira  “Ultra-light dark matter”  In _Astron. Astrophys. Rev._ 29.1, 2021, pp. 7  DOI: 10.1007/s00159-021-00135-6
  * [67] Mustafa A. Amin, Richard Easther and Hal Finkel  “Inflaton Fragmentation and Oscillon Formation in Three Dimensions”  In _JCAP_ 12, 2010, pp. 001  DOI: 10.1088/1475-7516/2010/12/001
  * [68] Mustafa A. Amin and David Shirokoff  “Flat-top oscillons in an expanding universe”  In _Phys. Rev. D_ 81, 2010, pp. 085045  DOI: 10.1103/PhysRevD.81.085045
  * [69] Mustafa A. Amin  “Inflaton fragmentation: Emergence of pseudo-stable inflaton lumps (oscillons) after inflation”, 2010  arXiv:1006.3075 [astro-ph.CO]
  * [70] Yu Sang and Qing-Guo Huang  “Stochastic Gravitational-Wave Background from Axion-Monodromy Oscillons in String Theory During Preheating”  In _Phys. Rev. D_ 100.6, 2019, pp. 063516  DOI: 10.1103/PhysRevD.100.063516
  * [71] Stefan Antusch, Francesco Cefala and Stefano Orani  “What can we learn from the stochastic gravitational wave background produced by oscillons?”  In _JCAP_ 03, 2018, pp. 032  DOI: 10.1088/1475-7516/2018/03/032
  * [72] Stefan Antusch, Francesco Cefalà and Francisco Torrentı́  “Properties of Oscillons in Hilltop Potentials: energies, shapes, and lifetimes”  In _JCAP_ 10, 2019, pp. 002  DOI: 10.1088/1475-7516/2019/10/002
  * [73] Daniel Baumann and Liam McAllister  “Inflation and String Theory”, Cambridge Monographs on Mathematical Physics  Cambridge University Press, 2015  DOI: 10.1017/CBO9781316105733
  * [74] Michele Cicoli et al.  “String cosmology: From the early universe to today”  In _Phys. Rept._ 1059, 2024, pp. 1–155  DOI: 10.1016/j.physrep.2024.01.002
  * [75] Fien Apers et al.  “String Theory and the First Half of the Universe”, 2024  arXiv:2401.04064 [hep-th]
  * [76] Swagat S. Mishra et al.  “Formation and decay of oscillons after inflation in the presence of an external coupling, Part-II: Analytical Results (In preparation)”
  * [77] Richard L. Arnowitt, Stanley Deser and Charles W. Misner  “The Dynamics of general relativity”  In _Gen. Rel. Grav._ 40, 2008, pp. 1997–2027  DOI: 10.1007/s10714-008-0661-1
  * [78] Juan Martin Maldacena  “Non-Gaussian features of primordial fluctuations in single field inflationary models”  In _JHEP_ 05, 2003, pp. 013  DOI: 10.1088/1126-6708/2003/05/013
  * [79] Daniel Baumann  “Primordial Cosmology”  In _PoS_ TASI2017, 2018, pp. 009  DOI: 10.22323/1.305.0009
  * [80] Andrew R Liddle and Samuel M Leach  “How long before the end of inflation were observable perturbations produced?”  In _Phys. Rev. D_ 68, 2003, pp. 103503  DOI: 10.1103/PhysRevD.68.103503
  * [81] Swagat S. Mishra, Varun Sahni and Alexei A. Starobinsky  “Curing inflationary degeneracies using reheating predictions and relic gravitational waves”  In _JCAP_ 05, 2021, pp. 075  DOI: 10.1088/1475-7516/2021/05/075
  * [82] Andrew R. Liddle, Paul Parsons and John D. Barrow  “Formalizing the slow roll approximation in inflation”  In _Phys. Rev. D_ 50, 1994, pp. 7222–7232  DOI: 10.1103/PhysRevD.50.7222
  * [83] Robert Brandenberger  “Initial conditions for inflation — A short review”  In _Int. J. Mod. Phys. D_ 26.01, 2016, pp. 1740002  DOI: 10.1142/S0218271817400028
  * [84] Swagat S. Mishra, Varun Sahni and Alexey V. Toporensky  “Initial conditions for Inflation in an FRW Universe”  In _Phys. Rev. D_ 98.8, 2018, pp. 083538  DOI: 10.1103/PhysRevD.98.083538
  * [85] Y. Akrami  “Planck 2018 results. X. Constraints on inflation”  In _Astron. Astrophys._ 641, 2020, pp. A10  DOI: 10.1051/0004-6361/201833887
  * [86] Renata Kallosh and Andrei Linde  “Universality Class in Conformal Inflation”  In _JCAP_ 07, 2013, pp. 002  DOI: 10.1088/1475-7516/2013/07/002
  * [87] Renata Kallosh, Andrei Linde and Diederik Roest  “Superconformal Inflationary $\alpha$-Attractors”  In _JHEP_ 11, 2013, pp. 198  DOI: 10.1007/JHEP11(2013)198
  * [88] Stefan Antusch, Daniel G. Figueroa, Kenneth Marschall and Francisco Torrenti  “Characterizing the postinflationary reheating history: Single daughter field with quadratic-quadratic interaction”  In _Phys. Rev. D_ 105.4, 2022, pp. 043532  DOI: 10.1103/PhysRevD.105.043532
  * [89] F. Finelli and Robert H. Brandenberger  “Parametric amplification of gravitational fluctuations during reheating”  In _Phys. Rev. Lett._ 82, 1999, pp. 1362–1365  DOI: 10.1103/PhysRevLett.82.1362
  * [90] Daniel G. Figueroa and Francisco Torrenti  “Parametric resonance in the early Universe—a fitting analysis”  In _JCAP_ 02, 2017, pp. 001  DOI: 10.1088/1475-7516/2017/02/001
  * [91] Daniel G. Figueroa, Adrien Florio, Francisco Torrenti and Wessel Valkenburg  “The art of simulating the early Universe – Part I”  In _JCAP_ 04, 2021, pp. 035  DOI: 10.1088/1475-7516/2021/04/035
  * [92] Daniel G. Figueroa, Adrien Florio, Francisco Torrenti and Wessel Valkenburg  “CosmoLattice: A modern code for lattice simulations of scalar and gauge field dynamics in an expanding universe”  In _Comput. Phys. Commun._ 283, 2023, pp. 108586  DOI: 10.1016/j.cpc.2022.108586
  * [93] G. Floquet  “Sur les équations différentielles linéaires à coefficients périodiques”  In _Annales scientifiques de l’École Normale Supérieure_ 2e série, 12 Elsevier, 1883, pp. 47–88  DOI: 10.24033/asens.220
  * [94] N.. McLachlan  “Theory and Application of Mathieu Functions”  Clarendon Press, 1947
  * [95] Milton Abramowitz and Irene A. Stegun  “Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables”  New York City: Dover, 1964
  * [96] “NIST Digital Library of Mathematical Functions” F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds., https://dlmf.nist.gov/, Release 1.1.12 of 2023-12-15  URL: https://dlmf.nist.gov/
  * [97] Swagat S. Mishra  “Cosmic Inflation: Background dynamics, Quantum fluctuations and Reheating”, 2024  arXiv:2403.10606 [gr-qc]
  * [98] Gary N. Felder et al.  “Dynamics of symmetry breaking and tachyonic preheating”  In _Phys. Rev. Lett._ 87, 2001, pp. 011601  DOI: 10.1103/PhysRevLett.87.011601
  * [99] Gary N. Felder, Lev Kofman and Andrei D. Linde  “Tachyonic instability and dynamics of spontaneous symmetry breaking”  In _Phys. Rev. D_ 64, 2001, pp. 123517  DOI: 10.1103/PhysRevD.64.123517
  * [100] Masahiro Ibe, Masahiro Kawasaki, Wakutaka Nakano and Eisuke Sonomoto  “Decay of I-ball/Oscillon in Classical Field Theory”  In _JHEP_ 04, 2019, pp. 030  DOI: 10.1007/JHEP04(2019)030
  * [101] Hong-Yi Zhang  “Gravitational effects on oscillon lifetimes”  In _JCAP_ 03, 2021, pp. 102  DOI: 10.1088/1475-7516/2021/03/102
  * [102] David Cyncynates and Tudor Giurgica-Tiron  “Structure of the oscillon: The dynamics of attractive self-interaction”  In _Phys. Rev. D_ 103.11, 2021, pp. 116011  DOI: 10.1103/PhysRevD.103.116011
  * [103] Laura Iacconi, Matteo Fasiello, Jussi Väliviita and David Wands  “Novel CMB constraints on the $\alpha$ parameter in alpha-attractor models”  In _JCAP_ 10, 2023, pp. 015  DOI: 10.1088/1475-7516/2023/10/015
  * [104] Thibault Damour and Viatcheslav F. Mukhanov  “Inflation without slow roll”  In _Phys. Rev. Lett._ 80, 1998, pp. 3440–3443  DOI: 10.1103/PhysRevLett.80.3440
  * [105] Andrew R. Liddle and Anupam Mazumdar  “Inflation during oscillations of the inflaton”  In _Phys. Rev. D_ 58, 1998, pp. 083508  DOI: 10.1103/PhysRevD.58.083508
  * [106] William T. Emond, Peter Millington and Paul M. Saffin  “Boltzmann equations for preheating”  In _JCAP_ 09, 2018, pp. 041  DOI: 10.1088/1475-7516/2018/09/041
  * [107] Marcos A.. Garcia et al.  “Effects of fragmentation on post-inflationary reheating”  In _JCAP_ 12, 2023, pp. 028  DOI: 10.1088/1475-7516/2023/12/028
  * [108] Marcos A.. Garcia and Mathias Pierre  “Reheating after inflaton fragmentation”  In _JCAP_ 11, 2023, pp. 004  DOI: 10.1088/1475-7516/2023/11/004
  * [109] Marcos A.. Garcia and Mathias Pierre  “Gravitational wave signatures of post-fragmentation reheating”, 2024  arXiv:2404.16932 [hep-ph]
  * [110] Sugumi Kanno, Jiro Soda and Akira Taniguchi  “Search for high-frequency gravitational waves with Rydberg atoms”, 2023  arXiv:2311.03890 [gr-qc]
  * [111] Nancy Aggarwal  “Challenges and opportunities of gravitational-wave searches at MHz to GHz frequencies”  In _Living Rev. Rel._ 24.1, 2021, pp. 4  DOI: 10.1007/s41114-021-00032-5
  * [112] Asher Berlin et al.  “Detecting high-frequency gravitational waves with microwave cavities”  In _Phys. Rev. D_ 105.11, 2022, pp. 116011  DOI: 10.1103/PhysRevD.105.116011
  * [113] Asuka Ito and Jiro Soda  “Exploring high-frequency gravitational waves with magnons”  In _Eur. Phys. J. C_ 83.8, 2023, pp. 766  DOI: 10.1140/epjc/s10052-023-11876-2
  * [114] Nicolas Herman, Léonard Lehoucq and André Fúzfa  “Electromagnetic antennas for the resonant detection of the stochastic gravitational wave background”  In _Phys. Rev. D_ 108.12, 2023, pp. 124009  DOI: 10.1103/PhysRevD.108.124009
  * [115] Kaloian D. Lozanov and Volodymyr Takhistov  “Enhanced Gravitational Waves from Inflaton Oscillons”  In _Phys. Rev. Lett._ 130.18, 2023, pp. 181002  DOI: 10.1103/PhysRevLett.130.181002
  * [116] Karsten Jedamzik, Martin Lemoine and Jerome Martin  “Collapse of Small-Scale Density Perturbations during Preheating in Single Field Inflation”  In _JCAP_ 09, 2010, pp. 034  DOI: 10.1088/1475-7516/2010/09/034
  * [117] Richard Easther, Raphael Flauger and James B. Gilmore  “Delayed Reheating and the Breakdown of Coherent Oscillations”  In _JCAP_ 04, 2011, pp. 027  DOI: 10.1088/1475-7516/2011/04/027
  * [118] Jérôme Martin, Theodoros Papanikolaou, Lucas Pinol and Vincent Vennin  “Metric preheating and radiative decay in single-field inflation”  In _JCAP_ 05, 2020, pp. 003  DOI: 10.1088/1475-7516/2020/05/003
  * [119] Tomas Andrade  “GRChombo: An adaptable numerical relativity code for fundamental physics”  In _J. Open Source Softw._ 6.68, 2021, pp. 3703  DOI: 10.21105/joss.03703
  * [120] Xiao-Xiao Kou, Chi Tian and Shuang-Yong Zhou  “Oscillon Preheating in Full General Relativity”  In _Class. Quant. Grav._ 38.4, 2021, pp. 045005  DOI: 10.1088/1361-6382/abd09f
  * [121] Josu C. Aurrekoetxea, Katy Clough and Francesco Muia  “Oscillon formation during inflationary preheating with general relativity”  In _Phys. Rev. D_ 108.2, 2023, pp. 023501  DOI: 10.1103/PhysRevD.108.023501
  * [122] Daniel del-Corral, Paolo Gondolo, K. Kumar and João Marto  “Revisiting primordial black holes formation from preheating instabilities: the case of Starobinsky inflation”, 2023  arXiv:2311.02754 [astro-ph.CO]
  * [123] Eemeli Tomberg and Hardi Veermäe  “Tachyonic preheating in plateau inflation”  In _JCAP_ 12.12, 2021, pp. 035  DOI: 10.1088/1475-7516/2021/12/035
  * [124] Niko Koivunen, Eemeli Tomberg and Hardi Veermäe  “The linear regime of tachyonic preheating”  In _JCAP_ 07.07, 2022, pp. 028  DOI: 10.1088/1475-7516/2022/07/028
  * [125] Stefan Antusch, Kenneth Marschall and Francisco Torrenti  “Characterizing the post-inflationary reheating history. Part II. Multiple interacting daughter fields”  In _JCAP_ 02, 2023, pp. 019  DOI: 10.1088/1475-7516/2023/02/019


# Efficient Subgraph Isomorphism using Graph Topology

Arpan Kusari and Wenbo Sun††footnotemark:
University of Michigan Transportation Research Institute,
University of Michigan
Ann Arbor
Authors contributed equally. Names listed in alphabetical order

###### Abstract

Subgraph isomorphism or subgraph matching is generally considered as an NP-
complete problem, made more complex in practical applications where the edge
weights take real values and are subject to measurement noise and possible
anomalies. To the best of our knowledge, almost all subgraph matching methods
utilize node labels to perform node-node matching. In the absence of such
labels (in applications such as image matching and map matching among others),
these subgraph matching methods do not work. We propose a method for
identifying the node correspondence between a subgraph and a full graph in the
inexact case without node labels in two steps - (a) extract the minimal unique
topology preserving subset from the subgraph and find its feasible matching in
the full graph, and (b) implement a consensus-based algorithm to expand the
matched node set by pairing unique paths based on boundary commutativity.
Going beyond the existing subgraph matching approaches, the proposed method is
shown to have realistically sub-linear computational efficiency, robustness to
random measurement noise, and good statistical properties. Our method is also
readily applicable to the exact matching case without loss of generality. To
demonstrate the effectiveness of the proposed method, a simulation and a case
study is performed on the Erdos-Renyi random graphs and the image-based affine
covariant features dataset respectively.

##  1 Introduction

In recent years, the number of applications employing graphs have increased
exponentially, with the scientific community using graphs to denote
relationships between the features of interest (Goyal and Ferrara,, 2018). A
fundamental problem arising in these applications is the subgraph isomorphism
problem, also known as subgraph matching, whereby an optimal correspondence is
sought between nodes of two given graphs, where one graph is derived from the
other. We note here that graph matching is a specialized case of the subgraph
matching problem when two graphs are of the same size. Optimality differs
according to problem description and domain, but generally refers to the
alignment of global graph structures, along with other node features when
available.

Ullman in his seminal paper (Ullmann,, 1976) showed that isomorphism can
always be estimated using brute-force enumeration as a default strategy which
amounts to a depth-first search. The backtracking strategy devised in the
paper is used to significantly reduce the size of the search space. Cordella
et al., (2004) proposed an improvement of Ullman’s algorithm, the VF2
algorithm, by utilizing a data structure during exploration to significantly
reduce memory requirements along with five feasibility rules for pruning the
search tree. Both the above approaches solve the exact matching problem where
the corresponding nodes and edges between the two graphs share same weights.
Conversely, in many practical applications, the nodes can be corrupted by
measurement noise and edges can be absent which renders exact matching
ineffective. In this case, the practitioners rely not on the exact
correspondence but rather some similarity measures. The first set of inexact
matching algorithms perform a series of editing operations such as node and
edge insertion, deletion or substitution, in order to attain an exact matching
(Cordella et al.,, 1998; Tsai and Fu,, 1983; Shapiro and Haralick,, 1985). As
is understood, in the case of subgraphs having a much lesser nodes than the
full graphs, these methods are really inefficient in solving the matching.

The second set of inexact matching algorithms formulate and solve an
optimization algorithm to obtain the node correspondence. Graduated assignment
algorithm (Gold and Rangarajan,, 1996) combines graduated nonconvexity, two-
way (assignment) constraints,and sparsity to match graphs even in the presence
of noise. Egozi et al., (2012) provides a probabilistic approach to spectral
matching which provides a computationally efficient approach towards graph
matching. Suh et al., (2015) augments the integer quadratic problem (IQP)
objective with a compactness prior to reduce the number of outliers matched.
They utilize a Markov chain Monte Carlo (MCMC) to solve the optimization
problem. They show that using these modifications can help improve the
matches. These approaches require construction of a weight matrix between the
vertices of the two graphs and are thus, constrained by small graph sizes and
pairwise interactions.

There have been some recent work in geometry preserving linear node
assignments. Zaslavskiy et al., (2008) proposed an approximate method as a
quadratic assignment problem over the set of permutation matrices, where a
permutation matrix is a binary matrix describing whether the particular nodes
are matched or not. There have been other approaches involving constructing
permutation matrices under specific real-valued approximations (Cho et al.,,
2010; Zhou and De la Torre,, 2012) but they have problems where there are poor
geometric alignments with high deformations. Byrne, (2013) proposed a graph
matching technique which preserves the global topological structure by finding
an optimal simplicial chain map. Although the author provides an explanation
of how this technique can be extended for higher topologies, the formulation
shows a matching between the nodes and the edges between the nodes to
constrain the homology. For these approaches, even though the topology is
constrained as a result of the optimization, the topology constraint is not
inherently employed in the matching.

In this work, we specifically focus on graph matching problems with
considering measurement noise on edge weights without the presence of node
labels, which are demanded in a variety of applications (Bunke,, 2000) such as
character recognition (Lu et al.,, 1991; Rocha and Pavlidis,, 1994), shape
analysis (Cantoni et al.,, 1998; Lourens,, 1998), chemical structure analysis
(Balaban,, 1985), etc. The idea of matching topology between subgraph and full
graph is inspired by the random sample consensus (RANSAC) literature (Fischler
and Bolles,, 1981; Derpanis,, 2010). Graph topology provides a description of
graph connectivity as the basic structure of a graph. Our intuition is that
matching the unique minimal unit that preserves a pre-specified graph topology
between the subgraph and full graph improves the matching efficiency and
robustness. Post to the topology-based matching, the remaining unmatched edges
in the subgraph are further matched to the corresponding edges in the full
graph. However, there are three fundamental questions that arise - (a) what is
the minimal topology-preserving unit? (b) how do we guarantee the matching
uniqueness? (c) how do we match the remaining edges after the initial match
has occurred?

The primary contributions of this work are:

  * •

Inexact matching in absence of node labels Almost all of the subgraph matching
methods rely on node labels to provide some information regarding the nodes.
In absence of such information, these methods fail to find feasible subgraph
matches. Our method therefore, fills this gap and provides a matching method
which does not need to incorporate labels for matching.

  * •

Globally consistent matching in the presence of large amounts of
outliers/deformation. Almost all subgraph matching algorithms rely on a single
optimization step to search for the node correspondence minimizing the
matching loss, which may lead to invalid matching since the true matching may
not result in the minimal matching loss under measurement noises. Our proposed
approach provides robust global matching under large noise by defining a
graph-topology-based structure in the subgraph and searching for feasible
matches of nodes sharing the same graph-topology-based structure between the
full graph and subgraph. Given the feasible initial match, the consensus based
searching expands the search for new matched nodes in a breadth-first tree
approach to collect the number of nodes which can be matched uniquely. This
formulation takes into account the noise in the edges by statistically
modeling them and proving theoretically that this can converge to the true
matching.

  * •

Providing efficiency in subgraph matching. Using a hypothesis based searching
leads to a fast convergence whereby a wrong initial matching leads to a weaker
final match set. Parallelizing the consensus based searching step can yield
further efficiencies without sacrificing accuracy.

The rest of the paper is organized as follows: Section 2 presents the
mathematical formulation of the graph matching problem and our proposed
approach. Section 3 provides a robust simulation study to compare the
strengths of our algorithm and in section 4, we evaluate performance on
homology preserving subgraph matching for image registration applications.
Section 5 provides the conclusions.

##  2 Method

###  2.1 Problem formulation

We start with a formal definition of graphs. A graph is denoted by
$G=\left(V,E,w\right)$, where $V$, $E$ and $w$ represent the set of nodes, the
set of edges and edge weights respectively. Let $|\cdot|$ denote the number of
elements in a set. Each node is denoted by a node index
$v\in\left\\{1,...,|V|\right\\}$. Each edge is expressed as $e=(i,j)$
indicating that nodes $i$ and $j$ are connected by an edge. The edge weight
$w:(i,j)\rightarrow\mathbb{R}$ is a function mapping an edge to the real-
valued weight.

In our problem setting, let $G_{f}=\left(V_{f},E_{f},w\right)$ and
$G_{s}=\left(V_{s},E_{s},w\right)$ denote the full graph and subgraph of
interest, respectively. Here the subscript “f” represents “full” while the
subscript “s” represents “sub”. The objective of the subgraph matching problem
is to seek a matching policy $\phi$, that projects each node $v\in V_{s}$ to
the corresponding node $\phi(v)\in V_{f}$. In some graph matching literature,
the policy function can be also expressed as a $|V_{f}|\times|V_{s}|$
permutation matrix $X$, whose $(i,j)$-th element equals to $1$ if $\phi(j)=i$.
Without ambiguity, for any $e=(i,j)\in E_{s}$, we define
$\phi(e)=\left(\phi(i),\phi(j)\right)\in E_{f}$ as the matched edge of $e$. In
this paper, we consider the situations where weights in the full graph are the
ground truth, and the weights in the subgraph are noisy observations from the
full graph, that is to assume:

$w(e_{s})=w\left(\phi(e_{s})\right)+\epsilon,$ (1)

where $\epsilon$ is the random noise following normal distribution
$N(0,\sigma^{2})$ with unknown variance $\sigma^{2}$. We assume $\epsilon$’s
are independent among different pairs of matched edges.

Subgraph matching can be formulated as an optimization problem. Under a pre-
specified distance measure $d(\cdot,\cdot)$ of edges, the overall matching
loss can be quantified as a $|V_{f}||V_{s}|\times|V_{f}||V_{s}|$ matrix $Q$
whose $(i+|V_{s}|j,i^{\prime}+|V_{s}|j^{\prime})$-th element equals to
$d\left((i,i^{\prime}),(j,j^{\prime})\right)$. Let $X_{v}$ denote the
vectorization of $X$ by stacking its columns on top of one another. Let
$X_{i}$ and $X_{(j)}$ denote the $i$-th row and $j$-th column vectors of $X$,
respectively. Let $\mathbb{I}$ be the vector of ones. The objective function
is written as:

$\widehat{X}_{v}=arg\,\min_{X_{v}}X_{v}^{T}QX_{v}$ (2)

subject to

$\displaystyle\mathbb{I}^{T}X_{i}$ $\displaystyle=$ $\displaystyle 1$
$\displaystyle X_{(j)}^{T}\mathbb{I}$ $\displaystyle=$ $\displaystyle 1,$

The exact solution to the optimization problem in Eq.(2) can be obtained
through an NP-complete non-convex integer programming. The problem can be
solved via a searching algorithm (e.g. the Hungarian algorithm) or
approximately solved via an integer programming solver (e.g. interior point
optimizer). There are two main limitations of these existing approaches. (i)
They require a computational time of $O(|V_{f}|^{3})$, which often exceeds the
computation budget for large-scale graphs. (ii) The existing methods are
designed based on the loss of edge-to-edge matching. Simply minimizing the
overall matching loss in Eq.(2) may lead to erroneous results due to the
measurement noise.

To overcome the above limitations, we propose a consensus-based approach that
performs graph matching on a set of nodes with topological constraints, which
introduces few key innovative components. First, topological constraints are
employed for matching a set of nodes. Compared to the edge weights, the
connectivity between nodes imposed by the topological constraints are
invariant to the measurement noise with probability one (in the case a
positive-valued edge weight is disturbed to zero). Moreover, simultaneously
matching multiple nodes which satisfy a topological constraint can reduce
computational time while shrinking the probability of mismatching due to the
measurement noise. Second, instead of searching for the unique solution to
Eq.(2) in a single step, we would like to bind the feasible matching policies
in a feasible set based on the revised version of Eq.(2) subject to some
topological constraints and then shrink the feasible set with accumulating
more nodes to the topology structure until a unique matching policy is
reached. In this way, the proposed method will be able to cover the true
matching policy in the feasible region with a high probability and identify it
as the unique solution when a sufficient number of nodes are introduced to the
topology structure. It is required to define the topological constraint and
feasible set to complete the proposed method. While we do not formally show
the proof, the exact matching case can be easily shown to be a special case of
our inexact matching algorithm. We will elaborate the technical details in the
following subsections.

###  2.2 Topology-based matching

In the first step, we aim to match a subset of $V_{s}$ to $V_{f}$, which is
defined as a topology-preserving unit. Different from the node-to-node
matching in Eq.(2), the graph matching is conducted based on multiple mutually
exclusive subsets of $V_{s}$ following a topology concept - p-simplex, which
is defined as the convex hull of $p+1$ affinely independent nodes. For
example, a $2$-simplex represents a triangle, and a $3$-simplex represents a
tetrahedron in a graph. The topology-preserving unit in the initial matching
is defined as a pair of connected p-simplexes in the graph along with the
connecting path. To guarantee the validity of the topology-preserving unit,
two basic assumptions are made while formulating our approach.

  * •

$G_{s}$ is a connected graph i.e. there do not exist two or more separate
clusters of nodes without edges in between, and

  * •

there exists more than two p-simplexes in $G_{s}$.

An illustrative example of a topology-preserving unit is displayed in Figure
1, where the two 2-simplexes in the subgraph on the left panel (nodes
$\\{3,8,9\\}$ and nodes $\\{6,11,12\\}$) and the shortest path in between
(nodes $\\{8,13,12\\}$) are considered as the topology-preserving unit and
matched with the corresponding nodes in the full graph. Let $\delta_{s,c}$
denote the topology-preserving unit consisting of $c$ nodes, the next step is
to seek a matching policy $\phi$ or equivalently a vectorized matching matrix
giving the $c$ matching nodes in $G_{f}$.

Figure 1: Illustration of the graph topology in initial matching. Left panel:
subgraph $G_{s}$ with two $2$-simplexes and the path highlighted in red. Right
panel: full graph $G_{f}$ with the unique matching simplexes highlighted in
orange.

Here we define a general feasible set of matching policies for any given
subset of $V_{s}$, which can be naturally applied to $\delta_{s,c}$. Let
$V_{c}$ denote any subset of $V_{s}$ with $c$ nodes. Let $Q(V_{c})$ denote the
distance measures between $V_{c}$ and the subsets in $V_{f}$ with the same
node size. The feasible set of vectorized matching matrix
$\widetilde{X}_{v,V_{c}}$ is defined by revising Eq.(2) as:

$\widetilde{X}_{v}(V_{c})=\left\\{X_{v,V_{c}}\big{|}X_{v,V_{c}}^{T}Q(V_{c})X_{v,V_{c}}\leq\tau_{c}\right\\},$
(3)

where $\tau_{c}$ is a pre-specified threshold determined by the set size $c$.
Note that $Q(V_{c})$ is a vector matching the subset $V_{c}$ towards all the
subsets of the same node size in $V_{f}$. Revising the optimization problem
into a feasible set searching problem enhances the robustness of the matching
under measurement noise. Under the statistical property developed in the later
subsection, the number of elements in $\widetilde{X}_{v}(V_{c})$ will shrink
as the set size $c$ increases. For a $c$ large enough, a unique element will
be eventually identified.

To improve the efficiency in calculating $\widetilde{X}_{v}(V_{c})$, we impose
the topology constraints on the elements in $\widetilde{X}_{v}(V_{c})$ to have
the same topology structure as in $V_{c}$. Let $E_{V_{c}}$ denote the edges in
$V_{c}$. Let $\phi_{X_{v,V_{c}}}$ denote the matching policy corresponding to
$X_{v,V_{c}}$. The topology constraint on $X_{v,V_{c}}$ is expressed in:

$\mathcal{T}_{v}(V_{c})=\left\\{X_{v,V_{c}}\big{|}\left(\phi_{X_{v,V_{c}}}(v_{1}),\phi_{X_{v,V_{c}}}(v_{2})\right)\in
E_{f}\iff\left(v1,v2\right)\in E_{V_{c}}\right\\}.$ (4)

Combining Eq.(3) and Eq.(4) yields the feasible set:

$\widetilde{X}_{v}(V_{c})=\left\\{X_{v,V_{c}}\in\mathcal{T}_{v}(V_{c})\big{|}X_{v,V_{c}}^{T}Q(V_{c})X_{v,V_{c}}\leq\tau_{c}\right\\},$
(5)

which can be used to search for $\widetilde{X}_{v}(\delta_{s,c})$.

The feasible set is well defined once $Q(V_{c})$ and $\tau_{c}$ are defined.
Here $Q(V_{c})$ is calculated based on the distance measure between $V_{c}$
and any given subset $v\subset V_{f}$. The distance measure is defined on the
differences in the averaged edge weights, that is

$d\left(V_{c},v\right)=\left|\frac{1}{c}\sum_{i\neq j\in
V_{c}}w_{s}(i,j)-\frac{1}{c}\sum_{i\neq j\in v}w_{f}(i,j)\right|.$ (6)

The averaged edge weight is invariant to possible rotation and scaling in real
graphs. Note that when all the edge weights in Eq.(6) are pairwise matched,
the distance measure $d$ is the arithmetic mean of $c$ independent normal
random variables, which shrinks to zero as $c$ grows large according to the
central limit theorem and thus guarantees the identification of the true
matching policy. Under this statistical property, we propose to set the
threshold $\tau_{c}$ as:

$\tau_{c}=\frac{\Phi^{-1}(1-\alpha/2)}{\sqrt{c}}\sigma,$ (7)

where $\Phi^{-1}(1-\alpha/2)$ represents the $1-\alpha/2$ quantile of the
standard normal distribution. Here $1-\alpha$ is interpreted as the coverage
rate of the true matching by the designed feasible set, and can be specified
by the users based on the requirement on the coverage rate. An estimate of
$\sigma$ should be provided when it is not available to the user, which will
be discussed in Section 3.

The first step is implemented via a depth-first searching of feasible elements
based on Eq.(5), as described in Algorithm 1. The topology structure and
shortest path information are used to prune the search. In particular, the
following searching algorithm is designed to find all feasible matching
policies for the given topology-preserving units, which will serve as a
foundation for the consensus-based searching elaborated in the next
subsection.

input : The subgraph $G_{s}$ and the full graph $G_{f}$, initial parameters
$\alpha$ and $\sigma$, number of tries $n$

initialize  _$\delta_{s,c}=\left\\{\right\\}$
,$\widetilde{X}_{v}(\delta_{s,c})=\left\\{\right\\}$_

for  _$i\in 1\dots n$_ do

Randomly draw a pair of p-simplexes $\nu_{1}\neq\nu_{2}$ from $G_{s}$

Calculate the shortest path $p(\nu_{1},\nu_{2})$ between $\nu_{1}$ and
$\nu_{2}$

$V_{c}\leftarrow\nu_{1}\cup p(\nu_{1},\nu_{2})\cup\nu_{2}$

$\tilde{\nu}_{1,f}\leftarrow$ all feasible matches to $\nu_{1}$ from $G_{f}$
based on Eq.(5)

$\tilde{\nu}_{2,f}\leftarrow$ all feasible matches to $\nu_{2}$ from $G_{f}$
based on Eq.(5)

for  _$(\nu_{1,f},\nu_{2,f})\in\tilde{\nu}_{1,f}\times\tilde{\nu}_{2,f}$_ do

for  _Path $p(\nu_{1,f},\nu_{2,f})$ between $\nu_{1,f}$ and $\nu_{2,f}$ having
the same length as $p(\nu_{1},\nu_{2})$_ do

$v=\nu_{1,f}\cup p(\nu_{1,f},\nu_{2,f})\cup\nu_{2,f}$

if  _$d(v,V_{c})\leq\tau_{c}$_ then

$\delta_{s,c}\leftarrow\delta_{s,c}\cup V_{c}$

$\phi_{X_{v,V_{c}}}(V_{c})\leftarrow v$

Calculate $X_{v,V_{c}}$ based on $\phi_{X_{v,V_{c}}}(V_{c})$

for  _$X_{v}(\delta_{s,c})\in\widetilde{X}_{v}(\delta_{s,c})$_ do

$X_{v}(\delta_{s,c})\leftarrow X_{v}(\delta_{s,c})\cup X_{v,V_{c}}$

end for

end if

end for

end for

end for

output : Topology-preserving unit $\delta_{s,c}$ and the feasible matching
$\widetilde{X}_{v}(\delta_{s,c})$

Algorithm 1 Topology-based matching

###  2.3 Consensus-based searching

After the initial matching of the topology preserving nodes in $\delta_{s,c}$,
a consensus-based algorithm is developed to sequentially match the remaining
nodes in $V_{s}\backslash\delta_{s,c}$. The basic concept is to seek for the
feasible matching of paths whose starting nodes are in $\delta_{s,c}$ and the
rest nodes are in $V_{s}\backslash\delta_{s,c}$. Instead of matching the
topology-preserving unit in the previous step, we aim to match paths of length
$q$ based on the feasible set defined in Eq.(5). The feasible set can be
narrowed down as the path length $q$ grows large. We start with an empty node
set $\eta_{s,v}$ with $v=0$ and expand the set by searching for feasible
matching of paths in $V_{s}\backslash\delta_{s,c}$ until it contains a unique
element. The algorithm is terminated when no new nodes can be matched between
$G_{s}$ and $G_{f}$. Specifically, we propose the following consensus-based
searching algorithm, Algorithm 2.

input : Topology-preserving unit $\delta_{s,c}$, unmatched node set
$V_{s}\backslash\delta_{s,c}$, initial parameters $\alpha$ and $\sigma$

initialize  _$v=0$ , $\eta_{s,v}=\left\\{\right\\}$,
$\widetilde{X_{v}}(\eta_{s,v})=\left\\{\right\\}$_

while  _$\eta_{s,v}\neq V_{s}\backslash\delta_{s,c}$_ do

Sample an edge $e_{s}\in E_{s}$ incidence to a node in $\eta_{s,v}$ and a node
$a_{s}\in V_{s}\backslash\eta_{s,v}$

$\eta_{s,v}\leftarrow\eta_{s,v}\cup\\{v:v\in e_{s}\\}$

$v\leftarrow v+1$

$v^{\prime}\leftarrow 0$

$\tilde{\eta}_{s,v^{\prime}}=\left\\{\right\\}$

while  _$\exists$ an edge $e_{s}^{\prime}\in E_{s}$ incidence to $a_{s}$ and a
node $b_{s}\in V_{s}\backslash\eta_{s,v}$_ do

$\tilde{\eta}_{s,v^{\prime}}\leftarrow\tilde{\eta}_{s,v^{\prime}}\cup
e_{s}^{\prime}$

$v^{\prime}\leftarrow v^{\prime}+1$

$a_{s}\leftarrow b_{s}$

Compute feasible matching for $\tilde{\eta}_{s,v^{\prime}}$ in $V_{f}$
according to Eq.(5), store it in
$\widetilde{X_{v}^{\prime}}(\eta_{s,v^{\prime}})$

if  _$\widetilde{X_{v}^{\prime}}(\eta_{s,v^{\prime}})$ is unique_ then

$\eta_{s,v}\leftarrow\eta_{s,v}\cup\tilde{\eta}_{s,v^{\prime}}$

$v\leftarrow v+v^{\prime}$

$\widetilde{X_{v}}(\eta_{s,v})\leftarrow\widetilde{X_{v}}(\eta_{s,v})\cup\widetilde{X_{v}^{\prime}}(\eta_{s,v^{\prime}})$

break

end if

end while

end while

output : Node correspondence for $\widetilde{X_{v}}(\eta_{s,v})$

Algorithm 2 Consensus-based searching

###  2.4 Theoretical Guarantees

Theorems are developed to evaluate the false positive and false negative rates
in Subsection 2.4. The first theorem shows that the feasible sets in Eq.(5)
will cover the true matching nodes with a pre-specified probability
$1-\alpha$, which guarantees the identification of the true matching nodes by
applying the proposed method.

######  Theorem 1.

Let $X_{v}^{*}(V_{c})$ denote the true vectorized matching matrix
corresponding for the node set $V_{c}$. The feasible set covers the true
matching nodes with probability $1-\alpha$, i.e.,

$\mathbb{P}\left[X_{v}^{*}(V_{c})\in\widetilde{X}_{v}(V_{c})\right]=1-\alpha.$
(8)

Proof: Let $V_{c,f}^{*}$ denote the true matching node set for $V_{c}$.
according to Eq.(5),

$\displaystyle\mathbb{P}\left[X_{v}^{*}(V_{c})\in\widetilde{X}_{v}(V_{c})\right]$
$\displaystyle=$
$\displaystyle\mathbb{P}\left[d(V_{c},V_{c,f}^{*})\leq\tau_{c}\right]$ (9)
$\displaystyle=$
$\displaystyle\mathbb{P}\left[\frac{1}{c}\left|\frac{1}{c}\sum_{i\neq j\in
V_{c}}w_{s}(i,j)-\frac{1}{c}\sum_{i\neq j\in
V_{c,f}^{*}}w_{f}(i,j)\right|\leq\tau_{c}\right]$ $\displaystyle=$
$\displaystyle\mathbb{P}\left[\left|\frac{1}{c}\sum_{j=1}^{c}\epsilon_{j}\right|\leq\tau_{c}\right],$

where $\epsilon_{j}$’s are $c$ independent and normally distributed random
variables defined in Eq.(1). Note that $\frac{1}{c}\sum_{j=1}^{c}\epsilon_{j}$
also follows a normal distribution $N(0,\sigma^{2}/c)$. Thus we have

$\displaystyle\mathbb{P}\left[X_{v}^{*}(V_{c})\in\widetilde{X}_{v}(V_{c})\right]$
$\displaystyle=$ $\displaystyle 2\Phi(\sqrt{c}\tau_{c}/\sigma)-1$ (10)
$\displaystyle=$ $\displaystyle 2(1-\alpha/2)-1$ $\displaystyle=$
$\displaystyle 1-\alpha.$

On the other hand, we would like the proposed feasible set to exclude any
irrelevant node set, especially when the number of nodes $c$ grows large. The
theorem is formally established as follows.

######  Theorem 2.

Let $\mu=\inf_{(e_{s},e_{f})\in E_{s}\times
E_{f}}\left|\mathbb{E}[w(e_{s})]-\mathbb{E}[w(e_{f})]\right|/\sigma>0$, where
$e_{s}$ and $e_{f}$ are not matched. The probability of excluding an
irrelevant vectorized matching matrix $X_{v}$ from $\widetilde{X}_{v}(V_{c})$,

$\mathbb{P}\left[X_{v}\notin\widetilde{X}_{v}(V_{c})\right]>\Phi\left(\frac{-\Phi^{-1}(1-\alpha/2)-\mu\sqrt{c}}{\sqrt{2c}}\right)+1-\Phi\left(\frac{\Phi^{-1}(1-\alpha/2)-\mu\sqrt{c}}{\sqrt{2c}}\right),$
(11)

Proof: given $X_{v}$ is not a true vectorized matching matrix for $V_{c}$,
there exists at least a pair of edges that are not matched. Let
$(e_{1},e_{2})$ denote the unmatched pair of edges, then we have

$\displaystyle\mathbb{P}\left[X_{v}\notin\widetilde{X}_{v}(V_{c})\right]$
$\displaystyle\geq$
$\displaystyle\mathbb{P}\left[\left|w(e_{1})-w(e_{2})\right|>\tau_{c}\right]$
(12) $\displaystyle\geq$
$\displaystyle\mathbb{P}\left[\left|Z_{1,2}\right|>\tau_{c}\right],$

where $Z_{1,2}$ is a normally distributed random variable with mean
$\mu\sigma$ and variance $2\sigma^{2}$. Therefore we have

$\displaystyle\mathbb{P}\left[\left|Z_{1,2}\right|>\tau_{c}\right]$
$\displaystyle=$
$\displaystyle\Phi\left(\frac{-\tau_{c}-\mu\sigma}{\sqrt{2}\sigma}\right)+1-\Phi\left(\frac{\tau_{c}-\mu\sigma}{\sqrt{2}\sigma}\right)$
(13) $\displaystyle=$
$\displaystyle\Phi\left(\frac{-\Phi^{-1}(1-\alpha/2)-\mu\sqrt{c}}{\sqrt{2c}}\right)+1-\Phi\left(\frac{\Phi^{-1}(1-\alpha/2)-\mu\sqrt{c}}{\sqrt{2c}}\right).$

Note that $\mathbb{P}\left[\left|Z_{1,2}\right|>\tau_{c}\right]$ is a
monotonic increasing function with $c$, implying that the proposed feasible
can exclude the irrelevant nodes when $c$ grows large.

###  2.5 Computational complexity

The computational complexity of the proposed method can be estimated as
follows. Suppose $k$ $p$-simplexes are initially matched. Let $\bar{c}$ denote
the average number of nodes in each $\delta_{s,c}$, $m_{f}$ denote the total
number of $p$ simplexes in $G_{f}$, and $\bar{d}_{f}$ denote the average
degree in $G_{f}$. The computational cost in the initial matching step is
$\mathcal{O}\left(m_{f}k\bar{d}_{f}\left(\bar{c}-2(p+1)\right)\right)$, where
the multiplication of the first two terms evaluates the computation time for
matching the $p$-simplexes while the last two terms quantifies the computation
time for matching the shortest path in between. The computation time for the
consensus step is dominated by the efforts in path-wise matching of the
remaining nodes after the initial matching, which is
$\mathcal{O}\left(\left(n_{s}-k\left(p+\bar{c}-1\right)\right)\bar{d}_{f}\right)$.

The total computation complexity of the proposed algorithm is expressed in

$\mathcal{O}\left(m_{f}k\bar{d}_{f}\left(\bar{c}-2(p+1)\right)+\left(n_{s}-k\left(p+\bar{c}-1\right)\right)\bar{d}_{f}\right).$
(14)

It is worth noting that Eq(14) is upper bounded by
$\mathcal{O}\left((n_{f}+m_{f})\bar{d}_{f}\right)=\mathcal{O}\left(n_{f}^{p}\right)$
in the worst case scenario where the subgraph and full graph have the same
size and both of them are fully connected. In this case, the number of
p-simplexes is $m_{f}=\mathcal{O}(n_{f}^{p})$, which greatly increase the
computational load. However, this case is not commonly seen in real
applications where the graphs are relatively sparse. For a sparsely connected
graph with a smaller $m_{f}$ and $\bar{d}_{f}$ and $n_{s}\ll n_{f}$, the
computational time is sub-linear with respect to $n_{f}$.

##  3 Simulation study

The proposed algorithm is evaluated in a simulation dataset generated from the
Erdős-Rényi model (Erdős et al.,, 1960). Specifically, the full graph is
generated as an Erdős-Rényi graph with $100$ nodes. Edges are connected
between each pair of nodes with a probability of $0.1$ to simulate graphs in
real applications where edges are sparsely presented. Each edge in the full
graph is assigned with a random weight sampled from the uniform distribution
$U(0,1)$. The subgraph is consisted of $20$ nodes randomly sampled from the
full graph. To ensure some edge connectivity in the subgraph, neighboring
nodes whose pair-wise distance is smaller than $0.5$ are included. Each edge
in the subgraph is then assigned with a normally distributed noise with zero
mean and standard deviation of $\sigma$. The algorithm performance will be
tested under different noise level $\sigma$.

A Monte-Carlo simulation study of $100$ iterations is conducted for each given
$\sigma\in\left\\{0.001,0.002,\dots,0.01\right\\}$. In each iteration, both a
full graph and a subgraph is generated according to the procedure described
above. The standard deviation $\sigma$ is presumed to be known, and $\alpha$
is set as $0.025$. The proposed algorithm takes the full graph and subgraph as
the input, and returns the node correspondence. When implementing the proposed
method, we set $p=2$ to focus the topology matching on triangles. This is
because a smaller $p$ provides a small searching space for all the feasible
$p-$simplex and hence improves the algorithm efficiency. A discussion on the
trade-off between the robustness of the method and the algorithm efficiency
can be found in Section 5.

The resultant node correspondence is compared with the true node
correspondence, and the percentage of the correctly matched nodes is reported
to evaluate the algorithm accuracy. The computational time in each iteration
is also recorded to evaluate the algorithm efficiency. The distributions of
the accuracy and computational time under the Monte-Carlo simulation are
summarized in Figure 2. We observe an decreasing trend of the matching
accuracy as the noise level increases as expected. The averaged accuracy is
above $95\%$ when the signal-to-noise ratio is above $20$ (corresponding to
$\sigma=0.05$). Regarding the computational time, it takes $4$ seconds on
average to seek a node correspondence between a full graph of $100$ nodes and
a subgraph of $20$ nodes, implying the relatively high efficiency of the
proposed method. We would like to mention that we tested four different
subgraph matching algorithms provided in the OpenGraphMatching toolbox
(Tianchang,, 2021)(GraphQL(He and Singh,, 2008), CECI (Bhattarai et al.,,
2019), Neural (Liu et al.,, 2020), QuickSI (Shang et al.,, 2008)) which
yielded no results without labels. Also VF2++ (Jüttner and Madarasi,, 2018), a
modification of VF2, provided no results with the addition of noise.

Figure 2: Performance of the proposed method under Monte-Carlo simulation. Top
panel: box plot of the matching accuracy under different noise levels. Bottom
panel: box plot of the computational time under different noise levels. The
upper and lower bounds of the boxes depict the $75\%$ and $25\%$ quantiles,
respectively. The top and bottom horizontal lines represent the $95\%$ and
$5\%$ quantiles, respectively. The scattered dots illustrate the outliers.

##  4 Case study

We investigate the performance of the proposed algorithm for image
correspondence using sampled images from the affine covariant features dataset
(Mikolajczyk and Schmid,, 2005). This dataset has been used as a traditional
dataset for local feature detector evaluation e.g. Scale Invariant Feature
Transform (SIFT) (Lowe,, 1999), Gradient location-orientation histogram (GLOH)
(Mikolajczyk and Schmid,, 2005) etc. Image correspondence or homography is
used to infer the relationship between features of two or more images of the
same planar object. As the first example, we create a subset of a given image
by cropping part of the original image as shown in the top left panel of Fig
3. Emulating other image graph extraction literature, we extracted keypoints
from the images using the SIFT feature extraction. SIFT is known to be robust
to changes in image scale, noise and illumination and thus can provide nodes
which have same relative positions for the same object irrespective of changes
in camera orientations. We utilize constrained Delaunay triangulations to
connect nodes in order to form the graphs. The edge weights are given by the
Euclidean distances between the nodes. Constraining the triangulations removes
spurious triangles which increase computational load without adding new
information. We show the graphs superimposed on the original and cropped image
in the top left panel of Fig. 3. The bottom left panel of Fig. 3 shows the
true matches between the nodes in the original and cropped image. We estimate
the node correspondence by getting the coordinates of the keypoints in the
cropped image and then offsetting them by the cropped origin coordinates. We
then search for keypoints in the original image which lie within certain
bounds and are not matched by any other keypoints.

The top right panel of Fig. 3 gives the estimated correspondences using the
proposed method. We choose a small set of triangles randomly in the subgraph
which have paths connecting them. Based on the matching performance in the
initial matching step, we tune the hyperparameter $\sigma$ and find that
$\sigma=1$ pixel results in the optimal matching performance. The feasible
sets are then constructed based on the constraints with the tuned $\sigma$.
Given there exists a feasible set, we move to the consensus step to
sequentially find nodes which are connected to the matched nodes and can be
uniquely matched to nodes in the full graph. From a visual inspection, we see
that even in presence of errors in SIFT keypoints, our method can robustly
match nodes in a topologically consistent manner. There are couple of node
matches at the top of the images which are missed due to the errors being
outside our bounds. In the bottom right panel of Fig. 3, we rotate the cropped
image by $30\degree$ clockwise and estimate the matching. Given the rotational
invariance in weights of the $2$-simplex, we see that the proposed method can
match the nodes even in this scenario. In this case though, some nodes in the
right side of the cropped image are unmatched due to keypoints being absent
and thus the connectivity assumption of the graph not being met completely.

Figure 3: Top left panel: The sampled graffiti image in affine covariant
features dataset and cropped image. The green lines depicts the constrained
Delaunay triangulations of SIFT features. Bottom left panel: The true image
correspondence between the two images. The red lines highlight the node
correspondence. Top right panel: The estimated image correspondences from the
proposed method. Bottom right panel: The estimated image correspondences
between the original image and rotated cropped image from the proposed method.

##  5 Conclusions

We have presented a subgraph matching algorithm, providing theoretical
guarantees about the accuracy of the method while looking at matching graph
topology as a way to reduce the search space. Through experiments, both in
simulation and real image datasets, we have presented the robustness and
accuracy of the algorithm while showcasing the efficiency in simulation. The
algorithm can deal with large graphs without any inherent assumptions other
than weak assumptions about connectivity of the graph and presence of adequate
number of simplexes. A distinctive feature of our algorithm is the concept of
feasibility where an edge match is informed by the statistical threshold value
found as a function of the path that the edge is part of. This can lead to
substantial efficiency gains. A second distinctive feature is in the choice of
the p-simplexes with the shortest path in between them as the minimum
topology-preserving unit which again reduces the search space substantially. A
combination of these two features leads to an algorithm that can solve an
inexact subgraph matching problem with large outliers and/or deformation in
sub-linear time (realistically), a large improvement over previous algorithms
which are polynomial complexity, which restricts these algorithms from
matching large scale graphs.

However, our approach suffers from limitations which plague all randomized
algorithms. The initial p-simplex chosen in sub-graph might not be present in
the full graph due to a missing edge or might be deformed beyond the
thresholding level set which can lead to wrongful initial matches. Our future
work would be in understanding how to utilize topology to further remove
inaccurate matches - particularly, moving away from the current technique of
assuming that the initial match is the truth and constantly checking the
previous matches according to global topology. Furthermore, there is a known
trade-off between the algorithm’s robustness and efficiency under different
choices of $p$. Specifically, A large $p$ improves the algorithm robustness by
comparing more edges in each topology unit. However, it sacrifices the
algorithm efficiency by allowing for more feasible $p-$simplex. In this paper,
the method is implemented with $p=2$ (triangles) to improve the efficiency.
For more noisy graph matching problems, it is recommended to increase $p$ for
better matching performances, which will be done in the future work.

## References

  * Balaban, (1985) Balaban, A. T. (1985).  Applications of graph theory in chemistry.  Journal of chemical information and computer sciences, 25(3):334–343.
  * Bhattarai et al., (2019) Bhattarai, B., Liu, H., and Huang, H. H. (2019).  Ceci: Compact embedding cluster index for scalable subgraph matching.  In Proceedings of the 2019 International Conference on Management of Data, pages 1447–1462.
  * Bunke, (2000) Bunke, H. (2000).  Graph matching: Theoretical foundations, algorithms, and applications.  In Proc. Vision Interface, volume 2000, pages 82–88.
  * Byrne, (2013) Byrne, J. (2013).  Topology preserving graph matching.
  * Cantoni et al., (1998) Cantoni, V., Cinque, L., Guerra, C., Levialdi, S., and Lombardi, L. (1998).  2-d object recognition by multiscale tree matching.  Pattern Recognition, 31(10):1443–1454.
  * Cho et al., (2010) Cho, M., Lee, J., and Lee, K. M. (2010).  Reweighted random walks for graph matching.  In European conference on Computer vision, pages 492–505. Springer.
  * Cordella et al., (1998) Cordella, L. P., Foggia, P., Sansone, C., and Vento, M. (1998).  Subgraph transformations for the inexact matching of attributed relational graphs.  In Graph based representations in pattern recognition, pages 43–52. Springer.
  * Cordella et al., (2004) Cordella, L. P., Foggia, P., Sansone, C., and Vento, M. (2004).  A (sub) graph isomorphism algorithm for matching large graphs.  IEEE transactions on pattern analysis and machine intelligence, 26(10):1367–1372.
  * Derpanis, (2010) Derpanis, K. G. (2010).  Overview of the ransac algorithm.  Image Rochester NY, 4(1):2–3.
  * Egozi et al., (2012) Egozi, A., Keller, Y., and Guterman, H. (2012).  A probabilistic approach to spectral graph matching.  IEEE transactions on pattern analysis and machine intelligence, 35(1):18–27.
  * Erdős et al., (1960) Erdős, P., Rényi, A., et al. (1960).  On the evolution of random graphs.  Publ. Math. Inst. Hung. Acad. Sci, 5(1):17–60.
  * Fischler and Bolles, (1981) Fischler, M. A. and Bolles, R. C. (1981).  Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography.  Communications of the ACM, 24(6):381–395.
  * Gold and Rangarajan, (1996) Gold, S. and Rangarajan, A. (1996).  A graduated assignment algorithm for graph matching.  IEEE Transactions on pattern analysis and machine intelligence, 18(4):377–388.
  * Goyal and Ferrara, (2018) Goyal, P. and Ferrara, E. (2018).  Graph embedding techniques, applications, and performance: A survey.  Knowledge-Based Systems, 151:78–94.
  * He and Singh, (2008) He, H. and Singh, A. K. (2008).  Graphs-at-a-time: query language and access methods for graph databases.  In Proceedings of the 2008 ACM SIGMOD international conference on Management of data, pages 405–418.
  * Jüttner and Madarasi, (2018) Jüttner, A. and Madarasi, P. (2018).  Vf2++—an improved subgraph isomorphism algorithm.  Discrete Applied Mathematics, 242:69–81.
  * Liu et al., (2020) Liu, X., Pan, H., He, M., Song, Y., Jiang, X., and Shang, L. (2020).  Neural subgraph isomorphism counting.  In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 1959–1969.
  * Lourens, (1998) Lourens, T. (1998).  A biologically plausible model for corner-based object recognition from color images.  Citeseer.
  * Lowe, (1999) Lowe, D. G. (1999).  Object recognition from local scale-invariant features.  In Proceedings of the seventh IEEE international conference on computer vision, volume 2, pages 1150–1157. Ieee.
  * Lu et al., (1991) Lu, S. W., Ren, Y., and Suen, C. Y. (1991).  Hierarchical attributed graph representation and recognition of handwritten chinese characters.  Pattern Recognition, 24(7):617–632.
  * Mikolajczyk and Schmid, (2005) Mikolajczyk, K. and Schmid, C. (2005).  A performance evaluation of local descriptors.  IEEE transactions on pattern analysis and machine intelligence, 27(10):1615–1630.
  * Rocha and Pavlidis, (1994) Rocha, J. and Pavlidis, T. (1994).  A shape analysis model with applications to a character recognition system.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 16(4):393–404.
  * Shang et al., (2008) Shang, H., Zhang, Y., Lin, X., and Yu, J. X. (2008).  Taming verification hardness: an efficient algorithm for testing subgraph isomorphism.  Proceedings of the VLDB Endowment, 1(1):364–375.
  * Shapiro and Haralick, (1985) Shapiro, L. G. and Haralick, R. M. (1985).  A metric for comparing relational descriptions.  IEEE transactions on pattern analysis and machine intelligence, (1):90–94.
  * Suh et al., (2015) Suh, Y., Adamczewski, K., and Mu Lee, K. (2015).  Subgraph matching using compactness prior for robust feature correspondence.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5070–5078.
  * Tianchang, (2021) Tianchang, W. (2021).  Open graph matching toolbox.
  * Tsai and Fu, (1983) Tsai, W.-H. and Fu, K.-S. (1983).  Subgraph error-correcting isomorphisms for syntactic pattern recognition.  IEEE Transactions on Systems, man, and cybernetics, (1):48–62.
  * Ullmann, (1976) Ullmann, J. R. (1976).  An algorithm for subgraph isomorphism.  Journal of the ACM (JACM), 23(1):31–42.
  * Zaslavskiy et al., (2008) Zaslavskiy, M., Bach, F., and Vert, J.-P. (2008).  A path following algorithm for the graph matching problem.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 31(12):2227–2242.
  * Zhou and De la Torre, (2012) Zhou, F. and De la Torre, F. (2012).  Factorized graph matching.  In 2012 IEEE Conference on Computer Vision and Pattern Recognition, pages 127–134. IEEE.


# Thermal photon anisotropic flow serves as a quark-gluon plasma viscometer

Chun Shen  Ulrich Heinz  Jean-François Paquet  Charles Gale  Department of
Physics, The Ohio State University, Columbus, Ohio 43210-1117, USA  Department
of Physics, McGill University, 3600 University Street, Montreal, Quebec, H3A
2T8, Canada

###### Abstract

Photons are a penetrating probe of the hot and dense medium created in heavy-
ion collisions. We present state-of-the-art calculations of viscous photon
emission from nuclear collisions at RHIC and LHC. Thermal photons’ anisotropic
flow coefficients $v_{n}$ are computed, both with and without accounting for
viscous corrections to the standard thermal emission rates. These corrections
are found to have a larger effect on the $v_{n}$ coefficients than the viscous
suppression of hydrodynamic flow anisotropies. For thermal photons, the ratio
$v_{2}\\{\mathrm{SP}\\}/v_{3}\\{\mathrm{SP}\\}$ shows stronger sensitivity to
the quark-gluon plasma (QGP) shear viscosity than for hadrons, and it can thus
serve as a sensitive QGP viscometer.

###### keywords:

thermal photons, anisotropic flow, shear viscous corrections

##  1 Introduction

Heavy-ion collisions at the Relativistic Heavy-Ion Collider (RHIC) and the
Large Hadron Collider (LHC) allow to study the physics of strongly interacting
quark-gluon plasma (QGP) under conditions of extreme heat and density.
Electromagnetic radiation is a sensitive and direct probe of the space-time
evolution in heavy-ion collisions. Thermal photons provide information about
the early dynamics of heavy-ion collisions that is complementary to that
obtained from the wealth of hadronic observables. Recent measurements show a
surprisingly large second order azimuthal anisotropy of direct photons,
comparable with the elliptic flow of hadrons, both in Au+Au collisions at RHIC
[1, 2] and Pb+Pb collisions at the LHC [3]. These results challenge our
current theoretical understanding of thermal photon production [4]. In this
work, we present state-of-the-art calculations of thermal photon anisotropic
flows from event-by-event viscous hydrodynamic simulations [5]. We find that
thermal photon anisotropic flows are more sensitive to the shear viscosity of
the medium than their charged hadron analogues. The ratio of elliptic and
triangular flow of direct photons is shown to serve as a sensitive QGP
viscometer.

##  2 Methodology

We employ the boost-invariant viscous hydrodynamic code VISH2+1 which has been
very successful in describing soft hadron observables at RHIC and LHC energies
[6, 7]. Initial conditions are generated from the Monte-Carlo Glauber (MCGlb)
and Monte-Carlo KLN (MCKLN) models and then evolved on an event-by-event basis
through viscous hydrodynamics, using a lattice-based equation of state with
chemical freeze-out at $T_{\mathrm{chem}}=165$ MeV [8]. The hydrodynamic
parameters were determined in [6, 7] by fitting hadronic observables at both
RHIC and LHC energies. We start the hydrodynamic simulations at $\tau_{0}=0.6$
fm/$c$ and decouple at $T=120$ MeV. For MCGlb initial conditions, we use the
specific shear viscosity $\eta/s=0.08$ while MCKLN initial conditions require
$\eta/s=0.20$ [6]. The thermal photon momentum distribution can be calculated
by folding the thermal emission rate with the evolving hydrodynamic medium:

$E\frac{dN^{\gamma}}{d^{3}p}=\int\tau d\tau
dxdyd\eta\left(\Gamma_{0}(q,T)+\frac{q_{\mu}q_{\nu}\pi^{\mu\nu}(x)}{2(e(x){+}P(x))}a_{\alpha\beta}\Gamma^{\alpha\beta}(q,T)\right)\bigg{|}_{q=p\cdot
u(x);T(x)}.$ (1)

The photon emission rate includes a dominant equilibrium contribution,
$\Gamma_{0}(q,T)$, and a first order viscous correction proportional to the
shear stress tensor,
$\frac{\pi^{\mu\nu}q_{\mu}q_{\nu}}{2(e+p)}a_{\alpha\beta}\Gamma^{\alpha\beta}(q,T)$,
where $a_{\mu\nu}=\frac{3q_{\mu}q_{\nu}}{2(u\cdot
q)^{4}}+\frac{u_{\mu}u_{\nu}}{(u\cdot q)^{2}}+\frac{g_{\mu\nu}}{2(u\cdot
q)^{2}}-\frac{3(q_{\mu}u_{\nu}+q_{\nu}u_{\mu})}{2(u\cdot q)^{3}}$ [5]. The
rates $\Gamma_{0}(q,T)$ and $\Gamma^{\alpha\beta}(q,T)$ are calculated using
local distribution functions with viscous corrections linear in the shear
stress,
$f(x,p)=f_{0}(x,p)\left[1+(1{\pm}f_{0}(x,p))\frac{p_{\mu}p_{\nu}\,\pi^{\mu\nu}(x)}{2T^{2}(x)(e(x){+}P(x))}\right]$.
The final thermal photon anisotropic flow coefficients are calculated using
the photon multiplicity weighted scalar product method:

$v^{\gamma}_{n}\\{\mathrm{SP}\\}(p_{T})=\frac{\left\langle\frac{dN^{\gamma}}{dyp_{T}dp_{T}}(p_{T})v_{n}^{\gamma}(p_{T})v_{n}^{\mathrm{ch}}\cos(n(\Psi_{n}^{\gamma}(p_{T})-\Psi_{n}^{\mathrm{ch}}))\right\rangle}{\left\langle\frac{dN^{\gamma}}{dyp_{T}dp_{T}}(p_{T})\right\rangle
v_{n}^{\mathrm{ch}}\\{2\\}}.$ (2)

Here $\langle...\rangle$ represents an average over events, and
$v_{n}^{\mathrm{ch}}$ and $\Psi_{n}^{\mathrm{ch}}$ are the $p_{T}$-integrated
charged hadron anisotropic flow coefficients and their associated flow plane
angles in each event. $v_{n}^{\mathrm{ch}}\\{2\\}$ is the charged hadron rms
flow of order $n$, extracted from two-particle cumulants.

##  3 Results

|
---|---
Figure 1: Thermal photon elliptic flow as a function of $p_{T}$ at 0-40%
centrality Pb+Pb collisions at $\sqrt{s}=2.76$ $A$ TeV. Left panel:
Comparisons between thermal photons emitted from ideal and viscous
hydrodynamic medium. Right panel: Comparisons between thermal photons elliptic
flow from smooth event averaged and event-by-event hydrodynamic medium.

In the left panel of Fig. 1, we study the shear viscous correction to the
thermal photon elliptic flow. Compared to results from an ideal hydrodynamic
medium, the viscous hydrodynamic evolution alone results in slightly larger
$v_{2}$ at high $p_{T}$. This is because the initial temperature of the
viscous hydrodynamic medium is lower than in the ideal case, in order to
compensate for viscous entropy production during the evolution. This reduces
the proportion of high momentum photons coming from the early hot region of
the fireball which carries small elliptic flow. However, this increase in
elliptic flow is overwhelmed by a much larger suppression arising from the
viscous correction to the thermal photon emission rate. In contrast to the
elliptic flow of hadrons, the viscous correction to the photon emission rates
dominates the suppression for photons. In the right panel of Fig. 1, we
illustrate the difference between thermal photon $v_{2}$ from smooth event-
averaged and event-by-event hydrodynamic calculations. For both ideal and
viscous cases, the photon $v_{2}\\{\mathrm{SP}\\}$ from event-by-event
simulations are smaller than the mean $v_{2}$ from the smooth event-averaged
hydrodynamic medium. Although the initial state fluctuation increases the
final flow anisotropy, the photon multiplicity factor in Eq. (2) biases the
measurement towards central collisions with smaller elliptic flow, especially
for bumpy profiles with temperature hot spots.

|
---|---
Figure 2: Left panel: Centrality dependence of the $p_{T}$-integrated
$v_{n}\\{\mathrm{SP}\\}$ of thermal photons from MCGlb. and MCKLN initial
conditions. $p_{T}$ is integrated from 1 to 4 GeV. Right panel: The
corresponding ratio $v_{2}\\{\mathrm{SP}\\}/v_{3}\\{\mathrm{SP}\\}$ as a
function of centrality.

The centrality dependence of the $p_{T}$-integrated thermal photon
$v_{n}\\{\mathrm{SP}\\}$ is shown in the left panel of Fig. 2 for $n=2,3$. The
elliptic flow of thermal photons increases with centrality due to the
increasing geometric eccentricity of the nuclear overlap area. The triangular
flow has little centrality dependence because it is purely driven by event-by-
event fluctuations. The ratio between thermal photon elliptic and triangular
flow is shown in the right panel of Fig. 2. Since prompt photons carry a small
anisotropy their contribution is expected to be negligible in this ratio:

$\frac{v^{\gamma}_{2}\\{\mathrm{SP}\\}(p_{T})}{v^{\gamma}_{3}\\{\mathrm{SP}\\}(p_{T})}\simeq\frac{\left\langle\frac{dN^{\gamma,\mathrm{thermal}}}{dyp_{T}\mathrm{dp}_{T}}(p_{T})v_{2}^{\gamma,\mathrm{thermal}}(p_{T})v_{2}^{\mathrm{ch}}\cos(2(\Psi_{2}^{\gamma,\mathrm{thermal}}(p_{T})-\Psi_{2}^{\mathrm{ch}}))\right\rangle}{\left\langle\frac{dN^{\gamma,\mathrm{thermal}}}{dyp_{T}\mathrm{dp}_{T}}(p_{T})v_{3}^{\gamma,\mathrm{thermal}}(p_{T})v_{3}^{\mathrm{ch}}\cos(3(\Psi_{3}^{\gamma,\mathrm{thermal}}(p_{T})-\Psi_{3}^{\mathrm{ch}}))\right\rangle}\frac{v_{3}^{\mathrm{ch}}\\{2\\}}{v_{2}^{\mathrm{ch}}\\{2\\}}.$
(3)

The ratio $v^{\gamma}_{2}\\{\mathrm{SP}\\}/v^{\gamma}_{3}\\{\mathrm{SP}\\}$
increases with the specific shear viscosity $\eta/s$ of the medium, reflecting
the fact that higher order harmonic flows are suppressed more strongly by
viscous effects. Comparing MCGlb results for different $\eta/s$, we find that
for photons this ratio has a larger sensitivity to shear viscosity than the
same ratio for thermal pions [5]. For fixed $\eta/s=0.20$, the centrality
dependence of this ratio is stronger for MCKLN than for MCGlb initial
conditions. This is due to the stronger centrality dependence of
$v^{\gamma}_{2}\\{\mathrm{SP}\\}$ in the MCKLN model, shown in the left panel
of Fig. 2.

Proceeding to Fig. 3, we show in the left panel the $p_{T}$-differential
elliptic and triangular flows of thermal photons. The two initial condition
models (with their correspondingly adjusted choices $\eta/s$) produce very
similar thermal photon elliptic flows, $v^{\gamma}_{2}\\{\mathrm{SP}\\}$,
whereas the triangular flow from MCKLN initial conditions is only about half
as large as for the MCGlb model, due to the larger $\eta/s$ used in evolving
the MCKLN profiles. Compared to hadrons, the thermal photon flow coefficients
$v_{n}\\{\mathrm{SP}\\}(p_{T})$ have a richer structure. They decrease

|
---|---
Figure 3: Left panel: $p_{T}$-differential $v_{2.3}\\{\mathrm{SP}\\}$ of
thermal photons at 0-40% centrality in Pb + Pb collisions at $\sqrt{s}=2.76$
$A$ TeV. Right panel: The corresponding ratio
$v_{2}\\{\mathrm{SP}\\}/v_{3}\\{\mathrm{SP}\\}$ as a function of $p_{T}$
compared with the same ratio for thermal $\pi^{+}$.

at $p_{T}\gtrsim 2.5$ GeV due to the increasing weight of photons coming from
early times when hydrodynamic flow is still weak [4]. The dip at $p_{T}$
around the mass of the $\rho$ meson is due to the transition from
$\pi{\,+\,}\pi{\,\rightarrow\,}\rho{\,+\,}\gamma$ to
$\pi{\,+\,}\rho{\,\rightarrow\,}\pi{\,+\,}\gamma$ as the dominant photon
production channel as $p_{T}$ increases [4]. Since for the hadronic mesons,
who transfer their flow to the photons, $v_{n}(p_{T}){\,\sim\,}p^{n}_{T}$ at
small $p_{T}$, this dip shifts slightly towards larger $p_{T}$ as $n$
increases. This small shift generates an interesting rise-and-fall structure
in the low-$p_{T}$ region for the ratio
$v^{\gamma}_{2}\\{\mathrm{SP}\\}(p_{T})/v^{\gamma}_{3}\\{\mathrm{SP}\\}(p_{T})$,
shown in the right panel of Fig. 3. Compared to charged hadrons, the thermal
photon $v^{\gamma}_{2}\\{\mathrm{SP}\\}/v^{\gamma}_{3}\\{\mathrm{SP}\\}$ ratio
shows higher sensitivity to the shear viscosity of the radiating medium. We
found that this sensitivity is even stronger for the ratios
$v^{\gamma}_{2}\\{\mathrm{SP}\\}/v^{\gamma}_{n}\\{\mathrm{SP}\\}$ with $n\geq
4$ (not shown).

Equation (3) shows that the
$v^{\gamma}_{2}\\{\mathrm{SP}\\}/v^{\gamma}_{3}\\{\mathrm{SP}\\}$ ratio in
Figs. 2 and 3 is insensitive to the yields from direct photon sources that
carry zero anisotropic flow. A measurement of this ratio thus allows to focus
attention on those photons that reflect and probe the space-time history of
the evolving near-thermal medium. Pre-equilibrium photons generated by the
strong magnetic field of the passing nuclei through the QCD anomaly [9] are
predicted to carry sizable $v_{2}$ but are expected to have small triangular
flow $v_{3}$; they would generate much larger $v_{2}/v_{3}$ than thermal
photons. Measuring this ratio will allow to disentangle these two mechanisms
for the already observed (large) direct photon elliptic flow.

##  4 Conclusions

State-of-the-art calculations of thermal photon anisotropic flow,
$v_{n}\\{\mathrm{SP}\\}$ ($n=2,3$), use event-by-event viscous hydrodynamic
simulations to account for event-by-event quantum fluctuations in the initial
state. Shear viscosity suppresses photon $v_{n}\\{\mathrm{SP}\\}$, with
viscous corrections to the photon production rates dominating this
suppression. For both the $p_{T}$-integrated and $p_{T}$-differential
anisotropic flows, the ratio
$v^{\gamma}_{2}\\{\mathrm{SP}\\}/v^{\gamma}_{3}\\{\mathrm{SP}\\}$ shows
stronger sensitivity to the specific shear viscosity of the QGP for thermal
photons than for charged hadrons. This ratio increases with $\eta/s$ because
the viscous suppression of $v_{n}$ increases with the harmonic order $n$.
Since the ratio
$v^{\gamma}_{2}\\{\mathrm{SP}\\}/v^{\gamma}_{3}\\{\mathrm{SP}\\}$ is
insensitive to photon sources that carry zero anisotropic flow, such as prompt
photons, the experimental measurements of this ratio for direct photons will
shed new and more direct light on the specific shear viscosity of the thermal
medium formed after the end of prompt photon emission, but well before most of
finally emitted hadrons are set free.

Acknowledgments:  This work was supported by the U.S. Department of Energy
under Grants No. DE-SC0004286 and (within the framework of the JET
Collaboration) DE-SC0004104 and by the Natural Sciences and Engineering
Research Council of Canada.

## References

  * [1] A. Adare, et al., Observation of direct-photon collective flow in $\sqrt{s_{NN}}=200$ GeV Au+Au collisions, Phys.Rev.Lett. 109 (2012) 122302.  arXiv:1105.4126, doi:10.1103/PhysRevLett.109.122302.
  * [2] I. Tserruya, Photons and low-mass dileptons: results from PHENIX, Nucl.Phys. A904-905 (2013) 225c–232c.  arXiv:1211.6002, doi:10.1016/j.nuclphysa.2013.01.063.
  * [3] D. Lohner, Measurement of Direct-Photon Elliptic Flow in Pb-Pb Collisions at $\sqrt{s_{NN}}=2.76$ TeV, J.Phys.Conf.Ser. 446 (2013) 012028.  arXiv:1212.3995, doi:10.1088/1742-6596/446/1/012028.
  * [4] R. Chatterjee, E. S. Frodermann, U. W. Heinz, D. K. Srivastava, Elliptic flow of thermal photons in relativistic nuclear collisions, Phys.Rev.Lett. 96 (2006) 202302.  arXiv:nucl-th/0511079, doi:10.1103/PhysRevLett.96.202302.
  * [5] C. Shen, U. W. Heinz, J.-F. Paquet, I. Kozlov, C. Gale, Anisotropic flow of thermal photons as a quark-gluon plasma viscometerarXiv:1308.2111.
  * [6] C. Shen, U. Heinz, P. Huovinen, H. Song, Radial and elliptic flow in Pb+Pb collisions at the Large Hadron Collider from viscous hydrodynamic, Phys.Rev. C84 (2011) 044903.  arXiv:1105.3226, doi:10.1103/PhysRevC.84.044903.
  * [7] Z. Qiu, C. Shen, U. Heinz, Hydrodynamic elliptic and triangular flow in Pb-Pb collisions at $\sqrt{s}=2.76$ATeV, Phys.Lett. B707 (2012) 151–155.  arXiv:1110.3033, doi:10.1016/j.physletb.2011.12.041.
  * [8] P. Huovinen, P. Petreczky, QCD Equation of State and Hadron Resonance Gas, Nucl.Phys. A837 (2010) 26–53.  arXiv:0912.2541, doi:10.1016/j.nuclphysa.2010.02.015.
  * [9] G. Basar, D. Kharzeev, D. Kharzeev, V. Skokov, Conformal anomaly as a source of soft photons in heavy ion collisions, Phys.Rev.Lett. 109 (2012) 202303.  arXiv:1206.1334, doi:10.1103/PhysRevLett.109.202303.