text:
\section{Introduction}

In the Cold Dark Matter (CDM) cosmogony, galactic stellar haloes are
built up in large part from the debris of tidally disrupted satellites
(e.g. Searle \& Zinn 1978; White \& Springel 2000; Bullock \& Johnston 2005; Cooper {et~al.} 2010). Discovering
and quantifying halo structures around the Milky Way may provide a
useful diagnostic of the Galaxy's merger history
(e.g. Helmi \& de~Zeeuw 2000; Johnston {et~al.} 2008; G{\'{o}}mez \& Helmi 2010). Upcoming Milky Way
surveys (for example with PanSTARRS1, LAMOST, HERMES and the LSST)
will provide large datasets in which to search for structure. The
\textit{Gaia} mission will determine six-dimensional phase-space
coordinates for all stars brighter than $V\sim17$, from which it
should be possible to untangle even well-mixed streams in the nearby
halo (G{\'{o}}mez {et~al.} 2010).

Testing the CDM model by comparing these observations to numerical
simulations of stellar halo formation requires a straightforward
definition for the `abundance of substructure', one that can be
quantified with a method equally applicable to simulations and
observations. Algorithms already exist for identifying substructure in
huge multidimensional datasets (e.g. Sharma \& Johnston 2009), such as the
data expected from \textit{Gaia} supplemented by chemical abundance
measurements (Freeman \& Bland-Hawthorn 2002). These algorithms can
also be applied to simulations, although this is not
straightforward. One problem is that current (cosmological)
hydrodynamic simulations still fall short of the star-by-star
`resolution' of the \textit{Gaia} data, particularly in the Solar
neighbourhood (e.g. Brown, Vel{\'{a}}zquez \& Aguilar 2005).

In the outer halo, longer mixing times allow ancient structures to
remain coherent in configuration space for many gigayears. However, 6D
\textit{Gaia} data will be restricted to relatively bright stars. In
the near future, studies of the outer halo (beyond $\sim20$~kpc) will
continue to rely on a more modest number of `tracers' (giant and
horizontal branch stars). For these stars, typically only angular
positions and (more uncertain) estimates of distance and radial
velocity are available. In this regime, current simulations contain as
many particles as there are (rare) tracer stars in observational
samples, and can be compared directly with data that are already
available. Here we focus on quantifying the degree of structure in
rare tracer stars with a generic method, which we apply to
observational data and to simulations of Milky Way-like stellar
haloes.

Most studies of spatial and kinematic structure in the Milky Way halo
have given priority to the discovery of individual overdensities
(exceptions include Gould 2003, Bell {et~al.} 2008, Xue, Rix \& Zhao 2009,
Xue {et~al.} 2011 and Helmi {et~al.} 2011a). Relatively few
have investigated global statistical quantities for the entire stellar
halo, although several authors have suggested an approach based on
clustering statistics. Re~Fiorentin {et~al.} (2005) analysed the
velocity-space clustering of a small number of halo stars in the Solar
neighbourhood, using a correlation function statistic. Following early
work by Doinidis \& Beers (1989), Brown {et~al.} (2004) examined the angular
two-point correlation function of photometrically selected blue
horizontal branch (BHB) stars in the Two Micron All Sky Survey,
probing from $\sim2-9$~kpc. They detected no significant correlations
at latitudes $|b|\ga50\degr$, but did detect correlations on small
scales ($1\degr$, $\sim100$~pc) at lower latitudes, which they
attributed to structure in the thick disc. Lemon {et~al.} (2004) performed a
similar analysis for nearby F-type stars in the Millennium Galaxy
Catalogue and found no significant clustering.

Starkenburg {et~al.} (2009) used a correlation function in \textit{four}
dimensions to identify substructures in the Spaghetti pencil-beam
survey of the distant halo (Morrison {et~al.} 2000; Dohm-Palmer {et~al.} 2000).
With this method they obtained a significant detection of clustering
and set a lower limit on the number of halo stars in all
substructures. Similarly, Schlaufman {et~al.} (2009) constrained the mass
fraction of the halo in detectable substructure by estimating the
completeness of their overdensity detection
algorithm. Starkenburg {et~al.} and Schlaufman {et~al.}
concluded that $>10\%$ (by number of stars) and $\sim30\%$ (by volume)
of the Milky Way halo belongs to groups meeting their respective
definitions of phase space substructure. These methods were tested on
`mock catalogues' of tracer stars based on simplified models of the
stellar halo.

The work of Xue {et~al.} (2009, 2011) is particularly
  relevant. Xue {et~al.} (2009) considered the pairwise radial velocity
  separation of a sample of 2558 halo BHB stars as a function of their
  separation in space, but found no evidence of clustering. From
  comparisons to the simulations of Bullock \& Johnston (2005),
  Xue {et~al.} concluded that a pairwise velocity statistic was
  not capable of detecting structure against a more smoothly
  distributed background in phase space (made up from stars in
  phase-mixed streams). However, their observed signal was not
  compared to an expected signal from random realizations. More
    recently, Xue {et~al.} (2011) studied an enlarged catalogue of
    BHBs comprising more than 4000 stars. They quantified clustering
    in this larger sample using a four-dimensional metric
    (similar to that of Starkenburg {et~al.} 2009), finding a significant
    excess of clustering on small scales by comparison to smooth
    models. The conclusions of this more recent study by
    Xue {et~al.} agree with our own, as we discuss further
    in \mnsec{sec:conclusion}.

The statistic we develop in this paper also builds on the
  approach of Starkenburg {et~al.} (2009). We define a two-point correlation
  function based on a metric combining pairwise separations in four
  readily obtained phase space observables for halo stars (angular
  position, radial distance and radial velocity). We apply this
  statistic to the catalogue of BHB stars published by
  Xue {et~al.} (2008)\footnote{The larger sample used by Xue {et~al.} (2011)
    was not publicly available at the time this paper went to press.}
  and demonstrate that a significant clustering signal can be
  recovered from these data.

A clustering metric of the kind we propose can be tuned to probe a
specific scale in phase space by adjusting the weight given to each of
its components. In the stellar halo, however, many `components' may be
superimposed with a complex assortment of scales and morphologies in
phase space (Johnston {et~al.} 2008; Cooper {et~al.} 2010). For this reason, it is not
clear, a priori, what sort of signal to expect, or which scales are
most relevant. We find that we cannot identify an `optimal'
  metric. Instead, we make a fiducial choice based on the
self-consistent accreted halo models of Cooper {et~al.} (2010). These
incorporate an ab initio \lcdm{} galaxy formation model using
high-resolution cosmological N-body simulations from the Aquarius
project (Springel {et~al.} 2008).  We apply our fiducial metric to the
data and to these models. We find that even though both the metric and
our choice of scaling are simple, this approach has the power to
discriminate quantitatively between qualitatively different stellar
haloes.

We describe the basis of our method in \mnsec{sec:method} and the SDSS
DR6 BHB catalogue of Xue {et~al.} (2008) in \mnsec{sec:obsdata}. In
\mnsec{sec:sims} we describe our simulations
(\mnsec{sec:nbodygalform}) and our procedure for constructing mock
catalogues (\mnsec{sec:tracer_stars}). \mnsec{sec:fiducial} describes
our choice of a fiducial metric. In \mnsec{sec:segue} we apply our
method to quantify clustering in the SDSS data and compare this with
our mock catalogues. Our conclusions are given in
\mnsec{sec:conclusion}.

\section{A metric for phase-space distance}
\label{sec:method}

The most readily obtained phase-space observables for halo stars are
their Galactic angular coordinates, $l$ and $b$, heliocentric radial
distance, $r_{\mathrm{hel}}$, and heliocentric line-of-sight velocity,
$v_{\mathrm{hel}}$. Using its angular position and distance estimate,
each star can be assigned a three-dimensional position vector in
galactocentric Cartesian coordinates, $\bmath{r}\,(X,Y,Z)$, and a
radial velocity corrected for the Solar and Local Standard of Rest
motions, $v$. We begin by defining a scaling relation (metric),
$\Delta$, which combines these observables into a simple `phase-space
separation' between two stars:

\begin{equation}
{\Delta}_{ij}^2=
{|\bmath{r_i}-\bmath{r_j}|}^2+w_{v}^2({v_{i}-v_{j}})^2.
\label{eqn:delta_metric}
\end{equation}

Here, $|\bmath{r_i}-\bmath{r_j}|$ is the separation of a pair of stars
in coordinate space (in kiloparsecs), and $v_{i}-v_{j}$ is the
difference in their radial velocities (in kilometres per second). The
scaling factor $w_{v}$ has units of $\mathrm{kpc\,km^{{-}1}\,s}$, such
that $\Delta$ has units of kpc. The choice of $w_{v}$ is arbitrary
unless a particular `phase space scale' of interest can be
identified. This is not straightforward; we discuss some possible
choices below.

The aim of this paper is to explore $\xi({\Delta})$, the cumulative
two-point correlation function of halo stars in the metric defined by
Equation \ref{eqn:delta_metric}. Throughout, we use the
estimator \begin{equation} 1+\xi({\Delta}) =
  \frac{DD(<{\Delta})}{\left \langle RR(<{\Delta}) \right
    \rangle}.  \label{eqn:estimator}
\end{equation} Here $DD(<\Delta)$ counts the number of pairs in the
sample separated by less than $\Delta$, and $\left \langle RR(<\Delta)
\right \rangle$ is the equivalent count for pairs of random points
within the survey volume, averaged over several realizations. To
  generate these realizations we `shuffle' the data by randomly
  reassigning $(r_{\mathrm{hel}},v_{\mathrm{hel}})$ pairs to different
  $(l,b)$ coordinates drawn from the original catalogue\footnote{The
    same randomisation procedure was used by
    Starkenburg {et~al.} (2009). Xue {et~al.} (2011) use a similar
    procedure, but for each galaxy they re-assign $r_{\mathrm{hel}}$
    and $v_{\mathrm{hel}}$ \textit{separately} to different $(l,b)$
    coordinates.}.

Similar methods for quantifying the clustering of stars in a
  four-dimensional space are described by Starkenburg {et~al.} (2009, applied to a sample
    of giant stars from the Spaghetti survey) and
  Xue {et~al.} (2011, applied to a sample of BHB stars from
    SDSS). Our $\Delta$ metric is very similar to the
  $\delta_{4\mathrm{d}}$ metric of Starkenburg et al. in the limit of
  small angular separations\footnote{Starkenburg et al. developed
    their metric with the aim of identifying `true' pairs of stars
    with high confidence. In their definition (equation 1 of
    Starkenburg {et~al.} 2009), the $\delta_{4\mathrm{d}}$ separation
    between two stars depends not only on their actual phase-space
    coordinates, but also on how accurately those coordinates are
    measured. For example, moving two stars 10 kpc further apart and
    simultaneously decreasing the error on their distance measurements
    by a factor of 10 (relative to the average of the sample) results
    in the same $\delta_{4\mathrm{d}}$. Thus $DD/RR$ for separations
    in $\delta_{4\mathrm{d}}$ is not a straightforward measurement of
    physical clustering.}. We have verified that our metric gives
  similar results when we repeat the analysis of Starkenburg et
  al. using the Spaghetti sample of 101 halo RGB stars. For the rest
  of this paper we will focus on clustering in the SDSS DR6 BHB
  catalogue of Xue {et~al.} (2008).

\section{Observational data}
\label{sec:obsdata}

Xue {et~al.} (2008) have published a catalogue of 2558 stars from SDSS DR6,
which they identify as high-confidence halo BHBs (contamination
$<10\%$) using a combination of colour cuts and Balmer line
diagnostics. This sample ranges in distance from $4-60\,\mathrm{kpc}$;
a cut on distance error excluded more distant stars observed in
SDSS. Xue {et~al.} (2008) claim distance errors of $\sim10\%$ and radial
velocity errors of $5-20\,\mathrm{km\,s^{-1}}$. This catalogue is not
a complete sample of halo BHBs. In particular, Xue {et~al.} note
that the targeting of SDSS spectroscopy is biased against the
follow-up of more distant stars.

To isolate stars representing the kinematics of the halo (in order to
study the Galactic circular velocity profile) Xue {et~al.} (2008)
further restricted their sample to stars at heights
$|z|>4\,\mathrm{kpc}$ above the plane (avoiding the thick disc). We
also apply this cut (which excludes thick disc stars and low-latitude
halo stars alike), leaving 2401 stars in the sample. Finally we
exclude 9 stars in the sample identified by Xue {et~al.} (2009) as probable
globular cluster members. Thus, the sample against which we compare our
models contains 2392 stars from the 2558 stars in the Xue {et~al.} (2008)
catalogue. The effects of these refinements to the sample on the
recovered \deltacf{} signal are discussed in
\mnsec{sec:obsdata_clustering}.

\section{Stellar Halo Simulations}
\label{sec:sims}

\subsection{N-body and galaxy formation model}
\label{sec:nbodygalform}

The mock observations that we use to test the \deltacf{} correlation
function are derived from simulations of the accreted stellar halo
presented in Cooper {et~al.} (2010). These simulations approximate the
dynamics of stars in dwarf satellites of Milky Way-like galaxies by
`tagging' appropriate particles (i.e. those strongly bound within
subhaloes) in the Aquarius suite of high-resolution N-body simulations
(Springel {et~al.} 2008). Each `tag' associates a dark matter (DM)
particle with a particular stellar population of a given mass, age
and metallicity. This `tagging' technique is reasonable in the regime
of high mass-to-light ratios, which is supported in this case by
observations of stellar kinematics in dwarf galaxies
(e.g. Walker {et~al.} 2009).

The tagging method has a single free parameter, the fraction of the
most-bound particles chosen in each DM halo for each assignment of
newly-formed stars (for further details see Cooper {et~al.} 2010). The
value of this parameter (1 per cent) was fixed by requiring the
population of \textit{surviving} satellites (at the present day) to
have a distribution of half-light radius as a function of luminosity
consistent with Milky Way and M31 observations\footnote{The luminosity
  function of surviving satellites in these models also agrees with MW
  and M31 data. This agreement is mostly due to the underlying galaxy
  formation model.}. The Cooper {et~al.} models differ from the
earlier models of Bullock \& Johnston (2005) in that they treat the full
cosmological evolution of all satellites self-consistently in a single
N-body simulation, and use a comprehensive semi-analytic model of
galaxy formation (Bower {et~al.} 2006) constrained by data on large scales
and compatible with the observed MW satellite luminosity
function. Both the Cooper {et~al.} and the
Bullock \& Johnston simulations produce highly structured stellar
haloes built from the debris of disrupted dwarf galaxies. As we
discuss further below, other halo components formed {\it in situ} may
be present in real galaxies
(e.g. Abadi, Navarro \& Steinmetz 2006; Zolotov {et~al.} 2009; Font {et~al.} 2011) but these are expected
to be more smoothly distributed than the accreted component
(e.g. Helmi {et~al.} 2011a).

 \begin{figure*}
\includegraphics[height=60mmi,clip=True,trim=1cm 0cm 1cm 0 ]{fig1a.pdf}
\includegraphics[height=60mm,clip=True]{fig1b.pdf}
\includegraphics[height=60mm,clip=True,trim=1cm 0cm 1cm 0 ]{fig1c.pdf}
\includegraphics[height=60mm,clip=True, trim=0 0 0 0 ]{fig1d.pdf}

\caption{\textit{Left panels:} An example of the sky distribution of
  halo BHB stars in Aq-A from the perspective of a `Solar' observer,
  shown as Mollweide projections in galactic coordinates centred on
  $(l,b)=(0,0)$. Colours indicate the mean heliocentric distance of
  stars in each pixel, in kiloparsecs. Pixels outside the SDSS DR6
  footprint are shown with lower contrast. The upper panel includes
  all BHB stars from $6$--$60$~kpc, the lower panel includes only
  stars in the range $30$--$60$~kpc. Our fiducial choice of the
  Galactic $Z$ axis with respect to the dark halo has been applied,
  but distances in these panels are not convolved with observational
  errors. The normalization of BHB stars per unit stellar mass in the
  halo has been increased in these panels to emphasise the
  distribution of structure. \textit{Right panels:} Blue histograms
  show the distribution of heliocentric distances (above) and
  heliocentric velocities (below) for the fiducial observer
  corresponding to the left panels, after convolution with
  observational errors (see text). Orange histograms are the
  distributions for stars in the Xue {et~al.} (2008) catalogue. To compare
  the shape of the two distributions, the normalization of the mock
  distributions in these panels has been matched to that of the
  observations for $r<20$~kpc, where the observations are most likely
  to be complete.}
\label{fig:fiducial}
 \end{figure*}

As in Springel {et~al.} (2008) and Cooper {et~al.} (2010), we refer to our
six simulations as haloes
Aq-A, Aq-B, Aq-C, Aq-D, Aq-E and Aq-F. From these simulations, we
construct catalogues of tracer stars (BHB stars, for example) by
converting the stellar mass assigned to each dark matter particle into
an appropriate number of stars.

Each DM particle can give rise to many tracer stars if it is tagged
with sufficient stellar mass. We therefore interpolate the positions
and velocities of these tracer stars between nearby tagged DM
particles. To accomplish this, the 32 nearest phase space neighbours
of each tagged particle are identified using a procedure which we
describe below. The mean dispersion in each of the six phase-space
coordinates is then calculated for each particle as an average over
its neighbours. These dispersions define a 6D ellipsoidal Gaussian
kernel centred on the particle, from which the positions and
velocities of its tracer stars are drawn randomly. Each progenitor
galaxy (a set of tagged DM particles accreted by the main `Milky Way'
halo as members of a single subhalo) is treated individually in this
smoothing operation, i.e. particles are smoothed using only neighbours
from the same progenitor (so there is no `cross talk' between streams
from different progenitors). This procedure can be thought of as a
crude approximation to running our original simulation again including
each tracer star as a test particle.

The `distance in phase space' used to identify neighbours in the
interpolation scheme is defined by a scaling relation between
distances in configuration space and velocity space\footnote{In this
  part of the calculation, we are only interested in finding
  neighbours, so the absolute values of these distances are not
  important. This scaling of velocity space to configuration space for
  the purpose of resampling the simulations should not be confused
  with the $\Delta$ metric we define for our analysis of
  clustering.}. For each progenitor, we adopt an individual scaling
which corresponds to making the median pairwise interparticle
separation of its particles in configuration space (at $z=0$) equal to
their median separation in velocity space. In practice, the value of
this scaling parameter makes very little difference to the results we
present, when compared to the extreme choice of selecting only 32
velocity or position neighbours (disregarding the other three
coordinates in each case). Giving more weight to configuration-space
neighbours smears out velocity substructure within the debris of a
progenitor (for example, where two wraps of a stream pass near one
another). Giving more weight to velocity neighbours has the opposite
effect -- stars can be interpolated over arbitrarily large separations
in configuration space, but coherent velocity structures are
preserved. Therefore, the `optimal' choice is the scaling which
balances smoothing in configuration space against smoothing in
velocity space.

To quantify this balance between smoothing in configuration and
velocity space, we compute six smoothing lengths for each particle,
$\epsilon_{x,i}$ and $\epsilon_{v,i}$, where $i$ represents a single
dimension in space or velocity. To compute these dispersions,
  we use the 32 nearest phase-space neighbours of each particle. We
define the `optimum' choice of scaling for \textit{each} progenitor
galaxy as that which minimises the quantity

\begin{equation}
\sigma_{\epsilon}^{2} = \left (\frac{1}{\bar{\epsilon}_{x,\mathrm{min}}} \sum_{i=0}^{3}{\epsilon_{x,i}} \right )^{2} + \left ( \frac{1}{\bar{\epsilon}_{v,\mathrm{min}}} \sum_{i=0}^{3}{\epsilon_{v,i}}\right )^{2}.
\end{equation} This is the sum in quadrature of the mean smoothing lengths in
configuration and velocity space, normalized respectively by
$\bar{\epsilon}_{x,\mathrm{min}}$, the `minimal' mean smoothing length
in configuration space (obtained from the 32 nearest configuration space
neighbours) and $\bar{\epsilon}_{v,\mathrm{min}}$, the `minimal' mean
smoothing length in velocity space (obtained from the 32 nearest
velocity space neighbours). We find that the scaling obtained by
matching the median interparticle separations in position and velocity
as described above is typically a good approximation to this optimal
value -- a similar result is discussed in more detail by
Maciejewski {et~al.} (2009).

In the Cooper {et~al.} model, when a given stellar population is
formed, the most bound 1\% of DM particles in the corresponding dark
halo at that time are chosen as dynamical tracers of that
population. Hence, each DM \textit{particle} to which stars are
assigned has an individual mass-to-light ratio, M/L, which can be as
high as $\sim1$ (i.e. $M_{\mathrm{stellar}}\sim10^{4}\,M_{\sun}$) and
as low as $\sim10^{-6}$. This will affect the density of stars seeded
by a DM particle independently of the density of its neighbours in
phase space (i.e. a low M/L particle will create a denser `cluster' of
tracers relative to a high M/L particle with the same neighbouring
positions and velocities). We have tested an alternative approach in
which the M/L of each particle in a given progenitor is resampled by
distributing the total stellar mass of the progenitor evenly amongst
its tagged particles\footnote{This is almost equivalent to choosing
  M/L only once, at the time in the simulation when the progenitor
  falls into the main halo (similar to the lower-resolution model of
  De~Lucia \& Helmi 2008).}. We find that the extra clustering due to a
few `hot' particles in our default approach makes no difference to our
results.

\subsection{Tracer stars and mock observations}
\label{sec:tracer_stars}

Each N-body dark matter particle in our simulation contributes a
number of tracer stars to our mock observations, based on the stellar
population with which it has been `tagged'. In the specific case of
our comparisons to SDSS, these tracers are BHB stars meeting the selection
criteria of Xue {et~al.} (2008). Here we assume a global scaling between the
stellar mass associated with each N-body particle, $M_{\star}$, and
the number of BHBs it contributes to our mock catalogues, i.e.
$N_{\mathrm{BHB}} = f_{\mathrm{BHB}}M_{\star}$ where
$f_{\mathrm{BHB}}$ is the number of tracer stars per unit mass of the
original stellar population\footnote{We do this because we prefer to
  make a straightforward comparison with the observational data in
  this paper. In principle, the age and metallicity information
  associated with each stellar population in our model could be used
  to populate an individual colour-magnitude diagram for each N-body
  particle, and make a detailed prediction for the appropriate number
  of tracers. The `bias' of BHBs relative to the total stellar mass
  distribution of the halo (Bell {et~al.} 2010) may affect the
  clustering statistic recovered (Xue {et~al.} 2011), but this effect
  is beyond the scope of the present paper. }. For each N-body
particle, the actual number of BHB stars generated is drawn from a
Poisson distribution with mean $N_{\mathrm{BHB}}$. The correlation
function results described below are not sensitive to the choice of
$f_{\mathrm{BHB}}$, provided that the underlying distribution is
well-sampled at a given scale. We have therefore selected a fiducial
value of $f_{\mathrm{BHB}}^{-1} = 3000\,M_{\sun}/\mathrm{star}$. In
creating the mock catalogue, we do not include any stars
gravitationally bound to satellites. However, we do include stars in
their tidal tails (which, by our definition, are part of the stellar
halo).

Using our simulated BHB catalogues, we create mock observations for
comparison to the Xue {et~al.} (2008) data as follows. First we located the
observer at a radius $r_{\sun}=8\,\mathrm{kpc}$ from the centre of the
halo. For our main comparison to the data, we restrict all observers
to the same `Galactic plane', with each random vantage point differing
only in its azimuthal angle in this plane and in the `polarity' of the
Galactic rotation axis (the $Z$ coordinate). However, the orientation
of the rotation axis cannot be directly constrained by the simulation,
which only models the \textit{accreted} component of the halo and the
bulge, and not the in situ formation of a stellar disc. As described
in Cooper {et~al.} (2010), the accreted `bulge' is prolate or mildly
triaxial. We define the minor axis of this bulge component
(conservatively defined by all accreted stars within $r<3$~kpc;
Cooper {et~al.} 2010) as the Galactic $Z$ axis. This axis is essentially
identical to the minor axis of the dark halo within $r<3$~kpc. There
are other plausible choices of Galactic plane (for example, relative
to the shape or spin vectors of the entire dark halo, rather than the
stars in its inner regions). However, any choice is somewhat arbitrary
without a self-consistent simulation of disc formation\footnote{In a
  full hydrodynamic simulation the effects of feedback and adiabatic
  contraction may also make the dark halo itself more spherical
  (e.g. Tissera {et~al.} 2010; Abadi {et~al.} 2010).}.

Having chosen a location for the observer, we select all tracer stars
within the spectroscopic footprint of SDSS DR6 having galactocentric
distance in the range 20--60~kpc (our principal comparison will focus
on the outer halo as defined by this distance range, although we also
study the ranges 5--60~kpc and 5--20~kpc below). Galactic longitude
and latitude are defined such that $(l,b)=(0,0)$ is the vector
directed from the observer to the centre of the halo. We set the
heliocentric velocity components of each star assuming a Solar motion
of $U,V,W = (10.0,5.2,7.2)\:\mathrm{km\,s^{{-}1}}$ (Dehnen \& Binney 1998) and
a velocity of the Local Standard of Rest about the Galactic centre
$v_{\mathrm{LSR}}=220\mathrm{\:km\,s^{{-}1}}$. We compute $(X,Y,Z)$
and $v_{\mathrm{los}}$ as described by Xue {et~al.} (2008). Finally,
distances and velocities are convolved with Gaussian observational
errors of $\sigma_{d}= 10\%$ and $\sigma_{v}=15\;\mathrm{km\,s^{-1}}$
respectively (Xue {et~al.} 2008).

In both the mock observations and the real data, the average random
pair count $\langle RR \rangle$ is calculated by reshuffling distances
and velocities among the positions on the sky of stars many times. We find that
by using 500 random catalogues to calculate $\langle RR \rangle$ for
each mock observation and performing 500 mock observations in each
halo, we obtain a sufficiently well-converged estimate of the
distribution of \deltacf{}.

\fig{fig:fiducial} illustrates the structure of one of our haloes and
verifies that our mock observations can result in distributions of
heliocentric distance and heliocentric radial velocity similar to the
SDSS data of Xue {et~al.} (2008). In this figure we have specifically chosen
an observer orientation in halo Aq-A such that the distributions of
distance and velocity we recover are close to those of the data, after
convolution with typical observational errors. We have aligned
the Galactic $Z$ axis of the mock observer with the minor axis of
the dark halo, as described above. This confines the most prominent
structures in the stellar halo to low Galactic latitudes, outside the
SDSS DR6 spectroscopic footprint. Of course, the simulated haloes are
inhomogeneous on large scales, and there are many choices of observer
in each halo that \textit{do not} resemble the SDSS data\footnote{As
  discussed by Xue {et~al.} 2008, the completeness of the data declines
  at larger distances ($r_{\sun}>20$~kpc). Mock catalogue distributions
  that match the observed distributions well at $r_{\sun}<20$~kpc
  typically show a flatter profile with identifiable overdensities
  (streams) at larger distances. As the SDSS spectroscopic selection
  function for the data we use is difficult to quantify (Xue {et~al.} 2008),
  we do not explore the effects of incompleteness in this paper.}.

\section{Distance - velocity scaling in the $\Delta$ metric}
\label{sec:fiducial}

\begin{figure}
\includegraphics[width=84mm,clip=True]{fig2.pdf}
\caption{Correlation functions in spatial separation (blue) and velocity
  separation (red) for stars in halo Aq-A. The velocity separation
  correlation function has been scaled by a factor
  $w_{v}=0.04\,\wvunits$ to match the turnover in the configuration
  space separation correlation function.}
\label{fig:accretedscale}
\end{figure}

Before \deltacf{} can be computed, a value must be chosen for the
  velocity-to-distance scaling $w_{v}$ in Equation
  \ref{eqn:delta_metric}. There is no clearly well-motivated way to
  choose this value, and in the absence of a physical justification,
it must be treated as a free parameter. The choice of $w_{v}$
determines the scale of substructure to which the correlation function
is most sensitive. Naively, we expect this to be the typical width and
transverse velocity dispersion of a `stream'. It is preferable to fix
this parameter in a universal manner that does not depend on the
  details of a particular survey. We make a fiducial choice of
$w_{v}$ as follows.

In each simulated halo we adopt the SDSS-like survey configuration
described in \mnsec{sec:tracer_stars} (without observational errors or
assumptions about the orientation of the Galactic plane). We
construct one dimensional distributions of the separation in radial
distance and velocity between stars. We generate many random
realizations of these distributions by first convolving each simulated
star with Gaussian smoothing kernels of width $8\,\kpc$ (distance) and
$80\,\kms$ (velocity), and then drawing randomly from these `smoothed'
distributions. The smoothing scales were chosen as a compromise
between signal (diminished by oversmoothing) and noise (increased by
undersmoothing). Using these random realizations we construct
one-dimensional correlation functions for each distribution. These two
correlation functions are shown for halo Aq-A in
\fig{fig:accretedscale}. Although the signals are intrinsically weak,
they have a very similar shape for both distributions, each with a
characteristic `turnover' scale. Matching this scale in the two
correlation functions corresponds to $w_{v}\sim0.04\pm0.01\,\wvunits$
for the six haloes, which we adopt as a fiducial value. We caution
that although the scales on which we match the one-dimensional
correlation functions are somewhat smaller than the smoothing scales
we adopt to create the random distributions, this does not guarantee
that our choice of $w_{v}$ is unaffected by our choice of smoothing.

Clearly, there are other ways of fixing $w_{v}$. In practice, however,
our conclusions are not highly sensitive to the value of this
parameter. Values of the order of $w_{v}\sim0.01$--$1.0\,\wvunits$
result in very similar $\xi({\Delta})$ correlation functions. Values
lower than $0.01\,\wvunits$ recover very little signal. Values above
$1\,\wvunits$ treat $1\,\kms$ velocity differences as equivalent to
$>1\,\kpc$ separations in space, and so make the cumulative
correlation function very noisy on small scales for only a marginal
increase in the overall signal. (This noise, in turn, increases the
scatter between signals measured by different observers.) We find that
our choice of $w_{v}\sim0.04\,\wvunits$ is a reasonable
compromise. Our method for choosing $w_{v}$ can be compared with that
of Starkenburg {et~al.} (2009), who determine the equivalent of $w_{v}$ in
their metric to be the ratio of the Spaghetti survey limits in radial
distance and velocity ($0.26\,\wvunits$). Either value is acceptable
to illustrate our approach and compare to simulations. We therefore
adopt $w_{v}\sim0.04\,\wvunits$.

\section{Clustering of SDSS BHB stars}
\label{sec:segue}

\subsection{Clustering in the Xue et al. sample}
\label{sec:obsdata_clustering}

\begin{figure}
\includegraphics[width=84mm,trim= 0.0cm 0cm 0cm 0.0cm,
  clip=True]{fig3.pdf}
\caption{\deltacf{} correlation function for the SDSS BHB sample of
  Xue {et~al.} (2008). Black points (with Poisson error bars) show
  \deltacf{} computed for all stars at galactocentric distances
  greater than 20~kpc. Grey points show the result for all stars in
  the main sample (galactocentric distances of 4--60~kpc). The blue
  squares include 9 stars suspected to belong to globular clusters,
  while red circles include stars at low Galactic $|z|$ heights
  (possible thick disc stars). Neither of these contributions are
  relevant for the more distant selection shown by the black points.}
\label{fig:sdss_real}
\end{figure}

\fig{fig:sdss_real} shows \deltacf{} computed for 2392 BHB stars in
the Xue {et~al.} (2008) sample (\mnsec{sec:obsdata}; grey points). Stars at
small separations in the metric of \eqn{eqn:delta_metric}
($\Delta<4\,\mathrm{kpc}$) show significant clustering. The amplitude
of the signal increases if we restrict the sample to larger
galactocentric distances, $r>20$~kpc (black points). At larger
distances substructure is expected to be dynamically young and to have
undergone less phase mixing. Our finding of stronger clustering for
more distant halo stars is in qualitative agreement with the results
of Xue {et~al.} (2011).

Although we appear to recover a significant clustering signal in
\fig{fig:sdss_real}, we have only one SDSS survey. The observed signal
may be an artifact of the particular structures covered by the SDSS
footprint. Other parts of the halo may be smoother or more structured,
or may appear so when viewed from different points around the Solar
circle. We will address this issue of sample variance in the following
section using our mock catalogues.

We show two further permutations of the Xue {et~al.} (2008) sample in
\fig{fig:sdss_real}. The first of these (red open circles) includes
stars close to the Galactic plane, $|Z|<4\,\mathrm{kpc}$. These were
excluded from the main sample of Xue {et~al.} (2008) to excise the thick
disc.

Although only $\sim150$ stars are excluded by the cut on
  $|Z|$, they make a substantial difference to the correlation
  function, suppressing the clustering signal on scales below
  $\Delta<\sim8\,\mathrm{kpc}$. In the SDSS data the majority of
  low-$|Z|$ stars are at small heliocentric radii. These stars
  constitute a foreground `screen' with a relatively smooth
  distribution, which may dilute the signal of correlated stars.

The final sample shown in \fig{fig:sdss_real} (blue open squares)
includes all stars from the main sample (grey points) and a further
nine BHB stars identified as globular cluster members by
Xue {et~al.} (2009). Two of these are from one cluster, and seven from
another. Including these stars marginally increases the clustering
signal in the smallest-separation bin. This shows that the technique
is sensitive to the clustering of stars on these scales, which
correspond to separations comparable to the distance and velocity
errors of the data.

\subsection{Comparison with Mock Catalogues}
\label{sec:mock_clustering}

\begin{figure}
\includegraphics[width=84mm, clip=True]{fig4.pdf}
\caption{\deltacf{} for halo stars of galactocentric distances
  $20<r<60$~kpc, computed for 500 observers in the six Aquarius
  simulations of Cooper {et~al.} (2010). Coloured regions (delineated by
dashed and solid lines) indicate the
  10--90\% range of \deltacf{} in each bin of $\Delta$, for each
  halo. All observers are restricted to the Galactic plane, and
  distances and velocities are convolved with observational errors (see
  text). Black points with error bars reproduce the observed BHB
  \deltacf{} shown in \fig{fig:sdss_real}.}
\label{fig:sdss_mock_pa}
\end{figure}

In \fig{fig:sdss_mock_pa} we show the range of \deltacf{} measured for
each of the six simulations described in \mnsec{sec:sims},
overplotting the data shown in \fig{fig:sdss_real}. As described in
\mnsec{sec:tracer_stars} these results correspond to mock BHB
catalogues for 500 Solar observers located in the `Galactic plane',
with typical observation errors in distance and velocity. The dashed
and solid lines bounding each coloured region in
\fig{fig:sdss_mock_pa} correspond to the 10th and 90th percentile
values we obtain in each bin of \deltacf{}. At a given scale in our
$\Delta$ metric, the strength of the clustering signal varies
considerably from halo to halo and between individual observers.

The overall trend of \deltacf{} is similar to the observations in all
haloes, although the clustering signal rises more steeply on small
scales in most of the simulations. On scales $\Delta < 4$~kpc, significant
clustering is detected by all observers in five of the six haloes. The
exception is halo Aq-E (dark blue). This stellar halo is highly
concentrated, and at $r>20$~kpc is dominated by a single radial stream
(see figures 6 and 7 of Cooper {et~al.} 2010).

We find that two haloes, Aq-E and Aq-F (red), are consistent with the
observed \deltacf{} on all scales. The structure of Aq-F is atypical
for the sample -- most of its stars are accreted in a late 3:1 merger
and its surface brightness at the Solar radius is substantially higher
than current estimates for the Milky Way halo. In projection, Aq-F
resembles the `shell'-dominated haloes observed in a number of nearby
elliptical galaxies. Meanwhile haloes Aq-A (black), Aq-B (cyan) and
Aq-D (green) are marginally inconsistent with the data: below
$\Delta\sim4$~kpc, $\sim90$ per cent of mock observations in these
haloes imply a greater degree of clustering than we find for the Milky
Way, particularly on small scales. Aq-C (purple) is entirely
inconsistent with the Milky Way observations on all scales, showing a
much higher degree of clustering. Beyond $20$~kpc, the sky of an
observer in Aq-C is dominated by two bright tidal streams on wide
($\sim100$~kpc) orbits. Although their orbital planes are
approximately coincident with our definition of the Galactic plane,
nevertheless sections of these streams intrude on the SDSS footprint
at low Galactic latitudes.

\begin{figure*}
\includegraphics[height=65mm,clip=True]{fig5a.pdf}
\includegraphics[height=65mm,trim=1.3cm 0 0 0,clip=True]{fig5b.pdf}
\includegraphics[height=65mm,trim=1.3cm 0 0 0,clip=True]{fig5c.pdf}
\caption{\deltacf{} for mock observations, as
  \fig{fig:sdss_mock_pa}. From left to right, we vary our modelling
  assumptions as follows: (a) no restriction on the alignment of the
  Galactic plane with respect to the dark halo (the observer is
  located randomly on a sphere of radius $r_{\sun}=8$~kpc) and no
convolution of the data with the expected observational errors; (b) the
  observer is restricted to the Galactic plane as in
  \fig{fig:sdss_mock_pa}, but the mock observations are \textit{not}
  convolved with expected observational errors; (c) as panel (a), but
  mock observations \textit{are} convolved with errors. Colours are as
  \fig{fig:sdss_mock_pa}}
\label{fig:sdss_mock_no}
\end{figure*}

The DR6 footprint and the cut on extra-planar height in the
Xue {et~al.} (2008) sample exclude stars near the Galactic plane from our
clustering analysis. \fig{fig:sdss_mock_no} illustrates how our
definition of the Galactic plane influences the halo clustering
signal. In panel (a) the orientation of the Galactic plane with
respect to the halo is chosen randomly for each of the 500 mock
observers (i.e. observers are distributed over a sphere of radius
$r_{\sun}=8$~kpc), whereas in panel (b) the galactic Z direction
is aligned with the minor axis of the halo for all observers as in
\fig{fig:sdss_mock_pa}. To focus on the effects of this alignment, the
distances and velocities of stars in these two panels have
\textit{not} been convolved with observational errors.

The systematically higher clustering signals in panel (b) of
\fig{fig:sdss_mock_no} suggest that the plane perpendicular to the
minor axis of the dark matter halo is special. In Cooper {et~al.} (2010) and
above, we have noted the strong correlation between the shape of the
dark halo and the inner regions of the stellar halo. This alignment of
halo structure also extends, more loosely, to other prominent stellar
halo structures at large distances. An overall flattening of the
stellar halo arises because our dark matter haloes are embedded in
long-lived filaments of the cosmic web. Typically one or two such
filaments dominate the infall directions of both satellite galaxies
and smoothly accreted dark matter, which also contributes to the shape
of the dark halo
(e.g. Libeskind {et~al.} 2005; Lovell {et~al.} 2011; Wang {et~al.} 2011; Vera-Ciro {et~al.} 2011). The
distribution of stars stripped from infalling satellites echoes the
large-scale correlation of their orbital planes.

Because of this flattened global structure, the distribution of halo
stars in our choice of Galactic plane tend to be more smoothly
distributed (i.e. this plane contains more diffuse phase-mixed
material as well as coherent substructure). Panel (b) demonstrates how
the `contrast' of small scale substructure in the outer halo is
enhanced when these smoother components are excluded from the
clustering analysis (through a combination of the SDSS footprint and
the cut on $|Z|$). This is particularly true in the case of Aq-F,
where the majority of the mass in the halo is contributed by one
extensive and relatively `smooth' component. By contrast, in Aq-C the
average clustering amplitude \textit{decreases} on large scales when
we fix the Galactic plane. As noted above, in this case the massive
coherent streams that dominate the clustering signal of this halo
mostly fall outside the SDSS footprint.

Panel (c) of \fig{fig:sdss_mock_no} shows the randomly aligned
observations of panel (a) convolved with observational errors in
distance and velocity. These errors `smooth out' the halo, suppress
the clustering signal overall and increase the variance between
observers on small scales. Again the effect is most pronounced for
Aq-F, where blurring of the dominant smooth component further
decreases the contrast of substructure. In most cases these two
effects (alignment and observational errors) counteract each other to
produce the distribution of signals shown in
\fig{fig:sdss_mock_pa}. In the case of Aq-E the signal suffers
disproportionately from errors in the aligned configuration, perhaps
because this signal is due to a small number of pairs at large
distances.

Finally, in \fig{fig:sdss_mock_distance} we examine differences
between nearby and more distant halo stars (also discussed by
Xue {et~al.} 2011). The left-hand panel corresponds to the full
range of the Xue {et~al.} (2008) sample ($5<r<60$~kpc), the central panel
corresponds to nearby stars ($5<r<20$~kpc) which we excluded in
\fig{fig:sdss_mock_pa}, and the right-hand panel corresponds to the
most distant stars in the sample ($30<r<60$~kpc). As shown in
\fig{fig:sdss_real}, including nearby stars considerably reduces the
amplitude of \deltacf{} for Milky Way halo BHBs, although the signal
in the lowest bins remains significantly above zero. The observational
data are dominated by stars at $r<20$~kpc and hence the Milky Way
signal is essentially identical in the left and central panels.

In the whole-halo sample (left-hand panel of
\fig{fig:sdss_mock_distance}) some of the simulations (particularly
Aq-A and Aq-B) show a reduction in \deltacf{} on small scales,
qualitatively similar to the observations. Nevertheless
\fig{fig:sdss_mock_distance} suggests that the contrast of structure
against smoothly distributed stars is too high in our models, and that
the discrepancy is most pronounced in the inner halo. Selecting an
inner-halo sample (central panel of \fig{fig:sdss_mock_distance})
emphasises this point. The inner halo shows a similar degree of
clustering in all six haloes, in all cases somewhat lower than the
outer halo signal (on intermediate scales), but much stronger than the
Milky Way data. This panel reinforces our comments above on the
structure of the haloes: Aq-C resembles Aq-A and Aq-D in the inner
halo, because this distance range excludes the massive coherent
streams that are the likely source of the strong signal seen in
\fig{fig:sdss_mock_pa}. Most of the coherent structure in Aq-E is
nearby and the inclusion of stars at greater distances
\textit{diminishes} the average signal. This is in contradiction to
the phase-mixing explanation for a more structured outer halo, but
Aq-E is unusual in that for $r<20$~kpc it is dominated by a highly
oblate, coherently rotating component (somewhat resembling a very
thick disc, but with a significant radial velocity dispersion). Aq-B
has a similar spherically averaged density profile to Aq-E, but here
the inner halo shows less correlation, as expected. In Aq-F \deltacf{}
is dominated by nearby stars, although these are highly clustered in
phase space.

\begin{figure*}
\includegraphics[height=65mm,clip=True]{fig6a.pdf}
\includegraphics[height=65mm,clip=True,trim=1.3cm 0 0 0]{fig6b.pdf}
\includegraphics[height=65mm,clip=True,trim=1.3cm 0 0 0]{fig6c.pdf}
\caption{\deltacf{} for mock observations, as
  \fig{fig:sdss_mock_pa}. Here we select three different
  galactocentric distance ranges in the halo: from $5<r<60$~kpc (left
  panel), from $5<r<20$~kpc (central panel) and from $30<r<60$~kpc
  (right panel). The Milky Way BHB clustering signal (black points) is
  recalculated for the distance range in each panel. Colours are as
  \fig{fig:sdss_mock_pa}.}
\label{fig:sdss_mock_distance}
\end{figure*}

The comparison with more distant stars shown in the right-hand
  panel of \fig{fig:sdss_mock_distance} suggest that most of the
  models are consistent with the data in this region of the halo. Only
  Aq-C is clearly discrepant (Aq-D is marginally inconsistent at
  scales below $\Delta\sim12$~kpc). That our simulations more closely
  match the outer regions of the halo (beyond $30$~kpc) supports the
  view that these regions are dominated by stars accreted from
  satellites (Xue {et~al.} 2011; Font {et~al.} 2011). The outer stellar halo
  of the Milky Way appears to be typical (or perhaps slightly more
  structured than average) with respect to our sample of six plausible
  realizations of dark matter haloes with similar mass.

\section{Conclusions}
\label{sec:conclusion}

We have analysed a correlation function for halo stars, \deltacf{},
defining their separation in four dimensions of phase space using a
metric (which we call $\Delta$) in readily-obtained observables
(angular and radial separation and radial velocity difference). A
statistic of this type usefully quantifies kinematic and spatial
substructure in the halo, and can be applied to observational data and
catalogues generated from theoretical models. This analysis is
particularly well suited to the distant halo -- other methods for
studying clustering in many dimensions may be more suitable for the
`fine grained' data on the nearby halo that will be obtained by the
\textit{Gaia} mission (e.g. G{\'{o}}mez {et~al.} 2010).

We have measured \deltacf{} for a large catalogue of spectroscopically
confirmed BHB stars from SDSS (Xue {et~al.} 2008). We find significant
clustering in these data, particularly when we restrict the sample to
stars with galactocentric distance $r>20$~kpc. This finding of
stronger phase space correlations between stars in the outer halo is
in agreement with that of Xue {et~al.} (2011).

To test models of the accreted components of stellar haloes and
understand the effects of sample variance, we have computed \deltacf{}
for mock observations constructed from the six $\Lambda$CDM
simulations of Cooper {et~al.} (2010) in which only the stellar haloes
produced by disrupted satellites are considered. Our statistic
distinguishes quantitatively between these six qualitatively different
halo realizations. When only stars with $r>20$~kpc are considered,
five of our six simulations show statistically significant
correlations on scales in our metric equivalent to $\sim 1-8$~kpc (for
all observers on the Solar circle). Most of the models are consistent
with the Milky Way data for the outer halo, $r>30$~kpc. For the inner
halo, however, particularly at galactocentric distances smaller than
20~kpc, the simulations tend to be significantly more strongly
clustered than the data. One possible explanation for this is a
deficiency of smoothly distributed halo stars in the models, perhaps
attributable to the absence of so-called `in situ' halo stars. These
stars may be scattered from the Galactic disc, or born on eccentric
orbits (in streams of accreted gas or an unstable cooling flow, for
example). Neither of these processes are included in our model of the
accreted halo.

Although it seems reasonable to expect that in situ haloes are
distributed with spherical or axial symmetry and have a low degree of
coherence in phase space, models of such components and predictions
for the fraction of stars they contain are not well constrained. Most
hypotheses for in situ halo formation limit these stars to an `inner'
halo and predict that the accreted component (which we simulate)
dominates at larger radii
(e.g. Abadi {et~al.} 2006; Zolotov {et~al.} 2009; Font {et~al.} 2011). However, the fraction
of the halo formed in situ and the boundary between `inner' and
`outer' halo are highly model-dependent. Detections of observable
`dichotomies' in the Milky Way halo (Carollo {et~al.} 2007) are still
debated (e.g. Sch{\"{o}}nrich, Asplund \&  Casagrande 2011; Beers {et~al.} 2011). It is
possible to place broad limits on the fraction of stars in a `missing'
smooth component, for example by comparing the RMS variation of
projected star-counts in our models (Helmi {et~al.} 2011b) to the Milky
Way (Bell {et~al.} 2008). However, the uncertainties involved are substantial.

Another factor in the discrepancy between the models and the data may
be the absence of a baryonic (disc) contribution to the gravitational
potential. A massive disc could alter the process of satellite
disruption in the inner halo and might make the potential within
$30$~kpc more spherical (Kazantzidis, Abadi \&  Navarro 2010), possibly distributing
more inner halo stars into the SDSS footprint (on the other hand, a
more axisymmetric or spherical dark halo might also result in fewer
chaotic orbits, hence more coherent streams). Because of these
modelling uncertainties, our application of the \deltacf{} statistic
can presently serve only as a basic test for the abundance of
substructure in the simulations.

Several aspects of our approach could be improved.
It seems desirable to use well-measured radial velocity data to
enhance clustering signals such as our correlation function relative
to those based on configuration space data alone. However, so far, no
proposal for
including these velocity data is well-supported by theory.
Here, we have used a straightforward choice of
parametrised metric to illustrate the concept of scaling radial
velocity separations to `equivalent' configuration space separations,
and this is empirically useful in recovering a measurable
signal. Nevertheless, we have not found any compelling or generic way
to select the scaling parameter
($w_{v}$). Improving either the definition of the metric itself or the
means of fixing this parameter is a clear priority for extensions of
this approach. A similar issue affects the weighting of velocity
information in clustering algorithms (e.g. Sharma {et~al.} 2010).

Finally, further comparisons between stellar halo models and
observational data should also account for the selection effects such
spectroscopic incompleteness and the potential bias of BHB stars as a
tracer of the stellar halo (Bell {et~al.} 2010; Xue {et~al.} 2011). For
statistical analysis, there is a pressing requirement for
observational samples with well-understood selection functions, even
if they do not probe the most distant halo. The LAMOST Galactic survey
is likely to be the first to approach this goal.

In summary, we have taken a first step in adapting a well-studied
cosmological statistic, the two-point correlation function, to the
study of the Milky Way halo. Our comparisons highlight the complexity
of statistical analysis in the stellar halo, and the importance of
interpreting observational results in the context of realistic models
of halo assembly. We have compared the SDSS data with the stellar
halos produced by disrupted satellites in {\em ab initio} galaxy
formation models constructed from the Aquarius N-body simulations of
galactic dark halos in the $\Lambda$CDM cosmology. With further
refinements and more data, our statistical approach to quantifying the
smoothness of the halo can provide a practical and productive way to
study the structure of the Milky Way halo and compare with theoretical
expectations.

\section*{Acknowledgements}

The authors thank Heather Morrison and the Spaghetti Survey team for
making their data available prior to publication, and Xiangxiang Xue
and Sergey Koposov for their assistance. We thank the referee for
their helpful comments which greatly improved the structure of this
paper. APC acknowledges an STFC studentship and thanks Else
Starkenburg for useful discussions. SMC acknowledges the support of a
Leverhulme Research Fellowship. CSF acknowledges a Royal Society
Wolfson Research Merit Award and ERC Advanced Investigator grant
267291- COSMIWAY. AH acknowledges funding support from the European
Research Council under ERC-StG grant GALACTICA-24027. This work was
supported in part by an STFC rolling grant to the ICC. The
calculations for this paper were performed on the ICC Cosmology
Machine, which is part of the DiRAC Facility jointly funded by STFC,
the Large Facilities Capital Fund of BIS, and Durham University.Fig. 1
was produced with the {\tt HEALPy} implementation of {\tt HEALPix}
[http://healpix.jpl.nasa.gov, G{\'{o}}rski {et~al.} 2005].

{}


text:
\section{introduction}
Iron selenide based superconductors have come to attract particular attention after the discovery that the binary material FeSe avoids the formation of the magnetic order but becomes  superconducting at $T_c$~=~8~K. \cite{hsu2008superconductivity} The value of $T_c$ can  be enhanced both by applying pressure, reaching a maximum value of about 37~K,\cite{medvedev2009electronic,margadonna2009pressure} or in  single layer films of FeSe grown on an SrTiO$_3$ substrate, \cite{qing2012interface,FeSe2013,FeSe2015} and it is suggested that large enhancements of $T_c$ arise due to an increased  electron carrier density. \cite{FeSeCarr} Another method of enhancing the superconductivity is by  intercalating alkali metals between  FeSe layers, as in the case of $A_x$Fe$_{2-y}$Se$_2$ ($A$~=~K, Rb, Cs,Tl/Rb, Tl/K) \cite{guo2010superconductivity,2011CsSCRep,wang2011superconductivity,2011TlRbSCRep,2011TlKSCRep}. A notable difference from both the iron arsenide based materials and the bulk binary compound FeSe is the  absence of the hole pocket in the Fermi surface at the Brillouin zone center.\cite{qian2011absence} This appears to contradict the  s$_{\pm}$ superconducting state often applied to iron arsenide superconductors, where there is nesting and a sign change of the order parameter between the hole and electron pockets.\cite{mazin2008unconventional,kuroki2008unconventional} While there have been a variety of alternative proposals for the pairing state, \cite{mazin2011symmetry,khodas2012interpocket,nica2015orbital} studying the intrinsic superconducting properties is greatly complicated by the clear evidence for phase separation between non-superconducting regions with an ordered arrangement of Fe vacancies and vacancy free superconducting regions. \cite{li2012phase,PhaseSep2011DL,PhaseSepNMR}.

Recently a new iron selenide based superconductor   (Li$_{1-x}$Fe$_x$)OHFeSe ($x\approx0.2$) was discovered, \cite{lu2015coexistence} with a high transition temperature of $T_c\approx$~40~K. It has a  quasi-two-dimensional crystal structure, with layers  of both (Li$_{1-x}$Fe$_x$)OH and superconducting FeSe. The   material displays   coexistence between superconductivity and  antiferromagnetism, \cite{lu2015coexistence} while in the phase diagram superconductivity occurs in close proximity to spin-density wave order. \cite{dong2014phase}  In common with $A_x$Fe$_{2-y}$Se$_2$, the hole pocket is also absent,\cite{niu2015surface,zhao2016common} but  the samples are much more homogeneous, indicating that (Li$_{1-x}$Fe$_x$)OHFeSe is a good candidate for probing iron based superconductors without a hole pocket at the zone center. Angle resolved photoemission spectroscopy (ARPES) measurements are consistent with the presence of nodeless superconductivity with a single isotropic energy gap, but  disagree over the gap magnitudes.\cite{niu2015surface,zhao2016common} However,  scanning tunneling spectroscopy (STS) studies show two distinct features in the conductance spectra, suggesting the presence of multiple gaps. \cite{du2016scrutinizing,yan2016surface} Meanwhile the in-plane superfluid density obtained from muon-spin rotation ($\mu$SR) measurements is consistent with either one or two gaps, but very different behavior is seen in the out-of-plane component, which shows a much more rapid drop of the superfluid density with temperature. \cite{khasanov2016proximity} Furthermore,  inelastic neutron scattering (INS) measurements give evidence for the presence of a spin resonance peak, \cite{Davies2016spin,pan2016structure} consistent with a sign change of the order parameter across the Fermi surface, while a lack of a sign change was suggested from an STS study, on the basis  of quasi-particle interference (QPI) results and the effect of impurities.\cite{yan2016surface}

To further characterize the superconducting gap structure,  we report London penetration depth measurements of (Li$_{1-x}$Fe$_x$)OHFeSe single crystals  using a tunnel-diode-oscillator (TDO) based technique, from which we obtain the temperature dependence of the London penetration depth shift $\Delta\lambda(T)$. The low temperature $\Delta\lambda(T)$ gives clear evidence for nodeless superconductivity, while a single-gap isotropic $s$-wave model is unable to account for the superfluid density.  The superfluid density is well fitted by a two-gap $s$-wave model, as well as models with anisotropic gaps.

\begin{figure}[t]
\begin{center}
  \includegraphics[width=\columnwidth]{Fig1.eps}
\end{center}
	\caption{Magnetic susceptibility $\chi$ of  (Li$_{1-x}$Fe$_x$)OHFeSe samples from two batches ($\#$A and $\#$B). Both field-cooled (FC) and zero-field-cooled (ZFC) measurements were performed in an applied magnetic field of 1~mT. The inset shows the magnetization as a function of applied field at 2~K. The solid line shows the linear fit at low fields while the arrow marks the position of the lower critical field $H_{c1}$, where there is deviation from linear behavior.}
   \label{samplecharacterization}
\end{figure}
\section{Experimental details}
Single crystals of (Li$_{1-x}$Fe$_x$)OHFeSe (with $x\approx0.2$) from two batches prepared by two different groups were measured ($\#$A and $\#$B), where the crystals were synthesized following Ref.~\onlinecite{DongSynth}. Using   the parameters from   Ref.~\onlinecite{DongSynth} ($\rho_0=0.1~$m$\Omega$-cm and a carrier density of $n=1.04\times10^{21}$cm$^{-3}$), we estimate a mean free path of $l=12.4$~nm using $l=\hbar(3\pi^2)^{\frac{1}{3}}/e^2n^{\frac{2}{3}}\rho_0$. This is considerably larger than the Ginzburg-Landau coherence length $\xi_{GL}=2$~nm calculated from an upper critical field of $H_{c2}(0)=79$~T, \cite{DongSynth} and therefore the material is in the clean limit. Magnetization measurements of samples from both batches were performed utilizing a SQUID magnetometer (MPMS-5T). The temperature dependence of the London penetration depth shift $\Delta\lambda(T)~=~\lambda(T)-\lambda(0)$ was measured in a ${^3}$He cryostat from 42~K down to a base temperature of about 0.5~K using a tunnel-diode-oscillator based method, with an operating frequency of about 7~MHz. The samples were cut into a regular shape and mounted on a sapphire rod which was inserted into the coil without any contact. A very small ac field of about 2~$\mu T$ is applied to the sample along the $c$~axis, which is much smaller than the lower critical field $\mu_0H_{c1}$ and therefore the sample is always in the Meissner state. As such the   shift in the resonant frequency from zero temperature $\Delta f(T)$ is related to the penetration depth shift in the $ab$~plane $\Delta\lambda(T)$ via $\Delta\lambda(T)$~=~G$\Delta f(T)$, where the calibration constant $G$  is calculated using the geometry of the coil and sample \cite{Gfactor}.

\begin{figure}[t]
\begin{center}
  \includegraphics[width=0.8\columnwidth]{Fig2.eps}
\end{center}
	\caption{Temperature dependence of the London penetration depth $\Delta\lambda(T)$ of two samples from (a) batch $\#$A and (b) batch $\#$B. The main panels show the low temperature behavior, where the data are fitted below $T_c/3$ with  a fully gapped model and a power law dependence. The insets display the frequency shift $\Delta f(T)$ from the lowest temperature up to above the superconducting transition temperature, normalized by the value at 40~K.}
   \label{penetrationdepth}
\end{figure}

\section{Results and discussion}

 Both the field-cooled (FC) and zero-field-cooled (ZFC) magnetic susceptibility [$\chi(T)$] measurements are shown in Fig.~\ref{samplecharacterization}, from above the superconducting transition temperature down to 2~K under a small applied field of 10~Oe, where corrections for the demagnetization effect were applied. The ZFC $\chi(T)$ of both samples show sharp superconducting transitions onsetting at  around 40~K
and 37~K for samples $\#$A and $\#$B, respectively. At low temperatures the data for $\#$A  flattens, while
for  $\#$B there continues to be a slight decrease with decreasing temperature. This indicates the high quality of the single crystals from batch $\#$A, whereas those from batch $\#$B show evidence for some inhomogeneity. Meanwhile the superconducting shielding fraction is around 100$\%$ in both samples. The inset displays the field dependence of the magnetization at 2~K, for a field applied in-plane so that the demagnetization effect is very small. There is a deviation of the low field linear behavior at $\mu_0H_{c1}\approx 4.5$~mT, confirming that $H_{c1}$ is significantly larger than the ac field applied in the  penetration depth measurements.

Figure~\ref{penetrationdepth} displays $\Delta\lambda(T)$ for single crystal samples from two batches, with samples from $\#$A and $\#$B displayed in (a) and (b) respectively. The insets of both panels display the temperature dependence of the frequency shift $\Delta f(T)$ from above the superconducting transition at 42~K  down to 0.5~K. The superconducting transition onsets at respective temperatures of 40~K and 39~K in samples $\#$A and $\#$B, while the corresponding endpoints of the transition are around 34~K and 35~K, and the latter values of $T_c$ are used in the subsequent analysis of the superfluid density. The main panels of Fig.~\ref{penetrationdepth} display the low temperature behavior of $\Delta\lambda(T)$. It can be seen in Fig.~\ref{penetrationdepth}(a) that $\Delta\lambda(T)$ for the  $\#$A samples decreases with decreasing temperature before flattening below about 0.1T$_c$, indicating a nodeless gap structure in (Li$_{1-x}$Fe$_x$)OHFeSe with a lack of low energy excitations, which is consistent with previous results. \cite{niu2015surface,zhao2016common,du2016scrutinizing,yan2016surface,khasanov2016proximity} Furthermore when fitted with a power law dependence  $\Delta\lambda(T)\sim T^n$, exponents of $n=2.83$ and $n=2.46$ are obtained for samples  $\#$A-1 and  $\#$B-1 respectively. In the case of nodal superconductors, $n=1$ for line nodes and $n=2$ for point nodes is generally anticipated. While impurity scattering, non-local effects and quantum fluctuations can all lead to  $n\approx2$ at low temperatures for $d$-wave superconductors with line nodes,\cite{Hirschfeld1993,Kosztin1997,Benfatto2001} our observation that $n$ is significantly larger than two gives further evidence for fully gapped behavior. For a fully-gapped superconductor at $T\ll T_c$, the penetration depth can be described by $\Delta\lambda(T)=\lambda_{eff}(0)\sqrt{\pi\Delta(0)/2k_BT}\textrm{exp}[-\frac{\Delta(0)}{k_BT}]$, where $\Delta(0)$ is the superconducting gap magnitude at zero temperature and $\lambda_{eff}(0)$ is an effective zero temperature penetration depth. The low temperature data  for sample $\#$A-1 is well fitted with a gap magnitude of  $\Delta(0)=~0.87k_BT_c$. A similar  value of $\Delta(0)=~0.78k_BT_c$ is obtained from the fitting for  sample $\#$B-1, although there is a small deviation of the  fitted curve for this sample at the lowest temperatures, where a weak anomaly is observed in the data, the origin of which  is not clear. The values of the fitted gaps are much smaller than the value from BCS theory of $1.76k_BT_c$ for weakly coupled isotropic $s$-wave superconductors, suggesting  multi-gap superconductivity or gap anisotropy in (Li$_{1-x}$Fe$_x$)OHFeSe. The fitted values of  $\lambda_{eff}(0)$ of 636~\AA\ and 1228~\AA\ for samples  $\#$A-1 and  $\#$B-1 respectively, are significantly different from the value of $\lambda(0)\approx2800$~\AA\ from  $\mu$SR measurements. \cite{khasanov2016proximity}. Such a difference is also expected for multi-gap or anisotropic superconductors. \cite{Malone2009}

In order to obtain further information about the superconducting pairing state, the normalized superfluid density [$\rho_s(T)$] was calculated  from the penetration depth using  $\rho_s(T)$~=~$[\lambda(0)$/$\lambda(T)]^2$, where $\lambda(0)\approx$~2800~\AA\ was estimated from  $\mu$SR measurements.\cite{khasanov2016proximity} Since the samples from batch  $\#$A are  higher quality and show a sharper superconducting transition, measurements from this batch were used in the subsequent analysis. The $\rho_s(T)$ of sample $\#$A-1 is displayed in Fig.~\ref{superfluiddensity}.  The superfluid density was modelled using
\begin{equation}
\rho_{\rm s}(T) = 1 + 2 \left\langle\int_{\Delta_k}^{\infty}\frac{E{\rm d}E}{\sqrt{E^2-\Delta_k^2}}\frac{\partial f}{\partial E}\right\rangle_{\rm FS},
\label{equation2}
\end{equation}
\noindent where $f(E, T)=[1+{\rm exp}(E/k_BT)]^{-1}$ is the Fermi function and $\left\langle\ldots\right\rangle_{\rm FS}$ denotes an average over the Fermi surface. The temperature dependence of the gap function $\Delta_k$ is approximated by \cite{carrington2003magnetic}
\begin{equation}
\delta(T)={\rm tanh}\left\{1.82\left[1.018\left(T_c/T-1\right)\right]^{0.51}\right\},
\label{equation3}
\end{equation}

\begin{figure}[t]
\begin{center}
  \includegraphics[width=\columnwidth]{Fig3.eps}
\end{center}
	\caption{Normalized superfluid density $\rho_s$ as a function of the reduced temperature $T/T_c$ for sample $\#$A-1 of (Li$_{1-x}$Fe$_x$)OHFeSe. The dashed, dashed-dotted and solid lines show fits to models with a single isotropic $s$-wave gap, two isotropic gaps and the $s\times\tau_3$ state respectively.   The inset shows  $\rho_s$ upon varying $\lambda(0)$ by $\pm30\%$, along with fits to a two-gap model.}
   \label{superfluiddensity}
\end{figure}

\noindent As discussed previously, the behavior of $\Delta\lambda(T)$ at low temperatures and previous experimental results indicate nodeless superconductivity in (Li$_{1-x}$Fe$_x$)OHFeSe and therefore we fitted $\rho_s(T)$ with various fully-gapped models. The simplest model is to assume an isotropic superconducting gap, with $\Delta_k(T)~=~\Delta(0)\delta(T)$. The fitted curve for such an isotropic single band $s$-wave model with $\Delta(0)~=~1.72k_BT_c$ is  shown by the dashed line in Fig.~\ref{superfluiddensity}, where although there is reasonable agreement above 0.5$T_c$, there is a significant discrepancy in the intermediate temperature range between 0.2$T_c$ and 0.5$T_c$. This difference arises due to the data dropping more quickly than the calculated $\rho_s(T)$, indicating the presence of an additional lower energy scale in the gap structure, which is consistent with the smaller gap value obtained from fitting  $\Delta\lambda(T)$ at low temperatures. Therefore the data are fitted with a two-gap model, as has been applied to many iron based superconductors. For this model the total superfluid density is given by the weighted sum of two components with different gaps,
\begin{equation}
\rho_{\rm s}(T) = \alpha\rho{_{\rm s}^1}(\Delta_k^1, T) + (1-\alpha)\rho{_{\rm s}^2}(\Delta_k^2, T),
\label{equation4}
\end{equation}

\noindent where $\rho_s^i (i=1, 2)$ is the normalized superfluid density with a gap function $\Delta_k^i (i=1, 2)$  and $\alpha$ is the relative weight for the component $\rho_s^1$ ($0\leq\alpha\leq1$).  The fitting result is shown in Fig.~\ref{superfluiddensity}  by the dashed-dotted  line, which agrees well with the data across the whole temperature range.  The fitted parameters are $\Delta_1(0)=0.8k_BT_c$, $\Delta_2(0)=1.9k_BT_c$ and $\alpha=0.13$, where the value of the small gap is close to that found from fitting $\Delta\lambda(T)$, which is further consistent with two-gap superconductivity. In order to consider possible uncertainties in the calibration constant $G$ or $\lambda(0)$, in the inset we show $\rho_s(T)$  upon varying $\lambda(0)$ by $\pm30\%$. The data can still be fitted by a two gap model with slightly different parameters of $\Delta_1(0)=0.8k_BT_c$, $\Delta_2(0)=1.6k_BT_c$ and $\alpha=0.15$ for  $\lambda(0)$~=~2000~\AA ~and $\Delta_1(0)=0.8k_BT_c$, $\Delta_2(0)=2.1k_BT_c$ and $\alpha=0.12$ for  $\lambda(0)$~=~3600~\AA. The data were also well fitted with an anisotropic single band model with $\Delta_k(T,\phi)=\Delta(0)(1+r{\rm cos}2\phi)\delta(T)$, using $\Delta(0)=1.32k_BT_c$ and $r=0.65$, which is not shown for the sake of clarity.

ARPES measurements indicate that the Fermi surface consists of electron pockets at the Brillouin zone corners, without the presence of hole pockets at the zone center \cite{niu2015surface,zhao2016common}. From a recent STM study it was proposed on the basis of QPI measurements, as well as examining the effects of impurities, that there is no sign change of the superconducting gap across the Fermi surface \cite{yan2016surface}, in which case such a two gap $s$-wave model readily explains the data. However INS measurements show evidence for a spin resonance peak,\cite{Davies2016spin,pan2016structure} which indicates that there is a sign change of the order parameter between regions of the Fermi surface connected by the resonance wave vector. This would be incompatible with two-gap $s$-wave superconductivity with no sign change and are also difficult to account for with the $s_{\pm}$ state proposed for  many iron arsenide superconductors, due to both the lack of a hole pocket at the zone center and a different nesting wave vector for the spin resonance. \cite{mazin2008unconventional,kuroki2008unconventional} On the other hand a different sign changing $s_{\pm}$ state has been suggested for $A_x$Fe$_{2-y}$Se$_2$  (A~=~K, Rb, Cs), \cite{khodas2012interpocket,mazin2011symmetry} with a very similar Fermi surface to  (Li$_{1-x}$Fe$_x$)OHFeSe.

Another proposed pairing state for nodeless sign changing superconductivity in iron based superconductors with only electron pockets is an orbital selective $s\times\tau_3$ state. \cite{nica2015orbital} In this scenario, intraband pairing has $d_{x^2-y^2}$ symmetry,  while the  interband pairing has $d_{xy}$ symmetry. As a result, the zeroes of the intraband and interband pairing  are offset by an angle of $\pi/4$ and therefore the gap remains nodeless. A simple model for the gap function of this state taking into account the Fermi surface is  $\Delta_k(T,\phi)$=[($\Delta_1(0)$)$^2$+($\Delta_2(0)$sin($\phi$))$^2$]$^{1/2}$$\delta(T)$ \cite{nica2015orbital,SiPrivate}. As shown in Fig.~\ref{superfluiddensity}, this model can also well fit the experimental data, with fitted parameters of $\Delta_1(0)=1.05k_BT_c$ and $\Delta_2(0)=3.2k_BT_c$. In this case the gap minimum $\Delta_1(0)$ is slightly larger than the smaller gap from the two-gap $s$-wave fit. Therefore the data can be well accounted for by fitting with models with either two-gaps or a strong gap anisotropy and our results are compatible with two-gap behavior, anisotropic $s$-wave superconductivity or an orbital selective $s\times\tau_3$ state. It is often very difficult to distinguish between scenarios with multiple gaps and those where there is one anisotropic gap. While the isotropic nature of the gap inferred from ARPES would favor the two-gap scenario over an anisotropic gap, this still requires further study. \cite{niu2015surface,zhao2016common} Furthermore, the superfluid density is only sensitive to the gap magnitude rather than the phase and  different measurements are required to clarify the presence of a sign change  and to determine the nature of the pairing state.

We note that there have been conflicting reports about the nature of the gap structure of (Li$_{1-x}$Fe$_x$)OHFeSe, with only a single gap being resolved from ARPES measurements \cite{niu2015surface,zhao2016common}, while two gaps are found from STS \cite{du2016scrutinizing,yan2016surface} and the in-plane superfluid density from $\mu$SR  measurements is compatible with both single-gap and two-gap models. \cite{khasanov2016proximity} The gap values we obtain from fitting the in-plane superfluid density are smaller than those reported from previous measurements, particularly in the case of the smaller gap. Evidence for this smaller gap is clearly observed from our measurements of $\Delta\lambda(T)$ and $\rho_s$ at low temperatures, which may be a result of the high sensitivity of the TDO-based technique. A further reason for the discrepancies could be due to the non-stoichiometric nature of (Li$_{1-x}$Fe$_x$)OHFeSe, where the exact composition and homogeneity may influence the doping level, $T_c$ and the gap magnitudes. In addition, the doping level can be different between the surface and in the bulk, \cite{niu2015surface} and therefore probes which primarily measure the surface properties may give different results. Indeed different results have been found from different measurements of the gap structure of other iron-based superconductors, \cite{evtushinsky2009momentum} and significant sample dependence has been suggested for Ba$_{1-x}$K$_x$Fe$_2$As$_2$.\cite{hashimoto2009microwave}

\section{summary}

We have measured the temperature dependence of the London penetration depth $\Delta\lambda(T)$ and the derived superfluid density $\rho_s(T)$ of the recently discovered high-$T_c$ iron-based superconductor (Li$_{1-x}$Fe$_x$)OHFeSe. The behavior of $\Delta\lambda(T)$ at low temperatures gives clear evidence for nodeless superconductivity, with a relatively small energy gap, while the analysis of $\rho_s(T)$ is consistent with both two-gap superconductivity and models with significant gap anisotropy such as an orbital selective $s\times\tau_3$ or anisotropic $s$-wave  state.

We thank Q. Si, E. M. Nica, R. Yu and X. Lu for helpful discussions and valuable suggestions. This work was supported by the National Natural Science Foundation of China (No.11574370), National Key Research and Development Program of China (No. 2016YFA0300202), and the Science Challenge Project of China.

redpajama_set_name: RedPajamaArXiv


text:
\section{Kronecker tensor structures for complex fields}
\label{App:ProjA}
In this appendix we demonstrate how a product of two operators can be decomposed onto irreps of $U(m)\times U(n)$ in the picture where we work with complex fields. For simplicity we will start with just $U(n)$. The generalization to $U(m)\times U(n)$ is trivial. The main utility of working with complex fields is that the projectors take their simplest form possible, namely as combinations of Kronecker deltas with the least number of indices possible. The real field picture can also have its projectors expressed in terms of Kronecker deltas, albeit at the cost of more indices. We hope that presenting our projectors in three different pictures will make our work more intuitive. In order to capture all irreps that can appear in the real field notation, and hence not miss any information in the bootstrap algorithm, we need to consider two OPEs, namely $\Phi_i^\dagger \times \Phi_j$ and $\Phi_i \times \Phi_j$. Below we show how they may be decomposed onto irreps:
\begin{equation}
\begin{split}
  \Phi_i^\dagger \times \Phi_j &\sim \Big(\Phi_i^\dagger \Phi_j - \frac{1}{n}\delta_{ij}  \Phi_k^\dagger \Phi_k \Big) +\frac{1}{n}\delta_{ij} \Phi_k^\dagger \Phi_k\,, \\
  \Phi_i \times \Phi_j &\sim (\Phi_i  \Phi_j + \Phi_j \Phi_i) + (\Phi_i  \Phi_j - \Phi_j \Phi_i)\,.
\end{split}
\label{complexOPE}
\end{equation}
The first line in \eqref{complexOPE} shows the decomposition into the adjoint ($R$) and singlet ($S$) representations, where as the second line shows the decomposition into the symmetric ($T$) and antisymmetric representations ($A$). To read off the projectors from \eqref{complexOPE} it is useful to remember
\begin{equation}
\begin{aligned}
  O^X_{ij}  \sim P^X_{ijkl} \Phi_k \Phi_l\,,\qquad
  O^X_{ij}  \sim P^X_{ijkl} \Phi^\dagger_k \Phi_l\,,
\end{aligned}
\label{projdefinitionscomplex}
\end{equation}
where $O^X$ is the exchanged field in some irrep, e.g.\ $O^T_{12} \sim \Phi_1 \Phi_2 + \Phi_2 \Phi_1$. The first relation in \eqref{projdefinitionscomplex} can be used to read off the $T$ and $A$ projectors, whereas the second can be used for $S$ and $R$. Notice that we have implicitly assumed that fields are inserted at different positions in order for antisymmetric irreps to not vanish identically. The projectors can be read off as\footnote{Note that we take the correlator to be $\langle \Phi_i^\dagger \Phi_j \Phi_k \Phi_l^\dagger \rangle$ which is why the projector of the adjoint representation is equal to $P^R_{ijkl} =\delta_{ik}\delta_{jl}-\frac{1}{n}\delta_{ij}\delta_{kl}$ instead of $P^R_{ijkl} =\delta_{il}\delta_{jk}-\frac{1}{n}\delta_{ij}\delta_{kl}$.}
\begin{align}
  P^S_{ijkl} &=\frac{1}{n}\delta_{ij}\delta_{kl}\,,\qquad
  P^R_{ijkl} =\delta_{ik}\delta_{jl}-\frac{1}{n}\delta_{ij}\delta_{kl}\,, \\
  P^T_{ijkl} &=\frac{1}{2}(\delta_{ik}\delta_{jl} +\delta_{il}\delta_{jk})\,,\qquad
  P^A_{ijkl}= \frac{1}{2}(\delta_{ik}\delta_{jl} -\delta_{il}\delta_{jk})\,.
\label{complexprojectorsUn}
\end{align}
The dimensions of the corresponding irreps are $(1,(n-1)(n+1),\frac12n(n+1),\frac12n(n-1))$. As the reader may have observed from the main text, or the next appendix, when going from the complex field picture to the real field picture the above dimensions get multiplied by a factor of $2$. For example $\frac{n(n+1)}{2}$ becomes $n(n+1)$. This is because each of the initial $\frac{n(n+1)}{2}$ complex elements contains two real elements. The last step now is to write down the projectors for $U(m) \times U(n)$. This is trivial, in the sense that they are just products of $U(m)$ with $U(n)$ projectors. We have
\begin{equation}
\begin{aligned}
  P^{S}_{ijklmnop} &=P^S_{ijkl}P^S_{mnop} \,,\quad
  P^{RS}_{ijklmnop} =P^R_{ijkl}P^S_{mnop} \,,\quad
  P^{SR}_{ijklmnop} =P^S_{ijkl}P^R_{mnop} \,,\\
  P^{RR}_{ijklmnop} &=P^R_{ijkl}P^R_{mnop} \,,\quad
  P^{TT}_{ijklmnop} =P^T_{ijkl}P^T_{mnop} \,,\quad
  P^{TA}_{ijklmnop} =P^T_{ijkl}P^A_{mnop} \,,\\
  P^{AT}_{ijklmnop} &=P^A_{ijkl}P^T_{mnop} \,,\quad
  P^{AA}_{ijklmnop} =P^A_{ijkl}P^A_{mnop}\,.
\end{aligned}
\label{complexprojectorsUnII}
\end{equation}
The sum rules that can be derived with the above projectors can be checked to be completely equivalent to those derived from the projectors of real fields outlined in the main text. Another observation is that, compared to real fields, we do not need separate projectors for even and odd spins.

\section{Kronecker tensor structures for real fields}
\label{App:ProjB}
The projectors that correspond to a four-point function of real fields can be intuitively presented in terms of Kronecker deltas if we add an additional index. This form is useful since one may directly extract the form of exchanged operators, as we will show. The form of exchanged operators is useful to know since it can guide us with respect to assumptions we may impose. Also, we expect it to be easier to work with in a mixed correlator system. We start by labeling the real and complex parts of an operator $\Phi_i$ with an index (we start with $U(n)$ for simplicity)
\begin{equation}
\Phi_i = \phi^1_{i} + i\hspace{1pt}\phi^2_{i}\,,
\label{complextoreal}
\end{equation}
where the upper case $\Phi$ denotes the complex operator and the lower case $\phi$ denote real fields. We must now simply plug in \eqref{complextoreal} to the expressions for the representations of the previous appendix. For simplicity we will do this for the singlet representation, and then quote the results for rest of the representations. Note that implicitly we consider the two external fields of the OPE at different positions, for otherwise the antisymmetric combinations would vanish identically. We have
\begin{equation}
  \Phi_i^\dagger \Phi_i = (\phi^1_i \phi^1_i + \phi^2_i \phi^2_i ) + i (\phi^1_i \phi^2_i - \phi^2_i \phi^1_i)\,,
\end{equation}
where the first parenthesis corresponds to what was called $S_{\text{even}}$ in the main text, and the second parenthesis corresponds to what was called $S_{\text{odd}}$. As expected $S_{\text{odd}}$ vanishes identically if we don't insert powers of derivatives between the operators. The projectors are now very straightforward to write down by recalling the relation
\begin{equation}
  O^X_{ij;ab} = P^X_{ijkl;abcd}\hspace{1pt} \phi^a_i \phi^b_j\,,
\label{projrelation}
\end{equation}
where $X$ stands for some specific irrep and indices from the beginning of the latin alphabet take the values $1,2$. Notice that \eqref{projrelation} is simply the statement that projectors must project products of operators onto irreps. We have
\begin{equation}
\begin{split}
  P^{S_{\text{even}}}_{ijkl;abcd} &= \frac{1}{2n} \delta_{ab}\delta_{cd}\delta_{ij}\delta_{kl}\,, \\
  P^{S_{\text{odd}}}_{ijkl;abcd}  &= \frac{1}{2n}(\delta_{ac}\delta_{bd}-\delta_{ad}\delta_{bc})\delta_{ij}\delta_{kl}\,.
\end{split}
\end{equation}
Indeed, one may confirm that, for example,
\begin{equation}
  O^{S_{\text{even}}}_{11;11}\sim(\phi^1_i \phi^1_i + \phi^2_i \phi^2_i ) \sim  P^{S_{\text{even}}}_{11kl;11cd} \hspace{1pt}\phi^{c}_k \phi^{d}_l\,.
\end{equation}
This procedure can be repeated for the rest of the irreps. The resulting projectors are
\begin{equation}
\begin{split}
  P^{S_{\text{even}}}_{ijkl;abcd} &= \frac{1}{2n} \delta_{ab}\delta_{cd}\delta_{ij}\delta_{kl}\,, \\
  P^{S_{\text{odd}}}_{ijkl;abcd}  &= \frac{1}{2n}(\delta_{ac}\delta_{bd}-\delta_{ad}\delta_{bc})\delta_{ij}\delta_{kl}\,,\\
  P^{R_{\text{even}}}_{ijkl;abcd}  &=  \tfrac{1}{2}\delta_{ab}\delta_{cd}\Big(\delta_{ik}\delta_{jl}-\frac{1}{n}\delta_{ij}\delta_{kl}\Big)\,,\\
  P^{R_{\text{odd}}}_{ijkl;abcd}  &= \tfrac{1}{2}(\delta_{ac}\delta_{bd}-\delta_{ad}\delta_{bc})\Big(\delta_{ik}\delta_{jl}-\frac{1}{n}\delta_{ij}\delta_{kl}\Big)\,,\\
  P^{T_{\text{even}}}_{ijkl;abcd}  &= \tfrac{1}{4}(\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc}-\delta_{ab}\delta_{cd})(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jl})\,,\\
  P^{A_{\text{odd}}}_{ijkl;abcd}  &= \tfrac{1}{4}(\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc}-\delta_{ab}\delta_{cd})(\delta_{ik}\delta_{jl}-\delta_{il}\delta_{jl})\,.
\end{split}
\end{equation}
The dimensions of the corresponding irreps are $(1,1,(n-1)(n+1),(n-1)(n+1),n(n+1),n(n-1))$.

Using the above expressions, it is now trivial to write down the $U(m)\times U(n)$ projectors:
\begin{align}
  P^{S_{\text{even}}}_{ijklmnop;abcdefgh} &= P^{S_{\text{even}}}_{ijkl;abcd}\hspace{1pt} P^{S_{\text{even}}}_{mnop;efgh}+ P^{S_{\text{odd}}}_{ijkl;abcd}\hspace{1pt} P^{S_{\text{odd}}}_{mnop;efgh}\,,\nonumber\\
  P^{S_{\text{odd}}}_{ijklmnop;abcdefgh} &= P^{S_{\text{even}}}_{ijkl;abcd}\hspace{1pt} P^{S_{\text{odd}}}_{mnop;efgh}+ P^{S_{\text{odd}}}_{ijkl;abcd}\hspace{1pt} P^{S_{\text{even}}}_{mnop;efgh}\,,\nonumber\\
  P^{RS_{\text{even}}}_{ijklmnop;abcdefgh} &= P^{R_{\text{even}}}_{ijkl;abcd}P^{S_{\text{even}}}_{mnop;efgh} + P^{R_{\text{odd}}}_{ijkl;abcd}P^{S_{\text{odd}}}_{mnop;efgh}\,,\nonumber\\
  P^{RS_{\text{odd}}}_{ijklmnop;abcdefgh} &= P^{R_{\text{even}}}_{ijkl;abcd}P^{S_{\text{odd}}}_{mnop;efgh} + P^{R_{\text{odd}}}_{ijkl;abcd}P^{S_{\text{even}}}_{mnop;efgh}\,,\nonumber\\
  P^{SR_{\text{even}}}_{ijklmnop;abcdefgh} &= P^{S_{\text{even}}}_{ijkl;abcd}P^{R_{\text{even}}}_{mnop;efgh} + P^{S_{\text{odd}}}_{ijkl;abcd}P^{R_{\text{odd}}}_{mnop;efgh}\,,\nonumber\\
  P^{SR_{\text{odd}}}_{ijklmnop;abcdefgh} &= P^{S_{\text{even}}}_{ijkl;abcd}P^{R_{\text{odd}}}_{mnop;efgh} + P^{S_{\text{odd}}}_{ijkl;abcd}P^{R_{\text{even}}}_{mnop;efgh}\,,\nonumber\\
  P^{RR_{\text{even}}}_{ijklmnop;abcdefgh} &= P^{R_{\text{even}}}_{ijkl;abcd}P^{R_{\text{even}}}_{mnop;efgh} + P^{R_{\text{odd}}}_{ijkl;abcd}P^{R_{\text{odd}}}_{mnop;efgh}\,,\nonumber\\
  P^{RR_{\text{odd}}}_{ijklmnop;abcdefgh} &= P^{R_{\text{even}}}_{ijkl;abcd}P^{R_{\text{odd}}}_{mnop;efgh} + P^{R_{\text{odd}}}_{ijkl;abcd}P^{R_{\text{even}}}_{mnop;efgh}\,,\nonumber\\
  P^{TT_{\text{even}}}_{ijklmnop;abcdefgh} &= P^{T_{\text{even}}}_{ijkl;abcd}P^{T_{\text{even}}}_{mnop;efgh}\,,\nonumber\\
  P^{TA_{\text{odd}}}_{ijklmnop;abcdefgh} &= P^{T_{\text{even}}}_{ijkl;abcd}P^{A_{\text{odd}}}_{mnop;efgh}\,,\nonumber\\
  P^{AT_{\text{odd}}}_{ijklmnop;abcdefgh} &= P^{A_{\text{odd}}}_{ijkl;abcd}P^{T_{\text{even}}}_{mnop;efgh}\,,\nonumber\\
  P^{AA_{\text{even}}}_{ijklmnop;abcdefgh} &= P^{A_{\text{odd}}}_{ijkl;abcd}P^{A_{\text{odd}}}_{mnop;efgh}\,.
\end{align}
From these expressions we can also see explicitly that when $m=n$, if we choose to consider the two $U(n)$ symmetries as indistinguishable (which we remind the reader is not strictly necessary), the $RS$ irreps become the same as the $SR$ irreps. The same also happens for $TA$ and $AT$.

\section{Numerical parameters}
\label{App:C}
For most of our plots, the bounds are obtained with the use of
\texttt{PyCFTBoot}~\cite{Behan:2016dtz} and \texttt{SDPB}~\cite{Landry:2019qug}. We use the numerical parameters $\texttt{m\_max}=6, \texttt{n\_max}=9, \texttt{k\_max}=36$ in \texttt{PyCFTBoot}, and we include spins up to $\texttt{l\_max}=26$. The binary precision for the produced \texttt{xml} files is 896 digits.  \texttt{SDPB} is run with the options \texttt{-}\texttt{-precision=896}, \texttt{-}\texttt{-detectPrimalFeasibleJump}, \texttt{-}\texttt{-detectDualFeasibleJump} and default values for other parameters. We refer to this set of parameters as ``\!\emph{A}''. Unless otherwise stated, our plots are ran with parameters ``\!\emph{A}''.

For some of the plots we used  $\texttt{m\_max}=5, \texttt{n\_max}=7, \texttt{l\_max}=36,  \texttt{k\_max}=42$ and
$\texttt{m\_max}=6, \texttt{n\_max}=9, \texttt{l\_max}=36,  \texttt{k\_max}=42$, referred to as ``\!\emph{B}'' and ``\hspace{-1pt}\emph{C}\hspace{1pt}'' respectively. Lastly, we also used \texttt{qboot} \cite{Go:2020ahx}, with $\Lambda=15$, $\texttt{n\_max}=500$, $\nu\texttt{\_max}=25$ and $l=\{0\text{--}49,55,56, 59, 60, 64, 65, 69, 70, 74, 75,\\ 79, 80, 84, 85, 89, 90\}$ referred to as ``\hspace{-1pt}\emph{D}\hspace{0.5pt}''.

\end{appendices}