# Incommensurate Magnetic Order in the $\mathbb{Z}_{2}$ Kagome Metal GdV6Sn6

Zach Porter  Materials Department, University of California, Santa Barbara,
California 93106, USA  Linac Coherent Light Source and Stanford Institute for
Materials and Energy Sciences, SLAC National Accelerator Laboratory, Menlo
Park, California, 94025, USA  Ganesh Pokharel  Materials Department,
University of California, Santa Barbara, California 93106, USA  Jong-Woo Kim
Phillip J. Ryan  X-ray Science Division, Argonne National Laboratory, Argonne,
Illinois 60439, USA  Stephen D. Wilson<EMAIL_ADDRESS>Materials
Department, University of California, Santa Barbara, California 93106, USA
California Nanosystems Institute, University of California, Santa Barbara,
California 93106, USA

###### Abstract

We characterize the magnetic ground state of the topological kagome metal
GdV6Sn6 via resonant X-ray diffraction. Previous magnetoentropic studies of
GdV6Sn6 suggested the presence of a modulated magnetic order distinct from the
ferromagnetism that is easily polarized by the application of a magnetic
field. Diffraction data near the Gd-$L_{2}$ edge directly resolve a $c$-axis
modulated spin structure order on the Gd sublattice with an incommensurate
wave vector that evolves upon cooling toward a partial lock-in transition.
While equal moment (spiral) and amplitude (sine) modulated spin states can not
be unambiguously discerned from the scattering data, the overall phenomenology
suggests an amplitude modulated state with moments predominantly oriented in
the $ab$-plane. Comparisons to the “double-flat” spiral state observed in Mn-
based $R$Mn6Sn6 kagome compounds of the same structure type are discussed.

Itinerant kagome lattices have drawn considerable recent attention due to the
possibility of realizing correlated topological states and unconventional
forms of electronic order. Recent examples of materials range from $A$V3Sb5
compounds hosting unconventional charge density wave order and
superconductivity [1, 2, 3, 4, 5], to compounds such as Fe3Sn2 [6, 7] and
Mn3Sn [8, 9] that manifest large anomalous/topological Hall effects, to
systems such as FeGe where charge order becomes intertwined with magnetic
order [10].

One versatile class of kagome compounds are materials of the form $R$V6Sn6
where $R$ is a rare earth ion that forms a triangular lattice network in close
proximity to a nonmagnetic kagome network of V-site ions [11]. Choice of the
$R$-site allows for control of the character of the magnetic order
proximitized with the nonmagnetic kagome network, which dominates the bands at
the Fermi level. These compounds are also interesting comparators to their Mn-
based $R$Mn6Sn6 cousins, where magnetic order directly forms within the kagome
network [12, 13, 14].

Of particular interest is the impact of magnetic order on the Dirac points and
saddle points native to the kagome network and very close to the Fermi level
in $R$V6Sn6 compounds [15, 16]. Strong spin-orbit coupling derived primarily
from the heavy Sn-sites intrinsically gaps Dirac points and can create a
tunable Chern gap with a large anomalous Hall effect [17], and while
ferromagnetism can modify the magnitude of this gap, other spin states
proximitized with the topological bands are also of interest.

A number of recent studies have reported differing magnetic anisotropies
manifest in $R$V6Sn6 via the choice of $R$-site [15, 18, 19, 20, 21, 22].
Strongly uniaxial ferromagnetism was reported in TbV6Sn6 [20, 19] and easy-
plane anistotropy was reported in compounds such as (Sm,Er)V6Sn6 [21, 22]. Of
particular note, nearly isotropic magnetic order was observed in spin-only
$S=7/2$ GdV6Sn6 below 5 K [15, 23]. Near the transition however, an unusual
anomaly was observed in magnetoentropy maps, hinting that the zero field state
was not trivially ferromagnetic. This raises the interesting possibility of
realizing an incommensurate, modulated spin state interfaced with the
conducting V-based kagome network, perhaps mimicking the modulation observed
in magnetic kagome lattices of the same structure type (e.g. YMn6Sn6) and key
to stabilizing the topological Hall effect in those materials [24, 12, 25]. To
address this possibility, momentum-resolved studies of the magnetic order in
GdV6Sn6 are needed.

In this paper, we present low-temperature resonant elastic X-ray scattering
(REXS) measurements characterizing the magnetic ground state of GdV6Sn6. Our
data reveal that the zero-field magnetic ground state is an inter-layer
modulated, incommensurate state and with the moments likely predominantly
oriented within the kagome planes. The incommensurate ordering wavevector
evolves upon cooling below the first T${}_{AF1}=5.2(1)$ K ordering transition
and develops an additional, coexisting commensurate component below
T${}_{AF2}=3.8(2)$ K, suggestive of a lower temperature lock-in transition.
Our results demonstrate that an incommensurate magnetic state can be
stabilized in $R$V6Sn6 via choice of an nearly isotropic $R$-site ion, and the
appearance of a low-temperature commensurate harmonic suggests that the nature
of the order is amplitude-modulated and driven by extended RKKY interactions.

For our studies, single crystal samples of GdV6Sn6 were synthesized via a
previously reported flux method [15]. The unit cell is $P6/mmm$ with lattice
parameters $a=b=5.532$ Å and $c=9.162$ Å, and Gd is in Wyckoff site $1b$ at
(0, 0, 0.5); see Figure 1. REXS measurements were performed on beamline 6-ID-B
at the Advanced Photon Source at Argonne National Laboratory. The
diffractometer endstation (Huber Psi-circle) was equipped with a Joule-Thomson
stage displex cryostat capable of reaching a base temperature of 1.9 K.
Measurements were performed near the Gd-$L_{2}$ resonance at $E_{i}=7.933$ keV
using an area detector (Dectris Pilatus 100K) directly after the sample.
Polarization analysis of the scattered beam was performed for select scans
using a pyrolytic graphite (0,0,6) flat single crystal analyzer placed between
the sample and a scintillation point detector. Magnetization was measured
using a Quantum Design Magnetic Property Measurement System (MPMS3) with a
single crystal attached to a quartz paddle with GE varnish.

Figure 1: Incommensurate magnetism in GdV6Sn6. Gd forms hexagonal planes with
one site per unit cell. Two possible Gd orders are shown: spiral and
amplitude-modulated. The propagation direction is along $c$ but the moment
orientation is uncertain.

Figure 2 shows thermodynamic measurements characterizing the onset of magnetic
order in GdV6Sn6. Prior magnetoentropic studies identified the ground state as
possibly noncollinear or modulated due to the presence of a low-field entropy
barrier in the ordered state [15], since inflections in the temperature-
dependent susceptibility evolve with increasing field [23]. Illustrating these
low-field inflections, dc susceptibility $\chi$(T) data are plotted in Figure
2(a) with 10 mT applied within the basal $ab$-plane. Two features appear: the
first is a cusp in the susceptibility at $T_{AF1}=5.2$ K, followed by a
second, more subtle cusp near $T_{AF2}=3.8$ K.

Figure 2: Thermodynamics of GdV6Sn6. Panel (a) shows the dc magnetic
susceptibility $\chi$ measured on (zero) field cooling, denoted (Z)FC. Panel
(b) shows the magnetic contribution to the heat capacity $C_{P,mag}/T$ (red
dots) and entropy $S_{mag}$ (black line). Vertical dashed lines mark two phase
transitions.

Heat capacity data also capture these two temperature scales as plotted in
Figure 2(b). Previously published $C_{P}$($T$) data for GdV6Sn6 [15] were
analyzed to extract the magnetic entropy via removal of the phonon and charge
contributions via a scaled subtraction of the nonmagnetic reference YV6Sn6
[26]. The resulting magnetic $C_{P,mag}$($T$)/$T$ data are plotted along with
the integrated entropy $S_{mag}$($T$). Both $T_{AF1}$ and $T_{AF2}$ can be
identified in $C_{P,mag}$($T$)/$T$ via the peak and inflection point marked as
dashed lines in Figure 2 (b). We note here that the Schottky anomaly expected
due to the mean-field splitting of the $J=7/2$ ground state multiplet is
expected to occur at $T_{AF1}/4\approx 1.3$ K and is likely hidden in the
shoulder of the broad magnetic entropy peak. There is substantial entropy that
extends to temperatures far above the initial ordering temperature, and the
entropy integration up to 30 K nearly reaches the expected $R$ln(8); however
uncertainties in the lattice subtraction and an imperfect lattice standard
begin to dominate at these temperatures. Notably, there is substantial entropy
that continues to be released at temperatures far below $T_{AF1}$, and as we
will show next, scattering data reflect this via a continued staging of wave
vectors as the low temperature limit is approached.

Figure 3: Overview of the resonant X-ray diffraction. Both panels show scans
along (0,0,$L$) at the Gd-$L_{2}$ resonance. Panel (a) shows scans at several
temperatures without an analyzer crystal. Panel (b) shows scans using an
analyzer crystal in the $\sigma-\pi^{\prime}$ and $\sigma-\sigma^{\prime}$
polarization channels.

To explore the origins of these two features in $\chi$($T$) and
$C_{P,mag}$($T$)/$T$ data, resonant X-ray scattering data were collected near
the Gd-$L_{2}$ absorption edge. Magnetic superlattice reflections appear with
a q=(0, 0, $l$) wave vector below the ordering temperatures $T_{AF1}$ and
$T_{AF2}$ in a manner which mirrors the staged transition observed in
$\chi$(T). To further illustrate this staging, $L$-scans were performed about
the Q=(0, 0, 3.5) position in three temperature regimes: $T>T_{AF1}$,
$T_{AF2}>T>T_{AF1}$, and $T_{AF1}>T$. The resulting data are plotted in Figure
3(a).

In Figure 3(a), upon cooling below $T_{AF1}$ an incommensurate set of
reflections appear near Q=(0, 0, 3.47) and Q=(0, 0, 3.53) with qIC ${\approx}$
$(0,0,0.47)$ below 5.2(1) K. Upon further cooling below $T_{AF2}$, an
additional, weaker reflection appears at the commensurate Q=(0, 0, 3.50)
position with qC${=}(0,0,0.50)$ below 3.8(2) K and coexists with the
incommensurate satellites. Both sets of peaks are long-range ordered, with
widths of 0.007 r.l.u. along $L$ and 0.003 r.l.u. along $K$ that mirror those
of nearby structural Bragg peaks. This indicates minimum correlation lengths
of 0.2 microns that are constrained by the crystallinity of the sample.

To further probe the origin of the incommensurate satellite peaks,
polarization analysis was performed to separate the magnetic
$\sigma-\pi^{\prime}$ scattering channel from the nonmagnetic
$\sigma-\sigma^{\prime}$ channel. The results are plotted in Figure 3(b),
where the low temperature incommensurate peaks appear only in the
$\sigma-\pi^{\prime}$ channel, confirming their magnetic origin. The weaker
commensurate peak is not clearly distinguishable in either polarization
channel below $T_{AF2}$. However, the analyzer greatly diminishes signal-to-
noise for these measurements. Perhaps the enhancement near Q=(0, 0, 3.50) in
the $\sigma-\pi^{\prime}$ channel suggests the peak is contributing. As will
be shown next, the commensurate peak’s temperature dependence correlates to
the lower temperature anomaly in $\chi$(T) and implies it also has a magnetic
origin.

Figure 4: Temperature dependence of the resonant X-ray diffraction. Panel (a)
is a false color map of the scattering near (0, 0, 3.5). Three features are
visible: two incommensurate reflections, and one half-integer commensurate
reflection. Panels (b,c) are Gaussian fit results to these data; qIC results
are the weighted mean of both peaks. Panel (b) shows the area (integrated
intensity) of qIC (orange dots) and qC (purple diamonds). Lines in (b) are
guides to the eye of the form $\sqrt{1-T/T_{AF}}$. Panel (c) shows the wave
vector $|$qIC$|$ as a function of temperature.

The detailed temperature dependence of both the qIC and qC peaks are plotted
in Figure 4(a). Peaks were parameterized via line scans at each temperature,
and fit with Gaussian line shapes that quantify the peak area and position
upon cooling toward base temperature. The resulting thermal evolution of both
commensurate and incommensurate order parameters was well as the
incommensurate wave vector qIC are shown in Figures 4(b,c).

Upon entering the incommensurate phase below $T_{AF1}$, qIC begins to decrease
with continued cooling and eventually reaches a minimum near $T_{AF2}$. Upon
further cooling below $T_{AF2}$, weak scattering weight appears at the
commensurate qC position and the wavevector qIC of the incommensurate peaks
begins to increase. This trend continues down to the lowest temperature probed
($\approx 2$ K). This interplay suggests that the discommensuration of the qIC
order is tied to the emergence of the qC order parameter, consistent with a
partial lock-in transition.

In the susceptibility $\chi$(T), the slight easy-plane anisotropy combined
with the fact that both $T_{AF1}$ and $T_{AF2}$ appear only under an in-plane
field suggest that the moments are predominantly in-plane. Magnetic scattering
surveyed in multiple zones reveals that the scattering is strongest in
$c$-axis aligned Q=(0, 0, $L$) type positions, also consistent with an in-
plane component 111See the Supplemental Material at [] for additional X-ray
scattering and susceptibility measurements.. However, in order to
unambiguously determine the moment orientation, future measurements such as
azimuthal scans at multiple wave vectors will be required.

The incommensurate state observed in our REXS measurements is commonly
indicative of either a flat spiral state or an amplitude-modulated state as
illustrated by cartoons in Fig. 1. A flat spiral is reminiscent of the “double
flat” spiral states reported in Mn-based $R$Mn6Sn6 compounds [28, 29]. A key
distinction, however, is that the magnetic order in compounds such as YMn6Sn6
and LuMn6Sn6 derives from their Mn-based bilayer kagome networks, which
naïvely host competing interlayer exchange interactions. Meanwhile magnetic
order in $R$V6Sn6 compounds is derived from their triangular lattice
lanthanide networks, where extended (beyond nearest neighbor) interlayer
exchange interactions are much less likely. This suggests that RKKY
interactions play a dominant role in the stabilization of the modulated
magnetic state.

A further distinction is the development of a commensurate harmonic qC at low
temperature in GdV6Sn6 versus the high temperature commensurate harmonic in
(Y, Lu)Mn6Sn6 that splits into incommensurate wavelengths upon cooling [29,
24]. The onset of the qC harmonic upon cooling below $T_{AF2}$ can be
envisioned as the local formation of phase-locked modulations of
antiferromagnetic planes that propagate within the longer wavelength,
incommensurate state. The fact that this commensurate harmonic strengthens
upon cooling toward the ground state points to the formation of a partial
lock-in transition [30]—one which may develop further at lower temperatures.
This transition toward an equal moment state at low temperature is suggestive
of an amplitude modulated character to the incommensurate order. The
substantial entropy remaining far below the onset of the incommensurate
$T_{AF1}$ [31, 32] and the associated minimization of entropy achieved as the
qC harmonic develops are also consistent with an amplitude modulated nature to
the transition.

While the incommensurate state that forms below $T_{AF1}$ mimics those of
(Y,Lu)Mn6Sn6, its origin should be distinct from the physics proposed in those
compounds. The evolution of qIC upon cooling below $T_{AF1}$ likely arises
from entropic factors associated with an amplitude modulated state rather than
lattice effects associated with unit cell changes tuning the balance of
interlayer exchange interactions in a flat spiral picture. The incommensurate
modulation wavevector qIC is roughly consistent with nesting between the small
Fermi surface pockets identified along the M-L line of the Brillouin zone
[23], again evoking an RKKY mechanism for the modulated state. Though future
measurements examining whether the ordered spin state is collinear or
noncollinear will be required, our current data demonstrate a complex
condensation of magnetic order in the topological kagome metal GdV6Sn6. This
provides an interesting avenue for engineering new magnetic states in
proximity to the kagome planes in $R$V6Sn6 compounds and for exploring their
impact on the topological band structures.

###### Acknowledgements.

This work was supported by the National Science Foundation (NSF) through
Enabling Quantum Leap: Convergent Accelerated Discovery Foundries for Quantum
Materials Science, Engineering and Information (Q-AMASE-i): Quantum Foundry at
UC Santa Barbara (Grant No. DMR-1906325). S.D.W. and Z.P. acknowledge support
from NSF Grant No. DMR-1905801. This research made use of the shared
facilities of the NSF Materials Research Science and Engineering Center at UC
Santa Barbara, Grant No. DMR-1720256. Z.P. acknowledges additional support
from the U.S. Department of Energy, Office of Science, Basic Energy Sciences,
Materials Sciences and Engineering Division, under Contract No. DE-
AC02-76SF00515. This research used resources of the Advanced Photon Source, a
U.S. Department of Energy (DOE) Office of Science User Facility operated for
the DOE Office of Science by Argonne National Laboratory under Contract No.
DEAC02-06CH11357.

## References

  * Ortiz  _et al._ [2021] B. R. Ortiz, P. M. Sarte, E. M. Kenney, M. J. Graf, S. M. L. Teicher, R. Seshadri, and S. D. Wilson, Phys. Rev. Mater. 5, 034801 (2021).
  * Ortiz  _et al._ [2020] B. R. Ortiz, S. M. L. Teicher, Y. Hu, J. L. Zuo, P. M. Sarte, E. C. Schueller, A. M. M. Abeykoon, M. J. Krogstad, S. Rosenkranz, R. Osborn, R. Seshadri, L. Balents, J. He, and S. D. Wilson, Phys. Rev. Lett. 125, 247002 (2020).
  * Yin  _et al._ [2021] Q. Yin, Z. Tu, C. Gong, Y. Fu, S. Yan, and H. Lei, Chinese Physics Letters 38, 037403 (2021).
  * Jiang  _et al._ [2021] Y.-X. Jiang, J.-X. Yin, M. M. Denner, N. Shumiya, B. R. Ortiz, G. Xu, Z. Guguchia, J. He, M. S. Hossain, X. Liu, J. Ruff, L. Kautzsch, S. S. Zhang, G. Chang, I. Belopolski, Q. Zhang, T. A. Cochran, D. Multer, M. Litskevich, Z.-J. Cheng, X. P. Yang, Z. Wang, R. Thomale, T. Neupert, S. D. Wilson, and M. Z. Hasan, Nature Materials 20, 1353 (2021).
  * Zhao  _et al._ [2021] H. Zhao, H. Li, B. R. Ortiz, S. M. L. Teicher, T. Park, M. Ye, Z. Wang, L. Balents, S. D. Wilson, and I. Zeljkovic, Nature 599, 216 (2021).
  * Li  _et al._ [2019] H. Li, B. Ding, J. Chen, Z. Li, Z. Hou, E. Liu, H. Zhang, X. Xi, G. Wu, and W. Wang, Applied Physics Letters 114, 192408 (2019), https://doi.org/10.1063/1.5088173 .
  * Ye  _et al._ [2018] L. Ye, M. Kang, J. Liu, F. von Cube, C. R. Wicker, T. Suzuki, C. Jozwiak, A. Bostwick, E. Rotenberg, D. C. Bell, L. Fu, R. Comin, and J. G. Checkelsky, Nature 555, 638 (2018).
  * Taylor  _et al._ [2020] J. M. Taylor, A. Markou, E. Lesne, P. K. Sivakumar, C. Luo, F. Radu, P. Werner, C. Felser, and S. S. P. Parkin, Phys. Rev. B 101, 094404 (2020).
  * Nakatsuji  _et al._ [2015] S. Nakatsuji, N. Kiyohara, and T. Higo, Nature 527, 212 (2015).
  * Teng  _et al._ [2022] X. Teng, L. Chen, F. Ye, E. Rosenberg, Z. Liu, J.-X. Yin, Y.-X. Jiang, J. S. Oh, M. Z. Hasan, K. J. Neubauer, B. Gao, Y. Xie, M. Hashimoto, D. Lu, C. Jozwiak, A. Bostwick, E. Rotenberg, R. J. Birgeneau, J.-H. Chu, M. Yi, and P. Dai, Nature 609, 490 (2022).
  * Romaka  _et al._ [2011] L. Romaka, Y. Stadnyk, V. Romaka, P. Demchenko, M. Stadnyshyn, and M. Konyk, Journal of Alloys and Compounds 509, 8862 (2011).
  * Ghimire  _et al._ [2020] N. J. Ghimire, R. L. Dally, L. Poudel, D. C. Jones, D. Michel, N. T. Magar, M. Bleuel, M. A. McGuire, J. S. Jiang, J. F. Mitchell, J. W. Lynn, and I. I. Mazin, Science Advances 6, eabe2680 (2020), https://www.science.org/doi/pdf/10.1126/sciadv.abe2680 .
  * Gorbunov  _et al._ [2012] D. Gorbunov, M. Kuz’min, K. Uhlířová, M. Žáček, M. Richter, Y. Skourski, and A. Andreev, Journal of Alloys and Compounds 519, 47 (2012).
  * Hu  _et al._ [1995] J. Hu, K.-Y. Wang, B.-P. Hu, Y.-Z. Wang, Z. Wang, F. Yang, N. Tang, R. Zhao, and W. Qin, Journal of Physics: Condensed Matter 7, 889 (1995).
  * Pokharel  _et al._ [2021] G. Pokharel, S. M. L. Teicher, B. R. Ortiz, P. M. Sarte, G. Wu, S. Peng, J. He, R. Seshadri, and S. D. Wilson, Phys. Rev. B 104, 235139 (2021).
  * Peng  _et al._ [2021] S. Peng, Y. Han, G. Pokharel, J. Shen, Z. Li, M. Hashimoto, D. Lu, B. R. Ortiz, Y. Luo, H. Li, M. Guo, B. Wang, S. Cui, Z. Sun, Z. Qiao, S. D. Wilson, and J. He, Phys. Rev. Lett. 127, 266401 (2021).
  * Yin  _et al._ [2020] J.-X. Yin, W. Ma, T. A. Cochran, X. Xu, S. S. Zhang, H.-J. Tien, N. Shumiya, G. Cheng, K. Jiang, B. Lian, Z. Song, G. Chang, I. Belopolski, D. Multer, M. Litskevich, Z.-J. Cheng, X. P. Yang, B. Swidler, H. Zhou, H. Lin, T. Neupert, Z. Wang, N. Yao, T.-R. Chang, S. Jia, and M. Zahid Hasan, Nature 583, 533 (2020).
  * Lee and Mun [2022] J. Lee and E. Mun, Phys. Rev. Mater. 6, 083401 (2022).
  * Rosenberg  _et al._ [2022] E. Rosenberg, J. M. DeStefano, Y. Guo, J. S. Oh, M. Hashimoto, D. Lu, R. J. Birgeneau, Y. Lee, L. Ke, M. Yi, and J.-H. Chu, Phys. Rev. B 106, 115139 (2022).
  * Pokharel  _et al._ [2022] G. Pokharel, B. Ortiz, J. Chamorro, P. Sarte, L. Kautzsch, G. Wu, J. Ruff, and S. D. Wilson, Phys. Rev. Mater. 6, 104202 (2022).
  * Huang  _et al._ [2023] X. Huang, Z. Cui, C. Huang, M. Huo, H. Liu, J. Li, F. Liang, L. Chen, H. Sun, B. Shen, Y. Zhang, and M. Wang, Anisotropic magnetism and electronic properties of the kagome metal smv6sn6 (2023), arXiv:2303.00627 [cond-mat.str-el] .
  * Zhang  _et al._ [2022] X. Zhang, Z. Liu, Q. Cui, Q. Guo, N. Wang, L. Shi, H. Zhang, W. Wang, X. Dong, J. Sun, Z. Dun, and J. Cheng, Phys. Rev. Mater. 6, 105001 (2022).
  * Ishikawa  _et al._ [2021] H. Ishikawa, T. Yajima, M. Kawamura, H. Mitamura, and K. Kindo, Journal of the Physical Society of Japan 90, 124704 (2021), https://doi.org/10.7566/JPSJ.90.124704 .
  * Dally  _et al._ [2021] R. L. Dally, J. W. Lynn, N. J. Ghimire, D. Michel, P. Siegfried, and I. I. Mazin, Phys. Rev. B 103, 094413 (2021).
  * Wang  _et al._ [2021] Q. Wang, K. J. Neubauer, C. Duan, Q. Yin, S. Fujitsu, H. Hosono, F. Ye, R. Zhang, S. Chi, K. Krycka, H. Lei, and P. Dai, Phys. Rev. B 103, 014416 (2021).
  * Bouvier  _et al._ [1991] M. Bouvier, P. Lethuillier, and D. Schmitt, Phys. Rev. B 43, 13137 (1991).
  * Note [1] See the Supplemental Material at [] for additional X-ray scattering and susceptibility measurements.
  * Rosenfeld and Mushnikov [2008] E. Rosenfeld and N. Mushnikov, Physica B: Condensed Matter 403, 1898 (2008).
  * Venturini  _et al._ [1996] G. Venturini, D. Fruchart, and B. Malaman, Journal of Alloys and Compounds 236, 102 (1996).
  * Penc  _et al._ [1999] B. Penc, M. Hofmann, J. Leciejewicz, and A. Szytula, Journal of Physics: Condensed Matter 11, 7579 (1999).
  * Blanco  _et al._ [1991] J. A. Blanco, D. Gignoux, and D. Schmitt, Phys. Rev. B 43, 13145 (1991).
  * Rotter  _et al._ [2001] M. Rotter, M. Loewenhaupt, M. Doerr, A. Lindbaum, and H. Michor, Phys. Rev. B 64, 014402 (2001).


# ATLAS: An IoT Architecture and Secure Open-source Networking Stack for
Anonymous Localization and Tracking Using
Smartphones and Bluetooth Beacons

Bharath Srinivas Prabakaran1,3, Felix Fasching1,3, Juri Schreib1,3, Andreas
Steininger1, Muhammad Shafique2 3 These authors have contributed to this work
equally. 1Institute of Computer Engineering, Technische Universität Wien (TU
Wien), Austria
{bharath.prabakaran<EMAIL_ADDRESS>{e11712208,
<EMAIL_ADDRESS>2eBrain Lab, Division of Engineering, New
York University Abu Dhabi, UAE
<EMAIL_ADDRESS>
###### Abstract

Bluetooth (BT) has revolutionized close-range communication enabling smart
capabilities in everyday devices through wireless technology. One of the most
important sub-domains of Internet-of-Things (IoT) specializes in the usage of
BT technologies to develop smart homes and environments, which include
hospitals, buildings, shopping facilities, etc. to offer a wide-range of
features, like instantaneous and remote access to ventilation, lighting,
security, localization, and tracking. However, the deployment of such features
in smart infrastructures are typically unaccompanied by appropriate security
measures that safeguard the data and protect its users. Towards this, we
propose the ATLAS framework, which is composed of our novel IoT architecture
and secure networking stack that can be used to anonymously localize and track
smartphones and wearables by deploying multiple Bluetooth Low Energy (BLE)
beacons across the environment. The proposed networking stack enables varying
levels of encryption across all layers of the communication stack to ensure an
easy-to-adopt, secure-by-design network architecture. We also deploy a novel
data transformation and fingerprinting-based localization algorithm, which is
highly effective in localizing user devices within a given area. The ATLAS
framework is open-sourced at https://atlas-tuw.sourceforge.io to enable wide-
spread adoption and further research and development.

###### Index Terms:

IoT, Smart, Environment, Security, Privacy, Anonymous, Tracing, BLE,
Localization, Tracking, Bluetooth

##  I Introduction

The next generation of automation-based quality of life improvements are
expected primarily due to the advancements in the domain of Internet-of-Things
(IoT) [1], by enables communication amongst physical devices. This
communication can be established via various wireless platforms, like Zigbee
[2], Wi-Fi [3], and Bluetooth (BT) [4], through which the user can also
interact, on-demand, with devices in their the surrounding environment.
Bluetooth has proven to be quite effective as a close-range communication
system, especially with the advent of Bluetooth Low Energy (BLE) [5], which
has been developed specifically for IoT platforms and use-cases, such as smart
homes, building, and environments.

Figure 1: An illustration of users interacting with the devices in smart
environments. Connecting everyday things to the Internet can enable on-demand
access to these devices via a single interface. Localization can further
enhance the experience to automatically activate such devices and features
when the user is in close proximity.

Smart environments tend to offer a wide-range of functionalities and services,
which can be controlled or operated using a single user device. These services
include remote control of heating, ventilation, and air conditioning (HVAC)
systems, lighting, security systems, robots, transportation systems, smart
grids, charging stations, indoor positioning, localization, tracking, etc. An
overview of such a smart environment is illustrated in fig. 1. For instance,
the smart television and sound system can be remotely enabled and controlled
by the user’s smartphone. Similarly, connecting video cameras and other
security systems to the internet enables the user to monitor their home and/or
their pet during their vacation and/or work, using just a smartphone.
Furthermore, indoor localization techniques can be highly effective in
automating features like heating, lighting, and ventilation, based on the
user’s proximity and preferred temperature settings.

Facilitating such localization techniques for larger environments, such as
public buildings and shopping malls, can be quite beneficial in reducing their
energy consumption and aid in emergency situations by tracking and tracing the
users, e.g., contact tracing. However, realizing such systems requires the
research and development of an IoT architecture with a prime focus on the
security, privacy, and anonymity of its users. Enabling the development of a
secure network infrastructure that can anonymously localize and track the user
devices requires us to address the following research challenges:

  1. (1)

How can Bluetooth low energy beacons be used to localize and track users in
smart environments?

  2. (2)

What kind of IoT architecture would ensure maximum coverage area with minimum
resource overhead?

  3. (3)

How do we design, develop, and deploy an easy-to-adopt secure networking stack
to ensure user privacy?

To address these challenges, we propose the ATLAS framework with the following
novel contributions:

  * •

We investigate the capability of BT beacons that can be deployed as part of
our novel IoT architecture. These beacons interact with BT-enabled devices
that can be used to locate and track its users anonymously.

  * •

A secure-by-design architecture that prevents user-to-user communication,
thereby enabling a higher level of security and privacy-preserving contact
tracing capabilities.

  * •

A fully open-source easy-to-adopt networking stack with varying levels of
encryption at each communication interface to ensure the privacy and security
of its users.

  * •

A novel data transformation technique and localization algorithm that can be
used to anonymously determine the location of the user and track them through
the environment.

  * •

Our framework is capable of localizing all user devices in range, at room-
level precision, while incorporating security measures that introduce minimal
communication overheads.

##  II Related Work

Plenty of works have focused on device localization using various
communication technologies and indoor positioning techniques, based on metrics
like angle-of-arrival, time-of-flight, and signal strength. Bahl et al. [6]
proposed a radio-frequency based system that relies on received signal
strength indicators (RSSIs) to localize and track users inside buildings with
the help of multiple base stations. Lim et al. [7] proposed a localization
algorithm that relies on the RSSIs  (1) between the WiFi access points and
(2) the WiFi access points and the user  to estimate the user location. Diaz
et al. [8] proposed a signal coverage density method that creates a virtual
matrix of the environment, which are composed of smaller blocks called cells.
The proposed approach uses multilateration to determine the most-likely cell
location of the user. Shirehjini et al. [9] proposed a low-range passive Radio
Frequency Identification-based (RFID) indoor positioning system that consists
of carpets and peripherals embedded with RFID tags for effective data
interpretation and device localization. Kotaru et al. [10] presented an
accurate indoor localization system that incorporates super-resolution
algorithms, which accurately calculate the angle-of-arrival to determine the
location of the WiFi users to decimeter-level accuracy. Ashraf et al. [11]
proposed the use of several convolutional neural networks to analyze
geomagnetic field patterns and determine the exact user location.

Figure 2: An overview of the proposed ATLAS framework. Note, the novel IoT
architecture depicted on top illustrates the lack of communication among
different user nodes, which provides a higher layer of security against
malicious devices, which interface solely with the beacons. After the initial
setup stage, the user device interacts with the BLE beacons, periodically,
which subsequently transmit data to the hub that is processed to generate
location information that can be used to track them. Each layer of the
networking stack is encrypted and secured to ensure the privacy, security, and
anonymity of its users.

However, more recently, BT-based localization has proven to be quite
effective, especially with the advent of BLE. For instance, Zafari et al. [12]
proposed a BT-based indoor localization system that relies on novel filtering
and averaging approaches, which are used to reduce the fluctuations in RSSIs
to improve the localization accuracy. Martin et al. [13] also presented a BLE-
based localization technique around the iBeacon technology protocol invented
by Apple. This technology is currently in use in the Apple AirTag [14]
devices, which can be leveraged by its users to locate and track things with
the help of BT and the network of apple smart phones around the world. Google
has had similar technology ventures in indoor localization via the Eddystone
project [15]. However, these devices act as passive beacons without any trust
verification for the beacons, require internet access, and user device
identification. To ensure the privacy of the users and the environmental BLE
beacons, Fawaz et al. [16] proposed a device-agnostic system that can enable
the users and system administrators to control who can discover, scan, and
connect to their devices. On the other hand, Kolias et al. [17] analyzed a
wide plethora of BLE nodes to identify vulnerabilities can track the users,
their behavior, and spoof the network.

These types of techniques primarily focus on a single aspect, which is either
localization and indoor positioning or security and privacy. To the best of
our knowledge, our ATLAS framework is the first to combine both aspects to
ensure that the privacy and security of its users are not compromised while
ensuring that the user devices are localized and tracked periodically, using
our localization algorithm. Note, the use of other localization algorithms is
orthogonal to the use of the proposed localization technique. Similarly, the
use of other secure networking stack protocols can be easily incorporated into
the ATLAS framework based on the system requirements and specifications.

##  III Our ATLAS Framework

Fig. 2 illustrates an overview of our ATLAS framework including the proposed
novel IoT architecture and the secure, open-source, and easy-to-adopt
networking stack, which can be used to anonymously track users using our novel
localization algorithm. During the initial setup stage, network administrators
and system designers are required to deploy the BLE beacons, interface them
with the hub, and generate a set of reference points known as the fingerprint
map, which can be used to subsequently localize the users to specific “areas”
of the environment. The entire networking stack is fully encrypted using
various encryption techniques to ensure the privacy, security, and anonymity
of its users.

###  III-A IoT Architecture

According to a $2018$ mobile consumer survey, smartphone penetration in
developed nations and regions, such as the EU and the US, averages at $91\%$,
i.e., on average, $91\%$ of the population owns a smartphone [18]. Due to
their wide-spread adoption and the low energy requirements of BT, we propose
to use these devices to actively localize and track the users by communicating
with them periodically, especially since the BT module on such devices are
typically active communicating with other wearables such as smart-watches and
wireless headphones. We do this by establishing an initial communication via a
quick response (QR) code (generated by the system integrator; block 0 in Fig.
2), which the users scan through the custom application to interface with the
smart environment (block 1). This one-time user interaction establishes device
communication with nearby BLE beacons [19] that are deployed as part of the
environment, as shown in Fig. 2. The beacons advertise and communicate with
the user devices, periodically, to estimate the received signal strength
indicator (RSSI) for each user in the network (block 2). The beacons
subsequently transmit this data to the hub (block 3), which connects to all
beacons in the “environment”, to process and determine the live fingerprints
of the users, along with a timestamp regarding their communication, to enable
tracking the user through the environment (see Section III-C). The system
integrator (block 0), which is active only during the setup stage, is
responsible for determining the initial fingerprint map, which serves as a
reference for localizing user devices. The data collected across multiple
environments, by their respective hub nodes, are transmitted to a persistent
cloud-based storage system (block 4), which is stored for a period of $28$
days as per the privacy regulations of the EU [20] in order to enable
anonymous tracking and contact tracing capabilities.

###  III-B Secure Networking Stack

To ensure the privacy, security, and anonymity of the users, we propose to
secure each communication layer in the network stack with varying levels of
encryption. The proposed secure networking stack is open-source and available
online at https://atlas-tuw.sourceforge.io, which can be easily adopted for
similar IoT architectures to ensure wide-spread adoption and further research
in this domain. The initial pairing of the user device with the BLE beacons
takes place using out-of-band (OOB) pairing with the OOB secret hidden in the
QR code, which needs to be scanned using the custom application. The
transmitted data packets are BLE link-layer encrypted, to ensure that the
packet sniffers are unable to extract user information. Next, the data
transmission between the hub and the beacons is completed using the
constrained application protocol (CoAP) [21] combined with the datagram
transport layer security (DTLS) [22] protocol to provide a higher layer of
security in order to prevent eavesdropping, tampering, and message forgery.
During beacon deployment and system setup, symmetric keys are pre-shared among
the beacons and the hub, to establish a secure line of DTLS-PSK communication,
followed by which the public key infrastructure (DTLS-PKI) is available to
provide reliable beacon authentication to the hub [23]. The hub interfaces
with the cloud-based storage device using the Bolt protocol [24] with a TLS
encryption to ensure the security requirements of our system. Similarly, the
database stored in the cloud is also encrypted. Investigating the security
aspects and vulnerabilities of the implemented communication and encryption
protocols lie beyond the scope of this work.

###  III-C Localization and Tracking

Figure 3: An overview of the localization technique used in the ATLAS
framework. The proposed technique works in two stages; the setup stage, during
which an initial fingerprint map is calculated to aid in user device
localization, and the live stage, which continuously executes at run-time to
actively locate and track users.

We propose to use the fingerprinting technique, shown in Fig. 3, for
localizing user devices, which works in two stages:

####  III-C1 Setup stage

Our technique relies on the existence of a predefined set of reference points,
which requires our system to undergo an initial setup stage composed of active
and passive phases. During the active phase, a smartphone with our customized
application and administrator privileges is used to generate the set of
reference points, using the RSSI values from each beacon and the device’s
coordinate values, at different locations in the environment (block 0 in Fig.
2). These reference points should be, typically, distributed evenly throughout
the environment in order to get optimal results; a non-evenly distributed set
of references may lead to incorrect localizations. In the passive phase, the
data collected during the active phase is sorted, based on their timestamps,
and filtered based on their RSSI values per beacon, using a Kalman filter to
smoothen potential outliers. The processed data act as individual fingerprints
that contain the RSSI per beacon at the reference coordinate. These
fingerprints are used to generate the fingerprint map, which is stored as a
database on the cloud. Fig. 4 illustrates an overview of this approach. Note,
any significant changes to the environment, i.e., changes to the existing
infrastructure, which could negatively affect the localization, can be
addressed by re-executing the initial setup stage to generate a new set of
reference points.

Figure 4: An overview of the data transformation technique used to generate
fingerprints and the fingerprint map. After collating the data in the required
format, the data is filtered and processed using weighted averaging to
generate the fingerprint map. Algorithm 1 The Localization Algorithm

1:$\text{User Fingerprint Array}(FP_{U})$

2:$\text{Fingerprint Map}(FM)$

3:$\text{Array of User Locations}(L_{U})$

4:$L_{U}=Dict(\\{\\});$

5:$boolArray=[];$

6:for $U$ in $FP_{U}$ do

7: $validFP=Dict(\\{\\});$

8: for $fp$ in $FM$ do

9: $fpB_{A}=fp.beaconArray()$

10: $UB_{A}=U.beaconArray()$

11: $boolArray=fpB_{A}$ $\&$ $UB_{A}$

12: if $size(boolArray)<2$ then

13: continue

14: end if

15: $validFP[fp]=fp.beaconArray()$

16: end for

17: $wtfp=Dict(\\{\\});$

18: for $vfp$ in $validFP$ do

19: $fpw=0$

20: $wt=0$

21: for $beacon$ in $U.fingerprint$ do

22: if $beacon.uuid$ in $vfp.beaconArray()$ then

23: $val1=beacon.rssi$

24: $val2=vfp.beaconArray(beacon.uuid).rssi$

25: $wt=wt+(val1-val2)^{2}$

26: end if

27: end for

28: $wt=wt/size(boolArray)$

29: $wt=sqrt(wt)$

30: $fpw=fpw+wt$

31: $wtfp[vfp]=fpw$

32: end for

33: $L_{U}[U]=max(wtfp)$

34:end for

35:return $L_{U}$

Figure 5: Timing analysis of our ATLAS framework. Collecting the reference
points depends on the environment and the number of beacons deployed; for the
purpose of simplicity, we consider $1$ clock-cycle duration. The beacons
connect to the user device every $15$s and since a beacon-user interaction
requires up to $3$s, theoretically, a total of up to $40$ user devices can
connect to a single BLE beacon.

####  III-C2 Live stage

After the setup stage, the live stage of our localization technique is
executed continuously at run-time to enable interactions with new users and
track them through the smart environment. Like in the previous stage, data is
collected during the active phase of this approach, and processed during the
passive phase to perform localization. During the active phase, the beacons
periodically advertise to the connected user devices and collect the RSSI
information for each device in the environment. The data processing, on the
other hand, is completed during the passive phase of the live stage, wherein
the data collected during the active phase is used to construct fingerprints,
similar to the approach presented in the passive phase of the setup stage.
Next, we fetch the fingerprint map from the database to evaluate the location
of the user by comparing the fingerprints obtained from the user to the
fingerprints in the map. We do this by iteratively comparing each generated
user fingerprint with fingerprints in the fingerprint map. First, we determine
the set of common beacons between the two compared fingerprints, which if less
than 2 are automatically discarded, and store them in an array of valid
fingerprints ($validFP$) (lines 8-16). Next, we determine the proximity of a
user to a given fingerprint by estimating the euclidean distance between the
RSSI value for each beacon in the user fingerprint array and the RSSI value of
the beacon from the valid fingerprint array. These distances are accumulated
per beacon to generate a weight distribution for each fingerprint in the
fingerprint map, thereby determining a proximal location for each user, based
on their fingerprint (lines 18-31). After iterating through each generated
fingerprint per user, an array of the proximal locations of each user, along
with their hashed user ID, is returned by the algorithm (line 35), which is
subsequently stored in the cloud, as discussed in Section III-A. Algorithm 1
discusses the pseudo-code for the proposed localization technique.

##  IV Experimental Setup

We have used the following list of equipment for the illustrating the results
presented in Section V: $9$ ESP-$32$ beacons, which are WiFi and BLE enabled,
$1$ Hub (Dell XPS $13$ with an Intel i$5-6200$U processor with $8$GB RAM and
$256$GB SSD, running Ubuntu $21.04$), and a WiFi-Access Point (Huawei $4$G
Router B$528$). The beacons are deployed in various configurations across two
floors of a building, an overview of which is provided in Section V. We use a
$2019$ iPad Pro as our system integrator (block 0 in Fig. 2) during the active
phase of the setup stage, which subsequently transmits the collected data to
the hub to generate the fingerprint map during the passive stage. The final
deployed environment is tested during the live stage using $4$ smartphones and
tablets, an iPhone $8$ Plus, a $2016$ iPhone SE, a $2019$ iPad Pro, and a
$2019$ iPad (all of them running iOS). The mobile application is programmed in
Swift utilizing the ARKit, a software development kit (SDK) that associates
beacon RSSI values with user device coordinates, and the Bluetooth core, which
is essential for communicating with our BLE beacons. Similarly, for android
platforms, Google offers an SDK called ARCore, which can be used to design
similar applications. The firmware for the ESP-$32$ BLE beacons are compiled
in C with the NimBLE and LibCoAP libraries for Bluetooth and WiFi
communication. Python is used as the default programming language on the hub
to realize the Django web service, the CoAPS, HTTPS, and Bolt communication
protocols. For simplicity, the persistent cloud-based storage system is
realized as a Virtual Private Server (VPS) in the hub, hosted by Netcup, to
store the encrypted fingerprint map database. The VPS is allocated $1$ core,
$2$GB RAM, and $20$GB of SSD space to run a Debian bullseye $11$ operating
system, which executes the Java-based Neo4J graph database that is ideal for
storing graph-based information, such as the timestamped fingerprint maps
generated by our framework. The Django-neomodel Python library is used as an
object graph mapper for communicating with the Neo4J database.

##  V Results and Discussion

Figure 8: Localization and tracking analysis of four users in a multi-
environment scenario. Our ATLAS framework is capable of recording and
processing the data of multiple users across multiple environments, which are
interconnected.

Fig. 5 presents an overview of the end-to-end timing analysis of our ATLAS
framework. After the initial beacon configuration and deployment, the system
integrator is also responsible for collecting the reference points, the
duration of which heavily relies on the environment and the number of beacons
deployed. For the sake of simplicity, we consider this duration to be $1$s in
our analysis. Once the beacons are setup and deployed, the hub configures them
to immediately start advertising, but the users can only connect to them once
the setup stage is complete as shown by the clock at $t_{11}$. During this
time, the beacon information is stored in the database and the fingerprint map
is calculated, which is also transferred to the database in the next
subsequent clock-cycle. Since the round-trip communication between the user
and the beacon requires $3$s, primarily due to the encryption layer, the BLE
advertisements and the active/passive phases of the live stage in the hub are
synchronized accordingly, as shown. Theoretically, every beacon can connect to
up to $8$ user devices in parallel using BLE and since we repeatedly ping the
discovered devices every $15$s, we can locate up to $40$ devices in an area.
Similarly, we can also increase the localization frequency of the discovered
devices by trading-off on the number of locatable devices in the environment.
The data collected by the beacon is transmitted to the hub for processing and
comparison with the fingerprint map during the next passive phase of the hub.
This process executes continuously until the system is manually shutdown or
paused by the administrator.

Figure 6: The floor-plan and the beacons deployed in Environment 1. With the
user generated data for fingerprints and their timestamps, we can localize and
track their devices across our environment. Figure 7: Per-area fingerprint
weights analysis for 4 different areas with varying number of beacons in Base
$B$. When the number of beacons in $B$ is 2, we observe the optimal
differences in fingerprint weights, thereby obtaining more accurate user
localizations.

Next, we study the number of beacons required for a given area by deploying
varying number of beacons ({$1$, $2$, $3$}) in a room on our test-floor ($B$
in Environment $1$ of floor-plan), to study the effect of their interaction on
the fingerprint weights, which determine the location of the user. Fig. 6
illustrates the floor-plan and the beacon deployment used for this experiment,
the results of which are illustrated in Fig. 7. With just $1$ beacon in $B$,
although clear localizations can be achieved when the user is in one of the
other three rooms, when the user is in $B$, the difference between the
fingerprint weights of the highest and the second highest values is quite
less, as illustrated by label A. This distinction becomes clearer when the
number of beacons in $B$ is increased to $2$, as shown by B. The fingerprint
weight difference between the first and second highest values also increase in
other rooms, making the localizations more accurate. When the number of
beacons are further increased to $3$, it causes localization problems in other
rooms as well, as illustrated by label C. Therefore, the optimal number of
beacons for $B$ is 2, which is the configuration used in all future
experiments. From this, we can estimate that a beacon is required for every
$40m^{2}-50m^{2}$ area, to achieve effective localization. With the timestamps
of the data packets and the fingerprint generated for each user, it is
possible to successfully localize and track them, with room-level precision,
across the environment, given its floor-plan at all time-instances. As
illustrated in Fig. 6, the user device moves to each location, in $15$-second
intervals.

Next, we evaluate the multi-user tracking capabilities of our framework by
deploying $6$ beacons (colored areas), in a multi-environment scenario as
illustrated in Fig. 8. We determine the location of the users every
$15$-seconds from $T_{0}$ to $T_{3}$ and track them through the multi-floor
multi-environment scenario by localizing their device at each time instant. As
illustrated, we successfully and simultaneously localize multiple user devices
across multiple floors in order to track them. Furthermore, we have also
evaluated the performance overheads associated with the deployed encryption
approaches, which have been presented in Table I. As anticipated, the
encryption techniques increase the latency of communication between each layer
of the network. However, the required data for the proposed localization
algorithm (Algorithm 1) is still collected and processed within the time
constraint, as discussed in Fig. 5, due to the IoT architecture proposed in
the ATLAS framework.

TABLE I: Overview of the performance overhead associated with the deployment of communication protocols with and without encryption. Communication Layer | Time (in seconds)
---|---
w/o Encryption | w/ Encryption
User Localization | 1.24 | 2.11
Database Store | 2.13 | 4.73
First Connection | 8.01 | 10.94

##  VI Conclusion and Future Work

In this work, we presented the ATLAS framework, which uses Bluetooth Low
Energy beacons in our proposed IoT architecture to anonymously localize and
track the users of the smart environment. Furthermore, to ensure the privacy
and the security of its users, we presented an easy-to-adopt secure networking
stack, which enables varying layers of encryption across the communication
stack. The ATLAS framework analyzes a combination of signal strength, hashed
user IDs, beacon UUIDs, and timestamps to effectively determine the proximal
location of each user in the environment, thereby enabling privacy-preserving
contact tracing capabilities. Our framework is open-sourced at https://atlas-
tuw.sourceforge.io to enable wide-spread adoption and reproducibility, which
can further research and development in this field.

In the next phase of our work, we propose to evaluate the scalability of our
ATLAS framework using a larger number of bluetooth beacons and smartphones,
followed by the deployment of these beacons across multiple buildings
and/other city infrastructures. This includes the evaluation of bluetooth
beacons and smartphones from other manufacturers.

## References

  * [1] L. Atzori _et al._ , “The internet of things: A survey,” _Computer networks_ , vol. 54, no. 15, pp. 2787–2805, 2010.
  * [2] S. C. Ergen, “Zigbee/ieee 802.15. 4 summary,” _UC Berkeley_ , 2004.
  * [3] G. R. Hiertz _et al._ , “The ieee 802.11 universe,” _IEEE Communications Magazine_ , vol. 48, no. 1, pp. 62–70, 2010.
  * [4] C. Bisdikian, “An overview of the bluetooth wireless technology,” _IEEE Communications magazine_ , vol. 39, no. 12, pp. 86–94, 2001.
  * [5] R. Heydon _et al._ , “Bluetooth low energy,” _Bluetooth SIG_ , 2012.
  * [6] P. Bahl _et al._ , “Radar: An in-building rf-based user location and tracking system,” in _Proceedings IEEE INFOCOM 2000._, vol. 2.
  * [7] H. Lim _et al._ , “Zero-configuration, robust indoor localization: Theory and experimentation,” Tech. Rep., 2005.
  * [8] J. J. Diaz _et al._ , “Bluepass: An indoor bluetooth-based localization system for mobile applications,” in _ISCC_. IEEE, 2010, pp. 778–783.
  * [9] A. A. N. Shirehjini _et al._ , “An rfid-based position and orientation measurement system for mobile objects in intelligent environments,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 61, 2012.
  * [10] M. Kotaru _et al._ , “Spotfi: Decimeter level localization using wifi,” in _ACM SIGCOMM_, 2015, pp. 269–282.
  * [11] I. Ashraf _et al._ , “Minloc: Magnetic field patterns-based indoor localization using convolutional neural networks,” _IEEE Access_ , 2020.
  * [12] F. Zafari _et al._ , “An ibeacon based proximity and indoor localization system,” _arXiv preprint arXiv:1703.07876_ , 2017.
  * [13] P. Martin _et al._ , “An ibeacon primer for indoor localization: demo abstract,” in _ACM BuildSys_, 2014, pp. 190–191.
  * [14] Apple AirTag, 2021. [Online]. Available: https://www.apple.com/airtag/
  * [15] Google Eddystone, 2015. [Online]. Available: https://git.io/vB9KU
  * [16] K. Fawaz _et al._ , “Protecting privacy of ble device users,” in _25th USENIX Security Symposium (USENIX Security 16)_ , 2016, pp. 1205–1221.
  * [17] C. Kolias _et al._ , “Breaking ble beacons for fun but mostly profit,” in _Proceedings of the 10th European Workshop on Systems Security_ , 2017.
  * [18] Deloitte, “Global mobile consumer survey,” 2018.
  * [19] A. Maier _et al._ , “Comparative analysis and practical implementation of the esp32 microcontroller module for the internet of things,” in _2017 Internet Technologies and Applications (ITA)_. IEEE.
  * [20] P. Voigt _et al._ , “The eu general data protection regulation (gdpr),” _A Practical Guide: Springer International Publishing_ , vol. 10, 2017.
  * [21] C. Bormann _et al._ , “CoAP: An application protocol for billions of tiny internet nodes,” _IEEE Internet Computing_ , vol. 16, no. 2, 2012.
  * [22] T. Kothmayr _et al._ , “Dtls based security and two-way authentication for the internet of things,” _Ad Hoc Networks_ , vol. 11, 2013.
  * [23] V. Seoane _et al._ , “Performance evaluation of the coap protocol with security support for iot environments,” in _ACM PE-WASUN_, 2020\.
  * [24] Neo Technology Inc., “Bolt Protocol,” 2016. [Online]. Available: https://boltprotocol.org/


Towards a study of the effects of dynamical factorization breaking at LHCb

J. D. Roth1$\star$ on behalf of the LHCb Collaboration

1 University of Michigan, Ann Arbor, USA

<EMAIL_ADDRESS>
Proceedings for the XXVIII International Workshop
on Deep-Inelastic Scattering and Related Subjects,
Stony Brook University, New York, USA, 12-16 April 2021
10.21468/SciPostPhysProc.?

## Abstract

The factorization of short-distance partonic cross-sections from universal
long-distance kinematic distributions is fundamental to phenomenology at
hadron colliders. It has been predicted however that observables sensitive to
momenta transverse to the direction of an energetic parton cannot be
factorized in the usual way, even at high energies. It should be possible to
study this factorization breaking using Z+jet production in high-energy
proton-proton collisions by studying azimuthal correlations between a Z boson
and associated charged hadrons. A plan to perform this measurement with data
collected by LHCb will be discussed, along with related work.

##  1 Introduction

Factorization theorems are useful because they allow the computation of cross-
sections for processes with non-perturbative components. Beyond their
phenomenological utility, however, factorization theorems make a formally
rigorous connection between the partonic picture and the hadronic picture of
Quantum Chromodynamics (QCD). This rigorous connection can be contrasted with
the more ad-hoc connection established by the hadron formation models included
in many Monte Carlo event generators. These models are equally useful in a
phenomenological sense, but they are not formally derived from QCD: they are
often inspired by specific features of QCD, but it is difficult or impossible
to derive higher-order corrections to these models or to clearly identify the
kinematic regions in which they fail. In regions where standard factorization
theorems start to fail, certain techniques have already been developed to
correct for effects that are typically neglected [1, 2, 3]. Hadron mass
corrections in particular have sometimes allowed analysts to use universal
parton distribution functions and fragmentation functions to fit cross-
sections that are characterized by hard scales of only a few ${GeV}$ at
colliders and fixed-target experiments with center-of-mass energies of a few
tens of ${GeV}$s. It has also been possible to develop special-purpose
factorization frameworks for certain processes, like the decay of heavy
hadrons [4], that allow for connections to perturbation theory in specific
kinematic limits.

In Section 2 we will introduce the idea of transverse-momentum-dependent (TMD)
factorization and provide a motivation for the study of TMD factorization
breaking. In Section 3 we detail the measurement that we plan to make, and in
Section 4 we present similar measurements that have already been made. We
present our conclusions in Section 5.

##  2 TMD factorization breaking

The TMD factorization framework [5] is used to compute cross-sections that
depend on the component of a hadron momentum that is transverse to the
direction of one of its constituent partons, in the limit where this
transverse momentum is much smaller than the largest energy scale that is
relevant to the scattering process. This transverse momentum is often called
$q_{T}$, and the large energy scale is called $Q$: the TMD framework was
established to compute $q_{T}$ spectra in the limit that $q_{T}\ll Q$. In
certain scattering channels, however, TMD factorization is expected to break
down in certain kinematic regions, even at very high energies. In particular:
a proof [6] has been written to show that it is not possible to use the TMD
framework to factorize the cross-section for the production of a pair of
hadrons from a proton-proton collision in the kinematic region where the two
final-state hadrons are produced nearly back-to-back in azimuth, which is a
region where $q_{T}$ tends to be very small. This breakage of factorization is
also expected to apply in back-to-back dijet production and $Z$+jet
production, where a colored parton coming out of the hard process can interact
with the beam remnants and absorb virtual emissions from initial-state
partons.

Unlike most well-known effects that complicate the factorized picture, TMD
factorization breaking is not suppressed at high energies. Therefore, at a
high-energy collider, it should be easy to isolate TMD factorization breaking
effects from any other type of factorization breaking. This breakage of TMD
factorization is also interesting because it does not generalize easily: that
is, there are certain factorizable observables that look very similar to
observables that break factorization. For example: proofs have shown that TMD
factorization can be used to compute $q_{T}$ spectra in back-to-back hadron
pair production in electron-positron annihilation and also back-to-back lepton
pair production in proton-proton collisions [5]. These processes look
superficially similar to back-to-back hadron pair production in proton-proton
collisions. In these processes, it may also be possible to use TMD
factorization to compute a wide variety of single-differential cross-sections,
in addition to $q_{T}$ spectra [7]. It is also expected to be possible to use
the collinear factorization framework to compute $q_{T}$ spectra for hadron
pair production in proton-proton collisions in the wide-angle kinematic region
where $q_{T}$ is similar in size to $Q$ [8], in which case a change of only
the kinematic region would distinguish a factorizable observable and a
factorization-breaking observable.

Cross-sections that do not factorize often tend to share certain
characteristics: it seems to become more difficult to factorize a cross-
section as the number of hadrons involved in the measurement increases, or as
the observable becomes less inclusive or more differential. But, there is not
yet any set of rules that is both strict and generally applicable that can
describe which observables factorize under which conditions [9]. For now, a
new proof must be written for more or less each observable that needs to be
factorized: factorization is handled on a case-by-case basis. Because back-to-
back production of a hadron pair in proton-proton collisions breaks
factorization and is also similar to processes that do not break
factorization, it might be used to bring attention to specific criteria that
prohibit factorization.

##  3 Plan for measurement at LHCb

We plan to measure a differential cross-section for unidentified charged
hadrons produced in association with a $Z$ boson and a jet. The same reasons
that motivate the study of dihadron production also motivate the study of
$Z$+hadron production, which is also expected to break TMD factorization. It
is also experimentally easier to extract clean $Z$-hadron correlations because
it is easy to reconstruct a $Z$ boson via its decay to $\mu^{+}\mu^{-}$: LHCb
in particular has a proven ability to precisely measure $Z$+jet cross sections
[10, 11], and the spectra of hadrons associated to a $Z$+jet pair [12]. If the
measurement is precise enough then a $Z$+hadron measurement can be compared to
a dihadron measurement to test if the number of colored partons coming out of
the hard vertex has an effect on the size of the factorization breaking.

###  3.1 Data sample and detector description

To make this measurement, we will use data that was collected by LHCb during
2016 with a $p+p$ center-of-mass energy of $13~{}{TeV}$. The LHCb detector
[13] is a single-arm forward spectrometer covering the pseudorapidity range
$2<\eta<5$. The detector includes a high-precision tracking system consisting
of a silicon-strip vertex detector around the $p+p$ interaction region, a
large area silicon-strip detector located upstream of a dipole magnet, and
three stations of silicon-strip detectors and straw drift tubes placed
downstream of the magnet. The tracking system provides for momentum
measurements with relative uncertainties that vary from $0.5\%$ at low
momentum to $1\%$ at $200~{}{GeV}/c$. We will take hadron candidates from
tracks that pass through all layers of the tracking system. The detector also
has electromagnetic and hadronic calorimeters: both tracks and calorimeter
clusters are used as input to a particle flow algorithm that is used as part
of the jet reconstruction. Muons are identified by a system composed of
alternating layers of iron and multiwire proportional chambers, and $Z$ boson
candidates will be reconstructed from high-mass muon pairs. We will require
both the jet and the $Z$ boson to have large transverse momenta: above about
$15$ or $20~{}{GeV}/c$. Only one jet will be used from each event: the jet
with largest transverse momentum. Hadrons and jets associated with the $Z$
boson should generally recoil against the $Z$ boson, so an azimuthal window of
width $\pi/3$ or $\pi/4$ on the away-side of the $Z$ boson will be established
as a signal region: the rest of azimuth will be used for background
estimation.

###  3.2 Kinematic variables

The hadronic cross-section will be binned in three kinematic variables:

$p_{out}\equiv p_{T}^{h^{\pm}}\sin\left(\phi^{h^{\pm}}-\phi^{jet}\right),$ (1)
$Q\equiv\left(p^{Z}+p^{jet}\right)^{2},$ (2)

and

$z\equiv{\mathbf{p}^{h^{\pm}}\cdot\mathbf{p}^{jet}\over\mathbf{p}^{jet}\cdot\mathbf{p}^{jet}}~{}.$
(3)

Here, $p_{T}^{h^{\pm}}$ is the transverse momentum of the charged hadron, and
$\phi^{h^{\pm}}$ and $\phi^{Z}$ are the azimuthal coordinates of the hadron
and the $Z$ boson. The four-momenta of the $Z$ boson and the jet are written
as $p^{Z}$ and $p^{jet}$, and the three-momenta of the jet and the hadron are
written as $\mathbf{p}^{jet}$ and $\mathbf{p}^{h^{\pm}}$. In order to mitigate
uncertainies that are associated with the beam luminosity and $Z$+jet
reconstruction, we normalize the hadronic cross-section by the cross-section
for $Z$+jet production:

$\left.{{d}N^{h^{\pm}}\over{d}p_{out}\,{d}z\,{d}Q}\middle/{{d}N^{Z{+jet}}\over{d}Q}\right.~{}.$
(4)

The variable $p_{out}$ is used to probe transverse momenta generated by the
parton shower and by long-range dynamics in both the initial state and final
state. To a first approximation, none of $p_{out}$ comes from transverse
momentum generated at the hardest scales: in the partonic center-of-mass
frame, the $Z$ boson and the outgoing parton come out exactly back-to-back.
Any transverse momentum imbalance is generated from processes characterized by
smaller energy scales. Some of the imbalance comes from transverse motion of
the initial-state partons inside the protons, much of the imbalance is
generated by splittings in the parton showers in the initial and final state,
and some comes from the hadron formation process where partons that are
separated in azimuth exert forces on each other. Note that $p_{out}$ ignores
components of the transverse momentum imbalance along the axis determined by
the $Z$ boson direction of motion, which is roughly the same as the jet axis.
The imbalance along this axis should be determined mostly by collinear aspects
of the parton splitting and fragmentation process: with a focus on the off-
axis component of the imbalance, we hope that we can improve our sensitivity
to uniquely transverse-momentum-dependent effects. This is important because
factorization is not expected to break in the collinear framework.

The $Z$+jet mass $Q$ is a proxy for the hard scale that characterizes the
scattering. The way that the $p_{out}$ distribution changes with the hard
scale is described by Collins-Soper-Sterman (CSS) evolution [14, 15]. CSS
evolution is similar to a TMD variant of Dokshitzer-Gribov-Lipatov-Altarelli-
Parisi (DGLAP) evolution [16, 17, 18], with the notable exception that the CSS
evolution kernels have both perturbative and universal non-perturbative
components while DGLAP’s evolution kernels can be computed entirely in
perturbation theory. CSS evolution takes as input the set of kinematic
distributions defined in the TMD framework, probed at one hard scale, and
describes how those distributions look at another hard scale. In particular:
the results of CSS evolution are only valid when the technique is applied to a
factorizable distribution defined in the TMD framework. Hence: if we can test
whether or not CSS evolution can correctly model the relationship between
$p_{out}$ distributions measured at different values of the hard scale $Q$,
then we can test whether or not TMD factorization holds. The alternative is to
compute $p_{out}$ distributions using TMD distributions that have already been
extracted from fits to $e^{+}e^{-}$ and $ep$ scatterings: but, as of now, no
such fits have been extracted with very good precision.

We also want to bin the hadronic cross-section in the fragmentation variable
$z$ in order to allow calculations that compare to this measurement to exclude
the low-$z$ and high-$z$ regions, where calculation is difficult. The $z$
variable will also provide a more complete or differential picture of the
hadron formation process, which might improve the power of the analysis if
factorization breaks most strongly in a sub-region of the full phase space.

In addition to measuring the differential cross-section, we plan to fit the
$p_{out}$ distribution in each bin of $Q\otimes z$ in order to show clearly
how the shapes of the distributions evolve with $Q$ in each bin of $z$. These
fits will hopefully make it easier to make a computation that determines
whether or not this measurement is consistent with the predictions of CSS
evolution.

##  4 Prior measurements

Qualitatively, we have a good idea of what to expect from this measurement
because similar measurements have been made at PHENIX [19, 20, 21]. Because
PHENIX uses $p+p$ collisions with lower beam energies than LHCb, which allow
less phase space for $Z$ boson production, they measured $\pi^{0}$-hadron
correlations and direct photon-hadron correlations instead of $Z$-hadron
correlations. Those measurements also used a slightly different set of
kinematic variables, since jet reconstruction was not feasible at PHENIX due
to a limited acceptance. In order to characterize the hard scale of the
scatterings, PHENIX used the transverse momentum of the $\pi^{0}$ or photon
instead of an invariant mass-type variable like $Q$, and in order to estimate
a fragmentation variable like $z$ they used $x_{E}\equiv
p_{T}^{h^{\pm}}/p_{T}^{\pi^{0},\,\gamma}\cdot\cos\left(\phi^{\pi^{0},\,\gamma}-\phi^{h^{\pm}}\right)$.
Some of their results are shown in Figure 1. The cores of their $p_{out}$
distributions can be fit to Gaussian shapes, but further out in their wings
the distributions fall off too slowly to fit an exponential: this is
consistent with expectations from perturbation theory for the power-law fall-
off of energetic gluon radiation densities. Qualitatively, the fitted Gaussian
widths of the $p_{out}$ distributions increased with the hard scale. This
result matches the predictions of CSS evolution: as more energy becomes
available to the particles involved in the scattering, more transverse
momentum is generated in the parton showers. No calculation has yet been
published that compares the PHENIX measurement to the quantitative predictions
of CSS evolution.

Figure 1: On the left are shown some of the $p_{out}$ distributions extracted
by the PHENIX collaboration. The solid lines show Gaussian fits to the data
points. On the right, the Gaussian widths are plotted against the hard scale.
Qualitatively, the widths increase with the hard scale. Both plots are from
Reference [20].

##  5 Conclusion

We want to quantify the breakdown of transverse-momentum-dependent
factorization, which is a fundamental prediction of QCD that has not yet been
verified. We hope that investigation into the regions where factorization
fails might inspire techniques that extend a type of factorization to
observables that cannot currently be factorized. In addition, these
investigations might help to develop a more general set of rules to determine
which processes do and do not factorize. Measurements that might be sensitive
to the breakdown of factorization have already been made by the PHENIX
collaboration, and we plan to make another set of measurements with data from
LHCb. These measurements allow us to test for the breakdown of factorization
via the breakdown of CSS evolution. We hope that the increased availability of
measurements from a variety of energy and rapidity ranges will encourage the
calculation of a quantitative comparison between these measurements and the
predictions of CSS evolution, especially with the improved array of kinematic
variables that the LHCb measurement will use to parameterize the scattering
and fragmentation process.

#### Funding information

This material is based upon work supported by the National Science Foundation
under Grant No. 2012926.

## References

  * [1] J. V. Guerrero, A. Accardi, Gauge Invariance and Kaon Production in Deep Inelastic Scattering at Low Scales, Phys. Rev. D 97 114012 (2018), 10.1103/PhysRevD.97.114012.
  * [2] T. B. Liu, J. W. Qiu, Power Corrections in Semi-Inclusive Deep Inelastic Scatterings at Fixed-Target Energies, Phys. Rev. D 101 014008 (2020), 10.1103/PhysRevD.101.014008.
  * [3] B. Z. Kopeliovich, R. Pasechnik, I. K. Potashnikova, Diffractive Dijet Production: Breakdown of Factorization, Phys. Rev. D 98 114021 (2018), 10.1103/PhysRevD.98.114021.
  * [4] M. Neubert, Heavy Quark Symmetry, Phys. Rept. 245 259 (1994), 10.1016/0370-1573(94)90091-4.
  * [5] J. Collins, Foundations of Perturbative QCD, Cambridge Monographs on Particle Physics, Nuclear Physics and Cosmology; Cambridge University Press: Cambridge (2011), 10.1017/CBO9780511975592.
  * [6] T. Rogers, P. Mulders, No Generalized TMD-Factorization in the Hadro-Production of High Transverse Momentum Hadrons, Phys. Rev. D 81 094006 (2010), 10.1103/PhysRevD.81.094006.
  * [7] M. Schwartz, K. Yan, H. X. Zhu, Factorization Violation and Scale Invariance, Phys. Rev. D 97 096017 (2018), 10.1103/PhysRevD.97.096017.
  * [8] J. Collins, J. W. Qiu, $k_{T}$-Factorization is Violated in Production of High-Transverse-Momentum Particles in Hadron-Hadron Collisions, Phys. Rev. D 75 114014 (2007), 10.1103/PhysRevD.75.114014.
  * [9] S. Catani, D. de Florian, G. Rodrigo, Space-Like (vs. Time-Like) Collinear Limits in QCD: is Factorization Violated?, JHEP 07 026 (2012), 10.1007/JHEP07(2012)026.
  * [10] LHCb Collaboration, Study of Forward $Z$+Jet Production in $pp$ Collisions at $\sqrt{s}=7$ TeV, JHEP 01 33 (2014), 10.1007/JHEP01(2014)033.
  * [11] LHCb Collaboration, Measurement of Forward $W$ and $Z$ Boson Production in Association with Jets in Proton-Proton Collisions at $\sqrt{s}=8$ TeV, JHEP 05 131 (2016), 10.1007/JHEP05(2016)131.
  * [12] LHCb Collaboration, Measurement of Charged Hadron Production in $Z$-Tagged Jets in Proton-Proton Collisions at $\sqrt{s}=8$ TeV, Phys. Rev. Lett. 123 232001 (2019), 10.1103/PhysRevLett.123.232001.
  * [13] LHCb Collaboration, The LHCb Detector at the LHC, JINST 3 S08005 (2008), 10.1088/1748-0221/3/08/S08005.
  * [14] J. C. Collins, D. E. Soper, Back-to-Back Jets in QCD, Nucl. Phys. B 193 381 (1981), Erratum: Nucl. Phys. B 213 545 (1983), 10.1016/0550-3213(81)90339-4.
  * [15] J. C. Collins, D. E. Soper, Parton Distribution and Decay Functions, Nucl. Phys. B 194 445 (1982), 10.1016/0550-3213(82)90021-9.
  * [16] V. N. Gribov, L. N. Lipatov, Deep Inelastic $e$ $p$ Scattering in Perturbation Theory, Sov. J. Nucl. Phys. 15 438 (1972), Yad. Fiz. 15 781 (1972).
  * [17] G. Altarelli, G. Parisi, Asymptotic Freedom in Parton Language, Nucl. Phys. B 126 298 (1977), 10.1016/0550-3213(77)90384-4.
  * [18] Y. L. Dokshitzer, Calculation of the Structure Functions for Deep Inelastic Scattering and $e^{+}e^{-}$ Annihilation by Perturbation Theory in Quantum Chromodynamics, Sov. Phys. JETP 46 641 (1977), Zh. Eksp. Teor. Fiz. 73 1216 (1977).
  * [19] PHENIX Collaboration, Nonperturbative-Transverse-Momentum Effects and Evolution in Dihadron and Direct Photon-Hadron Angular Correlations in $p+p$ Collisions at $\sqrt{s}=510~{}{GeV}$, Phys. Rev. D 95 072002 (2017), 10.1103/PhysRevD.95.072002.
  * [20] PHENIX Collaboration, Nonperturbative Transverse-Momentum-Dependent Effects in Dihadron and Direct Photon-Hadron Angular Correlations in $p+p$ Collisions at $\sqrt{s}=200~{}{GeV}$, Phys. Rev. D 98 072004 (2018), 10.1103/PhysRevD.98.072004.
  * [21] PHENIX Collaboration, Nonperturbative Transverse Momentum Broadening in Dihadron Angular Correlations in $\sqrt{s}=200~{}{GeV}$ Proton-Nucleus Collisions, Phys. Rev. C 99 044912 (2019), 10.1103/PhysRevC.99.044912.


$\begin{split}D^{\prime}_{\alpha}=D_{\alpha}-U_{\alpha}\frac{\frac{\partial
p}{\partial T}\frac{\partial n}{\partial\mu}-\frac{\partial
p}{\partial\mu}\frac{\partial n}{\partial T}}{\frac{\partial\rho}{\partial
T}\frac{\partial
n}{\partial\mu}-\frac{\partial\rho}{\partial\mu}\frac{\partial n}{\partial
T}}-N_{\alpha}\frac{\frac{\partial p}{\partial
T}\frac{\partial\rho}{\partial\mu}-\frac{\partial
p}{\partial\mu}\frac{\partial\rho}{\partial T}}{\frac{\partial\rho}{\partial
T}\frac{\partial
n}{\partial\mu}-\frac{\partial\rho}{\partial\mu}\frac{\partial n}{\partial
T}},\\\ D^{\prime}_{w}=D_{w}-U_{w}\frac{\frac{\partial p}{\partial
T}\frac{\partial n}{\partial\mu}-\frac{\partial p}{\partial\mu}\frac{\partial
n}{\partial T}}{\frac{\partial\rho}{\partial T}\frac{\partial
n}{\partial\mu}-\frac{\partial\rho}{\partial\mu}\frac{\partial n}{\partial
T}}-N_{w}\frac{\frac{\partial p}{\partial
T}\frac{\partial\rho}{\partial\mu}-\frac{\partial
p}{\partial\mu}\frac{\partial\rho}{\partial T}}{\frac{\partial\rho}{\partial
T}\frac{\partial
n}{\partial\mu}-\frac{\partial\rho}{\partial\mu}\frac{\partial n}{\partial
T}}.\end{split}$

##  Appendix B Relations between coefficients

Herein we derive the relations between coefficients (31) enforcing the
continuity equation for the mean value of the stress-energy tensor
$T_{\mu\nu}$ at second order in the thermal vorticity $\varpi$

$\begin{split}T_{\mu\nu}=&(\rho-\alpha^{2}U_{\alpha}-w^{2}U_{w})u_{\mu}u_{\nu}-(p-\alpha^{2}D_{\alpha}-w^{2}D_{w})\Delta_{\mu\nu}+A\alpha_{\mu}\alpha_{\nu}+Ww_{\mu}w_{\nu}+G(u_{\mu}\gamma_{\nu}+u_{\nu}\gamma_{\mu}).\end{split}$
(88)

Firstly, we observe that scalar thermodynamic functions depend on spacetime
coordinates only through the magnitude of the four-temperature
$F(|\beta|)=F(\sqrt{\beta^{2}(x)})$, thus:

$\partial_{\nu}F(|\beta|)=\frac{\partial|\beta|}{\partial
x^{\nu}}\frac{\partial
F(|\beta|)}{\partial|\beta|}=\frac{1}{2|\beta|}\partial_{\nu}\big{(}\beta^{\lambda}\beta_{\lambda}\big{)}\frac{\partial
F(|\beta|)}{\partial|\beta|}=\frac{1}{|\beta|}\beta^{\lambda}\partial_{\nu}\beta_{\lambda}\frac{F(|\beta|)}{\partial|\beta|}=-u^{\lambda}\varpi_{\nu\lambda}\frac{\partial
F(|\beta|)}{\partial|\beta|}=-\alpha_{\nu}\frac{\partial
F(|\beta|)}{\partial|\beta|}$ (89)

where we have used the definition of $u$ in eq. (3), the eqs. (6) and (20).

We start by reckoning the gradients of the four-vectors
$\\{u,\alpha,w,\gamma\\}$ which are needed to calculate the stress-energy
tensor divergence. First, we observe that, because of the Killing equation
(4), and using the eqs. (3), (6) and (20):

$0=\beta^{\mu}(\partial_{\mu}\beta_{\nu}+\partial_{\nu}\beta_{\mu})=|\beta|u^{\mu}\varpi_{\nu\mu}+\frac{1}{2}\partial_{\nu}{\beta^{2}}=|\beta|\alpha_{\nu}+\frac{1}{2}\partial_{\nu}{\beta^{2}}$

so that:

$\alpha_{\nu}=-\frac{1}{2|\beta|}\partial_{\nu}{\beta^{2}}$ (90)

and, contracting with $u$ again:

$u^{\nu}\partial_{\nu}\beta^{2}=D\beta^{2}=0$ (91)

what we already saw in section (IV). This equation implies that any scalar
function, whose argument are $\beta^{2}$ and $\xi$, has a vanishing derivative
along the flow, that is:

$DF(\beta^{2},\xi)=0.$ (92)

Now, we can find the gradient of $u$ as defined by the eq. (3):

$\partial_{\nu}u_{\mu}=\frac{\partial_{\nu}\beta_{\mu}}{|\beta|}+\beta_{\mu}\partial_{\nu}\Big{(}\frac{1}{|\beta|}\Big{)}=\frac{1}{|\beta|}\big{(}\varpi_{\mu\nu}+\alpha_{\nu}u_{\mu}\big{)},$
(93)

where we have used again the eq. (3) and the eq. (90). The divergence of $u$
then is:

$\partial_{\mu}u^{\mu}=\partial_{\mu}\frac{\beta^{\mu}}{|\beta|}=0$ (94)

because of the Killing vector equation (4), which obviously imply
$\partial_{\mu}\beta^{\mu}=0$ and the (91). Instead the derivative along its
direction is

$u_{\rho}\partial^{\rho}u_{\mu}=\frac{u_{\rho}}{|\beta|}\big{(}\varpi_{\mu}^{\,\rho}+\alpha^{\rho}u_{\beta}\big{)}=\frac{\varpi_{\mu\rho}u^{\rho}}{|\beta|}=\frac{\alpha_{\mu}}{|\beta|}.$
(95)

Then, let us calculate the derivative of $\alpha$, keeping in mind that
$\varpi$ is constant:

$\partial_{\mu}\alpha_{\nu}=\varpi_{\nu\rho}\partial_{\mu}u^{\rho}=\frac{\varpi_{\nu\rho}}{|\beta|}\big{(}\varpi^{\rho}_{\phantom{\mu}\mu}+\alpha_{\mu}u^{\rho}\big{)}=\frac{1}{|\beta|}\big{(}\varpi_{\nu\rho}\varpi^{\rho}_{\phantom{\mu}\mu}+\alpha_{\mu}\alpha_{\nu}\big{)};$
(96)

using the thermal vorticity decomposition (19) and the Levi-Civita tensor
properties we can express the previous formula in terms of the tetrad vectors

$\partial_{\mu}\alpha_{\nu}=\frac{1}{|\beta|}\left(-w_{\mu}w_{\nu}+\gamma_{\mu}u_{\nu}-u_{\mu}\gamma_{\nu}-(w^{2}+\alpha^{2})u_{\mu}u_{\nu}+g_{\mu\nu}w^{2}\right),$
(97)

whence we obtain the divergence:

$\partial\cdot\alpha=\frac{1}{|\beta|}\left(w^{2}-(w^{2}+\alpha^{2})+4w^{2}\right)=\frac{1}{|\beta|}\big{(}2w^{2}-\alpha^{2}\big{)},$

as well as the gradient of $\alpha^{2}$:

$\partial_{\mu}\alpha^{2}=2\alpha^{\nu}\partial_{\mu}\alpha_{\nu}=\frac{2}{|\beta|}(-\alpha\cdot
w\,w_{\mu}+w^{2}\alpha_{\mu})$ (98)

taking into account that $\alpha\cdot u=\alpha\cdot\gamma=0$.

Likewise, we can calculate the derivative of $w$ by using its definition (20)
and the (93):

$\partial_{\mu}w_{\nu}=-\frac{1}{2}\epsilon_{\nu\rho\sigma\lambda}\varpi^{\rho\sigma}\partial_{\mu}u^{\lambda}=\frac{-1}{2|\beta|}\epsilon_{\nu\rho\sigma\lambda}\varpi^{\rho\sigma}\varpi^{\lambda}_{\phantom{\mu}\mu}+\frac{\alpha_{\mu}w_{\nu}}{|\beta|}\,.$

Replacing the eq. (19) in the first term of the right hand side and using the
properties of the Levi-Civita tensor, it can be shown that:

$\partial_{\mu}w_{\nu}=\frac{1}{|\beta|}\left(-g_{\mu\nu}\alpha\cdot
w+\alpha_{\mu}w_{\nu}\right)$ (99)

so that its divergence is

$\partial_{\mu}w^{\mu}=-\frac{3}{|\beta|}(w\cdot\alpha)$ (100)

and the gradient of $w^{2}$:

$\partial_{\mu}w^{2}=2w^{\nu}\partial_{\mu}w_{\nu}=\frac{2}{|\beta|}(-\alpha\cdot
w\,w_{\mu}+w^{2}\alpha_{\mu})=\partial_{\mu}\alpha^{2}.$ (101)

In order to calculate the derivative of the last relevant vector field
$\gamma$ (21), one can first show that, by using
$\Delta_{\lambda\mu}=g_{\lambda\mu}-u_{\lambda}u_{\mu}$ and (19), it can be
expressed as:

$\gamma^{\mu}=\alpha^{\rho}\varpi_{\rho\lambda}\Delta^{\lambda\mu}=(\alpha\cdot\varpi)^{\mu}-\varpi_{\rho\lambda}\alpha^{\rho}u^{\lambda}u^{\mu}=(\alpha\cdot\varpi)^{\mu}-(\alpha_{\rho}u_{\lambda}-\alpha_{\lambda}u_{\rho})\alpha^{\rho}u^{\lambda}u^{\mu}=(\alpha\cdot\varpi)^{\mu}-\alpha^{2}u^{\mu}$

so that its divergence vanishes

$\partial_{\mu}\gamma^{\mu}=\partial_{\mu}\big{[}(\alpha\cdot\varpi)^{\mu}-\alpha^{2}u^{\mu}\big{]}=\varpi^{\rho\mu}\partial_{\mu}\alpha_{\rho}=\frac{1}{|\beta|}\underbrace{\varpi^{\rho\mu}}_{\textup{Antisym}}\big{(}\underbrace{\varpi_{\rho\lambda}\varpi^{\lambda}_{\phantom{\mu}\mu}}_{\textup{Sym}}+\underbrace{\alpha_{\mu}\alpha_{\rho}}_{\textup{Sym}}\big{)}=0,$

where we have used the eq. (96). Another useful relation involving $\gamma$ is
the contraction:

$\gamma^{\rho}\varpi_{\rho\kappa}=\epsilon^{\rho\mu\nu\sigma}w_{\mu}\alpha_{\nu}u_{\sigma}\varpi_{\rho\kappa}=\epsilon^{\rho\mu\nu\sigma}\epsilon_{\rho\kappa\lambda\tau}w^{\lambda}u^{\tau}w_{\mu}\alpha_{\nu}u_{\sigma}+\epsilon^{\rho\mu\nu\sigma}\alpha_{\rho}u_{\kappa}w_{\mu}\alpha_{\nu}u_{\sigma}-\epsilon^{\rho\mu\nu\sigma}\alpha_{\kappa}u_{\rho}w_{\mu}\alpha_{\nu}u_{\sigma},$

where we have used the decomposition (19). Then, because of the Levi-Civita
tensor properties we find:

$\gamma^{\rho}\varpi_{\rho\kappa}=w^{2}\alpha_{\kappa}-(\alpha\cdot
w)w_{\kappa}.$ (102)

We also need the derivative of the transverse projector $\Delta$:

$\partial_{\mu}\Delta_{\rho\sigma}=\partial_{\mu}\big{(}g_{\rho\sigma}-u_{\rho}u_{\sigma}\big{)}=-u_{\sigma}\partial_{\mu}u_{\rho}-u_{\rho}\partial_{\mu}u_{\sigma}$

whence:

$\partial^{\rho}\Delta_{\rho\sigma}=-u_{\sigma}\partial^{\rho}u_{\rho}-u_{\rho}\partial^{\rho}u_{\sigma}=-\frac{\alpha_{\sigma}}{|\beta|}.$
(103)

Finally, we observe that

$D\alpha^{2}=Dw^{2}=0$

as they are scalar functions.

To calculate the divergence of the stress-energy tensor we need to work out
some intermediate relations involving the derivatives of the four-vectors
$u,\alpha,w,\gamma$. The first relation can be obtained by using the (103) and
(98):

$\partial^{\mu}\left(\alpha^{2}\Delta_{\mu\nu}\right)=\alpha^{2}\partial^{\mu}\Delta_{\mu\nu}+\Delta_{\mu\nu}\partial^{\mu}\alpha^{2}=-\frac{\alpha^{2}}{|\beta|}\alpha_{\nu}+\partial_{\nu}\alpha^{2}=\frac{1}{|\beta|}\left((2w^{2}-\alpha^{2})\alpha_{\nu}-2(\alpha\cdot
w)w_{\nu}\right).$ (104)

The second by using the (103) and (101):

$\partial^{\mu}\left(w^{2}\Delta_{\mu\nu}\right)=\Delta_{\mu\nu}\partial^{\mu}w^{2}+w^{2}\partial^{\mu}\Delta_{\mu\nu}=\frac{2}{|\beta|}\left(w^{2}\alpha_{\nu}-(\alpha\cdot
w)w_{\nu}\right)-\frac{w^{2}\alpha_{\nu}}{|\beta|}=\frac{1}{|\beta|}\left(w^{2}\alpha_{\nu}-2(\alpha\cdot
w)w_{\nu}\right).$ (105)

Moreover with (97) and orthogonality properties:

$\alpha_{\mu}\partial^{\mu}\alpha_{\nu}=\frac{\alpha_{\mu}}{|\beta|}\left(-w^{\mu}w_{\nu}+\gamma^{\mu}u_{\nu}-u^{\mu}\gamma_{\nu}-(w^{2}+\alpha^{2})u^{\mu}u_{\nu}+g^{\mu}_{\,\nu}w^{2}\right)=\frac{w^{2}\alpha_{\nu}-(\alpha\cdot
w)w_{\nu}}{|\beta|};$ (106)

and, using (99):

$w_{\mu}\partial^{\mu}w_{\nu}=\frac{w_{\mu}}{|\beta|}\left(-g^{\mu}_{\hphantom{\mu}\nu}(\alpha\cdot
w)+\alpha^{\mu}w_{\nu}\right)=\frac{1}{|\beta|}\left(-w_{\nu}(\alpha\cdot
w)+(\alpha\cdot w)w_{\nu}\right)=0.$ (107)

Futhermore, taking advantage of (93) and (102):

$\gamma_{\mu}\partial^{\mu}u_{\nu}=\frac{\gamma_{\mu}}{|\beta|}(\varpi_{\nu}^{\,\mu}+\alpha^{\mu}u_{\nu})=-\frac{\gamma^{\rho}\varpi_{\rho\nu}}{|\beta|}=\frac{(\alpha\cdot
w)w_{\nu}-w^{2}\alpha_{\nu}}{|\beta|}.$ (108)

The last needed relation involves the gradients of $\gamma$ in eq. (21). By
expanding its definition and using the previous eqs. (93,97,99):

$\begin{split}u_{\mu}\partial^{\mu}\gamma_{\nu}&=\epsilon_{\nu\rho\sigma\tau}u_{\mu}\partial^{\mu}\left(w^{\rho}\alpha^{\sigma}u^{\tau}\right)\\\
&=\epsilon_{\nu\rho\sigma\tau}u_{\mu}(\partial^{\mu}w^{\rho})\alpha^{\sigma}u^{\tau}+\epsilon_{\nu\rho\sigma\tau}u_{\mu}w^{\rho}(\partial^{\mu}\alpha^{\sigma})u^{\tau}+\epsilon_{\nu\rho\sigma\tau}u_{\mu}w^{\rho}\alpha^{\sigma}\partial^{\mu}u^{\tau}\\\
&=-\epsilon_{\nu\rho\sigma\tau}u_{\rho}\alpha^{\sigma}u^{\tau}-\epsilon_{\nu\rho\sigma\tau}w^{\rho}(\gamma^{\sigma}-\alpha^{2}u^{\sigma})u^{\tau}+\epsilon_{\nu\rho\sigma\tau}w^{\rho}\alpha^{\sigma}\alpha^{\tau}\\\
&=-\epsilon_{\nu\rho\sigma\tau}w^{\rho}\gamma^{\sigma}u^{\tau}=-\epsilon_{\nu\rho\sigma\tau}\epsilon^{\sigma\kappa\lambda\pi}w^{\rho}u^{\tau}w_{\kappa}u_{\pi}=\frac{1}{|\beta|}(w_{\nu}(\alpha\cdot
w)-\alpha_{\nu}w^{2}-\alpha_{\nu}\alpha^{2}),\end{split}$ (109)

where the known contractions of Levi-Civita tensor have been employed in the
last equality.

We are now in a position to enforce the continuity equation
$\partial^{\mu}T_{\mu\nu}=0$, to the expression of the stress energy tensor
(88). The first term involving energy density gives rise to:

$\partial^{\mu}\big{(}\rho\,u_{\mu}u_{\nu}\big{)}=u_{\mu}\big{(}\partial^{\mu}\rho\big{)}u_{\nu}+\rho\,u_{\nu}\,\partial^{\mu}u_{\mu}+\rho\,u_{\mu}\partial^{\mu}u_{\nu}=\rho\frac{\alpha_{\nu}}{|\beta|}$
(110)

using eqs. (89), (94) and (95). Then, keeping in mind the (92) which implies
that $DU_{\alpha}=DU_{w}=\ldots=0$ for all coefficients, and using (89), (94),
(95), as well as $D\alpha^{2}=Dw^{2}=0$:

$\displaystyle\partial^{\mu}\big{(}U_{\alpha}\,\alpha^{2}u_{\mu}u_{\nu}\big{)}=DU_{\alpha}\,\alpha^{2}u_{\nu}+U_{\alpha}\Big{[}\alpha^{2}u_{\nu}\partial^{\mu}u_{\mu}+u_{\nu}D\alpha^{2}+\alpha^{2}u_{\mu}\partial^{\mu}u_{\nu}\Big{]}=U_{\alpha}\frac{\alpha^{2}\,\alpha_{\nu}}{|\beta|};$
(111)
$\displaystyle\partial^{\mu}\big{(}U_{w}\,w^{2}u_{\mu}u_{\nu}\big{)}=DU_{w}\,w^{2}u_{\nu}+U_{w}\Big{[}w^{2}u_{\nu}\partial^{\mu}u_{\mu}+u_{\nu}Dw^{2}+w^{2}u_{\mu}\partial^{\mu}u_{\nu}\Big{]}=U_{w}\frac{w^{2}\,\alpha_{\nu}}{|\beta|}.$
(112)

Now, using eqs. (89) and (103):

$\partial^{\mu}\big{(}p\,\Delta_{\mu\nu}\big{)}=\Delta_{\mu\nu}\partial^{\mu}p+p\partial^{\mu}\Delta_{\mu\nu}=-\Big{(}p+|\beta|\frac{\partial
p}{\partial|\beta|}\Big{)}\frac{\alpha_{\nu}}{|\beta|};$ (113)

and, thanks to eqs. (89) and (104):

$\partial^{\mu}\big{(}D_{\alpha}\,\alpha^{2}\Delta_{\mu\nu}\big{)}=\Delta_{\mu\nu}\alpha^{2}\partial^{\mu}D_{\alpha}+D_{\alpha}\partial^{\mu}\big{(}\alpha^{2}\Delta_{\mu\nu}\big{)}=-\Big{(}|\beta|\frac{\partial
D_{\alpha}}{\partial|\beta|}+D_{\alpha}\Big{)}\frac{\alpha^{2}\alpha_{\nu}}{|\beta|}+2D_{\alpha}\frac{w^{2}\alpha_{\nu}}{|\beta|}-2D_{\alpha}\frac{(\alpha\cdot
w)w_{\nu}}{|\beta|};$ (114)

likewise, because of (89) and (105):

$\partial^{\mu}\big{(}D_{w}w^{2}\Delta_{\mu\nu}\big{)}=w^{2}\Delta_{\mu\nu}\partial^{\mu}D_{w}+D_{w}\partial^{\mu}\big{(}w^{2}\Delta_{\mu\nu}\big{)}=\Big{(}D_{w}-|\beta|\frac{\partial
D_{w}}{\partial|\beta|}\Big{)}\frac{w^{2}\alpha_{\nu}}{|\beta|}-2D_{w}\frac{(\alpha\cdot
w)w_{\nu}}{|\beta|}.$ (115)

Furthermore, using (89), (94) and (106):

$\partial^{\mu}\big{(}A\,\alpha_{\mu}\alpha_{\nu}\big{)}=\alpha_{\mu}\alpha_{\nu}\big{(}\partial^{\mu}A\big{)}+A\,\alpha_{\nu}\partial^{\mu}\alpha_{\mu}+A\,\alpha_{\mu}\partial^{\mu}\alpha_{\nu}=-\big{(}|\beta|\frac{\partial
A}{\partial|\beta|}+A\big{)}\frac{\alpha^{2}\alpha_{\nu}}{|\beta|}+3A\frac{w^{2}\alpha_{\nu}}{|\beta|}-A\frac{(\alpha\cdot
w)\alpha_{\nu}}{|\beta|}$ (116)

and, similarly, by means of (89), (100) and (107):

$\partial^{\mu}\big{(}Ww_{\mu}w_{\nu}\big{)}=(\partial^{\mu}W)w_{\mu}w_{\nu}+Ww_{\nu}\partial^{\mu}w_{\mu}+Ww_{\mu}\partial^{\mu}w_{\nu}=-\Big{(}|\beta|\frac{\partial
W}{\partial|\beta|}+3W\Big{)}\frac{(\alpha\cdot w)w_{\nu}}{|\beta|}.$ (117)

Finally, using again (89) and the relations (108) and (109) we can determine
the term involving the gradients of $\gamma$:

$\partial^{\mu}\left[G(u_{\mu}\gamma_{\nu}+u_{\nu}\gamma_{\mu})\right]=(u_{\mu}\gamma_{\nu}+u_{\nu}\gamma_{\mu})\partial^{\mu}G+Gu_{\mu}\partial^{\mu}\gamma_{\nu}+G\gamma_{\mu}\partial^{\mu}u_{\nu}=2G\frac{(\alpha\cdot
w)w_{\nu}}{|\beta|}-2G\frac{w^{2}\alpha_{\nu}}{|\beta|}.$ (118)

The divergence of the stress-energy tensor is obtained by summing all right-
hand sides of the eqs. (110-118). As it can be readily checked, the resulting
expression involves the sum of terms multiplying
$\alpha_{\nu},\,\alpha^{2}\alpha_{\nu},\,w^{2}\alpha_{\nu}$ and $(\alpha\cdot
w)w_{\nu}$. As the divergence should vanish independently of the vectors $w$
and $\alpha$, which have as many degrees of freedom as the thermal vorticity
(that is 6), the conclusion is that each coefficient of the above combinations
must be zero. As a result, four equations are obtained. The coefficient of
$\alpha_{\nu}$ gives rise to:

$\rho+p+|\beta|\frac{\partial p}{\partial|\beta|}\Big{|}_{\zeta}=0,$

which is but the well known thermodynamic relation between energy density and
pressure. The other three coefficients yield three conditions on the second-
order coefficients:

$U_{\alpha}=-|\beta|\frac{\partial}{\partial|\beta|}\big{(}D_{\alpha}+A\big{)}-\big{(}D_{\alpha}+A\big{)};$
$U_{w}=-|\beta|\frac{\partial}{\partial|\beta|}D_{w}+(2D_{\alpha}+D_{w}+3A)-2G;$
(119)

and:

$2G=2\big{(}D_{\alpha}+D_{w}\big{)}+A+|\beta|\frac{\partial}{\partial|\beta|}W+3W.$
(120)

Subtracting eq. (120) from eq. (119) we can simplify some terms and obtain:

$U_{w}=-|\beta|\frac{\partial}{\partial|\beta|}\big{(}D_{w}+W\big{)}-D_{w}+2A-3W,$

which completes the proof of the relations in eq. (31).

## References

  * (1) R. Baier, P. Romatschke and U. A. Wiedemann, Phys. Rev. C 73 (2006) 064903.
  * (2) P. Huovinen and P. V. Ruuskanen, Ann. Rev. Nucl. Part. Sci. 56 (2006) 163.
  * (3) W. Florkowksi, Phenomenology of Ultra-Relativistic heavy ion collisions, World Scientific (2010).
  * (4) C. Gale, S. Jeon and B. Schenke, Int. J. Mod. Phys. A 28 (2013) 1340011.
  * (5) U. Heinz and R. Snellings, Ann. Rev. Nucl. Part. Sci. 63 (2013) 123.
  * (6) R. Derradi de Souza, T. Koide and T. Kodama, Prog. Part. Nucl. Phys. 86 (2016) 35.
  * (7) A. Muronga, Phys. Rev. C 76 (2007) 014910.
  * (8) K. Tsumura, T. Kunihiro and K. Ohnishi, [Phys. Lett. B 646 (2007) 134] Erratum: Phys. Lett. B 656 (2007) 274.
  * (9) B. Betz, D. Henkel and D. H. Rischke, Prog. Part. Nucl. Phys. 62 (2009) 556.
  * (10) M. A. York and G. D. Moore, Phys. Rev. D 79 (2009) 054011.
  * (11) A. Monnai and T. Hirano, Nucl. Phys. A 847 (2010) 283.
  * (12) G. S. Denicol, H. Niemi, E. Molnar and D. H. Rischke, Phys. Rev. D 85 (2012) 114047 [Erratum-ibid. D 91 (2015) 039902].
  * (13) P. Van and T. S. Biro, Phys. Lett. B 709 (2012) 106.
  * (14) A. Jaiswal, R. S. Bhalerao and S. Pal, Phys. Lett. B 720 (2013) 347.
  * (15) R. Baier, P. Romatschke, D. T. Son, A. O. Starinets and M. A. Stephanov, JHEP 0804 (2008) 100.
  * (16) S. Bhattacharyya, V. E. Hubeny, S. Minwalla and M. Rangamani, JHEP 0802 (2008) 045.
  * (17) M. Natsuume and T. Okamura, Phys. Rev. D 77 (2008) 066014 Erratum: [Phys. Rev. D 78 (2008) 089902].
  * (18) V. E. Hubeny, S. Minwalla and M. Rangamani, arXiv:1107.5780 [hep-th].
  * (19) T. Koide, G. S. Denicol, P. Mota and T. Kodama, Phys. Rev. C 75 (2007) 034909.
  * (20) M. Fukuma and Y. Sakatani, Phys. Rev. E 84 (2011) 026316.
  * (21) T. Koide and T. Kodama, Phys. Rev. E 78 (2008) 051107.
  * (22) Y. Minami and Y. Hidaka, Phys. Rev. E 87 (2013) 023007.
  * (23) T. Hayata, Y. Hidaka, M. Hongo and T. Noumi, Phys. Rev. D 92 (2015) 065008.
  * (24) N. Banerjee, J. Bhattacharya, S. Bhattacharyya, S. Jain, S. Minwalla and T. Sharma, JHEP 1209 (2012) 046.
  * (25) K. Jensen, M. Kaminski, P. Kovtun, R. Meyer, A. Ritz and A. Yarom, Phys. Rev. Lett. 109 (2012) 101601.
  * (26) P. Romatschke, Eur. Phys. J. C 77 (2017) 21.
  * (27) P. Romatschke, Class. Quant. Grav. 27 (2010) 025006.
  * (28) G. D. Moore and K. A. Sohrabi, Phys. Rev. Lett. 106 (2011) 122302.
  * (29) G. D. Moore and K. A. Sohrabi, JHEP 1211 (2012) 148.
  * (30) F. Becattini and E. Grossi, Phys. Rev. D 92 (2015) 045037.
  * (31) R. Panerai, Phys. Rev. D 93 (2016) no.10, 104021.
  * (32) L. Landau, L. Lifshitz, Statistical Physics, Pergamon Press (1980).
  * (33) A. Vilenkin, Phys. Rev. D 21 (1980) 2260.
  * (34) L. Adamczyk et al. [STAR Collaboration], arXiv:1701.06657 [nucl-ex].
  * (35) S. Bhattacharyya, JHEP 1207 (2012) 104.
  * (36) D. N. Zubarev, A. V. Prozorkevich, S. A. Smolyanskii, Theoret. and Math. Phys. 40 (1979) 821.
  * (37) Ch. G. Van Weert, Ann. Phys. 140 (1982) 133.
  * (38) H. A. Weldon, Phys. Rev. D 26 (1982) 1394.
  * (39) F. Becattini, Phys. Rev. Lett. 108 (2012) 244502.
  * (40) F. Becattini, L. Bucciantini, E. Grossi and L. Tinti, Eur. Phys. J. C 75 (2015) no.5, 191
  * (41) M. Hongo, arXiv:1611.07074 [hep-th].
  * (42) P. Ván and T. S. Biró, Proceedings of the 12th Joint European Thermodynamics Conference, Cartolibreria SNOOPY, 2013, ed. M. Pilotelli and G. P. Beretta, p546–551
  * (43) V. E. Ambrus and E. Winstanley, Phys. Rev. D 93 (2016) 104014.
  * (44) M. N. Chernodub and S. Gongyo, JHEP 1701 (2017) 136.
  * (45) C. G. Callan, Jr., S. R. Coleman and R. Jackiw, Ann. Phys. 59 (1970) 42.
  * (46) J. I. Kapusta and C. Gale, ``Finite-temperature field theory: Principles and applications,'' Cambridge University Press (2006).
  * (47) M. Laine and A. Vuorinen, Lect. Notes Phys. 925 (2016) pp.1.
  * (48) N. P. Landsman and C. G. van Weert, Phys. Rept. 145 (1987) 141.
  * (49) F. Becattini and F. Piccinini, Annals Phys. 323 (2008) 2452.
  * (50) D. E. Kharzeev, J. Liao, S. A. Voloshin and G. Wang, Prog. Part. Nucl. Phys. 88 (2016) 1.
  * (51) T. Kalaydzhyan, Phys. Rev. D 89 (2014) 105012.
  * (52) A. Vilenkin, Phys. Rev. D 20 (1979) 1807.
  * (53) D. T. Son and A. R. Zhitnitsky, Phys. Rev. D 70 (2004) 074018.
  * (54) D. T. Son and P. Surowka, Phys. Rev. Lett. 103 (2009) 191601.
  * (55) A. V. Sadofyev and M. V. Isachenkov, Phys. Lett. B 697 (2011) 404
  * (56) M. Torabian and H. Yee, JHEP 08 (2009) 020.
  * (57) J. Erdmenger, M. Haack, M. Kaminski and A. Yarom, JHEP 0901 (2009) 055.
  * (58) N. Banerjee, J. Bhattacharya, S. Bhattacharyya, S. Dutta, R. Loganayagam and P. Surowka, JHEP 1101 (2011) 094.
  * (59) I. Amado, K. Landsteiner and F. Pena-Benitez, JHEP 1105 (2011) 081.
  * (60) Y. Neiman and Y. Oz, JHEP 1103 (2011) 023.
  * (61) K. Landsteiner, E. Megias, L. Melgar and F. Pena-Benitez, JHEP 1109 (2011) 121.
  * (62) J. H. Gao, Z. T. Liang, S. Pu, Q. Wang and X. N. Wang, Phys. Rev. Lett. 109 (2012) 232301.
  * (63) K. Jensen, R. Loganayagam and A. Yarom, JHEP 1302 (2013) 088.
  * (64) K. Landsteiner, E. Megias and F. Pena-Benitez, Phys. Rev. Lett. 107 (2011) 021601.
  * (65) A. Avkhadiev and A. V. Sadofyev, arXiv:1702.07340 [hep-th].
  * (66) S. Golkar and S. Sethi, JHEP 1605 (2016) 105
  * (67) S. Golkar and D. T. Son, JHEP 1502 (2015) 169.
  * (68) K. Jensen, P. Kovtun and A. Ritz, JHEP 1310 (2013) 186.
  * (69) V. Braguta, M. N. Chernodub, V. A. Goy, K. Landsteiner, A. V. Molochkov and M. I. Polikarpov, Phys. Rev. D 89 (2014) 074510.
  * (70) F. Becattini, V. Chandra, L. Del Zanna and E. Grossi,
  * (71) A. Flachi and K. Fukushima, arXiv:1702.04753 [hep-th].
  * (72) M. Abramowitz and I. A. Stegun, _Handbook of Mathematical Functions_ , Applied Mathematics Series (Dover Publications, 1964).
  * (73) E. Megias and M. Valle, JHEP 1411 (2014) 005.
  * (74) J. Bhattacharya, S. Bhattacharyya, S. Minwalla and A. Yarom, JHEP 1405 (2014) 147.