2022

[1,2,3]Soroush Arabi

[2,8]Ruslan Temirov

[1,2,4]Markus Ternes

1]Institute of Physics IIB, RWTH Aachen University, 52074 Aachen, Germany

2]Peter-Grünberg-Institute for Quantum Nanoscience, Research Center Jülich,
52425 Jülich, Germany

3]Max Planck Institute for Solid State Research, Heisenbergstraße 1, 70569
Stuttgart, Germany 4]Jülich Aachen Research Alliance, Fundamentals of Future
Information Technology, 52425 Jülich, Germany

5]Institute of Physics IV, RWTH Aachen University, 52074 Aachen, Germany

6]King’s College London, Theory and Simulation of Condensed Matter (TSCM), The
Strand, London WC2R 2LS, UK

7]Institut de Physique, École Polytechnique Fédérale de Lausanne, 1015
Lausanne, Switzerland

8]Institute of Physics II, University of Cologne, 50937 Cologne, Germany

# Portrait of locally driven quantum phase transition cascades in a molecular
monolayer

<EMAIL_ADDRESS>Taner Esat  Aizhan Sabitova  Yuqi Wang  Hovan Lee
Cédric Weber  Klaus Kern  F. Stefan Tautz<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [  [  [  [  [  [  [

###### Abstract

Strongly interacting electrons in layered materials give rise to a plethora of
emergent phenomena, such as unconventional superconductivity cao2018(1) ;
Yankowitz2019 ; Chen2019 , heavy fermions Vano2021 ; Yazdani2016 ; Hunt2013 ,
and spin textures with non-trivial topology sharpe2019 ; Liu2022 ; Du2009 .
Similar effects can also be observed in bulk materials, but the advantage of
two dimensional (2D) systems is the combination of local accessibility by
microscopic techniques and tuneability. In stacks of 2D materials, for
example, the twist angle can be employed to tune their properties
Novoselov2016 ; Andrei2021 . However, while material choice and twist angle
are global parameters, the full complexity and potential of such correlated 2D
electronic lattices will only reveal itself when tuning their parameters
becomes possible on the level of individual lattice sites. Here, we discover a
lattice of strongly correlated electrons in a perfectly ordered 2D
supramolecular network by driving this system through a cascade of quantum
phase transitions using a movable atomically sharp electrostatic gate. As the
gate field is increased, the molecular building blocks change from a Kondo-
screened to a paramagnetic phase one-by-one, enabling us to reconstruct their
complex interactions in detail. We anticipate that the supramolecular nature
of the system will in future allow to engineer quantum correlations in
arbitrary patterned structures Khajetoorians2019 ; Gomes2012 ; Drost2017 ;
Green2014 .

The prototypical building block in strongly correlated systems is a localised
magnetic moment immersed in a bath of itinerant electrons. Below a
characteristic temperature, the Kondo temperature $T_{\rm K}$, the
hybridisation between the local moment and the itinerant electrons forms a
many-body singlet state at the Fermi energy $E_{\rm F}$, the so-called Kondo
singlet, and results in the quenching of the local moment. Arranging such
building blocks periodically leads to the lattice version of this phenomenon,
i. e., a Kondo lattice. Despite their complexities, these systems are
typically understood in terms of the competition of two antagonistic effects:
The local antiferromagnetic Kondo interaction which screens the local moments
at every site of the lattice by the spin of the conduction electrons, and a
long-range interaction between different local moments that favours magnetic
ordering and breaking the Kondo singlets Doniach1977 . In a subtle balance of
these two effects, close to absolute zero temperature, the systems become
critical, and quantum fluctuations dominate the physics Coleman2005 ;
Lohneysen2007 ; Kirchner2020 .

Whereas the understanding of the Kondo lattice physics has been developed to a
high degree through the study of transition-metal and heavy-fermion crystals,
attempts to access the local information regarding the quantum degrees of
freedom (spins) and to coherently manipulate them Heinrich2021 ; ternes2017
are still intensively stimulating the field. Naturally, this striving has
redirected the approach to realising strongly correlated lattices from three-
dimensional stoichiometric compounds Aynajian2010 ; Kavai2021  towards
artificial spin structures in reduced dimensions and also promoted the
utilisation of local probes for measurements. In this regard, studies have
shown the formation of correlated lattices using moiré superlattices in 2D
transition metal dichalcogenides Vano2021 ; Kumar2022 , and alternatively, in
bottom-up approaches using molecular/atomic building blocks which have been
arranged either by atomic manipulation techniques Figgins2019 ; Moro-
Lagares2019  or by self-assembly Jiang2011 ; Ternes2004 . Yet, despite the
enormous progress, a detailed microscopic understanding, particularly from the
local real-space perspective, is still out of reach.

Therefore, to fill this gap, we study a self-assembled relaxed monolayer of
1,4,5,6-naphthalene tetracarboxylic acid dianhydride (NTCDA) on an Ag(111)
substrate Stahl1998 ; Fink1999 . In gas phase the molecule is charge neutral,
but upon adsorption to Ag(111) its (formerly) lowest unoccupied molecular
orbital shifts below $E_{\rm F}$ ($E_{\rm SOMO}\approx\\!-0.5$ eV) and becomes
a singly occupied molecular orbital (SOMO) with spin 1/2 Bendounan2007 ;
Schoell2010 ; Tonner2016 ; Jakob2016 . In low-temperature scanning tunnelling
microscopy (STM), the molecules are imaged at low bias as dumbbell-shaped
objects forming a brick wall lattice with alternating rows of topographically
higher (brighter) and lower (darker) appearance (Fig. 1a) Eickhoff2019 . The
superlattice unit cell contains, therefore, two slightly different molecules:
a bright one and a dark one.

Figure 1: The correlated molecular lattice: a, Constant-current gray-scale STM
image ($V=-48\,\mbox{mV},I=1\,\mbox{nA}$) of a monolayer of NTCDA molecules on
Ag(111) showing a brick wall like structure consisting of molecules with
higher (brighter) and lower (darker) appearance (yellow and blue dots,
respectively). The apparent height difference is $\approx 30\,\mbox{pm}$.
Superimposed are stick-and-ball models of the molecules and the lattice basis
vectors $\mathbf{a_{1}}$ and $\mathbf{a_{2}}$ with lengths $1.5$ and
$1.16\,\mbox{nm}$. b, Constant-current gray-scale differential conductance
($dI/dV$) image measured at higher bias
($V=1.4\,\mbox{V},I=6\,\mbox{nA},V_{\rm mod}=8\,\mbox{mV}$) revealing ring-
like structures centred at each bright molecule. c, $dI/dV$ spectrum
(setpoint: $V=2\,\mbox{V},I=4\,\mbox{nA},V_{\rm mod}=8$ meV) measured with the
tip of the STM positioned above the centre of a bright molecule shows two
prominent features: the interface state (IS) of the Ag(111) surface, and a
sharp transition peak (TP) which is at the origin of the ring-like structures
seen in panel (b). The inset shows an over the lattice averaged $dI/dV$
spectra close to the Fermi energy $E_{\rm F}$. It reveals molecular vibration
modes (vib), a Kondo resonance, and Hubbard bands (HB) introduced by the
lattice. d, Schematic phase diagram of the lattice versus applied electric
($\vec{E}$) and magnetic ($\vec{B}$) field scaled by temperature $T$ (for
$T\ll T_{\rm K}$ and $\mu_{\rm B}\lvert\vec{B}\rvert\ll k_{\rm B}T_{\rm K}$).
The Kondo screened phase (KS) exists at $\vec{E}<\vec{E}_{\rm C}$ and
undergoes a quantum phase transition to a paramagnetic metallic (PM) phase
constituting free spins. The width $\Gamma$ of the transition regime where
fluctuations dominate the system narrows with the $\vec{B}$-field. A quantum
critical point at $T=0$ is expected to exist where the phases merge.

Interestingly, imaging the same lattice at large positive bias $V>1.3$ V
reveals ring-like structures in the differential conductance $(dI/dV)$ maps,
centred on the bright molecules (Fig. 1b). These rings originate from a strong
transition peak (TP) which can also be seen in the single-point spectroscopy
above the bright molecules (Fig. 1c). In the following, we will show that the
TP indicates a quantum phase transition (QPT) in the molecular layer in which
the electric field of the STM tip drives the Kondo-screened molecules of the
lattice successively into a paramagnetic free spin phase (Fig. 1d).

We focus now on energies close to $E_{\rm F}$ where we find very complex and
strongly position-dependent spectra in the lattice (see Supplementary Note 1).
A spectrum spatially averaged over the lattice (Fig. 1c inset) reveals a broad
peak close to zero bias as a manifestation of Kondo singlets. Furthermore, the
spectrum contains conductance steps at $\approx\pm\,47$ meV, which are due to
inelastic vibrational excitations of the individual molecules Braatz2012 ;
Tonner2016 ; Eickhoff2019 . At the flanks of the Kondo peak, we see well-
resolved electron- and hole-like quasiparticle excitations, which will turn
out to be Hubbard bands (HB).

To disentangle the spectroscopic features, we analyse individual spectra taken
on a grid of $58\times 58$ points on the molecular lattice (see Methods). In
this dataset, we search for peak structures and plot their weighted occurrence
in an energy distribution histogram (Fig. 2a). This histogram shows two peaks
very close to $E_{\rm F}$ whose averaged spectra correspond to the Kondo
resonances of the bright and dark molecules, apparent in the occurrence maps
(Fig. 2b II and III). Compared to the bright molecules’ $T_{\rm K}=132\pm 5$
K, the averaged spectrum of the dark molecules’ Kondo resonance reveals, in
agreement with ref. Eickhoff2019 , a higher $T_{\rm K}=340\pm 20$ K due to a
stronger hybridisation with the Ag(111). It shows instead of a peak an
asymmetric Fano line-shape which is the result of interference between two
possible tunnelling channels in the STM junction, i. e., electrons can tunnel
from the tip either via the resonance or directly into the substrate Fano1961
; Ternes2009 . The Kondo resonances are strongly localised at the ends of the
dumbbell-like topographic shape of the molecules, which correspond to their
C–H sites Eickhoff2019 . At the centre of the molecule, however, the resonance
vanishes. Such distribution cannot originate from a symmetric orbital, but
points to a $\pi$-like molecular orbital with a nodal plane along the
molecular long axis Eickhoff2019 . This $\pi$-like orbital contains the
spatially distributed spin with a phase shift $\theta=\pi$ between the two
lobes and enables a significant intermolecular coupling Esat2016 .

Figure 2: Low energy features of the lattice. a, Top: Energy distribution
histogram of peak occurrences in a map of $58\times 58$ $dI/dV$ spectra
(setpoint for each spectrum: $V=70\,\mbox{mV},I=1\,\mbox{nA},V_{\rm
mod}=0.2\,\mbox{mV}$) taken on the lattice. The full (hatched) bars mark areas
at negative (positive) bias ranges where spectra have been averaged or which
are mapped (panel b). Middle: Averaged Kondo resonance measured at bright and
dark molecules and corresponding Fano fits (full lines)
($\Gamma=12.2\,(28.7)\,\mbox{meV},q>300\,(=1.2),\epsilon_{0}=2\,(21)\,\mbox{meV},T_{\rm
K}=132\pm 5\,(340\pm 20)\,\mbox{K}$; values in brackets belong to dark
molecules). Bottom: Averaged spectra of the LHB (blue dots) and UHB (black
dots). b, Color-coded maps of the peak feature intensity in the energy ranges
as marked in panel (a) and superimposed stick-and-ball model of the molecules.

Also, the states in the flanks of the Kondo peak are clearly distinguishable
in the histogram with the maximum occurrence at $-23$ mV and $+20$ mV. Due to
their low intensity, their averaged spectra are outshone by the vibrational
excitations and the Kondo resonance. Their mapping (see Methods), however,
clearly reveals that the occupied state at negative bias is most strongly
located between the molecules, while the unoccupied one is localised at each
dark and bright molecule’s centre (Fig. 2b I and IV). This distribution, which
is reminiscent of bonding and antibonding combinations between the $\pi$
states on the three neighbouring molecules indicated in Fig. 2b, reveals that
the flank peaks are a consequence of the formation of a correlated lattice and
hence can be naturally understood as lower and upper Hubbard bands (LHB/UHB)
coleman_2015 . They have a composite structure containing contributions from
the localised $\pi$-electrons of the molecules as well as from the conduction
electrons. Note that their $U^{\ast}\approx 43$ meV is much smaller than the
Coulomb repulsion $U$ of the SOMO. The additional bands energetically further
away from $E_{\rm F}$ also show spatial dependence (see Supplementary Note 2),
but are for the following not important.

We now turn to the prominent TP feature of Fig. 1c by recording $dI/dV(z)$
spectra above the centre of the bright molecules by changing the tip-sample
tunnelling distance $z$ at successively reduced constant bias $V$ (Fig. 3a).
We define $z=0$ as the extrapolated height at which quantum of conductance
would be reached. The $z_{\rm TP}$ at which the peak occurs is proportional to
$V_{\rm TP}$ (Fig. 3c), leading us to an electric field driven origin of the
TP. Simplifying the STM junction as a plate capacitor with an electrical field
between tip apex and sample of $\vec{E}=(V/z)\vec{e}_{z}$, with $\vec{e}_{z}$
as the unit vector along the $z$-direction (Fig. 3b), the linear relation
corresponds to a critical field $\lvert\vec{E}_{\rm C}\rvert=(1.77\pm 0.02)$
GVm-1 at which the TP occurs.

Figure 3: Signatures of the quantum phase transition. a, Logarithmic plot of
the $dI/dV(z)$ signal above a bright molecule at constant $V$ ranging from
$1.4-2.0$ V and $B_{z}=5$ T. The peaks mark the QPT between the KS and the PM
state. The inset shows exemplarily the $V=1.7\,$V data in a linear plot and a
least-squares fit to the sum of an exponential and a Gaussian. b, Model of the
tip-sample tunnelling barrier $z$ and its change to $z+\Delta z$. c, Extracted
Gauss centre $z_{\rm TP}$ from (a) versus $V_{\rm TP}$. d, Extracted $\Delta
z$ from simultaneously measured $I(z)$ data in (a). e, Logarithmic plot of the
$I(z)$ signal for $V=1.7\,$V with exponential fits (dashed lines) separating
the two regimes of the system. f, Half-width at half-maximum $\Gamma$ of the
phase transition at different $\vec{B}$-fields and $T$ (dots) measured at
$V=1.7V$, and least-square model ($R^{2}=0.91$) of a QPT to a paramagnetic
spin 1/2 (solid lines).

The TP in $dI/dV$ is equivalent to an increase in conductance, as the example
of a simultaneously recorded $I(z)$ trace in Fig. 3e illustrates. Until the TP
is reached the current increases exponentially with $I\propto\exp(-\kappa z)$
and a $\kappa=15.3$ nm-1. Past the TP, $I(z)$ increases further with identical
$\kappa$, excluding a modification of the tunnel barrier by any significant
charge accumulation on the molecule. Thus, the TP is tantamount to a change of
the effective tip-sample distance from $z$ to $z-\Delta z$. This $\Delta
z\approx 10$ pm is bias independent (Fig. 3d), meaning that its origin lies in
an increase of the molecular lattice conductivity.

Gating of a molecular charge state by exerting the $\vec{E}$-field of an STM
tip has been observed Fernandez_Torrente2012 ; Pham2019 ; Kumar2019 ;
Mohammed2020 ; Li2021 ; Yan2021  and is usually explained in the framework of
a double barrier tunnel junction model in which the electric field shifts a
single particle molecular state across $E_{\rm F}$, thereby changing the
number of electrons on the molecule. This is different here. The only charge-
like (and therefore $\vec{E}$ gateable) state which lies close to $E_{\rm F}$
is the LHB which is a composite state of the lattice with relatively weak
intensity. The prominent Kondo resonance is fixed to $E_{\rm F}$ and,
therefore, cannot be gated.

The width $\Gamma$ of the TP shows a peculiar behaviour: With increasing
temperature $T$, $\Gamma$ increases linearly as expected from the thermal
broadening of the Fermi edge of the tip (Fig. 3f red line). However, an
additionally applied $\vec{B}$-field perpendicular to the surface of strength
$g\mu_{\rm B}\lvert\vec{B}\rvert>k_{\rm B}T$ (with $\mu_{\rm B}$ as Bohr’s
magneton, and $k_{\rm B}$ as Boltzmann’s constant) decreases $\Gamma$
significantly, as shown in Fig. 3f. This behaviour is reminiscent of the
freezing out of spin-fluctuations in an ideal paramagnetic metal with $g=2$
and can be well described by assuming a broadening due to the sum of intrinsic
quantum and thermal fluctuations:

$\Gamma=\Gamma_{\rm inst}+\Gamma_{0}\left(\overline{n}+\gamma\frac{k_{\rm
B}T}{E_{\rm C}}\sqrt{\overline{n}}\right),$ (1)

with $\Gamma_{\rm inst}=1.6\pm 1.3$ pm as the finite resolution of our
experimental setup and $\Gamma_{0}=7.7\,\pm 0.9$ pm as the extrapolated
quantum fluctuations per (spin) degree of freedom at $T=0$ where all thermal
fluctuations are frozen out. The average number of available spin states,
$\overline{n}=\exp(S/k_{\rm B})$ is given by the entropy $S$ of the
paramagnetic spin 1/2 in the field $\vec{B}$ and at temperature $T$.
$\gamma=4.0\,\pm\,1.3$ is the slope of the $T$-dependent broadening of the TP.
Note, however, that this model neglects any magnetic interactions between
neighbouring spins. Significant magnetic ordering would reduce $\overline{n}$
by freezing out fluctuations similar to the $\vec{B}$-field here (for further
details see Methods).

The proportionality between $z_{\rm TP}$ and $V_{\rm TP}$ allows us to express
$\Gamma_{0}$ also as a width in $V$,
$\Gamma_{0}^{\prime}=\Gamma_{0}\times\lvert\vec{E}_{\rm C}\rvert=14\pm 2$ mV
(see Supplementary Note 3). From the applied bias $V$ only a small fraction
$\alpha$ shifts the LHB (see Supplementary Note 4), which we can estimate for
$V_{\rm TP}=1.3$ V as $\alpha\approx E_{\rm LHB}/eV_{\rm TP}\approx 1/55$.
This allows to estimate the width of the quantum fluctuation per degree of
freedom in energy as $\Delta_{0}=\alpha\times\Gamma^{\prime}_{0}\approx 0.3$
meV, and with this the timescale of the quantum fluctuation
$\tau_{0}=\hbar/\Delta_{0}\approx 5$ ps (with $\hbar$ as the reduced Planck’s
constant) of the QPT extrapolated for $T=0$.

This QPT is driven by the $\vec{E}$-field of the tip. For a radially symmetric
tip, the field decreases with the lateral distance $d=\sqrt{x^{2}+y^{2}}$ to
the tip apex as $\vec{E}(d)=\vec{E}_{d=0}/\left(1+(d/r)^{\beta}\right)$,
whereby $r$ and $\beta$ approximate the decay of the $\vec{E}$-field due to
the shape of the tip Battisti2017 . At biases and tip heights where
$\vec{E}_{d=0}>\vec{E}_{\rm C}$ the TP occurs at a certain lateral distance
from the molecular centre, leading to growing ring-like structures when
$dI/dV$ images are taken at successively increased $V$. Fig. 4a shows a set of
such images. At $V=1.6$ V similar rings as in Fig. 1b form. The small
deviation from perfect circularity can be assigned to some tip shape
asymmetries. When rings from neighbouring bright molecules cross along the
short $\mathbf{a}_{2}$ direction, they overlap. Contrarily, when they cross
along the long $\mathbf{a}_{1}$ direction a gap opens (at $V=1.85$ V),
signalling an avoided crossing due to intermolecular interaction Li2021 .

Figure 4:  A portrait of strongly correlated quasiparticle excitations. a,
Colour-coded $dI/dV$ plots measured in an area of approx. $1.5\times 1.5\,$nm2
with a bright molecule in the centre ($I=4.9$ nA, $B_{z}=10$ T). b, TB
simulation using a nearest-neighbour hybridisation $V_{\rm AB}e^{i\theta}$,
with $V_{\rm AB}=4$ meV and $\theta=0.2\pi$ between bright and dark molecules
(arrows; their direction indicates hopping with positive phase shift), a next-
nearest neighbour hoping $t_{\rm A}=0.25(0.05)$ $\mu$eV between bright
molecules in direction $\vec{a}_{1}$($\vec{a}_{2}$) (dashed lines), $t_{\rm
B}=0$, and on-site energies (yellow: bright $E_{\rm A}=-23$ meV, blue: dark
$E_{\rm B}=20$ meV). For the tip we used $r=16$ nm and $\beta=2$. The numbers
in the last sub-panel indicate the number of QPTs in the lattice. c,
Calculated $V_{\rm TP}$ for a tip at different positions along $\vec{a}_{1}$.
Insets illustrate the tip position above or in between two bright molecules
for the TPs occurring at $V_{\rm TP}^{0,+,-}$ whose wavefunctions
$\Psi_{0,+,-}$ (coloured curves) and on-site energies $E_{\rm A}+U^{\rm eff}$
(black bars) of neighbouring sites are shown in panel d-f, respectively.

To rationalise the ring formation in our bipartite lattice we use a single-
particle tight-binding model in which we include the effect of the
$\vec{E}$-field of the tip. We use a complex hybridisation amplitude $V_{\rm
AB}e^{i\theta}$ between nearest-neighbour bright (A) and dark (B) molecules
and different next-nearest-neighbour hopping between bright molecules along
the directions $\vec{a}_{1}$ and $\vec{a}_{2}$ as illustrated in Fig. 4b, left
panel (for details see Methods).

The use of a finite phase $\theta$ is inspired by the seminal work of Haldane
on graphene, in which he assumes a phase-dependent hopping between next-
nearest neighbours Haldane1988 , which is based on the observation that
electrons face different energy landscapes on their left and right while
hopping forward and backward Kane2005 . Similarly, but for nearest-neighbours,
this argument holds for the NTCDA lattice: electrons hopping from the marked
dark molecule in Fig. 2b to the lower-left bright molecule face on their right
side the broad Kondo resonance of the dark molecule, while on the left side,
they face the strong Kondo resonance of the bright molecule (see maps in Fig.
2). For electrons in opposite directions, the situation is reversed. This
symmetry breaking does not occur for hopping between two bright or two dark
molecules.

With this model, we can reproduce our observations as the comparison of Fig.
4a and 4b shows. When the tip is placed above a bright molecule at site $i$
and $V=V^{0}_{\rm TP}$ (Fig. 4c), the wavefunction $\Psi_{0}$ which crosses
$E_{\rm F}$ is mainly localised at the site $i$ underneath the tip (Fig. 4d).
However, when the tip is placed between two bright molecules along the
$\vec{a}_{1}$ direction, two transitions occur: At $V=V^{-}_{\rm TP}$ the
antisymmetric superposition
$\Psi_{-}\approx\sqrt{0.5}(\phi_{\text{A},i}-\phi_{\text{A},i+1})$ crosses
$E_{\rm F}$ (Fig. 4e), where $\phi_{\text{A},i}$ is the orbital wavefunction
of the individual bright molecule at site $i$. This is then followed by the
observed gap between the rings, and finally by the crossing of the symmetric
superposition
$\Psi_{+}\approx\sqrt{0.5}(\phi_{\text{A},i}+\phi_{\text{A},i+1})$ from
$E_{\rm F}$ at $V=V^{+}_{\rm TP}$ (Fig. 4f). Thus, at this tip position, each
QPT involves a coherent superposition between two neighbours.

When we increase the bias further, more and more rings cross each other, and
when the tip is placed on a dark molecule, as many as four lattice sites are
coherently brought into the paramagnetic state at $V=2.22$ V. This is
remarkable because even though the gated molecules have a lateral distance of
$0.95$ nm from the tip apex, we only see a minor reduction of the TP intensity
compared to the situation when the tip is on top of a bright molecule, clearly
revealing the coherent lattice nature of the effect. Importantly, without the
complex hybridisation between the sublattices, we cannot reproduce correctly
the experimentally observed intensity modulation (see Supplementary Note 5).

To summarise, taking advantage of sub-angstrom-resolution spectroscopic
imaging offered by STM, we unraveled the role of strong electronic
correlations in the formation of Hubbard bands and Kondo resonances in a self-
assembled molecular lattice. Build out of simple molecular building blocks,
the lattice can reach macroscopic dimensions. We demonstrated that the ground
state of the lattice can be controlled and manipulated by the local electric
field of the STM tip, which drives the charge localised in the Hubbard bands
across $E_{\rm F}$ and leads to a fundamental change of the lattice properties
from a KS lattice to a PM metal. We found clear signals of the quantum nature
of this transition as finite fluctuations drive the transition even at zero
temperature. For specific tip positions, interactions between neighbouring
building blocks lead to a closed sequence of QPTs involving several molecules.
These sequences manifest themselves in ring-like structures which coherently
interact with each other even nanometres away from the tip’s position and
create an interference pattern in real space. Finally, we showed that these
interference patterns in the $dI/dV$ maps can be effectively simulated with a
single-particle model Hamiltonian which has a complex hybridisation.

The observed QPTs shed light on the fragility of the correlations in
supramolecular Kondo lattices. Even though the single-ion Kondo effect in our
system has a high $T_{\rm K}$ originating from the strong interactions between
the substrate and the SOMO, which would require a very large $\vec{E}$-field
to be brought to $E_{\rm F}$, the lattice singlets are fragile objects that
can be broken by relatively moderate fields.

Our study not only expands the current understanding of strong correlations
and topological properties in 2D materials but also offers a new
supramolecular platform for further explorations. This can be done
particularly by engineering complex structures out of molecular building
blocks with the possibility of arbitrary pattering using molecular
manipulation techniques and application of long-range gate operations at
atomic scales.

## Methods

### Experimental procedure

The measurements were taken using a home-built combined STM and atomic force
microscope operating in ultrahigh vacuum ($\leq 10^{-10}$ mbar), at fields
perpendicular to the sample surface of up to $B_{z}=14$ T, and, if not
mentioned otherwise, at a base temperature of $T=1.2-1.4$ K. $dI/dV$ spectra
were acquired by modulating the bias voltage with a sinusoidal $V_{m}$ of
frequency $617$ Hz and employing a lock-in amplifier. Note, the number of
gated molecules is only limited by tip induced instabilities of the molecular
lattice that start to occur at biases $V\gtrsim 2.4$ V. The bare Pt tip was
conditioned by indentation into the Ag(111) surface. The Ag(111) crystal was
prepared by repeated cycles of Ar+ sputtering and annealing. Afterwards, NTCDA
molecules were deposited onto the clean Ag(111) surface held at room
temperature from a Knudsen cell molecule evaporator.

### Data analysis and fitting

The statistical analysis of the $dI/dV$ grid spectra in Fig. 2 was performed
by using feature-detection scanning tunnelling spectroscopy (FD STS)
Sabitova2018 . For that a second order Savitzky-Golay filter (SGF) with a
sliding window of $\pm 2\mbox{ mV}$ was applied at each data point. The SGF
performs a quadratic fit $(y=a_{i}\,x^{2}+b_{i}\,x+c_{i})$, enabling the
determination of peak positions when $b_{i}$ changes sign and $a_{i}>0$.
Subsequently, a weight $w$ is assigned to each peak by determining the
intensity difference $b_{i}-b_{i+1}$ between the neighbouring extrema in the
first derivative of the signal, i. e. at $a_{i}=0$. This $w$ is then mapped
producing the images in Fig. 2b, and used as weight for the averaged spectra
shown in Fig. 2a.

The fits of Fig. 2a are least-square fits of the Fano line-shape $g(eV)\propto
g_{0}+(q+\epsilon)^{2}/(1+\epsilon^{2})$, with $g$ as the differential
conductance ($dI/dV$), $g_{0}$ as the background conductance,
$\epsilon=(eV-\epsilon_{0})/\Gamma$, $\epsilon_{0}$ as the center of the
resonance, $\Gamma$ as its width, and $q$ as the Fano form-factor Fano1961 .
Because $\Gamma\ll k_{\rm B}T$ and $\Gamma\ll eV_{m}$, we neglected
temperature and modulation broadening effects.

The empirical Eq. 1 is motivated by the following assumptions: The quantum
fluctuations are linearly dependent on the number of available states which
can be reached when crossing the QPT, i. e., for the transition from the non-
magnetic Kondo singlet state to the paramagnetic state it depends on the ratio
between magnetic field energy $g\mu_{\rm B}\vec{B}$ and thermal energy $k_{\rm
B}T$. This changes from $\overline{n}$=1 when the spin degrees are frozen out,
to $\overline{n}$=2 when $\vec{B}=0$ resulting in a total broadening of
$\overline{n}\,\Gamma_{0}^{\prime}$. At $T>0$, an additional thermal
broadening $\propto k_{\rm B}T$ occurs due to the Fermi-Dirac distribution of
the tip’s electrons which smears out any measurement in STM spectroscopy.
Additionally, we have to account for the thermally activated broadening of the
paramagnetic spin when $\overline{n}>1$. Both thermal broadening can be
combined to $\propto k_{\rm B}T\sqrt{\overline{n}}$. Thus, for
$\overline{n}=1$ only the tip effect is accounted for, while at
$\overline{n}=2$, the two thermally activated broadening sources are
geometrically added to the prefactor $\sqrt{2}$.

### Tight-binding calculation

We used a single-particle tight-binding Hamiltonian:

$\mathcal{H}=\mathcal{H}_{\rm sub}+\mathcal{H}_{\rm hyb}.$ (2)

Here, $H_{\rm sub}$ accounts for the on-site energies and the hopping within
each sublattice (A$=$ bright, B$=$ dark molecules) and is given by:

$\mathcal{H}_{\rm sub}=\sum_{\nu=\rm
A,B}\left[\sum_{i}\left(E_{\nu}+U_{\nu,i}^{\rm
eff}\right)\hat{c}_{\nu,i}^{\dagger}\hat{c}_{\nu,i}-\sum_{\langle\langle
i,j\rangle\rangle}t_{\nu,ij}\left[\hat{c}_{\nu,i}^{\dagger}\hat{c}_{\nu,j}+h.c.\right]\right],$
(3)

and $\mathcal{H}_{\rm hyb}$ for the hybridisation between the two sublattices:

$\mathcal{H}_{\rm hyb}=-\sum_{\langle i,j\rangle}\left[V_{\rm
AB}\,e^{i\theta_{ij}}\,\hat{c}_{\text{A},i}^{\dagger}\hat{c}_{\text{B},j}+h.c.\right].$
(4)

The creation $\hat{c}^{\dagger}_{\rm\nu,i}$ and annihilation operators
$\hat{c}_{\rm\nu,i}$ act on the quasiparticle at site $i$ of the sublattice
$\nu$. The single bracket $\langle i,j\rangle$ and double bracket
$\langle\langle i,j\rangle\rangle$ in the sums stand for the nearest- and
next-nearest-neighbour sites, where we allowed the hopping amplitude
$t_{\nu,ij}$ to differ in each sub-lattice ($\nu$) and also along the short
and long unit vectors. We used the energies of the LHB and UHB for the on-site
energies $E_{\rm A}$ and $E_{\rm B}$. However, the simulations remain
indifferent to the exact value of $E_{\rm B}$ as long as $E_{\rm B}>0$. To
include the local $\vec{E}$-field of the tip we added at each lattice site
$U_{\nu,i}^{\rm eff}=\alpha z\,\vec{e}_{z}\cdot\vec{E}(d_{\nu,i})$, where
$d_{\nu,i}$ is the lateral distance between the tip apex and site $i$ at
sublattice $\nu$. The $dI/dV$ maps of Fig. 4b were then calculated by
diagonalization of Eq. 2 at each pixel and by taking the absolute square of
the coherent sum of all spatially extended wavefunctions $\Psi_{n}$ which
energy $\xi_{n}$ lie within the derivative of Fermi-Dirac distribution

$dI/dV\propto\sum_{n}\left\lvert\Psi_{n}\right\rvert^{2}\left.\frac{\partial
f_{\rm F}(\varepsilon)}{\partial\varepsilon}\right|_{\varepsilon=\xi_{n}}.$
(5)

Note, that each wavefuction
$\Psi=\sum_{\nu=\text{A,B}}\sum_{i}a_{\nu,i}\phi_{\nu,i}$ is a sum of
spatially extending Warnier-like orbitals $\phi_{\nu,i}$ which are centred at
the molecules at site $(\nu,i)$, with complex coefficients $a_{\nu,i}$. We
used $s$-like orbitals for the representation in Fig. 4d-f, whereby we mark
that the shape of the $\phi_{\nu,i}$ did not play a role in the calculation of
the $dI/dV$. Acknowledgments We thank Christian Ast, Johann Kroha, Christopher
Leon, and Nicolás Lorente for fruitful discussions. M.T. acknowledges support
from the Heisenberg Program (Grant No. TE 833/2-1) and the Priority Program
(SPP 2244) of the German Research Foundation. Author contributions R.T.,
F.S.T., and M.T. conceived the experiment. S.A., T.E., R.T., A.S., and Y.W.
performed the STM measurements. S.A., M.T., T.E., and R.T. analysed the
experimental data. S.A. and M.T. developed the theoretical modelling and
tight-binding calculations. S.A. and M.T. wrote the manuscript with
significant input from T.E., H.L., C.W., F.S.T., and R.T. All authors
discussed the results and reviewed the manuscript.

## References

  * (1) Cao, Y. _et al._ Unconventional superconductivity in magic-angle graphene superlattices.  _Nature_ 556, 43–50 (2018).
  * (2) Yankowitz, M. _et al._ Tuning superconductivity in twisted bilayer graphene.  _Science_ 363, 1059–1064 (2019).
  * (3) Chen, G. _et al._ Signatures of tunable superconductivity in a trilayer graphene moiré superlattice.  _Nature_ 572, 215–219 (2019).
  * (4) Vaňo, V. _et al._ Artificial heavy fermions in a van der Waals heterostructure.  _Nature_ 599, 582–586 (2021).
  * (5) Yazdani, A., da Silva Neto, E. H. & Aynajian, P.  Spectroscopic imaging of strongly correlated electronic states.  _Annual Review of Condensed Matter Physics_ 7, 11–33 (2016).
  * (6) Hunt, B. _et al._ Massive Dirac Fermions and Hofstadter butterfly in a van der Waals heterostructure.  _Science_ 340, 1427–1430 (2013).
  * (7) Sharpe, A. L. _et al._ Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene.  _Science_ 365, 605–608 (2019).
  * (8) Liu, X. _et al._ Visualizing broken symmetry and topological defects in a quantum Hall ferromagnet.  _Science_ 375, 321–326 (2022).
  * (9) Du, X., Skachko, I., Duerr, F., Luican, A. & Andrei, E. Y.  Fractional quantum Hall effect and insulating phase of Dirac electrons in graphene.  _Nature_ 462, 192–195 (2009).
  * (10) Novoselov, K. S., Mishchenko, A., Carvalho, A. & Neto, A. H. C.  2D materials and van der Waals heterostructures.  _Science_ 353, aac9439 (2016).
  * (11) Andrei, E. Y. _et al._ The marvels of moiré materials.  _Nature Reviews Materials_ 6, 201–206 (2021).
  * (12) Khajetoorians, A. A., Wegner, D., Otte, A. F. & Swart, I.  Creating designer quantum states of matter atom-by-atom.  _Nature Reviews Physics_ 1, 703–715 (2019).
  * (13) Gomes, K. K., Mar, W., Ko, W., Guinea, F. & Manoharan, H. C.  Designer Dirac fermions and topological phases in molecular graphene.  _Nature_ 483, 306–310 (2012).
  * (14) Drost, R., Ojanen, T., Harju, A. & Liljeroth, P.  Topological states in engineered atomic lattices.  _Nature Physics_ 13, 668–671 (2017).
  * (15) Green, M. F. B. _et al._ Patterning a hydrogen-bonded molecular monolayer with a hand-controlled scanning probe microscope.  _Beilstein Journal of Nanotechnology_ 5, 1926–1932 (2014).
  * (16) Doniach, S.  The Kondo lattice and weak antiferromagnetism.  _Physica B+C_ 91, 231–234 (1977).
  * (17) Coleman, P. & Schofield, A. J.  Quantum criticality.  _Nature_ 433, 226–229 (2005).
  * (18) Löhneysen, H. v., Rosch, A., Vojta, M. & Wölfle, P.  Fermi-liquid instabilities at magnetic quantum phase transitions.  _Rev. Mod. Phys._ 79, 1015–1075 (2007).
  * (19) Kirchner, S. _et al._ Colloquium: Heavy-electron quantum criticality and single-particle spectroscopy.  _Rev. Mod. Phys._ 92, 011002 (2020).
  * (20) Heinrich, A. J. _et al._ Quantum-coherent nanoscience.  _Nature Nanotechnology_ 16, 1318–1329 (2021).
  * (21) Ternes, M.  Probing magnetic excitations and correlations in single and coupled spin systems with scanning tunneling spectroscopy.  _Prog. Surf. Sci._ 92, 83–115 (2017).
  * (22) Aynajian, P. _et al._ Visualizing the formation of the Kondo lattice and the hidden order in URu2Si2.  _Proceedings of the National Academy of Sciences_ 107, 10383–10388 (2010).
  * (23) Kavai, M. _et al._ Inhomogeneous Kondo-lattice in geometrically frustrated Pr2Ir2O7.  _Nature Communications_ 12, 1377 (2021).
  * (24) Kumar, A., Hu, N. C., MacDonald, A. H. & Potter, A. C.  Gate-tunable heavy fermion quantum criticality in a moiré Kondo lattice.  _Phys. Rev. B_ 106, L041116 (2022).
  * (25) Figgins, J. _et al._ Quantum engineered Kondo lattices.  _Nature Communications_ 10, 55–88 (2019).
  * (26) Moro-Lagares, M. _et al._ Real space manifestations of coherent screening in atomic scale Kondo lattices.  _Nature Communications_ 10, 2211 (2019).
  * (27) Jiang, Y., Zhang, Y. N., Cao, J. X., Wu, R. Q. & Ho, W.  Real-space imaging of Kondo screening in a two-dimensional O2 lattice.  _Science_ 333, 324–328 (2011).
  * (28) Ternes, M. _et al._ Scanning tunneling spectroscopy of surface-state electrons scattered by a slightly disordered two-dimensional dilute “solid”: Ce on Ag(111).  _Phys. Rev. Lett._ 93, 146805 (2004).
  * (29) Stahl, U., Gador, D., Soukopp, A., Fink, R. & Umbach, E.  Coverage-dependent superstructures in chemisorbed NTCDA monolayers: a combined LEED and STM study.  _Surface Science_ 414, 423–434 (1998).
  * (30) Fink, R., Gador, D., Stahl, U., Zou, Y. & Umbach, E.  Substrate-dependent lateral order in naphthalene-tetracarboxylic-dianhydride monolayers.  _Phys. Rev. B_ 60, 2818–2826 (1999).
  * (31) Bendounan, A. _et al._ Electronic structure of 1ML NTCDA/Ag(111) studied by photoemission spectroscopy.  _Surface Science_ 601, 4013–4017 (2007).
  * (32) Schöll, A. _et al._ Disordering of an organic overlayer on a metal surface upon cooling.  _Science_ 329, 303–305 (2010).
  * (33) Tonner, R., Rosenow, P. & Jakob, P.  Molecular structure and vibrations of NTCDA monolayers on Ag(111) from density-functional theory and infrared absorption spectroscopy.  _Phys. Chem. Chem. Phys._ 18, 6316–6328 (2016).
  * (34) Jakob, P. _et al._ Adsorption geometry and interface states: Relaxed and compressed phases of NTCDA/Ag(111).  _Phys. Rev. B_ 94, 125–436 (2016).
  * (35) Eickhoff, F. _et al._ Inelastic electron tunneling spectroscopy for probing strongly correlated many-body systems by scanning tunneling microscopy.  _Phys. Rev. B_ 101, 125–405 (2020).
  * (36) Braatz, C. R., Öhl, G. & Jakob, P.  Vibrational properties of the compressed and the relaxed 1,4,5,8-naphthalene-tetracarboxylic dianhydride monolayer on Ag(111).  _The Journal of Chemical Physics_ 136, 134706 (2012).
  * (37) Fano, U.  Effects of configuration interaction on intensities and phase shifts.  _Phys. Rev._ 124, 1866–1878 (1961).
  * (38) Ternes, M., Heinrich, A. J. & Schneider, W.-D.  Spectroscopic manifestations of the Kondo effect on single adatoms.  _Journal of Physics: Condensed Matter_ 21, 053001 (2009).
  * (39) Esat, T. _et al._ A chemically driven quantum phase transition in a two-molecule kondo system.  _Nature Physics_ 12, 867–873 (2016).
  * (40) Coleman, P.  _Introduction to Many-Body Physics_ (Cambridge University Press, 2015).
  * (41) Fernández-Torrente, I., Kreikemeyer-Lorenzo, D., Stró$\dot{\text{z}}$ecka, A., Franke, K. J. & Pascual, J. I.  Gating the charge state of single molecules by local electric fields.  _Phys. Rev. Lett._ 108, 036–801 (2012).
  * (42) Pham, V. D. _et al._ Selective control of molecule charge state on graphene using tip-induced electric field and nitrogen doping.  _npj 2D Materials and Applications_ 3, 5 (2019).
  * (43) Kumar, D., Krull, C., Yin, Y., Medhekar, N. V. & Schiffrin, A.  Electric field control of molecular charge state in a single-component 2D organic nanoarray.  _ACS Nano_ 13, 11882–11890 (2019).
  * (44) Mohammed, M. S. G. _et al._ Electronic decoupling of polyacenes from the underlying metal substrate by $sp^{3}$ carbon atoms.  _Communications Physics_ 3, 159 (2020).
  * (45) Li, H. _et al._ Imaging local discharge cascades for correlated electrons in WS2/WSe2 moiré superlattices.  _Nature Physics_ 17, 1114–1119 (2021).
  * (46) Yan, L. _et al._ Synthesis and local probe gating of a monolayer metal-organic framework.  _Advanced Functional Materials_ 31, 2100519 (2021).
  * (47) Battisti, I. _et al._ Poor electronic screening in lightly doped Mott insulators observed with scanning tunneling microscopy.  _Phys. Rev. B_ 95, 235141 (2017).
  * (48) Haldane, F. D. M.  Model for a quantum Hall effect without Landau levels: Condensed-matter realization of the “parity anomaly”.  _Phys. Rev. Lett._ 61, 2015–2018 (1988).
  * (49) Kane, C. L. & Mele, E. J.  Quantum spin Hall effect in graphene.  _Phys. Rev. Lett._ 95, 226801 (2005).
  * (50) Sabitova, A., Temirov, R. & Tautz, F. S.  Lateral scattering potential of the PTCDA/Ag(111) interface state.  _Phys. Rev. B_ 98, 205429 (2018).


# New stellar velocity substructures from Gaia DR3 proper motions

Daniel Mikkola1, Paul J. McMillan1, David Hobbs1
1Lund Observatory, Lund University, Department of Astronomy and Theoretical
Physics, Box 43, SE-22100, Lund, Sweden E-mail<EMAIL_ADDRESS>E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Local stellar motions are expected, and have been shown, to include signatures
of the Galaxy’s past dynamical evolution. These are typically divided into the
disc, which shows the dynamical effects of spiral arms and the bar, and the
stellar halo, with structures thought to be debris from past mergers. We use
Gaia Data Release 3 to select large samples of these populations without
limiting them to sources with radial velocities. We apply a penalised maximum
likelihood method to these samples to determine the full 3D velocity
distribution in Cartesian $(U,V,W)$ or spherical $(v_{r},v_{\phi},v_{\theta})$
coordinates. We find that the disc population is dominated by four moving
groups and also detect a new moving group at $(U,V)=(-10,-15)$ km s-1 which we
call MMH-0. For the stellar halo, we isolate the accreted component with cuts
in transverse velocity and the colour-magnitude diagram. In this component we
find several known structures believed to be caused by past mergers,
particularly one around $(v_{r},v_{\phi},v_{\theta})=(-150,-300,-100)$ km s-1
appears more prominent than previously claimed. Furthermore we also identify
two new structures near $(v_{r},v_{\phi},v_{\theta})=(225,25,325)$ km s-1 and
$(0,150,-125)$ km s-1 which we refer to as MMH-1 and MMH-2 respectively. These
results give new insights into local stellar motions and shows the potential
of using samples that are not limited to stars with measured line-of-sight
velocities, which is key to providing large samples of stars, necessary for
future studies.

###### keywords:

methods: statistical - methods: data analysis - Galaxy: structure - Galaxy:
Solar neighbourhood - stars: kinematics and dynamics - Galaxy: kinematics and
dynamics.

††pubyear: 2022††pagerange: New stellar velocity substructures from Gaia DR3
proper motions–10

##  1 Introduction

As our Galaxy evolves, the kinematics of the stars that reside in it are
imprinted by the various external and internal processes that affect it. The
volume near the Sun is no different and contains footprints of possible
dynamical resonances and interactions with nearby dwarf galaxies. For this
reason, untangling the causes of the kinematic structure can give us vital
information about the history of the Milky Way and its interactions with its
nearest neighbours.

The approach to studying nearby kinematic space has bifurcated into the study
of the Galactic disc (e. g. Dehnen 1998; Kushniruk et al. 2017; Antoja et al.
2018; Lucchini et al. 2022; McMillan et al. 2022) and the study of the
Galactic stellar halo (e. g. Koppelman et al. 2019a, b; Koppelman & Helmi
2021a; Lövdal et al. 2022; Ruiz-Lara et al. 2022; Dodd et al. 2022). This
seems only natural with the two stellar components recording different
dynamical processes. The Galactic disc in the Solar neighbourhood shows
evidence of dynamical resonances from the spiral arms and the bar (Antoja et
al., 2010; Trick et al., 2021). The stellar halo, however, records evidence of
mergers between the Milky Way and its neighbours (e.g., Helmi, 2020) which in
the $\Lambda$ cold dark matter ($\Lambda$CDM) model is how galaxies build up
their halos. Finding the causes behind the structures that we can see in the
velocity distributions of the local Galaxy will be an important step towards
fully understanding its complex history.

As noted by Helmi (2020), large samples with accurate kinematics are required
if we are to detect each individual structure. This requirement is now
starting to met by the advent of Gaia (Gaia Collaboration et al., 2016) and
it’s subsequent data releases: DR2 (Gaia Collaboration et al., 2018a), EDR3
(Gaia Collaboration et al., 2021a), and DR3 (Gaia Collaboration et al., 2022).
Thanks to this, we now have over 1.4 billion sources measured with 5D phase-
space coordinates: positions and proper motions. As of DR3, ${\sim}$33 million
sources also have radial velocities, about 2% of all sources, increased from
0.5% in EDR3. Unfortunately as we start to look at more local samples, apply
quality cuts, and pick out specific populations, the number of useful sources
can rapidly decline. As an example consider the local ($\varpi>1/3$ mas)
stellar halo (defined as $v_{\mathrm{T}}>200\ \mathrm{kms}^{-1}$) with good
parallaxes ($\varpi/\sigma_{\varpi}>10$). This sample will contain 503 572
sources with 5 parameters, which is reduced to only 84 784 with measured
radial velocities. By working without radial velocities we are able access
significantly larger datasets and important discoveries can still be reached
as demonstrated by previous works using only proper motions (e.g., Dehnen &
Binney 1998; Dehnen 1998; Antoja et al. 2017; Koppelman & Helmi 2021b; Gaia
Collaboration et al. 2021b; McMillan et al. 2022).

In our previous paper Mikkola et al. (2022) (hereafter referred to as Paper
I), we implemented the penalized maximum-likelihood method of Dehnen (1998) to
infer the 3D velocity distribution of white dwarfs in Gaia EDR3. We apply the
same method here to an extended Solar neighbourhood sample, split into a more
local stellar disc sample and a stellar halo sample. This allows us to view
the velocity distribution in unprecedented velocity resolution. We
particularly focus on the stellar halo which has been shown to host many
structures which are likely due to the merger history of the Galaxy (see e. g.
Dodd et al. 2022 and references therein) with Naidu et al. (2020) even
suggesting the halo could be almost entirely comprised of substructure.

The paper is organised as follows: In Section 2 we describe the data selection
and the quality cuts that have been made to provide a local disc-dominated
sample and a stellar halo sample. We outline some of the differences and new
additions we have made to the method from Paper I in Section 3. Then in
Section 4 we present the velocity distribution for the stellar disc
population. The velocity distribution for the stellar halo is shown in Section
5 where we discuss each of the velocity features we see and compare to
literature, as well as present new features that we identify, and discuss our
findings. We finally summarise with our conclusions in Sections 6.

##  2 Sample selection

Figure 1: Colour-Magnitude diagram of our Solar neighbourhood sample. The colour shows the number density of sources. We exclude bins with N < 5 stars in them. Figure 2: Left: Colour-Magnitude diagram of our sample of stars with $V_{\mathrm{T}}>200\mathrm{kms}^{-1}$ corrected for extinction effects. The two sequences from Gaia Collaboration et al. (2018b) are clearly visible. right: The same plot but showing the regions used to isolate the left and right sequences in red and blue respectively. The color shows the number density of sources and again we exclude bins with N < 5 stars in them. Table 1: The names of the various samples used and the number of sources in them Name | $N$ |  Description
---|---|---
DISC | 1 171 846 |  Stars in the Solar neighbourhood disc defined as $\varpi>5$ mas
DISC_NORTH | 578 368 |  Same as DISC but with $b>0$
DISC_SOUTH | 593 478 |  Same as DISC but with $b<0$
DISC_RVS | 510 478 |  Same as DISC but with measured radial velocities
HALO_BLUE | 239 115 |  Blue halo sequence, see Section 2.
HALO_RED | 194 507 |  Red halo sequence, see Section 2
HALO_RVS | 69 820 |  Our halo sample but with measured radial velocities

The largest sample of stars with proper motions is available from Gaia DR3
through the Gaia Archive111https://gea.esac.esa.int/archive/. Our samples are
selected to be within an extended Solar neighbourhood (Snbh, $d\lesssim$ 3
kpc) in order to meet the criteria of homogenous velocity distribution over
the volume in question as outlined in Section 2.1.1 of Paper I. we perform a
series of quality cuts, most of which are visible in our ADQL queries (which
we provide in appendix A). We select stars with $\varpi/\sigma_{\varpi}>10$ as
a strong cut on parallax uncertainties lets us approximate distance as
$d=1/\varpi$. We also filter ruwe $<1.15$ to ensure good quality astrometric
solutions (Lindegren, 2018) after inspection of the ruwe distribution.
Additionally, Lindegren et al. (2018) explains that for five-parameter
solutions to be accepted at least six separate observations are used (called
visibility_periods_used in the archive). As in Gaia Collaboration et al.
(2018b) we use a stronger filter of at least 8 which removes outliers at the
fainter end. In addition, we use their criteria for astrometric_chi2_al and
astrometric_n_good_obs_al to remove artefacts due to excess astrometric noise.
Finally we use their quality filters for relative flux error on photometry:

$\displaystyle\textrm{{phot\\_g\\_mean\\_flux\\_over\\_error}}>50,$
$\displaystyle\textrm{{phot\\_rp\\_mean\\_flux\\_over\\_error}}>20,\textrm{\qquad
and}$ $\displaystyle\textrm{{phot\\_bp\\_mean\\_flux\\_over\\_error}}>20.$

Beyond this, we also calculate the flux excess in BP and RP, $C^{*}$,
following the procedure of Riello et al. (2021) directly in our query as in
Gaia Collaboration et al. (2021a). The scatter of $C^{*}$ with magnitude,
$\sigma_{C^{*}}(G)$ is fitted with a power law

$\sigma_{C^{*}}(G)=c_{0}+c_{1}G^{m},$ (1)

where $c_{0}=0.0059898$, $c_{1}=8.817481\times 10^{-12}$, and $m=7.618399$.
For each star, the value of $G$ and $C^{*}$ is input and we then select stars
such that $C^{*}<3|\sigma_{C^{*}}(G)|$.

Our disc sample contains stars which have $\varpi>5$ mas ($d\lesssim$ 200 pc).
This sample is further split into north and south Galactic hemisphere samples
with $b>0$ and $b<0$ respectively. The final result of our filters on this
sample can be seen in the colour-magnitude diagram (CMD) in Fig. 1. We also
create a disc sample from the sources with RVs for comparison.

Figure 3: Galactic velocity distributions. The colormap shows the probability
distribution, $f(v)$. The contour lines contain 99, 98, 96, 92, 84, 68, 35,
27, 18, 10, and 1 percent of all the stars, from outside going inward. The
white contours start at 27%. First row: The velocity distribution of the
DISC_FULL sample. In the first column, the location of the first 9 groups of
Antoja et al. (2012) are shown as red crosses. Second row: Same as the first
row, but for the DISC_SOUTH sample. Third row: Velocity distribution of
DISC_NORTH. Figure 4: The conditional probability on $f(U,V)$ on either $U$ or
$V$ for the top and bottom plot respectively. The density is scaled such that
$P(v)^{0.25}$ to reveal low probability structure. Some features discussed in
the text are highlighted.

To generate a sample of halo stars we select stars with a transverse velocity

$v_{\mathrm{T}}=\frac{4.7405}{\varpi}\sqrt{\mu_{l}^{2}+\mu_{b}^{2}}>200\
\mathrm{kms}^{-1},$ (2)

as in Gaia Collaboration et al. (2018b) where it was shown this reveals two
dual well-defined sequences, the standard sequence of the Solar neighbourhood
and a new, bluer, sequence to the left consistent with a more metal-poor
population. This cut on velocity will remove some low velocity halo stars, but
more importantly removes the majority of disc stars. This means our sample
contains slightly fewer stars but is purer, allowing us to identify the
sequences more easily. The smaller sample also makes computations less
demanding. We further require $\varpi>1/3$ mas ($d\lesssim$ 3 kpc) to ensure
we avoid orbits belonging to the bulge.

Our disc sample should be mostly unaffected by extinction due to its proximity
and it will not have any bearing its analysis either. For the halo sample
however, we need to isolate the two sequences in the CMD and therefore it is
important we treat it correctly. To correct our colours and magnitudes, we use
the extinction map of Capitanio et al. (2017)222https://stilism.obspm.fr/ to
determine the reddening for our stars . We pair this with extinction
coefficients for the Gaia bands $G,G_{\mathrm{BP}},G_{\mathrm{RP}}$ from
Sanders & Das (2018) to reduce the effects of extinction. Using the corrected
values, we show the CMD in Fig. 2 which clearly shows the two separate
sequences. Splitting the sample into a left and right sequence sample can be
done reliably by eye and the selection is shown with blue and red shaded
regions respectively.

We list the names of our sample as well as the number of sources contained
within them in Table 1.

##  3 Velocity distribution without RVs

As explained in Section 1, only a small minority of stars with the full
astrometric solution in Gaia DR3 also has radial velocity measurements. In the
following section, we outline how we are able to estimate the full 3D velocity
distribution from stars without measured radial velocities and therefore use
larger samples than would otherwise be possible.

We have used the same maximum penalized-likelihood algorithm from Dehnen
(1998) that was used for Paper I. We summarise the key elements of this
approach here, and a more detailed explanation can be found in our earlier
paper. The method makes use of the projection of transverse velocities on the
sky. The probability distribution of these transverse velocities in a specific
direction, $\hat{\bm{r}}$, we write as $\rho(\bm{q|\hat{\bm{r}}})$ where
$\bm{q}$ is the 2D transverse velocity. This relates to the 3D velocity
distribution, $f(\bm{v})$, as:

$\rho(\bm{q|\hat{\bm{r}}})=\int\mathrm{d}v_{r}f(\bm{v})=\int\mathrm{d}v_{r}f(\bm{p}+v_{r}\hat{\bm{r}}).$
(3)

Here, $\bm{p}$ is the tangential motion of a star projected into 3D.
Tangential motion is not sufficient to determine a true distribution and
instead we must estimate it with a log-likelihood of a model for it.
Numerically, we use the discrete velocity distribution

$f(\bm{v})=e^{\phi(\bm{v})},$ (4)

where $\phi(\bm{v})$ is the logarithm of the probability density, which we
discretize on a 3D grid of $L_{U}\times L_{V}\times L_{W}$ cells with widths
$h_{U}\times h_{V}\times h_{W}$333$U$, $V$ and $W$ are the usual heliocentric
velocity components in the direction towards the Galactic centre, towards
Galactic rotation, and towards the north Galactic pole, respectively. The
resulting function which we use in our maximum penalized-likelihood estimation
is:

$\tilde{\mathscr{Q}}_{\alpha}(\bm{\phi})=N^{-1}\sum_{k}\ln\left[\sum_{\bm{l}}e^{\phi_{\bm{l}}}K(k|\bm{l})\right]-\sum_{\bm{l}}e^{\phi_{\bm{l}}}\\\
-\frac{1}{2}\alpha
h_{U}h_{V}h_{W}\sum_{\bm{l}}\left(\sum_{\bm{n}}\phi_{\bm{n}}\Xi_{\bm{n}\bm{l}}\right)^{2}.$
(5)

The first term is the sum of the probability distribution function with $N$
being the sample size. For each star, $k$, the length in velocity space of the
line formed by its tangential velocity and all possible radial velocities
through a cell, $\bm{l}$, is $K(k|\bm{l})$. The second term is a normalising
term and the third is the penalizing term where
$\left(\sum_{\bm{n}}\phi_{\bm{n}}\Xi_{\bm{n}\bm{l}}\right)$ is a numerical
approximation for the second derivative of $\phi(\bm{v})$ for a given cell.
This term therefore penalises unsmooth solutions and is scaled by the
smoothing parameter, $\alpha$.

We use a similar method for determining the optimal value of $\alpha$ as in
Paper I. That is, we compare many estimations of an equally sized RVS sample
using different $\alpha$. The comparable RVS sample to the full disc sample is
only half as large as the sample without $v_{r}$ so we have to upscale our
sample. We create a copy of the RVS sample, DISC_RVS where the Galactic
positions are sampled randomly from the original. The velocities are taken
from the RVS, with each proper motion and radial velocity resampled from a
multivariate Gaussian with the measured values as mean, and with uncertainties
and correlation coefficients in the covariance matrix. With a randomly
selected Galactic position this then transforms to Galactic velocities,
$U,V,W$, and this copy is then projected back into on-sky motions with the
radial velocities discarded.

For the disc sample, the upscaling only needs to be done once as the RVS is
half the size of the full dataset. We then compare the estimated
$\tilde{f}^{\alpha}_{\bm{v}}$ with the real RVS $f_{\bm{v}}$ and select the
$\alpha$ that gives the smallest integrated square error (ISE)

$D(\tilde{f}^{\alpha}_{\bm{v}},f_{\bm{v}})=\int
d^{3}~\bm{v}(\tilde{f}^{\alpha}_{\bm{v}}-f_{\bm{v}})^{2}.$ (6)

This gives for the full disc sample an optimal smoothing of $\alpha=10^{-11}$
for a $\bm{n}=[304,304,192]$ grid over the ranges $U\in[-150,150]\
\mathrm{kms}^{-1}$, $V\in[-200,100]\ \mathrm{kms}^{-1}$, $W\in[-100,100]\
\mathrm{kms}^{-1}$, corresponding to a resolution of ${\sim}1\
\mathrm{kms}^{-1}$. The same setup is used for the north and south Galactic
hemisphere samples to ensure the differences are not by construction. This
will however result in a slightly under-smoothed distribution for the smaller
samples.

For the stellar halo, the HALO_RVS sample is about 3 times smaller than the
two samples HALO_BLUE and HALO_RED and we upscale it to three times its
original size. We also no longer use heliocentric Cartesian velocity
coordinates, $U,V,W$ but instead use Galactocentric spherical velocities,
$v_{r},v_{\phi},v_{\theta}$ defined such that for a star in the galactic plane
at the position of the Sun, $v_{\phi}$ and $v_{\theta}$ are in the same
direction as V and W, respectively, to make comparisons easier. To make the
transformation, we have assumed the Sun’s position to be $(R,z)=(8122,20.8)$
pc (GRAVITY Collaboration et al., 2018; Bennett & Bovy, 2019) and its velocity
$(U,V,W)=(12.9,245.6,7.78)$ km s-1 (Drimmel & Poggio, 2018) with respect to
the Galactic centre. On a grid of $\bm{n}=[240,240,240]$ with
$v_{r},v_{\phi}$, and $v_{\theta}$ all in the range $[-600,600]\
\mathrm{kms}^{-1}$ corresponding to a resolution of $5\ \mathrm{kms}^{-1}$,
the value of $\alpha$ that minimises the ISE is $4.64\times 10^{-13}$ and is
used for both halo samples.

One difference between our approach when handling the disc sample and when
handling the halo sample is how the velocity dispersion, $\bm{\sigma}$, and
average velocity, $\langle\bm{v}\rangle$, are determined. For the former
sample we can determine $\bm{\sigma}$ and $\langle\bm{v}\rangle$ of the sample
directly, following the procedure of Dehnen & Binney (1998) (as in Paper I),
but for the latter case we use spherical coordinates and it becomes
unnecessarily complex. Instead, we determine $\langle\bm{v}\rangle$ and
$\bm{\sigma}$ for the subset of the sample that has measured radial
velocities. However, since $\bm{\sigma}$ and $\langle\bm{v}\rangle$ are only
used for the initial guess of $\phi(\bm{v})$ and as scaling factors for
determining $\alpha$, we are free to find $\sigma$ however we wish without
significant consequences for the analysis.

##  4 The stellar disc

Table 2: Average velocities of moving groups mentioned in the text, based on reported values in Kushniruk et al. (2017). $(\langle U\rangle,\langle V\rangle)$ | Moving group
---|---
$[\mathrm{km\ s}^{-1}]$ |
$(8,3)$ | Sirius
$(-7,-8)$ | Coma Berenices
$(-37,-17)$ | Hyades
$(-15,-22)$ | Pleiades
$(-35,-48)$ | Hercules
$(46,-24)$ | Dehnen98
$(31,-24)$ | wolf630
$(-86,-46)$ |  $\epsilon$Ind
$(-95,-8)$ | bobylev16

The velocity distribution estimated for the Solar neighbourhood sample is seen
in Fig. 3. The first row shows distribution for the complete sample and the
second and third rows show the south and north Galactic hemisphere samples
respectively. We find that the distribution is mostly dominated by the common
features: Sirius, Coma Berenices, Hyades, Pleiades, Hercules. We can also
identify Dehnen98 and Wolf 630 to an extent. The first four of the major
groups occupy a region that contains roughly 35% of all the stars in the
sample, shown by the white contour lines. We can see the incomplete vertical
face-mixing of Coma Berenices (Quillen et al., 2018; Monari et al., 2018;
Bernet et al., 2022) as it is a much stronger feature in the southern Galactic
hemisphere. A curious feature that appears more strongly in the northern
hemisphere is the rather strong overdensity between Pleiades and the expected
position of Coma Berenices at roughly $(U,V)=-10,-15\ \mathrm{kms}^{-1}$. This
feature is clearly separate from Pleiades and we find no match for it in the
list of moving groups in Antoja et al. (2017), Kushniruk et al. (2017), or
Lucchini et al. (2022). The close proximity to Pleiades suggests that this
feature is now visible thanks to the improved velocity resolution. When
looking for this feature in the 6D sample, we find that pleiades appears to
stretch toward this region but there is no separate feature. For this reason
we tentatively names this a new velocity feature which we call MMH-0. Overall
we find limited detailed substructure in the direct velocity distribution
which shows the dominance of the major moving groups in this space.

###  4.1 Conditional $f(\bm{v})$ of the stellar disc

To unravel low-level structure that the representations of the velocity
distribution of Fig. 3 may have missed, in Fig. 4 we renormalize the plots so
that, rather than showing the full 2D probability density of $U$ and $V$, we
show the _conditional_ probabilities of $V$ or $U$ for each $U$ or $V$,
respectively. That is, the colour represents the probability of the star
having a specific $V$ given that it has certain $U$ velocity (or vice versa).

In addition to the structure we have seen above, we can see in $P(V|U)$ around
$(U,V)=-100,-50\ \mathrm{km\ s}^{-1}$ a structure that matches well with
estimates of $\epsilon$Ind (e.g., Antoja et al. 2012; Bobylev & Bajkova 2016;
Kushniruk et al. 2017). Above it, closer to $V=-10\ \mathrm{km\ s}^{-1}$ is
another feature that matches to a group identified by Bobylev & Bajkova (2016)
and Kushniruk et al. (2017). There are also features with positive $U$ sitting
at $V\approx-30$ km s-1 with $U=50$ km s-1 and $U=100$ km s-1, the first of
which is the Dehnen98 group from Antoja et al. (2012), which itself is from
Dehnen 1998. The second group is likely Antoja12 (see Kushniruk et al. 2017
and references therein). This demonstrates the strength of plotting
conditional probabilities for the inferred velocity distributions to gain
insight into low-level structures.

##  5 The local stellar halo

Figure 5: Velocity distributions in spherical coordinates. The colormap shows
the square root of the probability distribution, $\sqrt{f(v)}$, to enhance
fainter structure. The contour levels are the same as in Fig. 3. First row:
The velocity distribution of the HALO_BLUE sample. Five distinct features
thought belonging to the region occupied by GSE are labelled. Second row: Same
as the first row, but overlaid with rough expected positions of reported
substructures from literature in similar style to Naidu et al. (2020) and
Mardini et al. (2022) but in velocity space. Third row: Velocity distribution
of HALO_RED, with very little substructure. Because we are interested in the
substructure found in HALO_BLUE we focus exclusively on this sample in
subsequent figures. Figure 6: The conditional probability of the three
different spaces in each column and velocity space as Fig. 5 for the HALO_BLUE
sample. Top row is conditional probability on the x-axis coordinate and bottom
row on the y-axis coordinate. We show specifically $P(v)^{0.25}$ to reveal low
probability structure even further.

We present the velocity distribution of the stellar halo in the planes
$(v_{r},v_{\phi})$, $(v_{r},v_{\theta})$, and $(v_{\phi},v_{\theta})$ in Fig.
5 estimated through our penalized maximum-likelihood method. The figure very
clearly shows the separation into a blue accreted population and a red in-situ
halo or hot thick disc. The blue accreted population occupies phase space
almost symmetrically around $v_{\phi}=0$ and does not continue smoothly from
the disc population, with multiple substructures. The red, in-situ halo or hot
thick disc on the other hand has very little substructure and is mostly a
continuation of the disc phase space distribution.

Since our main interest is in the substructure of the accreted halo,
throughout the rest of this paper we only look at the HALO_BLUE sample.

For this sample we can see that there are a plethora of features. We overlay
the distribution with known features from literature, revealing which of these
we do or do not see from this perspective of the population. The extent and
distribution of these shapes has been adapted from visual inspection of the
features as seen in Koppelman et al. (2019a); Naidu et al. (2020); Dodd et al.
(2022); Oria et al. (2022); Tenachi et al. (2022) and are not exact but simply
to be used as a guiding tool. We will now go over and discuss the features we
can see in our distributions.

###  5.1 Gaia-Sausage-Enceladus

The strongest feature across all velocity space is the GSE and this feature is
not continuous but rather appears to be a composition of multiple different
features which we have labelled in the top-left of Fig. 5 as G1-5 in
$(v_{r},v_{\phi})$. Since our cut of $v_{\mathrm{T}}>200\ \mathrm{km\ s}^{-1}$
is rather generous, we expect there to be some contamination between the
samples. Thus, G2 and G4 are likely contaminants from the two features seen in
the HALO_RED sample in the third row of Fig. 5. G1 and G5 appear in the space
of what is typically associated with the GSE (e. g. Koppelman et al. 2019a;
Feuillet et al. 2021; Dodd et al. 2022) while the central group, G3, appears
likely to be the L-RL3 group in Dodd et al. (2022) who identifies it with
Cluster 3 in Lövdal et al. (2022). The GSE appears slightly asymmetric on
either side of $v_{\phi}=0$ here because of our cut on tangential velocity
removing some of its lower $v_{\phi}$ members and including some contaminants
from the disc.

###  5.2 Sequoia & Antaeus

Some noticeable features are the ‘horns’ sticking out at the bottom in
$(v_{r},v_{\phi})$ at $v_{\phi}=-250\ \mathrm{km\ s}^{-1}$ and around
$v_{r}=\pm 200\ \mathrm{km\ s}^{-1}$ and to the left in
$(v_{\phi},v_{\theta})$ around $v_{\theta}=\pm 150\ \mathrm{km\ s}^{-1}$.
These features are in the location associated with Sequoia (Myeong et al.,
2019). In Naidu et al. (2020) this region is mixed with groups Arjuna and
I’itoi which are distinguished from one another by their metallicities. In
Ruiz-Lara et al. (2022), three separate clusters are also identified as
belonging to Sequoia but share similar metallicity. Instead the kinematics
distinguish these three structures and they attribute one of the clusters to
the metal-poor end of the GSE based on its kinematics. Since we do not have
metallicity measurements here we will refer to the dynamical space occupied by
all of these features simply as Sequoia.

In the same space we can also see the Antaeus group from Oria et al. (2022)
which shares many of the same attributes as Sequoia. It is not clear to what
extent these features are separate but Oria et al. (2022) claims that the low
$J_{z}$ and position in the disc plane of Antaeus are unique.

###  5.3 Helmi streams

Two of the most prominent substructures that appears in all the distributions
are the Helmi streams (Helmi et al., 1999). They are also among the first
identified substructures. An updated view of the streams lets us narrow them
down in velocity space (Koppelman et al. 2019a, b; Koppelman & Helmi 2021a).
This feature is particularly strong and is bimodal in $v_{\theta}$, with the
lower $v_{\theta}$ group being far more represented as expected (e.g.,
Koppelman et al. 2019b). At slightly larger $v_{\phi}$ we would expect to find
Aleph, reported first in Naidu et al. (2020), but in our sample it appears to
be absent. Similarly in Lövdal et al. (2022) the absence of Aleph was noted
and argued to be caused by the velocity cut removing almost all of its stars
from the sample. The same is likely the cause behind its absence in our
sample.

###  5.4 Thamnos

The structure Thamnos was identified by Koppelman et al. (2019b) using the
Gaia DR2 RVS sample supplemented with line-of-sight velocities and abundances
from RAVE (Kunder et al., 2017), APOGEE (Abolfathi et al., 2018), and LAMOST
(Cui et al., 2012). Their sample was limited to 3 kpc like halo samples. In
our distributions there is no distinct separate feature corresponding to
Thamnos, but in both $(v_{r},v_{\phi})$ and $(v_{\phi},v_{\theta})$ the
distribution’s densest parts extends to slightly lower $v_{\phi}$ than for the
GSE structure, which may be due to the presence of Thamnos. It is suggested by
Naidu et al. (2020) that Thamnos may be more discernible at larger distances
where the GSE and disc-like stars contribute less to the distribution.

###  5.5 Other structures

In addition to the structures mentioned above we find several others in these
projections. One of these can be associated with a known velocity structure,
ED-2 from Dodd et al. (2022) and is marked with a cerise-shaded region. This
feature is close to Sequoia at $(v_{r},v_{\phi})=(-150,-300)\ \textrm{km
s}^{-1}$ and presents in $(v_{\phi},v_{\theta})=(-300,-100)\ \textrm{km
s}^{-1}$. To verify if both velocity representations are ED-2, we investigate
the full 3D velocity structure to find that the features overlap and are one
and the same. In Dodd et al. (2022) there is limited metallicity information
for the cluster and it has average metallicity
$\langle\mathrm{[Fe/H]}\rangle=-2.05$, which would match the metal poor part
of Sequoia. In Section 5.2 we discussed the different parts of Sequoia which
we are not able to distinguish. It is possible that this feature is yet
another part of the same velocity feature.

The cluster is assigned 33 members in Dodd et al. (2022) out of a sample of 72
274 stars (or ${\sim}$0.05%). We can look at the fraction of the probability
density that occupies the region. We define the region around the group with
$v_{r}\in[-175,-100]$ km s-1, $v_{\phi}\in[-325,-275]$ km s-1, and
$v_{\theta}\in[-100,-150]$ km s-1 and find that ${\sim}$0.074% of the sample
lies there, corresponding to ${\sim}$180 stars in the sample. This suggests
that the very dense feature is slightly more prominent than previous believed.

Above the high-$v_{\theta}$ part of the Helmi stream, at $v_{\theta}\approx
300\ \textrm{km s}^{-1}$, we find a new feature split across two different
values of $v_{r}$, one around $(v_{r},v_{\theta})=(-250,300)\ \textrm{km
s}^{-1}$ and the other at $(v_{r},v_{\theta})=(200,350)\ \textrm{km s}^{-1}$.
At such a large $v_{\theta}$ we also find a new structure in
$(v_{\phi},v_{\theta})$ at $(25,300)\ \textrm{km s}^{-1}$. We refer to these
groups as a single feature which we call MMH-1. In the full 3D probability
distribution the positions of MMH-1 overlap and we consider them the same
feature.

Lastly, at $(v_{\phi},v_{\theta})=(150,-100)\ \textrm{km s}^{-1}$, another
feature can be seen. At such values of $v_{\phi}$ and $v_{\theta}$ it is
difficult to discern any stronger feature in the spaces of $(v_{r},v_{\phi})$
and $(v_{r}v_{\theta})$ as the region is crowded particularly in
$(v_{r},v_{\phi})$ where it would lie close to the cutoff caused by our
tangential velocity limit. We refer to this feature as MMH-2.

###  5.6 Conditional $f(\bm{v})$ of the local halo

In much the same way as we did in Fig. 5, we again use conditional
probabilities to illustrate our halo velocity maps with respect to one of the
two velocity dimensions to investigate faint structure that otherwise may not
be visible. We show these conditional probability maps in Fig. 6 which reveal
more of the surrounding velocity structure with certain features become
strikingly visible. For example the two-pronged structure around Sequioa,
Antaeus, and ED-2 is much more readily apparent in $P(v_{r}|v_{\phi})$ and
$P(v_{\theta}|v_{\phi})$ than before. The separation of the GSE from the hot
thick disc is apparent in both $P(v_{r}|v_{\phi})$ and $P(v_{\phi}|v_{r})$,
and it is located primarily around two features around $v_{\phi}\pm[150,300]$
km s-1. The region occupied by Thamnos is now also readily apparent where it
was not before.

Our two novel groups, MMH-1 and MMH-2, appear in these figures as well. The
first, MMH-1, which has a large $v_{\theta}$ of around 300 km s-1 appears in
the conditional probabilities $P(v_{r}|v_{\theta})$ and
$P(v_{\phi}|v_{\theta})$ corresponding to the bottom row, middle and last
columns. It also stretches up towards much larger extents in $v_{\theta}$, but
it is unclear if this structure is physical given how far out in the edges of
the distribution it lies. Similarly there is a symmetric feature at large
negative $v_{\theta}$ near -400 km s-1 which, if real, could relate to MMH-1.
The group MMH-2 appears strongly in $(v_{\phi},v_{\theta})$ for both
$P(v_{\phi}|v_{\theta})$ and $P(v_{\theta}|v_{\phi})$.

A feature we have not discussed previously that is present in
$P(v_{\theta}|v_{r})$ and $P(v_{r}|v_{\theta})$ is the sloped feature around
$(v_{r},v_{\theta})=(0,-150)$ km s-1 which we associate with MMH-2. This
feature is also visible in Fig. 5 and does not appear in the other papers we
have reviewed and is very difficult to find in the other spaces, suggesting it
lies at rather low $v_{\phi}$ as this would place it close to the densest
parts of the distribution, thus obscuring it from detection. We have confirmed
this by limiting the $(v_{r},v_{\theta})$-space to separate bins of
$v_{\phi}$, which reveals that the feature only appears between
$v_{\phi}\in[0,100]$ km s-1 (this can be seen in our binned velocity
distributions in Appendix C). This means the feature is most likely the
representation of MMH-2 in this space.

###  5.7 Action space distribution

Figure 7: Action distribution associated with our HALO_BLUE sample. The
probability density is in units of
$\mathrm{kpc}^{-2}\;\mathrm{km}^{-2}\,\mathrm{s}^{2}$. While there are strong
selection effects that shape these distributions, we are still able to see
some of the substructure in this projection of the data, most strikingly the
Helmi stream in the lower panel near
$(J_{\phi},J_{\theta})=(1300,1200)\;\mathrm{kpc}\,\mathrm{km}\,\mathrm{s}^{-1}$.

A very common method of searching for accreted substructure in the Milky Way
halo is to work in terms of integrals of motion (e.g., Helmi et al., 1999). In
particular, the orbital actions have become commonly used since a convenient
approximation within galaxy potentials became available (Binney, 2012). When
stars in a sample are distributed over a large volume in the Milky Way it is
essential to use techniques like these to identified substructure accreted
long ago because, while we can expect stars accreted together to have very
similar integrals of motion, they will have very different velocities if they
are in very different parts of the Galaxy.

We have very deliberately limited our sample to a relatively local volume
within the stellar halo so that stars on similar orbits have similar
velocities. Nonetheless, to make consistent comparison with other studies
possible, it is still valuable to determine the distribution of the stars in
terms of their orbital actions. We use the agama software package (Vasiliev,
2019) to determine the actions in the Milky Way gravitational potential from
McMillan (2017) rescaled such that the Sun’s position and velocity are
consistent with the values given in Section 3. Fig. 7 shows this action
distribution in terms of $J_{r}$ or $J_{\theta}$ against
$J_{\phi}$,444$J_{\phi}$ is the conserved component of angular momentum, while
$J_{r}$ characterises the radial oscillation and $J_{\theta}$ that of
oscillation in the $\theta$ direction. where we have approximated that stars
are at the position of the Sun. Changing this assumed position within our
survey volume does not make any important qualitative difference.

The distributions in Fig. 7 are clearly shaped by the selection effects that
apply to them. At high $|J_{\phi}|$ there is a lower limit for $J_{r}$ below
which there is no orbits that pass the Sun’s position. At lower $|J_{\phi}|$,
there is a minimum $J_{r}$ that reaches the Sun’s position for $J_{\theta}=0$,
and this is where we find the bright maximum in the $(J_{\phi},J_{r})$
distribution that runs from approximately $(-1500,0)$ to
$(0,600)\;\mathrm{kpc}\,\mathrm{km}\,\mathrm{s}^{-1}$, which is a selection
effect. For $J_{\theta}>0$ there are orbits that reach the solar position at
these $J_{\phi}$ with lower $J_{r}$, so the density does not fall to zero
below this point. There is an upper envelope for $J_{\theta}$ at a given
$J_{\phi}$, above which the density becomes very low. This is a consequence of
the requirement that the orbit reaches the Sun’s position and is bound to the
Galaxy.

Nonetheless, there are features of the velocity distribution that stand out in
these plots too. The GSE stands out as a strikingly high density of stars over
a large range in $J_{r}$ around ${J_{\phi}}={0}$. ED-2 is clearly seen in both
panels at
$(J_{r},J_{\phi},J_{\theta})=(600,-2450,90)\;\mathrm{kpc}\,\mathrm{km}\,\mathrm{s}^{-1}$;
the Helmi Stream is clearly visible in the $(J_{\phi},J_{\theta})$ plane
around $(1300,1200)\;\mathrm{kpc}\,\mathrm{km}\,\mathrm{s}^{-1}$ and expected
at a $J_{r}$ of 80 kpc km s-1; the Sequoia group produces an overdensity that
can be seen around
$(J_{\phi},J_{r})=(-2000,1000)\;\mathrm{kpc}\,\mathrm{km}\,\mathrm{s}^{-1}$.
Otherwise, for our sample the substructure is substantially clearer in the
velocity distribution than in these action distributions.

We note, in the interests of finding members of our newly discovered
substructure in future study of samples beyond the Solar neighbourhood, that
in our assumed potential these have actions $(J_{r},J_{\phi},J_{\theta})$
around $(1450,200,2300)$ and
$(70,1200,200)\;\mathrm{kpc}\,\mathrm{km}\,\mathrm{s}^{-1}$ for MMH-1 and
MMH-2 respectively.

Finally, we can use angle-action modelling to justify an assumption underlying
our approach: that it is reasonable to approximate the velocity distribution
of the stellar halo as independent of position in our sample volume. To
provide a realistic example, we use the Torus Code (Binney & McMillan, 2016)
we sample the points from the phase-mixed orbit corresponding to the
approximate actions of the Helmi and ED-2 streams. For the Helmi stream we
focus only on component at $v_{\theta}<0$ (by symmetry the equivalent at
$v_{\theta}>0$ will have the same spread in velocity), while for ED-2 we also
limit it to $v_{r}<0$ to match the major component we observe. For the points
on the orbital torus sampled within our survey volume, we have a dispersion in
$v_{\phi}$ of $12\,\mathrm{km}\,\mathrm{s}^{-1}$ and in $v_{\theta}$ of
$17\,\mathrm{km}\,\mathrm{s}^{-1}$ for the Helmi Stream, and of
$(50,40,30)\,\mathrm{km}\,\mathrm{s}^{-1}$ in $(v_{r},v_{\phi},v_{\theta})$
for ED-2. This spread is comparable to that seen in for these groups in Fig.
5, which is quite small on the scale of the velocity distribution that we are
studying, and clearly does not prevent us from finding substructure. We note
also that this is likely to be an overestimate of the associated dispersion,
because the real sample has a smaller spread in position, being preferentially
near the Sun, and this smaller spread in position, for a given orbital torus,
corresponds to a smaller spread in velocity.

##  6 Conclusions

We use DR3 astrometry data without radial velocities, giving us access to a
significantly larger catalogue of stars. With the penalized maximum likelihood
algorithm implemented in Paper I we can then infer full 3D velocity
distributions to investigate the Solar neighbourhood for substructure. We
analyse the extended Solar neighbourhood in two separate stellar components:
the Galactic disc ($d<200$ pc) and the stellar halo ($d<3$ kpc). The disc is
also split into a north and southern hemisphere based on Galactic latitude as
in Monari et al. (2018) and we find that the overall velocity distribution is
dominated by the four major moving groups; Sirius, Coma Berenices, Hyades, and
Pleiades, with 35% of stars lying in and around them. However, we find some
degree of asymmetry with Galactic latitude with Coma Berenices being most
prominent in the southern hemisphere in agreement with previous results
(Quillen et al., 2018; Monari et al., 2018; Bernet et al., 2022). We also
identify a new structure at $(U,V)=-10,-15\ \mathrm{kms}^{-1}$ which does not
align with any known moving groups.

For the local stellar halo, we use the same approach as Gaia Collaboration et
al. (2018b) to reveal a double main sequence for stars with
$v_{\mathrm{T}}>-200\ \textrm{km s}^{-1}$ which we split into in an ‘in-situ’
and ‘accreted’ population to the right and left in the CMD respectively. These
samples are then used to infer the velocity distributions in spherical
Galactocentric coordinates, $v_{r},v_{\phi}$, and $v_{\theta}$. We see that we
can reliably make out several of the more well-known features of the stellar
halo in the ‘accreted’ sample: GSE, Sequoia, Helmi streams, and Thamnos are
all visible in our sample. We also find three additional structures, the first
of which is identified already as ED-2 in Dodd et al. (2022). We then have a
new structure, MMH-1, appear at large $v_{\theta}$, split into two locations
at $(v_{r},v_{\theta})=(-250,300)\ \textrm{km s}^{-1}$ and
$(v_{r},v_{\theta})=(200,350)\ \textrm{km s}^{-1}$. It also appears at
$(v_{\phi},v_{\theta})=(25,300)\ \textrm{km s}^{-1}$. By inspection of the
full 3D velocity space we confirm this as one feature with velocities
$(v_{r},v_{\phi},v_{\theta})=(\pm 225,25,325)$ km s-1. Lastly we also have the
new feature MMH-2 at $(v_{\phi},v_{\theta})=(150,-100)\ \textrm{km s}^{-1}$
which we trace into $v_{r}$, which gives it velocities
$(v_{r},v_{\phi},v_{\theta})=(0,150,-125)$ km s-1. These velocity
distributions gives us the appearance of the stellar halo at ‘face value’,
which provides a clear idea of what structure can be expected there and to
what extent.

In addition to this, we investigated the conditional velocity distributions
which provides further support for the existence of these structures and their
extent. This is where we also find a match for MMH-2 in $v_{r}$, which was not
as readily apparent in the standard distributions.

Furthermore, we transform our velocity distributions into action space
distributions and identify the location of several of the previous
substructures there as well. This further demonstrates the possibilities of
our approach and lets us connect between velocity space and orbital space.

A promising benefit of these velocity distributions is the ability to produce
stellar candidate lists of velocity features. While the method obviously does
not allow us to determine the $v_{r}$ of individual stars, it does allow for
determining a probability of having the necessary $v_{r}$ to belong to a
certain velocity feature. Consider again our 3D grid of velocities from
$\bm{v}_{m}in$ to $\bm{v}_{m}ax$ with some grid spacing. In this box, a star
with unknown $v_{r}$ forms a line $\bm{p}+v_{r}\bm{\hat{r}}$. If we make the
reasonable assumption that the true $v_{r}$ lies within our box, then we can
determine the integral of the probability distribution where the line crosses
the velocity feature. This is then normalised against the probability
distribution along the entire line:

$P_{k}({\bm{l^{*}}})=\frac{{\sum}\limits_{\bm{l^{*}}}e^{\phi_{\bm{l^{*}}}}K(k|{\bm{l^{*}}})}{{\sum}\limits_{\bm{l}}e^{\phi_{\bm{l}}}K(k|\bm{l})},$
(7)

where $\bm{l^{*}}$ are the cells assigned to a specific velocity feature
(where cells $\bm{l}$ are defined in Section 3). Once this is determined, a
candidate list is created by requiring that $P_{k}(\bm{l^{*}})$ be greater
than some threshold probability.

This work demonstrates what can be achieved without needing to rely on the
full 6D phase-space information. During the era of Gaia, there will be more
sources with astrometry alone than with added radial velocities due to the
inherent differences in the methods by which the measurements are obtained.
Currently around 2% of the data has radial velocities and with spectroscopic
follow-up this is likely to increase, but not match the amount of pure
astrometric sources. In the next era with a successor mission in the infrared
(Hobbs 2022, submitted)555Proceedings of the XXXI IAU General Assembly, to be
published in Cambridge University Press the amount of radial velocities could
increase significantly. This will mean that enough sources will be available
that more discoveries can be made directly with 6D data and in our case it
could allow us to determine the $\alpha$-parameter more finely, but the proper
motions will remain more numerous and so methods such as these will be
pivotal. Additionally, as the number of 5D sources increase, the velocity
resolution that can be achieved in the estimation of the velocity distribution
will also increase.

## Acknowledgements

This work has made use of data from the European Space Agency (ESA) mission
Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing
and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. This work made use of
Astropy:666http://www.astropy.org a community-developed core Python package
and an ecosystem of tools and resources for astronomy (Astropy Collaboration
et al., 2013, 2018, 2022).

We thank members of Lund Observatory for helpful comments and ideas.
Computations for this study were performed on equipment funded by a grant from
the Royal Physiographic Society in Lund. PM is supported by research project
grants from the Swedish Research Council (Vetenskaprådet Reg: 20170-03721 and
2021-04153). DH and PM gratefully acknowledge support from the Swedish
National Space Agency (SNSA Dnr 74/14 and SNSA Dnr 64/17).

## Data availability

All data analysed in this paper are publicly available from the Gaia archive
(http://gea.esac.esa.int/archive/). The 3D probability distributions used in
Figures 3 and 5 are available upon request to the corresponding author.

## References

  * Abolfathi et al. (2018) Abolfathi B., et al., 2018, ApJS, 235, 42
  * Antoja et al. (2010) Antoja T., Figueras F., Torra J., Valenzuela O., Pichardo B., 2010, in , Vol. 4, Lecture Notes and Essays in Astrophysics. Tórculo Press, pp 13–31
  * Antoja et al. (2012) Antoja T., et al., 2012, MNRAS, 426, L1
  * Antoja et al. (2017) Antoja T., de Bruijne J., Figueras F., Mor R., Prusti T., Roca-Fàbrega S., 2017, A&A, 602, L13
  * Antoja et al. (2018) Antoja T., et al., 2018, Nature, 561, 360
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2022) Astropy Collaboration et al., 2022, ApJ, 935, 167
  * Bennett & Bovy (2019) Bennett M., Bovy J., 2019, MNRAS, 482, 1417
  * Bernet et al. (2022) Bernet M., Ramos P., Antoja T., Famaey B., Monari G., Al Kazwini H., Romero-Gómez M., 2022, A&A, 667, A116
  * Binney (2012) Binney J., 2012, MNRAS, 426, 1324
  * Binney & McMillan (2016) Binney J., McMillan P. J., 2016, MNRAS, 456, 1982
  * Bobylev & Bajkova (2016) Bobylev V. V., Bajkova A. T., 2016, Astronomy Letters, 42, 90
  * Capitanio et al. (2017) Capitanio L., Lallement R., Vergely J. L., Elyajouri M., Monreal-Ibero A., 2017, A&A, 606, A65
  * Cui et al. (2012) Cui X.-Q., et al., 2012, Research in Astronomy and Astrophysics, 12, 1197
  * Dehnen (1998) Dehnen W., 1998, AJ, 115, 2384
  * Dehnen & Binney (1998) Dehnen W., Binney J. J., 1998, MNRAS, 298, 387
  * Dodd et al. (2022) Dodd E., Callingham T. M., Helmi A., Matsuno T., Ruiz-Lara T., Balbinot E., Lovdal S., 2022, arXiv e-prints, p. arXiv:2206.11248
  * Drimmel & Poggio (2018) Drimmel R., Poggio E., 2018, Research Notes of the American Astronomical Society, 2, 210
  * Feuillet et al. (2021) Feuillet D. K., Sahlholdt C. L., Feltzing S., Casagrande L., 2021, MNRAS, 508, 1489
  * GRAVITY Collaboration et al. (2018) GRAVITY Collaboration et al., 2018, A&A, 615, L15
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2018a) Gaia Collaboration et al., 2018a, A&A, 616, A1
  * Gaia Collaboration et al. (2018b) Gaia Collaboration et al., 2018b, A&A, 616, A10
  * Gaia Collaboration et al. (2021a) Gaia Collaboration et al., 2021a, A&A, 649, A1
  * Gaia Collaboration et al. (2021b) Gaia Collaboration et al., 2021b, A&A, 649, A8
  * Gaia Collaboration et al. (2022) Gaia Collaboration Vallenari A., Brown A. G. A., Prusti T., 2022, A&A
  * Helmi (2020) Helmi A., 2020, ARA&A, 58, 205
  * Helmi et al. (1999) Helmi A., White S. D. M., de Zeeuw P. T., Zhao H., 1999, Nature, 402, 53
  * Koppelman & Helmi (2021a) Koppelman H. H., Helmi A., 2021a, A&A, 645, A69
  * Koppelman & Helmi (2021b) Koppelman H. H., Helmi A., 2021b, A&A, 649, A136
  * Koppelman et al. (2019a) Koppelman H. H., Helmi A., Massari D., Roelenga S., Bastian U., 2019a, A&A, 625, A5
  * Koppelman et al. (2019b) Koppelman H. H., Helmi A., Massari D., Price-Whelan A. M., Starkenburg T. K., 2019b, A&A, 631, L9
  * Kunder et al. (2017) Kunder A., et al., 2017, AJ, 153, 75
  * Kushniruk et al. (2017) Kushniruk I., Schirmer T., Bensby T., 2017, A&A, 608, A73
  * Lindegren (2018) Lindegren L., 2018, Re-normalising the astrometric chi-square in Gaia DR2, GAIA-C3-TN-LU-LL-124, http://www.rssd.esa.int/doc_fetch.php?id=3757412
  * Lindegren et al. (2018) Lindegren L., et al., 2018, A&A, 616, A2
  * Lövdal et al. (2022) Lövdal S. S., Ruiz-Lara T., Koppelman H. H., Matsuno T., Dodd E., Helmi A., 2022, A&A, 665, A57
  * Lucchini et al. (2022) Lucchini S., Pellett E., D’Onghia E., Aguerri J. A. L., 2022, MNRAS,
  * Mardini et al. (2022) Mardini M. K., Frebel A., Chiti A., Meiron Y., Brauer K. V., Ou X., 2022, ApJ, 936, 78
  * McMillan (2017) McMillan P. J., 2017, MNRAS, 465, 76
  * McMillan et al. (2022) McMillan P. J., et al., 2022, MNRAS, 516, 4988
  * Mikkola et al. (2022) Mikkola D., McMillan P. J., Hobbs D., Wimarsson J., 2022, MNRAS, 512, 6201
  * Monari et al. (2018) Monari G., et al., 2018, Research Notes of the American Astronomical Society, 2, 32
  * Myeong et al. (2019) Myeong G. C., Vasiliev E., Iorio G., Evans N. W., Belokurov V., 2019, MNRAS, 488, 1235
  * Naidu et al. (2020) Naidu R. P., Conroy C., Bonaca A., Johnson B. D., Ting Y.-S., Caldwell N., Zaritsky D., Cargile P. A., 2020, ApJ, 901, 48
  * Oria et al. (2022) Oria P.-A., Tenachi W., Ibata R., Famaey B., Yuan Z., Arentsen A., Martin N., Viswanathan A., 2022, ApJ, 936, L3
  * Quillen et al. (2018) Quillen A. C., et al., 2018, MNRAS, 478, 228
  * Riello et al. (2021) Riello M., et al., 2021, A&A, 649, A3
  * Ruiz-Lara et al. (2022) Ruiz-Lara T., Matsuno T., Lövdal S. S., Helmi A., Dodd E., Koppelman H. H., 2022, A&A, 665, A58
  * Sanders & Das (2018) Sanders J. L., Das P., 2018, MNRAS, 481, 4093
  * Tenachi et al. (2022) Tenachi W., Oria P.-A., Ibata R., Famaey B., Yuan Z., Arentsen A., Martin N., Viswanathan A., 2022, ApJ, 935, L22
  * Trick et al. (2021) Trick W. H., Fragkoudi F., Hunt J. A. S., Mackereth J. T., White S. D. M., 2021, MNRAS, 500, 2645
  * Vasiliev (2019) Vasiliev E., 2019, MNRAS, 482, 1525

##  Appendix A Gaia archive query

The following query has been used on the Gaia
archive777https://gea.esac.esa.int/archive/ to obtain our DISC samples and is
detailed in Section 2:

⬇

select source_id, bp_rp, phot_g_mean_mag, phot_bp_rp_excess_factor, ruwe, ra,
dec,

parallax, pmra, pmdec, parallax_error, pmra_error, pmdec_error,
parallax_pmra_corr,

parallax_pmdec_corr, pmra_pmdec_corr, visibility_periods_used,
astrometric_chi2_al,

astrometric_n_good_obs_al, radial_velocity,

radial_velocity_error,

if_then_else(

bp_rp > -20,

to_real(case_condition(

phot_bp_rp_excess_factor - (1.162004 + 0.011464* bp_rp +
0.049255*power(bp_rp,2)

- 0.005879*power(bp_rp,3)),

bp_rp < 0.5,

phot_bp_rp_excess_factor - (1.154360 + 0.033772* bp_rp +
0.032277*power(bp_rp,2)),

bp_rp >= 4.0,

phot_bp_rp_excess_factor - (1.057572 + 0.140537*bp_rp)

)),

phot_bp_rp_excess_factor

) as excess_flux

from gaiadr3.gaia_source

where parallax_over_error > 10

and parallax > 5

and ruwe < 1.15

and phot_g_mean_flux_over_error > 50

and phot_rp_mean_flux_over_error > 20

and phot_bp_mean_flux_over_error > 20

and visibility_periods_used > 8

and astrometric_chi2_al/(astrometric_n_good_obs_al-5)

< 1.44*greatest(1,exp(-0.4*(phot_g_mean_mag-19.5)))

To create our HALO samples, we have used a similar query with a different
parallax cut:

⬇

select source_id, bp_rp, phot_g_mean_mag, phot_bp_rp_excess_factor, ruwe, ra,
dec,

parallax, pmra, pmdec, parallax_error, pmra_error, pmdec_error,
parallax_pmra_corr,

parallax_pmdec_corr, pmra_pmdec_corr, visibility_periods_used,
astrometric_chi2_al,

astrometric_n_good_obs_al, radial_velocity,

radial_velocity_error,

if_then_else(

bp_rp > -20,

to_real(case_condition(

phot_bp_rp_excess_factor - (1.162004 + 0.011464* bp_rp +
0.049255*power(bp_rp,2)

- 0.005879*power(bp_rp,3)),

bp_rp < 0.5,

phot_bp_rp_excess_factor - (1.154360 + 0.033772* bp_rp +
0.032277*power(bp_rp,2)),

bp_rp >= 4.0,

phot_bp_rp_excess_factor - (1.057572 + 0.140537*bp_rp)

)),

phot_bp_rp_excess_factor

) as excess_flux

from gaiadr3.gaia_source

where parallax_over_error > 10

and parallax > power(3, -1)

and ruwe < 1.15

and phot_g_mean_flux_over_error > 50

and phot_rp_mean_flux_over_error > 20

and phot_bp_mean_flux_over_error > 20

and visibility_periods_used > 8

and astrometric_chi2_al/(astrometric_n_good_obs_al-5)

< 1.44*greatest(1,exp(-0.4*(phot_g_mean_mag-19.5)))

and 4.74*sqrt(power(pmra, 2) + power(pmdec, 2))/parallax > 200

##  Appendix B Red and blue CMD sequence selection

Table 3 lists the vertices for the intersect between red and blue regions of
the CMD outlined in Section 2. In addition to this, the blue section starts
from the point $(G_{\mathrm{BP}}-G_{\mathrm{RP}},M_{G})=(-0.125,11)$ mag and
passes through the points (0.8, 1) mag and (-0.125, 1.95) mag after the
intersecting vertices. The red section starts from (3.125, 11) mag and after
the intersecting vertices ends at (3.125, 5.25) mag.

Table 3: Vertices for the regions in $M_{G}$ and $G_{\mathrm{BP}}-G_{\mathrm{RP}}$ that make up our red and blue halo sequences. $G_{\mathrm{BP}}-G_{\mathrm{RP}}$ | $M_{G}$
---|---
$\mathrm{mag}$ | $\mathrm{mag}$
2.4 | 11
2.110 | 10.087
2.012 | 9.522
1.789 | 8.630
1.637 | 8.178
1.452 | 7.565
1.262 | 7.021
1.142 | 6.652
1.034 | 6.239
0.936 | 5.826
0.843 | 5.326
0.756 | 4.782
0.729 | 4.478
0.708 | 4.000
0.713 | 3.695
0.756 | 3.413
0.849 | 3.326
0.925 | 3.260
0.963 | 3.239
0.990 | 2.412
1.066 | 1.173

##  Appendix C Halo velocity distributions in binned velocities

In Section 5 we discuss certain features in the concept of the full 3D
velocity structure. The way we visualise the full structure is by looking at
the different velocity spaces in bins of the third velocity component, here in
steps of 100 km s-1. These figures show how some of the features are related
across velocity spaces.

Figure 8: Velocity distributions in spherical coordinates $v_{r}$ and
$v_{\phi}$, binned by their $v_{\theta}$. Figure 9: Velocity distributions in
spherical coordinates $v_{r}$ and $v_{\theta}$, binned by their $v_{\phi}$.
Figure 10: Velocity distributions in spherical coordinates $v_{\phi}$ and
$v_{\theta}$, binned by their $v_{r}$.