# Coupling of the triple-Q state to the atomic lattice by anisotropic
symmetric exchange

Felix Nickel<EMAIL_ADDRESS>Institut für Theoretische Physik und
Astrophysik, Christian-Albrechts-Universität zu Kiel, D-24098 Kiel, Germany
André Kubetzka  Department of Physics, University of Hamburg, Jungiusstraße
11, 20355 Hamburg, Germany  Soumyajyoti Haldar  Institut für Theoretische
Physik und Astrophysik, Christian-Albrechts-Universität zu Kiel, D-24098 Kiel,
Germany  Roland Wiesendanger  Department of Physics, University of Hamburg,
Jungiusstraße 11, 20355 Hamburg, Germany  Stefan Heinze  Institut für
Theoretische Physik und Astrophysik, Christian-Albrechts-Universität zu Kiel,
D-24098 Kiel, Germany  Kiel Nano, Surface, and Interface Science (KiNSIS),
University of Kiel, Germany  Kirsten von Bergmann
<EMAIL_ADDRESS>Department of Physics, University
of Hamburg, Jungiusstraße 11, 20355 Hamburg, Germany

(August 28, 2024)

###### Abstract

We identify the triple-Q (3Q) state as magnetic ground state in Pd/Mn and
Rh/Mn bilayers on Re(0001) using spin-polarized scanning tunneling microscopy
and density functional theory. An atomistic model reveals that in general the
3Q state with tetrahedral magnetic order and zero net spin moment is coupled
to a hexagonal atomic lattice in a highly symmetric orientation via the
anisotropic symmetric exchange interaction, whereas other spin-orbit coupling
terms cancel due to symmetry. Our experiments are in agreement with the
predicted orientation of the 3Q state. A distortion from the ideal tetrahedral
angles leads to other orientations of the 3Q state which, however, results in
a reduced topological orbital magnetization compared to the ideal 3Q state.

Figure 1: Illustration of the 3Q state in different highly symmetric
orientations with respect to the hexagonal atomic lattice; the respective
direction of the topological orbital moment is indicated. Note that the 3Q
states in (a) and (b) are connected by an inversion of all spins, whereas the
states in (b-d) can be converted into each other by a continuous rotation of
the spin structure.

The triple-Q (3Q) state is a three-dimensional spin structure on a two-
dimensional hexagonal lattice. It can be understood as a superposition of
three symmetry equivalent spin spiral (1Q) states resulting in a non-coplanar
magnetic state with tetrahedral angles between all adjacent magnetic moments
and four atoms in the magnetic unit cell (Fig. 1). This fascinating magnetic
state was predicted more than 20 years ago [1, 2], but the first experimental
observation was reported only recently for an hcp-stacked Mn monolayer on
Re(0001) using spin-polarized scanning tunneling microscopy (SP-STM) [3]. The
ideal 3Q state does not exhibit a net spin moment, however, theoretical
investigations have shown a significant topological orbital magnetization and
a spontaneous topological Hall effect, without the necessity of spin-orbit
coupling [4, 5, 6, 7]. Recently, neutron scattering experiments on the layered
material Co1/3TaS2 were interpreted as validation of the 3Q state, and indeed,
transport measurements of this bulk system are in agreement with a magnetic-
field-induced switching of the orbital moment direction and reveal the
topological Hall effect [8, 9]. For magnetic 3Q states in a two-dimensional
system the orbital moment is always perpendicular to the layer and the two
different directions are related to the sign of the scalar spin chirality in
an _all-in_ versus _all-out_ configuration (see Fig. 1). Another emergent
phenomenon of the 3Q state is highlighted by a theoretical proposal that
demonstrates a topological superconducting phase induced by this non-coplanar
magnetic state when it is adjacent to a conventional superconductor [10].

The 3Q state can arise in frustrated antiferromagnets, for instance as a
ground state in a hexagonal lattice of spins with antiferromagnetic nearest
and next-nearest neighbor interactions ($1<J_{1}/J_{2}<8$) [2, 3]. To lower
the energy with respect to the otherwise degenerate 1Q state higher-order
interactions (HOI) are necessary. One prominent HOI term is the biquadratic
interaction, which is one of the four-spin interactions which arise in fourth
order in a perturbative expansion of the Hubbard model [11, 12, 13]. The
aforementioned topological orbital moments that can occur in these non-
coplanar magnetic states can interact with the emergent magnetic field leading
to the topological chiral magnetic interaction that constitutes a sixth order
term [7]. Due to this topological chiral magnetic interaction large
distortions from the perfect tetrahedron angle can occur as proposed for
Mn/Re(0001) [14].

The possibility of a coupling of the ideal 3Q state to the crystal lattice is
an intriguing question which has not been addressed so far. Figure 1 displays
several 3Q states which have different orientations of the spins with respect
to the magnetic lattice plane: a 3Q state with a given orbital moment can
occur in three symmetric spin orientations with respect to the plane, denoted
as 3Q1, 3Q2, and 3Q3 [3]. An energy variation between these states, i.e. a
preferred coupling to the lattice, must originate from spin-orbit interaction.
The most prominent spin-orbit interaction term is the magneto-crystalline
anisotropy energy (MAE). However, for the ideal 3Q state with its 4 spins
arranged in a tetrahedral fashion the leading MAE term vanishes. In non-
collinear magnetic states often the spin-orbit coupling induced
Dzyaloshinskii-Moriya interaction (DMI) contributes, but in the 3Q state for
symmetry reasons also this term cancels (because the two spins on opposite
sides of any reference spin have the same magnetization direction). A further
interaction that originates from spin-orbit coupling is the anisotropic
symmetric exchange interaction (ASE), also referred to as pseudo-dipolar
interaction or compass anisotropy [15, 16]. It has been shown to play a role
for the coupling of the spins to the lattice in a uniaxial antiferromagnetic
state [3], and it can induce skyrmion lattices in two-dimensional van der
Waals magnets [17] and in centroysymmetric crystals lacking DMI [18]. Recently
another spin-orbit coupling induced term has been proposed, namely the spin-
chiral interaction (SCI) which results from the coupling of a spin to the
topological orbital moment arising in non-coplanar states [19, 20, 21, 7].

Here, we use SP-STM to investigate the magnetic ground states of Pd/Mn and
Rh/Mn bilayers on Re(0001) in real space. We find that both systems exhibit a
hexagonal magnetic superstructure, indicative of the 3Q state with tetrahedral
arrangement of the magnetic moments. Based on density functional theory (DFT)
we confirm the 3Q magnetic ground state of these two systems and calculate the
resulting topological orbital moments. We show that the ASE couples the state
to the atomic lattice resulting for both systems in the 3Q1 configuration, in
which one of the four spins of the tetrahedron points perpendicular to the
surface, Fig. 1(a,b). In addition, the ASE also determines the direction of
the in-plane magnetization components within this 3Q1 configuration, selecting
a group of very specific states out of the large number of possible 3Q
configurations.

Figure 2: (a) Perspective view of a sample of about 0.3 atomic layers of Pd on
about 0.8 atomic layers of Mn on Re(0001); the dotted white line marks a
buried Re(0001) step edge. (b) Spin-resolved constant-current STM image of the
same area as shown in (a), with the height contrast adjusted separately for
the two Pd/Mn bilayer areas to $\delta z=\pm 10$ pm ($U=+10$ mV, $I=5$ nA,
$T=4.2$ K, Cr-bulk tip). (c) Selected area of (b) imaged with a different tip
magnetization direction, exhibiting a different pattern. (d) SP-STM images of
the 3Q1 ground state of Pdfcc/Mnfcc/Re(0001) calculated from DFT based on the
spin-polarized model of STM [22]. The images were calculated for a bias
voltage of $-100\,\text{meV}$ and at a distance of $6$ Å above the surface.
The tip magnetization direction is indicated by the orange arrow. The green
and red spheres with arrows denote the Pd and Mn atoms, respectively, and
their magnetic moment directions. The two images represent the same 3Q rotated
by $180^{\circ}$, as realized in two Pd/Mn bilayers on adjacent terraces. (e)
Map of differential tunnel conductance acquired simultaneously to (b).

The starting point of our SP-STM investigation is the Mn monolayer on
Re(0001), which preferentially grows in fcc stacking. The experimentally found
magnetic ground state is the row-wise antiferromagnetic state that is the
constituing 1Q state of the tetrahedral 3Q state [3], however, DFT predicts
that the 1Q and the 3Q states are nearly degenerate. To drive the system into
a 3Q ground state we try to tune the interactions by covering the extended
fcc-Mn monolayer by non-magnetic overlayers [23, 24], namely Pd and Rh. When
submonolayer amounts of Pd or Rh are deposited we observe pseudomorphic growth
of the overlayer on top of the fcc-stacked Mn monolayer and a decoration at
the Mn step edges, see perspective view of a Pd/Mn/Re(0001) sample in Fig.
2(a); note that due to the hcp crystal structure of Re(0001) the relative
layer positions of the Pd and Mn layers are interchanged on adjacent Re
terraces, cf. the two sketches of the bilayer structures, which are rotated by
$180^{\circ}$ with respect to each other. Analysis of the step edge directions
of the overlayer islands reveals that Pd grows in only one stacking on the
fcc-Mn monolayer, with small patches of Pd islands on top of the Pd/Mn
bilayer. In contrast, the Rh monolayer occurs in both possible stackings on
top of the fcc-Mn/Re(0001) (see Supplemental Material [25]).

In Fig. 2(b) the height contrast for the two Pd/Mn bilayer areas is adjusted
separately to the same value for better visibility. We find that the Pd/Mn
exhibits an atomic-scale pattern of magnetic origin, which slightly varies
over the islands but dominantly appears as vertical stripes on the left
island, and like a hexagonal pattern on the right island, see enlarged insets.
The length scale of the magnetic superstructure indicates that the 3Q state is
the magnetic ground state. Depending on the tip magnetization direction
relative to the 3Q state, different superstructure patterns are expected;
indeed, in a measurement with a different tip magnetization direction the left
Pd/Mn island can also exhibit a hexagonal pattern, see Fig. 2(c). Also the
Rh/Mn bilayers show the same magnetic superstructure unit cell for both Rh
stackings [25] albeit with very small contrast amplitude. Nevertheless, this
is an indication that also the Rh/Mn bilayer exhibits the 3Q state. With the
help of SP-STM simulations we can find a combination of 3Q state and tip
magnetization that reproduces the magnetic patterns observed on the two Pd/Mn
areas in the measurement of Fig. 2(b). We use the constraint that the
configuration of a bilayer, and thus also a given 3Q state, must be rotated by
$180^{\circ}$ on adjacent Re terraces and find a canted tip magnetization that
nicely reproduces the experimentally observed magnetic contrast, see the DFT-
based SP-STM simulations displayed in (d). Note that an inversion of all
sample spins (or the tip magnetization direction), only slightly changes the
observed patterns.

The d$I$/d$U$ map in Fig. 2(e) is acquired simultaneously to (b). It shows a
signal variation that is spatially correlated to the slight changes of the
magnetic pattern on the Pd/Mn islands, suggesting a contribution of an
electronic contrast mechanism due to local changes in the magnetic state [26,
27]. We interpret the darker areas as magnetic domains, and the brighter areas
as transition regions between them. Indeed, close analysis shows that there is
a phase shift between different areas of the stripe contrast in the left
island [25]. In contrast to the sharp domain walls previously observed for the
3Q state in hcp-stacked Mn on Re(0001) [3], the transitions between different
3Q areas in Pd/Mn are much more extended. It is worth noting that, despite the
presence of several independent domains, no rotational domains are observed.
For both the 3Q2 and 3Q3 state rotational domains are expected to coexist, and
with an arbitrary and unknown tip magnetization direction, as in our case,
they could be easily discriminated. This suggests that we have the 3Q1
configuration in our Pd/Mn bilayer.

To understand the experimental observations we have performed DFT calculations
for Pd/Mn and Rh/Mn bilayers on Re(0001) based on the full-potential
linearized augmented plane-wave method as implemented in the FLEUR code [28,
29, 30] and using the projector augmented wave method as implemented in the
VASP code [31, 32] (see Supplemental Material for computational details). In
order to scan a large part of the phase space of magnetic states we have first
calculated the total energy of spin spiral (1Q) states which are the
fundamental solution of the classical Heisenberg model and allow to calculate
the isotropic pair-wise exchange constants 111Note, that the DMI does not play
a role for the magnetic ground state in these films. The contribution of spin-
orbit coupling (SOC) to the energy dispersion of spin spirals, which allows
the determination of the DMI parameters, is displayed in the Supplemental
Material. [25]. These calculations show that out of all 1Q states the row-wise
antiferromagnetic (RW-AFM) state (Fig. 3(a)) and the Néel state are the lowest
1Q state for Pd/Mn/Re(0001) and Rh/Mn/Re(0001), respectively [25]. However, if
also superposition states that can be stabilized by HOI are considered our DFT
calculations show that the 3Q state (Fig. 1) is by about 20 meV/Mn atom lower
than the lowest 1Q state for both systems. This is in agreement with the
experimentally observed $p(2\times 2)$ magnetic unit cell and we conclude that
the 3Q state is the ground state for both the Pd/Mn and the Rh/Mn bilayer on
Re(0001).

Figure 3: (a) Sketch of the 2Q, the 3Q, a distorted 3Q and the 1Q (RW-AFM)
state. (b) Energy dispersion for spin states along the path 2Q-3Q-1Q according
to Eq. (1) for Pd/Mn/Re(0001), Rh/Mn/Re(0001) and hcp-stacked Mn/Re(0001). The
symbols represent DFT data, the lines show fits to the DFT values using HOI
contributions of 4th and 6th order. The insets show the contributions of the
4th and 6th order separately. (c) Absolute value of the topological orbital
moment of the Mn atoms in the films. Symbols show DFT values and lines a fit
to the scalar spin chirality. Data for Mn/Re(0001) are taken from Ref. [14].

The 1Q state and its corresponding 3Q superposition state are degenerate in
energy within the Heisenberg model [14], and the large total energy difference
obtained by DFT must stem from HOIs. In order to evaluate the role of the HOIs
for these systems we calculate the energy of different superposition states
along the path continuously connecting the 1Q and the 2Q state via the 3Q
state, see Fig. 3 [14].

The configuration of the superposition state is given by

$\mathbf{s}_{i}(\alpha)=\mathbf{s}_{i}^{\rm
2Q}\cos(\alpha)+\mathbf{s}_{i}^{\rm 1Q}\sin(\alpha),$ (1)

where $\mathbf{s}_{i}(\alpha)$ represents a spin at lattice site $i$ in the
state characterized by $\alpha$, and $\mathbf{s}_{i}^{\rm 2Q}$ and
$\mathbf{s}_{i}^{\rm 1Q}$ are the spin orientations for the 2Q and 1Q state,
respectively (Fig. 3(a), see also Ref. [25]). The angle $\alpha$ is varied
between $0^{\circ}$ (2Q state) and $90^{\circ}$ (1Q state). The 3Q state
occurs at an angle of $\arcsin(1/\sqrt{3})\approx 35.26^{\circ}$.

The energy difference between the 3Q and the 1Q state (RW-AFM), $\Delta E_{\rm
3Q-1Q}$, is on the order of 20 and 40 meV/Mn atom for the Pd/Mn and the Rh/Mn
bilayer, respectively. The DFT total energy differences along the path can be
fitted by HOI contributions of 4th and 6th order (Fig. 3(b) and insets). We
find that while $\Delta E_{\rm 3Q-1Q}$ for both systems is dominated by the
4th order interaction, also the 6th order interaction is significant with a
value of about $15\%$ compared to the 4th order contribution. As a reference,
we also show the energies for the hcp-Mn/Re(0001), which was found to exhibit
the 3Q state (note that in the bilayers studied here the Mn is in fcc-
stacking). We find that the energy variation along the 2Q-3Q-1Q path is
significantly smaller; here the contributions of 4th and 6th order terms are
of similar size, giving rise to a minimum at an angle of about $55^{\circ}$,
which can be interpreted as a distorted 3Q state [14]. One can also calculate
the HOI constants explicitly [13], see Supplemental Material for the values.

Figure 4: Energy contributions of (a) the magnetocrystalline anisotropy energy
(MAE) and (b) the anisotropic symmetric exchange (ASE) for spin states along
the path 2Q-3Q-1Q according to Eq. (1) with interaction strengths obtained via
DFT for Pd/Mn/Re(0001). The line indicates the minimal (green) and maximal
(red) energy contributions out of all possible global spin rotations. For
selected points along the path the spin configurations that minimize or
maximize the energy are displayed. (c) Sketches of the two 3Q states with
minimal and maximal ASE energy, the bonds are colorized according to their
individual contributions to the total ASE. In the 3Q${}^{1}_{\textrm{M}}$
state all in-plane components are along the $\overline{\Gamma\textrm{K}}$
direction, whereas in the 3Q${}^{1}_{\textrm{K}}$ state they are along
$\overline{\Gamma\textrm{K}}$, i.e. the close packed atomic rows.

The 6th order terms have the largest effective value for Rh/Mn/Re(0001). This
is accompanied by a large topological orbital moment (Fig. 3(c)), which
maximizes for the 3Q state, and decreases if the 3Q state is distorted towards
the 2Q or the 1Q state. The direction of the orbital moment is linked to the
specific spin configuration and inverts between the all-in and the all-out
states; however, the continuous rotation of all spin directions does not alter
the size or the sign of the orbital moment (cf. Fig. 1). The presence of
topological orbital moments, which occur due to the scalar spin chirality
$\chi_{ijk}=\mathbf{s}_{i}\cdot(\mathbf{s}_{j}\times\mathbf{s}_{k})$ even in
the absence of spin-orbit coupling (SOC), is an indicator for topological-
chiral interactions [7, 14].

We have shown that the ideal undistorted 3Q state is the ground state of Pd/Mn
and Rh/Mn bilayers on Re(0001). Now we want to analyze the impact of SOC onto
the specific configuration of the 3Q state. In Fig. 4 we show the energy due
to MAE and ASE along the 2Q-3Q-1Q path for the Pd/Mn bilayer 222Note that the
spin-chiral interaction which also arises due to SOC [7] does not play a role
for the coupling of the 3Q state to the atomic lattice as can be shown
analytically [25].. For each state along the path, we allowed global rotations
of the complete spin structure. Then we performed a global optimization to
find the spin orientations of lowest (green) and highest (red) energy. From
our DFT calculations we found that the MAE prefers an in-plane orientation of
the magnetic moments with a strength of $-1.70\,\text{meV/Mn atom}$ compared
to an out-of-plane alignment. The 2Q state and the 1Q state have
configurations with all spins in the easy plane and accordingly there is a
selection of the specific configuration by the MAE, with an energy difference
of 0.7 and 1.7 meV/Mn atom between different configurations for the 2Q and 1Q,
respectively. In contrast, an undistorted 3Q state with tetrahedral angles
between all nearest neighbor moments cannot gain MAE (Fig. 4(a)), as obvious
from symmetry considerations. Therefore, the 3Q always has the higher MAE
compared to the 2Q or the 1Q state, for both easy plane and easy axis systems
(by $0.7$ and $0.2-1$ meV/Mn atom, respectively). However, for Pd/Mn the large
HOI contribution favoring the 3Q state (see inset) more than overcompensates
the MAE between the 3Q and the 2Q/1Q states. In contrast, when a distorted 3Q
state is favored by the 6th order HOIs, as is the case for hcp-Mn/Re(0001)
with $\alpha\approx 55^{\circ}$, the MAE can select between different
orientations of spin configurations.

Figure 4(b) shows the same graph for the ASE energy, which can be calculated
for a given magnetic state via

$E_{\rm ASE}=-\sum_{i,j}J_{\rm
ASE}(\mathbf{s}_{i}\cdot\mathbf{d}_{ij})(\mathbf{s}_{j}\cdot\mathbf{d}_{ij})\quad.$
(2)

Here $\mathbf{d}_{ij}$ is the normalized connection vector between two lattice
sites defined by $i,j$ and $\mathbf{s}_{i},\mathbf{s}_{j}$ are the
corresponding normalized spin moments. This interaction favours either a
ferromagnetic or an antiferromagnetic alignment along the connection
$\mathbf{d}_{ij}$ depending on the sign of $J_{\rm ASE}$. If one or both
magnetic moments are perpendicular to $\mathbf{d}_{ij}$ the contribution
vanishes. For Pd/Mn the value of the ASE calculated by DFT is $J_{\rm
ASE}=-0.30\,\text{meV/Mn atom}$.

In contrast to the MAE, the ASE favors a specific orientation of the spin
structure with respect to the lattice for every $\alpha$, including the 3Q
state (Fig. 4(b)). For a given sign of the ASE, the energy difference between
the different superposition states, i.e. a variation of $\alpha$, is small
with only up to 0.2 meV/Mn atom. In contrast, the energy difference between
different spin configurations for the same $\alpha$ are much larger, on the
order of up to 1 meV/Mn atom. We have also studied the effect of magnetic
dipole-dipole interactions on the different states, however, the energy scale
is much lower with a maximum energy difference of 0.25 meV/Mn atom [25]. We
find that for a 3Q ground state the ASE always favors the 3Q1 orientation,
regardless of the sign of the ASE. The sign of the ASE selects the specific
orientation of the spins within the 3Q1 state.

The two configurations corresponding to lowest energy states for the different
signs of the ASE are displayed in Fig. 4(c) and differ only by a rotation
around the surface normal. For a negative sign of the ASE the
3Q${}^{1}_{\textrm{M}}$ state is favored and for a positive sign the
3Q${}^{1}_{\textrm{K}}$ state. For each of these 3Q1 states the energy is
invariant under the inversion of all spins (e.g. from all-out to all-in) and
under the inversion of the $z$-components. This leads to four degenerate
states [25], which is confirmed by self-consistent DFT calculations including
SOC for the four 3Q${}^{1}_{\textrm{M}}$ states and one
3Q${}^{1}_{\textrm{K}}$ state for Pd/Mn/Re(0001). The obtained total energy
difference between the 3Q${}^{1}_{\textrm{M}}$ states is negligible and the
energy difference $\Delta E$ between the 3Q${}^{1}_{\textrm{M}}$ and
3Q${}^{1}_{\textrm{K}}$ state amounts to $0.75\,\text{meV}$ per Mn atom [25],
in very good agreement with the results obtained by the atomistic model (Fig.
4(b,c)) with the ASE constant determined from DFT.

In conclusion, we have inspected the different energy terms that govern the
detailed spin configuration of the 3Q state. We find that in general the ASE
determines both the orientation and the relative spin alignments of the ideal
3Q state. In the case of Pd/Mn/Re(0001), it is the 3Q${}^{1}_{\textrm{M}}$
configuration which exhibits three-fold symmetry and therefore does not show
rotational domains, in agreement with the experimental data. Distorted 3Q
states are always of 3Q2 or 3Q3 type, depending on the sign of MAE, and can be
identified in experiments by the presence of rotational domains according to
their two-fold symmetry on a triangular lattice. In general, for distorted
states both the topological orbital moment and the topological Hall effect
will be reduced compared to the ideal 3Q state. Our work provides key
microscopic insight into the detailed spin configuration of the 3Q state. The
knowledge of the magnetic moment directions is important for an understanding
of the response of domains and domain walls in such topological orbital
ferromagnets to external magnetic fields, which in turn is crucial for their
transport properties.

###### Acknowledgements.

We thank Moritz Goerzen for valuable discussions. K.v.B. and S.H. gratefully
acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) via projects no. 402843438, no. 418425860, and no.
462602351, and computing time provided by the North-German Supercomputing
Alliance (HLRN).

## References

  * Momoi  _et al._ [1997] T. Momoi, K. Kubo, and K. Niki, Possible chiral phase transition in two-dimensional solid ${}^{3}\mathrm{He}$, Phys. Rev. Lett. 79, 2081 (1997).
  * Kurz  _et al._ [2001] P. Kurz, G. Bihlmayer, K. Hirai, and S. Blügel, Three-Dimensional Spin Structure on a Two-Dimensional Lattice: Mn $/$Cu(111), Phys. Rev. Lett. 86, 1106 (2001).
  * Spethmann  _et al._ [2020] J. Spethmann, S. Meyer, K. von Bergmann, R. Wiesendanger, S. Heinze, and A. Kubetzka, Discovery of magnetic single- and triple-$\mathbf{q}$ states in $\mathrm{Mn}/\mathrm{Re}(0001)$, Phys. Rev. Lett. 124, 227203 (2020).
  * Martin and Batista [2008] I. Martin and C. D. Batista, Itinerant electron-driven chiral magnetic ordering and spontaneous quantum hall effect in triangular lattice models, Phys. Rev. Lett. 101, 156402 (2008).
  * Hoffmann  _et al._ [2015] M. Hoffmann, J. Weischenberg, B. Dupé, F. Freimuth, P. Ferriani, Y. Mokrousov, and S. Heinze, Topological orbital magnetization and emergent hall effect of an atomic-scale spin lattice at a surface, Phys. Rev. B 92, 020401(R) (2015).
  * Hanke  _et al._ [2016] J.-P. Hanke, F. Freimuth, A. K. Nandy, H. Zhang, S. Blügel, and Y. Mokrousov, Role of Berry phase theory for describing orbital magnetism: From magnetic heterostructures to topological orbital ferromagnets, Phys. Rev. B 94, 121114(R) (2016).
  * Grytsiuk  _et al._ [2020] S. Grytsiuk, J.-P. Hanke, M. Hoffmann, J. Bouaziz, O. Gomonay, G. Bihlmayer, S. Lounis, Y. Mokrousov, and S. Blügel, Topological–chiral magnetic interactions driven by emergent orbital magnetism, Nat. Commun. 11, 511 (2020).
  * Park  _et al._ [2023] P. Park, W. Cho, C. Kim, Y. An, Y.-G. Kang, M. Avdeev, R. Sibille, K. Iida, R. Kajimoto, K. H. Lee, W. Ju, E.-J. Cho, H.-J. Noh, M. J. Han, S.-S. Zhang, C. D. Batista, and J.-G. Park, Tetrahedral triple-q ordering in the metallic triangular lattice antiferromagnet Co1/3TaS2, arXiv:2303.03760 (2023).
  * Takagi  _et al._ [2023] H. Takagi, R. Takagia, S. Minami, T. Nomoto, K. Ohishi, M.-T. Suzuki, Y. Yanagi, M. Hirayama, N. D. Khanh, K. Karube, H. Saito, D. Hashizume, R. Kiyanagi, Y. Tokura, R. Arita, T. Nakajima, and S. Seki, Spontaneous topological hall effect induced by non-coplanar antiferromagnetic order in intercalated van der waals materials, Nat. Phys. (2023).
  * Bedow  _et al._ [2020] J. Bedow, E. Mascot, T. Posske, G. S. Uhrig, R. Wiesendanger, S. Rachel, and D. K. Morr, Topological superconductivity induced by a triple-q magnetic structure, Phys. Rev. B 102, 180504(R) (2020).
  * Takahashi [1977] M. Takahashi, Half-filed Hubbard model at low temperature, J. Phys. C Solid State Phys. 10, 1289 (1977).
  * MacDonald  _et al._ [1988] A. H. MacDonald, S. M. Girvin, and D. Yoshioka, $\frac{t}{U}$ expansion for the Hubbard model, Phys. Rev. B 37, 9753 (1988).
  * Hoffmann and Blügel [2020] M. Hoffmann and S. Blügel, Systematic derivation of realistic spin models for beyond-heisenberg solids, Phys. Rev. B 101, 024418 (2020).
  * Haldar  _et al._ [2021] S. Haldar, S. Meyer, A. Kubetzka, and S. Heinze, Distorted $3q$ state driven by topological-chiral magnetic interactions, Phys. Rev. B 104, L180404 (2021).
  * Smith [1976] D. Smith, New mechanisms for magnetic anisotropy in localised s-state moment materials, Journal of Magnetism and Magnetic Materials 1, 214 (1976).
  * Staunton  _et al._ [1988] J. B. Staunton, B. L. Gyorffy, J. Poulter, and P. Strange, A relativistic RKKY interaction between two magnetic impurities-the origin of a magnetic anisotropic effect, Journal of Physics C: Solid State Physics 21, 1595 (1988).
  * Amoroso  _et al._ [2020] D. Amoroso, P. Barone, and S. Picozzi, Spontaneous skyrmionic lattice from anisotropic symmetric exchange in a Ni-halide monolayer, Nat. Commun. 11, 5784 (2020).
  * Hirschberger  _et al._ [2021] M. Hirschberger, S. Hayami, and Y. Tokura, Nanometric skyrmion lattice from anisotropic exchange interactions in a centrosymmetric host, New J. Phys. 23, 023039 (2021).
  * Brinker  _et al._ [2019] S. Brinker, M. dos Santos Dias, and S. Lounis, The chiral biquadratic pair interaction, New J. Phys. 21, 083015 (2019).
  * Lászlóffy  _et al._ [2019] A. Lászlóffy, L. Rózsa, K. Palotás, L. Udvardi, and L. Szunyogh, Magnetic structure of monatomic Fe chains on Re(0001): Emergence of chiral multispin interactions, Phys. Rev. B 99, 184430 (2019).
  * Mankovsky  _et al._ [2020] S. Mankovsky, S. Polesya, and H. Ebert, Extension of the standard heisenberg hamiltonian to multispin exchange interactions, Phys. Rev. B 101, 174401 (2020).
  * Wortmann  _et al._ [2001] D. Wortmann, S. Heinze, P. Kurz, G. Bihlmayer, and S. Blügel, Resolving complex atomic-scale spin structures by spin-polarized scanning tunneling microscopy, Phys. Rev. Lett. 86, 4132 (2001).
  * Romming  _et al._ [2013] N. Romming, C. Hanneken, M. Menzel, J. E. Bickel, B. Wolter, K. von Bergmann, A. Kubetzka, and R. Wiesendanger, Writing and deleting single magnetic skyrmions, Science 341, 636 (2013).
  * Romming  _et al._ [2018] N. Romming, H. Pralow, A. Kubetzka, M. Hoffmann, S. von Malottki, S. Meyer, B. Dupé, R. Wiesendanger, K. von Bergmann, and S. Heinze, Competition of Dzyaloshinskii-Moriya and higher-order exchange interactions in Rh/Fe atomic bilayers on Ir(111), Phys. Rev. Lett. 120, 207201 (2018).
  * [25] See Supplemental Material for details on the experiments, DFT calculations and the atomistic model.
  * von Bergmann  _et al._ [2012] K. von Bergmann, M. Menzel, D. Serrate, Y. Yoshida, S. Schröder, P. Ferriani, A. Kubetzka, R. Wiesendanger, and S. Heinze, Tunneling anisotropic magnetoresistance on the atomic scale, Phys. Rev. B 86, 134422 (2012).
  * Hanneken  _et al._ [2015] C. Hanneken, F. Otte, A. Kubetzka, B. Dupé, N. Romming, K. von Bergmann, R. Wiesendanger, and S. Heinze, Electrical detection of magnetic skyrmions by tunneling non-collinear magnetoresistance, Nat. Nanotechnol. 10, 1039 (2015).
  * [28] See https://www.flapw.de.
  * Kurz  _et al._ [2004] P. Kurz, F. Förster, L. Nordström, G. Bihlmayer, and S. Blügel, Ab initio treatment of noncollinear magnets with the full-potential linearized augmented plane wave method, Phys. Rev. B 69, 024415 (2004).
  * Heide  _et al._ [2009] M. Heide, G. Bihlmayer, and S. Blügel, Describing Dzyaloshinskii-Moriya spirals from first principles, Phys. B 404, 2678 (2009).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * [32] See https://www.vasp.at.
  * Note [1] Note, that the DMI does not play a role for the magnetic ground state in these films. The contribution of spin-orbit coupling (SOC) to the energy dispersion of spin spirals, which allows the determination of the DMI parameters, is displayed in the Supplemental Material. [25].
  * Note [2] Note that the spin-chiral interaction which also arises due to SOC [7] does not play a role for the coupling of the 3Q state to the atomic lattice as can be shown analytically [25].


How the Blast-Wave Model Describes PID Hadron Spectra
from 5 TeV p-Pb Collisions

Thomas A. Trainor1,

1 University of Washington, Seattle, USA

<EMAIL_ADDRESS>
51st International Symposium on Multiparticle Dynamics (ISMD2022)
Pitlochry, Scottish Highlands, 1-5 August 2022 10.21468/SciPostPhysProc.?

## Abstract

The blast-wave (BW) spectrum model has been applied extensively to nucleus-
nucleus collision data with the intention to demonstrate formation of a quark-
gluon plasma (QGP) in more-central A-A collisions. More recently the BW model
has been applied to p-p, d-Au and p-Pb collisions. Such results are
interpreted to indicate that “collectivity” (flows) and QGP appear in smaller
systems. In this talk I review BW analysis of identified-hadron spectra from 5
TeV p-Pb collisions and examine the shape evolution of model spectra with
collision centrality. I evaluate data-model fit quality using conventional
statistical measures. I conclude that the BW model is not a valid data model.

##  1 Introduction

The blast-wave (BW) spectrum model has been extended in recent years to
identified-hadron (PID) spectrum data from small collision systems (e.g.
$p$-$p$, $p$-A and $d$-A collisions at the RHIC and LHC). It is now
conventional to interpret the existence of such BW fit results as confirming
the presence of hydrodynamic flows in small systems and to infer quark-gluon
plasma (QGP) formation as well [1]. As a result, the intended role of small
systems as control experiments relative to QGP formation in more-central
nucleus-nucleus or A-A collisions is vacated. In response, several questions
emerge [2]: Can BW model fits actually demonstrate flows? Is there an
alternative spectrum model with more likely physical interpretation that
better describes spectrum data? What role do jets (nonflow) play in nuclear
collisions? How should spectrum models be evaluated as to fit quality? And, is
QGP actually formed in small collision systems?

##  2 Blast Wave (BW) Spectrum Model

The BW model applied in Ref. [3] to 5 TeV $p$-Pb spectrum data is taken from
Ref. [4] that introduced a hydrodynamics-based formula to describe pion
spectra from fixed-target $\sqrt{s_{NN}}\approx 19$ GeV S-S collisions at the
CERN SPS. The relevant formula is Eq. (7) of Ref. [4]

$\displaystyle dn_{ch}/m_{t}dm_{t}$ $\displaystyle\propto$ $\displaystyle
m_{t}\int_{0}^{R}rdrI_{0}[p_{t}\sinh(\rho)/T]K_{1}[m_{t}\cosh(\rho)/T],~{}~{}~{}~{}$
(1)

with source boost $\rho=\tanh^{-1}(\beta_{t})$, transverse speed
$\beta_{t}(r)$ and mean transverse speed $\langle\beta_{t}\rangle$. $I_{0}$
and $K_{1}$ are modified Bessel functions. Equation (1) represents a thermal
energy spectrum (Boltzmann exponential) in a boost (comoving) frame convoluted
with a source boost ($\sim$speed) distribution on source radius to describe a
particle spectrum in the lab frame.

Figure 1 shows BW fits (solid) based on Eq. (1) compared to PID spectrum data
(points) from 5 TeV $p$-Pb collisions. The spectra are plotted as densities on
$p_{t}$ (i.e. $p_{t}$ spectra as published) vs pion transverse rapidity
$y_{t}=\ln[(p_{t}+m_{t\pi})/m_{\pi}]$ that provides improved visual access at
low $p_{t}$. Model curves are generated using fitted BW parameters from Table
5 of Ref. [3]. The arrows indicate fit intervals for each hadron species.
Substantial data-model deviations are notable.

Figure 1:  Data spectra (points) and BW model fits (curves) for four
identified (PID) hadron species from 5 TeV $p$-Pb collisions [3]. Arrows
indicate fit intervals.

##  3 Two-component Spectrum Model (TCM)

The two-component (soft+hard) model (TCM) for hadron spectra was first
introduced in Ref. [5] for 200 GeV $p$-$p$ collisions. Given a $p$-$p$
spectrum TCM for unidentified-hadron spectra with soft and hard charge
densities $\bar{\rho}_{s}$ and $\bar{\rho}_{h}$ [6], a TCM for PID hadrons can
be generated by assuming that each hadron species $i$ comprises certain
fractions of soft and hard TCM components denoted by $z_{si}$ and $z_{hi}$
(both $\leq 1$). The PID spectrum TCM can then be written as [8, 9]

$\displaystyle\bar{\rho}_{0i}(y_{t},n_{s})$ $\displaystyle\approx$
$\displaystyle d^{2}n_{chi}/m_{t}dm_{t}dy_{z}\approx
z_{si}(n_{s})\bar{\rho}_{s}\hat{S}_{0i}(y_{t})+z_{hi}(n_{s})\bar{\rho}_{h}\hat{H}_{0i}(y_{t}),$
(2)

where $\hat{S}_{0i}(y_{t})$ and $\hat{H}_{0i}(y_{t})$ are unit-normal model
functions approximately independent of $n_{ch}$, total charge density
$\bar{\rho}_{0}\equiv n_{ch}/\Delta\eta=\bar{\rho}_{s}+\bar{\rho}_{h}$, and
$n_{s}=\Delta\eta\bar{\rho}_{s}$ serves as an event-class index.

Figure 2:  First: $y_{t}$ spectra for ten $n_{ch}$ classes of 200 GeV $p$-$p$
collisions. Second: Fragmentation functions for a range of jet energies from 2
TeV $p$-$\bar{p}$ collisions. Third: Jet energy spectra for several $p$-$p$
collision energies. Fourth: Predicted (curves) and measured (points) hard
components for several $p$-$p$ collision energies.

Figure 2 illustrates definition and physical interpretation of TCM model
functions. The first panel shows $p_{t}$ spectra from 200 GeV $p$-$p$
collisions in relation to fixed soft-component model $\hat{S}_{0}(y_{t})$
(curves) that describes the data asymptotic limit for $n_{ch}\rightarrow 0$
[5]. Data hard components are complementary to $\hat{S}_{0}(y_{t})$ and are
modeled by $\hat{H}_{0}(y_{t})$. TCM hard components are quantitatively
predicted by convoluting $p$-$p$ jet fragmentation functions (second panel)
with a minimum-bias jet spectrum appropriate to the $p$-$p$ collision energy
(third panel) [7]. Predictions (curves) compared to data hard components
(points) are shown in the fourth panel.

Figure 3:  PID spectrum data from Fig. 1 (points) compared to a TCM
description (curves) [9, 10]. The TCM is not fitted to individual spectra.

Figure 3 shows PID TCM spectra (curves) compared to data spectra (the same
points appearing in Fig. 1). Derivation of the PID spectrum TCM for 5 TeV
$p$-Pb collisions is described in Refs. [8, 9, 10]. All data are described
within their published statistical uncertainties. It is important to note that
the TCM does not result from fits to individual spectra.

##  4 Spectrum Shape Evolution

Whereas the TCM provides absolute predictions for hadron yields as well as
spectrum shapes, the BW model is not expected to provide such absolute
predictions: “…the normalization of the spectrum…we will always adjust for a
best fit to the data, because we are only interested in the shape of the
spectra to reveal the dynamics of the collision zone at freeze-out [emphasis
added]” [4]. “This [assumed collective hydrodynamic flow] results in a
characteristic dependence of the [spectrum] shape which can be described with
a common kinetic freeze-out temperature parameter $T_{kin}$ and a collective
average expansion velocity $\langle\beta_{t}\rangle$ [citing Ref. [4]]” [3].
Is a BW model shape physically interpretable? Given that limitation one may
elect to invoke model comparisons based on neutral shape measures so as to
provide unbiased comparisons. “Neutral measure” here means a statistical
measure motivated by effective statistical analysis via standard practice
rather than by a sought-after result.

One possibility is logarithmic derivatives to determine local spectrum
curvature. A logarithmic derivative is $(1/f)df/dx=d\ln f/dx$. The second
derivative of the logarithm of spectrum $\bar{\rho}_{0}(y_{t})$,
$-d^{2}\ln[\bar{\rho}_{0}(y_{t})]/dy_{t}^{2}$, approximates local curvature of
the spectrum plotted in a semilog format. Since Fig. 3 and Refs. [9, 10]
establish that the TCM is statistically equivalent to PID spectrum data,
elements of the TCM can be used to explain curvature trends. For the TCM hard
component, curvature is approximately a fixed value $1/\sigma_{y_{t}}^{2}$ (in
terms of the $\hat{H}_{0}(y_{t})$ Gaussian width) near its mode but falls to
zero for the power-law tail at higher $y_{t}$. For the TCM soft component,
curvature is approximated by $\propto\cosh(y_{t})$ at lower $y_{t}$ since
spectra (except for pions) are well-approximated there by a Boltzmann
exponential on $m_{t}=m_{i}\cosh(y_{t})$ [2].

Figure 4:  First: Local curvature trends for neutral kaons and for seven event
classes of 5 TeV $p$-Pb collisions represented by the TCM. Second: Local
curvatures for BW model fits to data demonstrating large deviations. Third: BW
model curves corresponding to fixed zero radial flow
$\langle\beta_{t}\rangle=0$ that closely approximate a Boltzmann distribution
(bold dotted curve). Fourth: BW model curves corresponding to fixed
temperature $T=145$ MeV compared to the soft component $\hat{S}_{0}(y_{t})$ of
the $p$-Pb TCM (bold dashed).

Figure 4 (first) shows local curvatures vs $y_{t}$ for seven event classes of
5 TeV $p$-Pb collisions. As noted, curvature goes as $\propto\cosh(y_{t})$ at
lower $y_{t}$, rises toward a saturation value (hatched band) near the hard-
component mode and then falls to zero for the power-law tail. Detailed
centrality dependence corresponds to varying relative amplitudes of hard and
soft components.

Figure 4 (second) shows the same procedure applied to BW model functions
(curves of several line styles) as in Fig. 1. Included for reference is the
TCM curve at left corresponding to the most-central event class (bold solid).
There is a dramatic difference between TCM (statistically equivalent to data
as noted) and BW model above $y_{t}$ = 2.7 ($p_{t}\approx 1$ GeV/c).
Ironically, the BW model has no sensitivity to detailed variation of data
spectrum shapes.

Figure 4 (third) shows BW model functions corresponding to fit parameters from
Ref. [3] except $\langle\beta_{t}\rangle\rightarrow 0$ (zero radial flow) in
which case the curves approximate a Boltzmann exponential (bold dotted) which
is consistent with the basic assumptions of Ref. [4]. The fourth panel shows
BW model functions corresponding to fit parameters from Ref. [3] except
temperature $T_{kin}$ is held fixed at 145 MeV. Comparison of third and fourth
panels relative to fixed model $\hat{S}_{0}(y_{t})$ demonstrates the main
effect of nonzero $\langle\beta_{t}\rangle$ within the BW model. Note that the
most-peripheral event class in the fourth panel (the lowest solid curve), with
$\langle\beta_{t}\rangle\approx 0.25$, is consistent with TCM soft-component
model $\hat{S}_{0}(y_{t})$ (bold dashed) that represents the
$n_{ch}\rightarrow 0$ limiting case (zero particle density). What mechanism
generates flows at zero particle density?

##  5 Spectrum Model Fit Quality

A standard measure of data-model fit quality is the Z-score (for an $i^{th}$
observation) [11]

$\displaystyle Z_{i}$ $\displaystyle=$
$\displaystyle\frac{O_{i}-M_{i}}{\sigma_{i}},$ (3)

where $O_{i}$ is an observation, $M_{i}$ is a model prediction and
$\sigma_{i}$ is the uncertainty for the observation. For an acceptable model
one expects $Z_{i}$ to have an r.m.s. value near 1.

Figure 5 shows BW model Z-scores for pions, neutral kaons, protons and
Lambdas, where uncertainties $\sigma_{i}$ are published statistical errors
from Ref. [3]. Since $\chi^{2}=\sum_{i}Z_{i}^{2}$ these results imply
$\chi^{2}/\text{ndf}\sim O(25-100)$. Even within imposed fit intervals
(arrows) chosen to favor the model the Z-scores are sufficiently large as to
falsify the BW model as applied to these data.

Figure 5: Z-scores for BW model fits to pion, neutral kaon, proton and Lambda
spectra as in Fig. 1. Statistical uncertainties are used. Arrows indicate fit
intervals.

Figure 6 shows TCM Z-scores for the same data with the same uncertainties. The
TCM is applied to all available data over full $y_{t}$ acceptances (no
restricted fit intervals) and is not fitted to individual spectra. Z-scores
are consistent with an acceptable model: random $O(1)$ fluctuations with the
exception of narrow excursions (for charged pions and protons) that are
statistically significant but may be local data anomalies associated with
$dE/dx$ PID [9].

Figure 6: Z-scores for TCM applied to four hadron species. The TCM is not
fitted to individual spectra. Note that kaon data are described accurately in
$p_{t}\in[0,7]$ GeV/c.

##  6 What is an Elementary Collision?

The high-energy heavy ion program initially assumed that QGP formation might
occur in more-central A-A collisions, but confirmation of that achievement
required control (null) experiments in the form of $p$-$p$, $p$-A and $d$-A
(small systems) data, i.e. hadron production from elementary collisions (e.g.
“cold nuclear matter”) where QCD is nominally well understood. Data
manifestations of QGP formation in A-A should contrast dramatically (?) with
data trends from elementary collisions. However, arguments based on certain
correlation features (ridges) and evolution of hydrodynamic theory to achieve
“good agreement with the data” assert that nominally “elementary” (small
system) collisions actually support hydrodynamic evolution manifested by
“flow-like features,” [1] which begs the question: what is an elementary
collision?

Figure 7 (first) shows a $\sqrt{s_{NN}}\approx 19$ GeV S-S pion spectrum
(dots), the data that motivated Ref. [4], compared to a Boltzmann exponential
(dash-dotted). Within the BW model context any deviations from a Boltzmann
reference curve in the lab frame must indicate a boosted particle source:
radial flow. Also plotted are data (open circles) from 17 GeV $p$-$p$
collisions and TCM soft-component model $\hat{S}_{0}(m_{t})$ (dashed, $T=145$
MeV) appropriate for 19 GeV [7]. The second panel shows an $m_{t}$ spectrum
from 91 GeV $e^{+}$-$e^{-}$ collisions with $q$-$\bar{q}$ dijet final state
[12] compared to the same $\hat{S}_{0}(m_{t})$. Does that mean there is radial
flow in $e^{+}$-$e^{-}$ collisions?

Figure 7:  First: TCM soft-component $\hat{S}_{0}(m_{t})$ compared to 19 GeV
S-S data (solid points) and Boltzmann exponential (dash-dotted). Second:
$m_{t}$ spectrum from 91 GeV $e^{+}$-$e^{-}$ collisions (points) [12] compared
to $\hat{S}_{0}(m_{t})$ with $T=145$ MeV. Third: The same $e^{+}$-$e^{-}$
spectrum but with $T=90$ MeV for $\hat{S}_{0}(m_{t})$. Fourth: $e^{+}$-$e^{-}$
longitudinal momentum spectrum on rapidity with pion mass assumed [12].

Figure 7 (third) replots the same $e^{+}$-$e^{-}$ $m_{t}$ spectrum vs pion
$y_{t}$ with $\hat{S}_{0}(m_{t})$ now based on $T=90$ MeV. $e^{+}$-$e^{-}$
data are described within uncertainties. The fourth panel presents an
$e^{+}$-$e^{-}$ dijet longitudinal momentum spectrum as a density on $y_{z}$
(pion mass assumed). The point of this comparison is that hadrons, even from
the most elementary $e^{+}$-$e^{-}$ collisions, follow a Boltzmann
distribution with power-law tail $\hat{S}_{0}(m_{t})$ on transverse momentum,
a basic characteristic of high-energy collisions and parton or nucleon
fragmentation. Deviations from a Boltzmann exponential cannot be used to claim
emission from a flowing particle source.

##  7 Conclusion

Model-independent shape measures and Z-scores (based on statistical
uncertainties) falsify the BW spectrum model from Ref. [4] as applied to 5 TeV
$p$-Pb PID spectra from Ref. [3]. Monolithic BW model parameters,
conventionally associated with radial flow, mimic TCM nonjet and jet
contributions over limited $p_{t}$ intervals. The soft components of $p$-Pb
$p_{t}$ spectra are consistent with $e^{+}$-$e^{-}$ dijet $p_{t}$ spectra. The
data hard components are predicted by measured jet properties. It is certainly
true of high-energy nuclear collisions that almost all high-$p_{t}$ hadrons
are jet fragments. But it is equally true that almost all jet fragments are
low-$p_{t}$ hadrons. Spectrum models with a single component (monolithic), and
especially with no jet description, cannot successfully model spectrum data,
especially the strong low-$p_{t}$ jet fragment contribution with peak near 1
GeV/c. BW model fits thus cannot provide evidence for flows, “collectivity” or
QGP in small collision systems. The parameter values may have no physical
significance.

## References

  * [1] J. L. Nagle and W. A. Zajc, Small System Collectivity in Relativistic Hadronic and Nuclear Collisions, Ann. Rev. Nucl. Part. Sci. 68, 211-235 (2018). 10.1146/annurev-nucl-101916-123209.
  * [2] T. A. Trainor, Evaluating the blast-wave model as a description of 5 TeV p-Pb $p_{t}$ spectra, arXiv:2206.07791.
  * [3] B. B. Abelev et al. (ALICE Collaboration), Multiplicity Dependence of Pion, Kaon, Proton and Lambda Production in p-Pb Collisions at $\sqrt{s_{NN}}$ = 5.02 TeV, Phys. Lett. B 728, 25 (2014). 10.1016/j.physletb.2013.11.020.
  * [4] E. Schnedermann, J. Sollfrank and U. W. Heinz, Thermal phenomenology of hadrons from 200-A GeV S+S collisions, Phys. Rev. C 48, 2462 (1993). 10.1103/PhysRevC.48.2462.
  * [5] J. Adams et al. (STAR Collaboration), The multiplicity dependence of inclusive $p_{t}$ spectra from $p$-$p$ collisions at $\sqrt{s}=200$ GeV, Phys. Rev. D 74, 032006 (2006). 10.1103/PhysRevD.74.032006.
  * [6] T. A. Trainor, A two-component model of hadron production applied to $p_{t}$ spectra from 5 TeV and 13 TeV $p$-$p$ collisions at the large hadron collider, arXiv:2104.08423.
  * [7] T. A. Trainor, Manifestations of minimum-bias dijets in high-energy nuclear collisions, arXiv:1701.07866.
  * [8] T. A. Trainor, A two-component model for identified-hadron $p_{t}$ spectra from 5 TeV p-Pb collisions, J. Phys. G 47, no.4, 045104 (2020). 10.1088/1361-6471/ab5831.
  * [9] T. A. Trainor, Precision identified-hadron spectrum analysis for 5 TeV $p$-Pb collisions – Part I, arXiv:2112.09790.
  * [10] T. A. Trainor, Precision identified-hadron spectrum analysis for 5 TeV $p$-Pb collisions – Part II, arXiv:2112.12330.
  * [11] E. Kreyszig (1979), Advanced Engineering Mathematics (Fourth ed.), Wiley, p. 880, eq. 5. ISBN 0-471-02140-7.
  * [12] D. Buskulic et al. (ALEPH Collaboration), Properties of hadronic Z decays and test of QCD generators, Z. Phys. C 55, 209 (1992). 10.1007/BF01482583.


$\displaystyle\tau(x)-\tau(y)$
$\displaystyle=\int_{y}^{x}\tau^{\prime}(u)\mathrm{d}u$ (287)
$\displaystyle\geq(x-y)\int_{0}^{1}(1-t)\tau^{\prime}(y)+t\tau^{\prime}(x)\,\mathrm{d}t$
(288)
$\displaystyle=\frac{x-y}{2}\left(\tau^{\prime}(x)+\tau^{\prime}(y)\right)\,.$
(289)

For the upper bound, concavity of $\tau^{\prime}$ implies the existence of an
affine linear function $h$ with $h(u)\geq\tau^{\prime}(u)$ for all
$u\in[0,\infty)$ and

$h\\!\left(\frac{x+y}{2}\right)=\tau^{\prime}\\!\left(\frac{x+y}{2}\right)\,.$
(290)

Thus,

$\displaystyle\tau(x)-\tau(y)$ $\displaystyle\leq\int_{y}^{x}h(u)\mathrm{d}u$
(291) $\displaystyle=\frac{x-y}{2}\left(h(x)+h(y)\right)$ (292)
$\displaystyle=(x-y)h\\!\left(\frac{x+y}{2}\right)\,.$ (293)

∎

###### Proof of Lemma 2.5.

First, consider the case $x\geq y$. Define
$f(x,y)=\tau(x-y)-\tau(x)-\tau(y)+2y\tau^{\prime}(x)$. We want to show
$f(x,y)\geq 0$. The derivative of $f$ with respect to $y$ is

$\partial_{y}f(x,y)=-\tau^{\prime}(x-y)-\tau^{\prime}(y)+2\tau^{\prime}(x)\,.$
(294)

As $\tau^{\prime}$ is nondecresing and $x\geq\max(y,x-y)$, we obtain
$\partial_{y}f(x,y)\geq 0$. Hence, $f(x,y)\geq f(x,0)=0$, as $\tau(0)=0$.

Now, consider the case $x\leq y$. Set
$g(x,y)=\tau(y-x)-\tau(x)-\tau(y)+2y\tau^{\prime}(x)$, which yields

$\partial_{y}g(x,y)=\tau^{\prime}(y-x)-\tau^{\prime}(y)+2\tau^{\prime}(x)\,.$
(295)

As $\tau^{\prime}$ is subsadditive (Lemma B.6), we obtain
$\partial_{y}g(x,y)\geq 0$. Thus, $g(x,y)\geq
g(x,x)=-2\tau(x)+2x\tau^{\prime}(x)$ as $\tau(0)=0$. As $\tau^{\prime}$ is
nondecreasing and $\tau(0)=0$, we have $\tau(x)\leq x\tau^{\prime}(x)$. Hence,
$g(x,y)\geq 0$. ∎

###  C.2 From Section 4

###### Proof of Proposition 4.2.

  1. (i)

The quadratic polynomial $t\mapsto t^{2}+at+b$ attains its minimum
$-\frac{a^{2}}{4}+b$ at $t_{0}=-\frac{a}{2}$. It is nonnegative if and only if
$4b\geq a^{2}$. A function of the form $t\mapsto(t-t_{0})^{2}+h^{2}$ can be
written as $t^{2}+\tilde{a}t+\tilde{b}$ with $\tilde{a}=-2t_{0}$ and
$\tilde{b}=t_{0}^{2}+h^{2}$. Thus, $4\tilde{b}\geq\tilde{a}^{2}$. Furthermore,
we can clearly choose $h\in[0,\infty)$ and $t_{0}\in\mathbb{R}$ to obtain any
$\tilde{a},\tilde{b}\in\mathbb{R}$ with $4\tilde{b}\geq\tilde{a}^{2}$.

  2. (ii)

Let $\|\cdot\|$ be the norm of the Hilbert space $\mathcal{Q}$. By the
Pythagorean theorem,

$\displaystyle\|y-\gamma(t)\|^{2}$
$\displaystyle=\|y-\gamma(t_{0})\|^{2}+\|\gamma(t_{0})-\gamma(t)\|^{2}$ (296)
$\displaystyle=\|y-\gamma(t_{0})\|^{2}+(t-t_{0})^{2}\,,$ (297)

where $\gamma(t_{0})$ is the orthogonal projection of $y$ onto $\gamma$. If
the dimension of $\mathcal{Q}$ is at least $2$, we can choose $y$ to obtain
any value for $t_{0}\in\mathbb{R}$ and $\|y-\gamma(t_{0})\|\in[0,\infty)$.

∎

###### Proof of Proposition 4.3.

  1. (i)

Trivial.

  2. (ii)

If $g_{1},g_{2}$ are parameterized by $t_{0,1},h_{1}$ and $t_{0,2},h_{2}$,
respectively, then squaring the two equations yields

$\displaystyle(s-t_{0,1})^{2}+h_{1}^{2}$
$\displaystyle=(s-t_{0,2})^{2}+h_{2}^{2}\,,$ (298)
$\displaystyle(t-t_{0,1})^{2}+h_{1}^{2}$
$\displaystyle=(t-t_{0,2})^{2}+h_{2}^{2}\,.$ (299)

The difference of these two equations yields

$(s-t_{0,1})^{2}-(t-t_{0,1})^{2}=(s-t_{0,2})^{2}+(t-t_{0,2})^{2}\,,$ (300)

which is equivalent to $(t-s)(t_{0,1}-t_{0,2})=0$. Using $t_{0,1}=t_{0,2}$ in
$f_{1}(s)=f_{2}(s)$ yields $h_{1}^{2}=h_{2}^{2}$.

  3. (iii)

Set $h(x)=g_{2}(x)^{2}-g_{1}(x)^{2}$. Then $h$ is an affine linear function as
the squared terms cancel. Furthermore, $h(r)\geq 0$, $h(s)\leq 0$, and
$h(t)\geq 0$. Thus $h\equiv 0$ and we have $g_{1}=g_{2}$.

  4. (iv)

Follows directly from Lemma B.2.

∎

###### Proof of Proposition 4.4.

  1. (i)

All functions in $\mathcal{G}$ are nonnegative.

  2. (ii)

All functions in $\mathcal{G}$ are $1$-Lipschitz. Let
$t_{1},t_{2}\in\mathbb{R}$. Without loss of generality, assume $f(t_{2})\geq
f(t_{1})$. Let $g\in\mathcal{G}$ be a $\mathcal{G}$-tangent of $f$ at $t_{2}$.
Then $f(t_{2})-f(t_{1})\leq g(t_{2})-g(t_{1})\leq\left|t_{2}-t_{1}\right|$.

  3. (iii)

All functions in $\mathcal{G}$ are convex. Thus, Proposition 4.6 (i) yields
convexity of $f$ (we do not use (iii) in the proof of Proposition 4.6 (i)).

  4. (iv)

Let $g\in\mathcal{G}$ be a $\mathcal{G}$-tangent of $f$ at $t_{1}$. Then
$f(t_{1})+f(t_{2})\geq g(t_{1})+g(t_{2})$. By Lemma B.2,
$g(t_{1})+g(t_{2})\geq\left|t_{1}-t_{2}\right|$.

  5. (v)

The only function $g\in\mathcal{G}$ with $g(t_{0})=0$ is
$\left|t-t_{0}\right|$. As $f$ is $1$-Lipschitz and
$t\mapsto\left|t-t_{0}\right|$ is a lower bound of $f$ with equality at
$t_{0}$, $f$ must also be the function $t\mapsto\left|t-t_{0}\right|$.

∎

###### Proof of Proposition 4.14.

  1. (i)

The functions $g,g^{2}\colon\mathbb{R}\to\mathbb{R}$ are convex. Thus, their
subdifferentials are the closed intervals between the respective left and
right derivative. Let $t_{0}\in\mathbb{R}$. By Proposition 4.6 (ii), for all
$s\in\mathbb{R}$,

$g(s)^{2}\geq
g(t_{0})^{2}+2(s-t_{0})g(t_{0})g^{\oplus}(t_{0})+(s-t_{0})^{2}\,.$ (301)

This inequality holds if we replace $g^{\oplus}(t_{0})$ by
$g^{\ominus}(t_{0})$ (follow the proof of Proposition 4.6 (ii)). Thus it is
true for every subderivative (subgradient in one dimension) in $\partial
g(t_{0})$. Furthermore,

$\partial g^{2}(t_{0})=\left\\{2g(t_{0})v\
\kern-1.2pt\mathopen{}\left|\vphantom{2g(t_{0})v}\ v\in\partial
g(t_{0})\right.\mathclose{}\kern-1.2pt\right\\}\,.$ (302)

Thus, by Proposition 4.13, $g^{2}$ is strongly convex with modulus $1$.

  2. (ii)

Set $g:=\sqrt{f}$. As $f$ is strongly convex with modulus $1$, by Proposition
4.13,

$g(s)^{2}\geq
g(t_{0})^{2}+(s-t_{0})(g^{2})^{\oplus}(t_{0})+(s-t_{0})^{2}:=\tilde{f}(s)\,.$
(303)

As the square root of a degree two polynomial with positive second order
coefficient is convex, $g$ is convex. Thus, its one-sided derivatives exist.
There is $v_{0}\in\partial^{\pm}g(t_{0})$ such that
$(g^{2})^{\oplus}(t_{0})=2g(t_{0})v_{0}$. As we assume $g$ to be
$1$-Lipschitz, we have $v_{0}^{2}\leq 1$ and $(1-v_{0}^{2})g(t_{0})^{2}\geq
0$. We use the first of these two inequalities to show $\tilde{f}\geq 0$.
Hence, we can define $\tilde{g}:=\sqrt{\tilde{f}}$. The second inequality
allows us to write

$\tilde{g}(s)=\sqrt{\left(s-s_{0}\right)^{2}+h^{2}}$ (304)

with $s_{0}=t_{0}+g(t_{0})v_{0}$ and $h^{2}=(1-v_{0}^{2})g(t_{0})^{2}$. Now we
have $g(s)\geq\tilde{g}(s)$ for all $s\in I$ by (303). Furthermore, we
calculate $g(t_{0})=\tilde{g}(t_{0})$. Clearly, $\tilde{g}\in\mathcal{G}$.
Hence, $\tilde{g}$ is a $\mathcal{G}$-tangent of $g$ at $t_{0}$.

∎

###  C.3 From Section 5.4

###### Proof of Remark 5.13.

Assume $\\#M>1$. Let $m,p\in M$ with $p\neq m$. Let $\eta\in(0,1]$. Then, by
Theorem 5,

$0=\mathbb{E}\left[\overline{Yp}-\overline{Ym}\right]\geq\frac{1}{2}\eta^{2}\,\overline{pm}^{2}\,\mathbb{E}\left[\max\\!\left(\overline{Ym},\overline{pm}\right)^{-1}\mathds{1}_{\\!A(m,p,\eta)}(Y)\right]\,.$
(305)

For $r\in(0,\infty)$, we obtain

$0\geq\max\\!\left(r,\overline{pm}\right)^{-1}\mathbb{P}\mathopen{}\left(\left\\{Y\in
A(m,p,\eta)\right\\}\cap\left\\{\overline{Ym}\leq
r\right\\}\right)\mathclose{}\,.$ (306)

Thus, $\mathbb{P}\mathopen{}\left(\left\\{Y\in
A(m,p,\eta)\right\\}\cap\left\\{\overline{Ym}\leq
r\right\\}\right)\mathclose{}=0$ and we can let $r$ be arbitrarily large. If
$\eta^{\prime}\geq\eta$, then $A(m,p,\eta^{\prime})\subseteq A(m,p,\eta)$.
Thus,

$\mathbb{P}\mathopen{}\left(Y\in\bigcup_{\eta\in(0,1]}A(m,p,\eta)\right)\mathclose{}=\lim_{\eta\to
0}\mathbb{P}\mathopen{}\left(Y\in A(m,p,\eta)\right)\mathclose{}=0\,.$ (307)

Thus, one can show $\mathbb{P}(Y\in\tilde{\mathcal{Y}})=1$, where

$\tilde{\mathcal{Y}}:=\mathcal{Q}\setminus\bigcup_{\begin{subarray}{c}p,m\in
M\\\ p\neq m\end{subarray}}\bigcup_{\eta\in(0,1]}A(m,p,\eta)\,.$ (308)

∎

## References

  * [ABY13] Marc Arnaudon, Frédéric Barbaresco and Le Yang  “Medians and means in Riemannian geometry: existence, uniqueness and computation”  In _Matrix information geometry_ Springer, Heidelberg, 2013, pp. 169–197  DOI: 10.1007/978-3-642-30232-9“˙8
  * [ALP20] A. Ahidar-Coutrix, T. Le Gouic and Q. Paris  “Convergence rates for empirical barycenters in metric spaces: curvature, convexity and extendable geodesics”  In _Probab. Theory Related Fields_ 177.1-2, 2020, pp. 323–368  DOI: 10.1007/s00440-019-00950-0
  * [Bač14] Miroslav Bačák  “Computing medians and means in Hadamard spaces”  In _SIAM J. Optim._ 24.3, 2014, pp. 1542–1566  DOI: 10.1137/140953393
  * [Bač14a] Miroslav Bačák  “Convex analysis and optimization in Hadamard spaces” 22, De Gruyter Series in Nonlinear Analysis and Applications  De Gruyter, Berlin, 2014, pp. viii+185  DOI: 10.1515/9783110361629
  * [Bar13] Frédéric Barbaresco  “Information geometry of covariance matrix: Cartan-Siegel homogeneous bounded domains, Mostow/Berger fibration and Fréchet median”  In _Matrix information geometry_ Springer, Heidelberg, 2013, pp. 199–255  DOI: 10.1007/978-3-642-30232-9“˙9
  * [BBI01] Dmitri Burago, Yuri Burago and Sergei Ivanov  “A course in metric geometry” 33, Graduate Studies in Mathematics  American Mathematical Society, Providence, RI, 2001, pp. xiv+415  DOI: 10.1090/gsm/033
  * [BH06] Rajendra Bhatia and John Holbrook  “Riemannian geometry and matrix geometric means”  In _Linear Algebra Appl._ 413.2-3, 2006, pp. 594–618  DOI: 10.1016/j.laa.2005.08.025
  * [BHV01] Louis J. Billera, Susan P. Holmes and Karen Vogtmann  “Geometry of the space of phylogenetic trees”  In _Adv. in Appl. Math._ 27.4, 2001, pp. 733–767  DOI: 10.1006/aama.2001.0759
  * [BJ22] Moise Blanchard and Adam Quinn Jaffe  “Fréchet Mean Set Estimation in the Hausdorff Metric, via Relaxation”  In _arXiv preprint_ , 2022  DOI: 10.48550/arXiv.2212.12057
  * [BP05] Rabi Bhattacharya and Vic Patrangenaru  “Large sample theory of intrinsic and extrinsic sample means on manifolds. II”  In _Ann. Statist._ 33.3, 2005, pp. 1225–1259  DOI: 10.1214/009053605000000093
  * [CCG17] Hervé Cardot, Peggy Cénac and Antoine Godichon-Baggioni  “Online estimation of the geometric median in Hilbert spaces: nonasymptotic confidence balls”  In _Ann. Statist._ 45.2, 2017, pp. 591–614  DOI: 10.1214/16-AOS1460
  * [Cha+94] P. Charbonnier, L. Blanc-Feraud, G. Aubert and M. Barlaud  “Two deterministic half-quadratic regularization algorithms for computed imaging”  In _Proceedings of 1st International Conference on Image Processing_ 2 Los Alamitos, CA, USA: IEEE Computer Society, 1994, pp. 168169170171172  DOI: 10.1109/ICIP.1994.413553
  * [EH19] Benjamin Eltzner and Stephan F. Huckemann  “A smeary central limit theorem for manifolds with application to high-dimensional spheres”  In _Ann. Statist._ 47.6, 2019, pp. 3360–3381  DOI: 10.1214/18-AOS1781
  * [EJ23] Steven N. Evans and Adam Q. Jaffe  “Limit Theorems for Fréchet Mean Sets”  In _arXiv preprint_ , 2023  DOI: 10.48550/arXiv.2012.12859
  * [Eva08] Steven N. Evans  “Probability and real trees” Lectures from the 35th Summer School on Probability Theory held in Saint-Flour, July 6–23, 2005 1920, Lecture Notes in Mathematics  Springer, Berlin, 2008, pp. xii+193  DOI: 10.1007/978-3-540-74798-7
  * [Fol99] Gerald B. Folland  “Real analysis” Modern techniques and their applications, A Wiley-Interscience Publication, Pure and Applied Mathematics (New York)  John Wiley & Sons, Inc., New York, 1999, pp. xvi+386
  * [Fré48] Maurice Fréchet  “Les éléments aléatoires de nature quelconque dans un espace distancié”  In _Ann. Inst. H. Poincaré_ 10, 1948, pp. 215–310  URL: http://www.numdam.org/item?id=AIHP_1948__10_4_215_0
  * [Gre90] Peter J. Green  “Bayesian reconstructions from emission tomography data using a modified EM algorithm.”  In _IEEE transactions on medical imaging_ 9 1, 1990, pp. 84–93  DOI: 10.1109/42.52985
  * [Hub64] Peter J. Huber  “Robust estimation of a location parameter”  In _Ann. Math. Statist._ 35, 1964, pp. 73–101  DOI: 10.1214/aoms/1177703732
  * [Huc11] Stephan F. Huckemann  “Intrinsic inference on the mean geodesic of planar shapes and tree discrimination by leaf growth”  In _Ann. Statist._ 39.2, 2011, pp. 1098–1124  DOI: 10.1214/10-AOS862
  * [Kem87] J… Kemperman  “The median of a finite measure on a Banach space”  In _Statistical data analysis based on the $L_{1}$-norm and related methods (Neuchâtel, 1987)_ North-Holland, Amsterdam, 1987, pp. 217–230
  * [Le ̵+23] Thibaut Le Gouic, Quentin Paris, Philippe Rigollet and Austin J. Stromme  “Fast convergence of empirical barycenters in Alexandrov spaces and the Wasserstein space”  In _J. Eur. Math. Soc. (JEMS)_ 25.6, 2023, pp. 2229–2250  DOI: 10.4171/jems/1234
  * [Mac67] J. MacQueen  “Some methods for classification and analysis of multivariate observations”  In _Proc. Fifth Berkeley Sympos. Math. Statist. and Probability (Berkeley, Calif., 1965/66), Vol. I: Statistics_ Univ. California Press, Berkeley, CA, 1967, pp. 281–297
  * [Min15] Stanislav Minsker  “Geometric median and robust estimation in Banach spaces”  In _Bernoulli_ 21.4, 2015, pp. 2308–2335  DOI: 10.3150/14-BEJ645
  * [MN10] Nelson Merentes and Kazimierz Nikodem  “Remarks on strongly convex functions”  In _Aequationes Math._ 80.1-2, 2010, pp. 193–199  DOI: 10.1007/s00010-010-0043-0
  * [MNO10] Jyrki Möttönen, Klaus Nordhausen and Hannu Oja  “Asymptotic theory of the spatial median”  In _Nonparametrics and robustness in modern statistical inference and time series analysis: a Festschrift in honor of Professor Jana Jurečková_ 7, Inst. Math. Stat. (IMS) Collect.  Inst. Math. Statist., Beachwood, OH, 2010, pp. 182–193
  * [MS23] Stanislav Minsker and Nate Strawn  “The Geometric Median and Applications to Robust Mean Estimation”  In _arXiv preprint_ , 2023  DOI: 10.48550/arXiv.2307.03111
  * [Nik14] Kazimierz Nikodem  “On Strongly Convex Functions and Related Classes of Functions”  In _Handbook of Functional Equations: Functional Inequalities_ New York, NY: Springer New York, 2014, pp. 365–405  DOI: 10.1007/978-1-4939-1246-9˙16
  * [NP18] Constantin P. Niculescu and Lars-Erik Persson  “Convex functions and their applications” A contemporary approach, Second edition of [ MR2178902], CMS Books in Mathematics/Ouvrages de Mathématiques de la SMC  Springer, Cham, 2018, pp. xvii+415  DOI: 10.1007/978-3-319-78337-6
  * [Oht12] Shin-ichi Ohta  “Barycenters in Alexandrov spaces of curvature bounded below”  In _Adv. Geom._ 12.4, 2012, pp. 571–587  DOI: 10.1515/advgeom-2011-058
  * [PM19] Alexander Petersen and Hans-Georg Müller  “Fréchet regression for random objects with Euclidean predictors”  In _Ann. Statist._ 47.2, 2019, pp. 691–719  DOI: 10.1214/17-AOS1624
  * [Pol66] B.. Poljak  “Existence theorems and convergence of minimizing sequences for extremal problems with constraints”  In _Dokl. Akad. Nauk SSSR_ 166, 1966, pp. 287–290
  * [Roc70] R. Rockafellar  “Convex Analysis”, Princeton Landmarks in Mathematics and Physics  Princeton University Press, 1970
  * [Sch19] Christof Schötz  “Convergence rates for the generalized Fréchet mean via the quadruple inequality”  In _Electron. J. Stat._ 13.2, 2019, pp. 4280–4345  DOI: 10.1214/19-EJS1618
  * [Sch22] Christof Schötz  “Strong laws of large numbers for generalizations of Fréchet mean sets”  In _Statistics_ 56.1, 2022, pp. 34–52  DOI: 10.1080/02331888.2022.2032063
  * [Sch23] Christof Schötz  “Quadruple Inequalities: Between Cauchy-Schwarz and Triangle”  In _arXiv preprint_ , 2023  DOI: 10.48550/arXiv.2307.01361
  * [Stu03] Karl-Theodor Sturm  “Probability measures on metric spaces of nonpositive curvature”  In _Heat kernels and analysis on manifolds, graphs, and metric spaces (Paris, 2002)_ 338, Contemp. Math.  Amer. Math. Soc., Providence, RI, 2003, pp. 357–390  DOI: 10.1090/conm/338/06080
  * [Stu05] Karl-Theodor Sturm  “A semigroup approach to harmonic maps”  In _Potential Anal._ 23.3, 2005, pp. 225–277  DOI: 10.1007/s11118-004-7740-z
  * [Stu06] Karl-Theodor Sturm  “On the geometry of metric measure spaces. I”  In _Acta Math._ 196.1, 2006, pp. 65–131  DOI: 10.1007/s11511-006-0002-8
  * [Sve81] Harald Sverdrup-Thygeson  “Strong law of large numbers for measures of central tendency and dispersion of random variables in compact metric spaces”  In _Ann. Statist._ 9.1, 1981, pp. 141–145  URL: http://links.jstor.org/sici?sici=0090-5364(198101)9:1%3C141:SLOLNF%3E2.0.CO;2-Y&origin=MSN
  * [Yan10] Le Yang  “Riemannian median and its estimation”  In _LMS J. Comput. Math._ 13, 2010, pp. 461–479  DOI: 10.1112/S1461157020090531
  * [YP23] Ho Yun and Byeong U. Park  “Exponential concentration for geometric-median-of-means in non-positive curvature spaces”  In _Bernoulli_ 29.4, 2023, pp. 2927–2960  DOI: 10.3150/22-bej1569
  * [Zho18] Xingyu Zhou  “On the Fenchel Duality between Strong Convexity and Lipschitz Continuous Gradient”  In _arXiv preprint_ , 2018  DOI: 10.48550/arXiv.1803.06573