\section{Introduction}
\label{sec:introduction}
Dwarf satellite galaxies provide a powerful laboratory for studying physical processes in galaxy formation and 
the assembly history of the Milky Way (MW) halo \citep{1999ApJ...522...82K,2017ARA&A..55..343B}, 
making a stringent test of one of the most important model predictions in the $\Lambda$CDM structure formation scenario on small scales, 
such as a universal dark matter density profile \citep{1994Natur.370..629M,1997ApJ...490..493N,1997ApJ...477L...9F}, 
and exploring the particle nature of dark matter via a search for their annihilation and decay signals 
\citep{1990Natur.346...39L,2013pss5.book.1039W,2015ApJ...801...74G,2016MNRAS.461.2914H}.

Sextans dwarf spheroidal 
(hereafter just Sextans) is one of the most luminous ``classical'' dwarf galaxy satellites of the MW, at a heliocentric distance of 86~kpc \citep{1990MNRAS.244P..16I,1991AJ....101..892M,1995MNRAS.277.1354I,1995AJ....110.2166M,2003AJ....126.2840L,10.1093/mnras/stx086}. 
Sextans displays a larger spatial extent 
\citep{1995MNRAS.277.1354I,10.1093/mnras/stw949,10.1093/mnras/stx086,2018ApJ...860...66M}, 
which might be a signature of tidal disruption, although any clear signs of tidal tails and S-shaped contours have not been found in Sextans.  
From studies on the internal kinematics,  
Sextans is found to be a dark matter dominated system; 
the total mass is about $10^8M_\odot$, with a very large mass-to-light ratio ($M/L$) in the range of $M/L\sim 500$--$ 900$ 
depending on the aperture radii
\citep{2007ApJ...667L..53W,2009AJ....137.3100W,2008Natur.454.1096S,2011MNRAS.411.1013B}.
Thus Sextans is one of the most important targets for exploring the nature of dark matter \citep[e.g.][]{2020ApJ...904...45H}.

On the other hand, the determination of the orbital motions of the dwarf galaxies allows us to infer the mass profile 
of the MW and to estimate the orbital histories of the satellites themselves
\citep{2018AA...619A.103F,McConnachie_2020}. 
To obtain the orbital information of each dwarf galaxy, 
we need both systemic radial and tangential velocities estimated from the velocities of member stars, 
in addition to their spatial positions with respect to the MW centre.
Among these phase information, proper motions are most expensive to be measured precisely, 
especially for faint stars out to large heliocentric distances that are indeed the case for Sextans at 86~kpc. 
It usually requires a long time baseline, accurate astrometry, and careful control of systematic errors. 
$Gaia$ \citep{2018A&A...616A...1G} has transformed this situation: 
precise proper motion measurements are now readily available for bright stars. 
For Sextans, many works used the $Gaia$ proper motion measurements for member star candidates 
on the red giant branch and blue horizontal branch (and blue stragglers) to infer the systemic proper motion of Sextans
\citep{2018AA...616A..12G,2018AA...619A.103F,McConnachie_2020,Li_2021,10.1093/mnras/stab1568}. 
Some of these studies, even using the same $Gaia$ datasets, show discrepancies depending on which selection of member stars to be used. 
However, these bright stars cover only a small portion of the member populations, 
e.g. compared to main sequence stars, so it is still important to measure 
the systemic proper motion of Sextans from a majority of members, i.e. fainter stars including main-sequence stars.

Hence the purpose of this paper is to study structural parameters and systemic proper motion of Sextans using
the Subaru Hyper Suprime-Cam (HSC) \citep{2018PASJ...70S...1M,2018PASJ...70S...2K} datasets, which were taken 
with multiple passbands at different epochs, about 2.66~years time baseline on average. 
Thanks to its superb image quality, large aperture,  and wide field-of-view, 
the Subaru HSC data allows a clean discrimination of stars from galaxies 
and a secure selection of Sextans member star candidates, based on the colour-magnitude cuts,  
down to faint magnitudes $i\sim 24$. 
We can then use an unprecedentedly large number of the selected stars to accurately estimate the structural parameters. 
Even if the Subaru HSC is a ground-based dataset, 
we will show that the data allows us to estimate the systemic proper motion of Sextans,
using the method developed in \citet{2021MNRAS.501.5149Q} for the HSC data. 
To do these we use stars sufficiently outside Sextans to estimate the number density and systemic proper motions of foreground MW halo stars, to minimize the contamination of foreground stars to our measurements. 
Our study gives a demonstration that a ground-based dataset can be used to study the proper motions of stars, globular clusters and 
dwarf galaxies if a longer time baseline and a control of systematic errors are available. 
For example, our study is relevant to  the upcoming Rubin Observatory's 
LSST\footnote{\url{https://www.lsst.org}}.   

This paper is structured as follows. 
We describe details of the Subaru HSC data used in this paper in Section~\ref{sec:data}.
In Section~\ref{sec:dataprocess}, we describe details of the data processing: star/galaxy separation, matching objects between the two datasets, 
and the recalibration of astrometric solutions for the proper motion measurement.  
In Section~\ref{sec:results}, we show the main results of this paper, estimations of the structural parameters and the systemic proper motions of Sextans. 
We conclude in Section~\ref{sec:conclusion}. 
\section{Data}
\label{sec:data}
\begin{table}
    \begin{center}
  \begin{tabular}[t]{l|cc} \hline\hline
  Type & Number of objects \\\hline
  HSC-SSP \\
  \hspace{1em} all objects $(i_\mathrm{PSF} < 24)$ &1,444,206 \\ 
  \hspace{1em} stars used in this work  & 210,506 \\
  \hspace{1em} galaxies used in this work & 406,220 \\ \hline 
  HSC-PI\\
  \hspace{1em} all objects $(i_\mathrm{PSF} < 24.5)$ & 2,743,551 \\ 
  \hspace{1em} stars used in this work & 397,309 \\ 
  \hspace{1em} galaxies used in this work  & 879,791 \\ \hline 
  Matched stars for proper motion measurements &121,015 \\
  Matched galaxies used for calibration & 215,325 \\ \hline\hline
        \end{tabular}  
\end{center}
\caption{
The number of objects used in this paper. 
``All objects'' denotes the total number of detected objects with $i_{\rm PSF}<24$ or $24.5$ 
for the HSC-SSP and HSC-PI catalogues, respectively. 
The deeper magnitude cut for the HSC-PI data is taken because 
the HSC-PI data is shallower than the HSC-SSP data in that we do not miss any object, 
which exists in the HSC-SSP catalogue, from the HSC-PI catalogue due to larger photometric errors.
``Stars'' and ``galaxies'' in each catalogue are defined by objects that are 
selected by the size threshold based on $i_{\rm cModel}-i_{\rm PSF}$ and the colour-magnitude cuts (see text for details).
The last rows denote the number of the matched stars used for proper motion measurements 
and that of the matched galaxies used for calibration.}
\label{tab:numberofobjects}
\end{table}
In this section, we describe the details of the data that we use in this paper.
To measure the proper motions of stars, we use the imaging data taken at two different epochs, 
which allow us to infer an angular displacement for each star over the time interval, therefore, its proper motion.
All the images we use were taken by the Subaru Hyper Suprime-Cam (HSC), as we describe below.
The HSC camera at the prime focus of the 8.2~m Subaru Telescope is a wide-field imaging camera with
1.77~deg$^2$ field-of-view \citep{2018PASJ...70S...1M,2018PASJ...70S...2K,2018PASJ...70S...3F,2018PASJ...70...66K}.
We summarise the number of objects used in this paper in Table \ref{tab:numberofobjects}. 
\subsection{HSC-SSP}
\label{sec:hsc-ssp}
\begin{figure}
    \centering
    \includegraphics[width=0.95\hsize]{img/visit_delta.pdf}
    \caption{
    The HSC-SSP co-add images are constructed from the data taken by multiple visits to each field. 
    The co-addition is done based on each HSC “tract” region.
    colour in each grid means the number of HSC tracts lying in a
    given value of the 2D space, parameterised by the number of visits of each tract ($y$-axis) and the separation time between the first and last visits of the tract ($x$-axis). Based on this information of epochs and tract, we compute the mean epoch of the co-add image for each object in HSC-SSP, which is used to measure proper motion.
    }
    \label{fig:visitdelta}
\end{figure}
The first HSC data we use is from the HSC Subaru Strategic Program (HSC-SSP) survey that started in 2014 
and was recently completed in the early 2022 \citep{10.1093/pasj/psx066}. 
The HSC-SSP survey used 330 Subaru nights to conduct a five-band ($grizy$) wide-area imaging survey over about $1,200$~deg$^2$.
The total exposure time for each field is 10~min for $gr$ and 20~min for $izy$, respectively.
Among the five filters, the $i$-band images were taken when both the weather and the seeing were good (typically $0\farcs6$ seeing FWHM) \citep{10.1093/pasj/psx066,10.1093/pasj/psx130}. 
We use the catalogue in about $20$~deg$^2$ field containing Sextans at the centre, which is taken
from the internal S21A data release product based on the data taken from March 2014 to January 2021.
The HSC-SSP catalogue consists of the primary photometric sources in $g,r,i$ filters with $i_\mathrm{PSF} < 24$ and contains about 1.78 million objects.
The sample selection query is provided in Appendix~\ref{ap:sspsql}.
We use the co-add images, which are based on different exposures taken at different epochs \citep{2018PASJ...70S...5B}. 
All single exposures overlapped with a particular ``tract'' are projected to an output plane centred on the tract. 
Then the co-add images are produced by the weighted average of the tract \citep{2018PASJ...70S...5B}.
For each co-add source, we calculate its mean modified Julian date (MJD) based on epochs of individual exposures (obviously, this is important for the HSC-SSP). 
Fig.~\ref{fig:visitdelta} shows the number distribution of the HSC-SSP co-add images lying in a 2D space, 
parameterised by the number of visits to each tract and the separation time between the first and last visits
(star/galaxy definition will be described in Section~\ref{sec:dataprocess}).
The HSC-SSP co-add images consist of 1 to 16 exposures, and most of them are taken within a few days or a year.
The details of the analysis pipeline used for data reduction are given in \cite{10.1093/pasj/psx080}.
The calibration of the photometry and astrometry is done against the Pan-STARRS1 \citep{2012ApJ...756..158S, 2012ApJ...750...99T,2013ApJS..205...20M, 10.1093/pasj/psx080, 2020ApJS..251....3M, 2020ApJS..251....4W, 2020ApJS..251....5M, 2020ApJS..251....6M,2020ApJS..251....7F}.
Detection is first performed in each band, and then the footprints and peaks of sources detected in different bands are merged to ensure consistency across different bands.
These peaks are separated and independent measurements of the positions and source parameters in each band which are obtained by the source measurement algorithm.
The reference band of detection is then defined for each object based on both the signal-to-noise ratio and the objective of maximising the number of objects in the reference band.
Finally, the source measurement is performed again with the position and shape parameters fixed to the values of the measurement in the reference band \citep{2018PASJ...70S...5B}.

To remove objects affected by artefacts around bright stars, such as halo, ghost, and blooming, we employ the ``bright'' star masks \citep[see Section~4 of ][]{10.1093/pasj/psab122} 
that are created using the $Gaia$ second data release \citep{Li_2021} \citep[also see][for discussion on the bright star masks]{2022PASJ...74..421L}. 
\subsection{The open-use HSC data of Sextans: HSC-PI}
The second HSC data we use is based on the PI-based HSC data of the Sextans field, taken through the Subaru open-use program and 
the Subaru/Keck Time-Exchange program, which is available from the Subaru archive system \citep[SMOKA][]{2002ASPC..281..298B}.
The data were taken on 2014-11-19, 2014-11-20, and 2016-04-03, with $g, i, i2$-band filters.
The typical exposure time for each field is about 210, 190, and 190~sec for $g, i$ and $i2$ filters, 
respectively, which is shallower than in the HSC-SSP data.
The typical seeing of HSC-PI data is about 
$0.8$--$1.2^{\prime\prime}$ 
for $g$, with some variations in the observational conditions, 
and 
$\sim 0.7^{\prime\prime}$ for the $i$ and $i2$ data, respectively.  
We re-analyzed the HSC-PI data using the same version (v8.5.3) of {\tt hscPipe} pipeline \citep{2010SPIE.7740E..15A} as that used in the HSC-SSP data.
We obtain about 2.74 million objects with $i_{\rm PSF}<24.5$ from the HSC-PI data.
In this paper, we use the primary photometric sources with $i_\mathrm{PSF} < 24.5$.
We employed the fainter magnitude cut than that for the SSP data ($i_{\rm PSF}<24$)
because the HSC-PI data is shallower than the HSC-SSP data 
so that we do not miss any object in the HSC-PI data, which exists in the HSC-SSP data, due to the larger photometric errors.
For the positions and the observed time epoch of objects, we use the centroid positions and the observed MJD in $i,i2$-band, respectively. If an object is observed in both $i,i2$-band, data of $i2$-band is given priority.
\section{Data Processing}
\label{sec:dataprocess}
\subsection{Star/galaxy separation}

\begin{figure}
    \centering
    \includegraphics[width=0.95\hsize]{img/galaxy_star_separation.pdf}
    \caption{
    The number distribution of all objects as a space of 
    the $i$-band PSF magnitudes ($x$-axis) and 
    the HSC $i$-band cModel minus PSF magnitudes ($y$-axis). 
    Colour denotes the number of objects in each grid of the two-dimensional space. 
    We select stars with the cut of $(i_\mathrm{cModel} - i_\mathrm{PSF})> -0.015$ as denoted 
    by the blue dashed line. 
    To select relatively compact galaxies that are well distinguishable from point sources, including stars, 
    we adopt a cut of $ -0.5 < (i_\mathrm{cModel} -  i_\mathrm{PSF})< -0.15$ as denoted by the green shaded region.
    We use galaxies for the recalibration of the astrometric solutions of stars.
    }
    \label{fig:galstarsep}
\end{figure}
To measure the proper motions of Sextans accurately, we first need to construct a secure sample of stars.
For this purpose, we use the high-resolution $i$-band data of HSC-SSP and identify point sources
on the basis of the ratio between PSF and cModel fluxes \citep[e.g.][]{2004AJ....128..502A}.
Fig.~\ref{fig:galstarsep} shows the distribution of objects in the master catalogue in the plane, 
$i_{\rm PSF}$ vs. $i_{\rm CModel}-i_{\rm PSF}$, where $i_{\rm PSF}$ and $i_{\rm CModel}$ are the PSF and cModel magnitudes, respectively.
In this work, we adopt the criterion $i_\mathrm{cModel} - i_\mathrm{PSF} > -0.015$ 
following the threshold in the HSC pipeline
to define point sources, i.e. candidates of stars. 

We also use a catalogue of galaxies to correct for systematic errors 
in astrometric solutions of stars between the HSC-SSP and HSC-PI catalogues. 
We define objects satisfying the condition of 
$ -0.5 < i_\mathrm{cModel} -  i_\mathrm{PSF} < -0.15$ as galaxies.
The upper limit $-0.15$ is sufficiently far from the threshold ($-0.015)$ of star separation, 
and therefore we believe we can safely select galaxies even in the presence of photometric errors 
thanks to the deep HSC photometry (a scatter of stars into this galaxy threshold due to photometric errors would be negligible). 
We also adopt the lower limit $-0.5$ because too extended galaxies are not helpful for the astrometry calibration 
because the accuracy of the centroid of such an extended galaxy is not so good.
\subsection{Matching objects in the HSC-SSP and HSC-PI catalogues}
Following the method in \citet{2021MNRAS.501.5149Q}, we measure the proper motion of each star from its angular offset between the HSC-SSP and 
HSC-PI catalogues over the time baseline relative to the reference frame defined by the galaxy positions.
To implement this method, we construct a matched catalogue of stars and galaxies between the HSC-SSP and HSC-PI catalogues, respectively.

For the matching, we start with each star or galaxy in the HSC-SSP catalogue, which 
has higher resolution and deeper depth than in the HSC-PI catalogue, 
and then search for a counterpart object in the HSC-PI within a radius of $1 \arcsec$.
Here the matching radius $1 \arcsec$ is large enough because most of the stars we are interested in (stars 
in the Sextans are at 86~kpc in the distance) should have much smaller angular offsets over our 
time baseline up to $\sim 4\,$yr\footnote{For comparison, 
if a star at 1~kpc distance has a velocity greater than 2300~km~s$^{-1}$ with respect to us, the star
has an angular offset larger than $1 \arcsec$ over 2~years. Most of the stars we are interested in are 
in greater distances than 1~kpc.}. 
To perform this matching, we used the software {\tt TOPCAT}\footnote{\url{http://www.star.bris.ac.uk/~mbt/topcat/}}.

\begin{figure*}
    \centering
    \includegraphics[width=0.95\hsize]{img/magnitude_hist.pdf}
    \caption{The magnitude histograms of star/galaxy catalogues. 
    The left panel is about stars and the right panel is about galaxies.
    The orange and green histograms denote objects of the HSC-PI catalogue and the HSC-SSP catalogue respectively. The blue histograms show the distribution of matched catalogues. The distribution of
    objects in HSC-SSP that matched multiple objects in HSC-PI is shown in the red histogram. 
    }
    \label{fig:galstarmaghist}
\end{figure*}
Fig.~\ref{fig:galstarmaghist} shows the number of star and galaxy objects in each bin of the $i$-band PSF magnitude 
for the HSC-SSP and HSC-PI catalogues, respectively.
To increase the purity of the matched objects, we discard objects that have a significant difference 
in their colours in between the HSC-SSP and HSC-PI catalogues.
To be more precise, we divide star or galaxy objects in each bin of the $i$-band PSF magnitudes, 
make a Gaussian fit to the distribution of colour for objects in the bin, and perform a $3\sigma$ clipping.
Here we use the magnitude difference between $g$-band and $i$-band as their colour, 
and use the difference between them, $(g-i)_{\rm SSP} - (g-i)_{\rm PI}$ for clipping.
We discard $9,083$ stars ($7\%$ of all stars) and $9,135$ galaxies ($4\%$ of all galaxies).
We found that the star objects for the HSC-PI catalogue are contaminated by galaxies at fainter magnitudes, 
which are properly removed by the better-seeing HSC-SSP catalogue. 
There is some contamination of the objects that have multiple matching, 
although the number of these objects is very small compared to that of the single-matched objects. 
We will not use the multiple-matched objects in the following analysis. 
In total, we have 121,015 matched stars and 215,325 matched galaxies, respectively. 
\subsection{Recalibration}
\label{sec:recalibration}

Quasars and galaxies are distant objects that should not have any proper motions. 
Hence, any apparent angular offsets of galaxy/quasar positions between the HSC-SSP and HSC-PI catalogues 
reflect residual systematic errors in the astrometric solutions of the two catalogues. 
Conversely, we use the observed angular offsets of galaxies to recalibrate the systematic errors \citep{2013ApJ...766...79K,2021MNRAS.501.5149Q}.
Following the method in \citet{2021MNRAS.501.5149Q}, we build a ``recalibration'' map as follows. 
The HSC CCD chip size is about $0.196\times 0.096$~deg$^2$ ($4,176\times 2,048$~pixels with mean pixel size of $0.168\arcsec$) \citep{2018PASJ...70S...1M}. 
Hence we split the sky footprint into small grids with sizes of $0.05\times0.05$~deg$^2$ each, 
which is chosen to be smaller than the CCD chip size and thus is able to capture spatial variations in the systematic errors within the CCD chip.
Also, It is notable that each grid contains a sufficient number of galaxies for the recalibration (each grid typically contains $\sim 25$ galaxies). 
For each grid, we compute the mean angular offset of matched galaxies between the HSC-SSP and HSC-PI 
after 3$\sigma$ clipping based on their angular offset distribution in the grid. 
Then we divide the angular offsets by the mean time separation, estimated based on the method in Section~\ref{sec:data}, 
to convert the angular offset into units of proper motions in R.A. and Dec. directions, i.e. mas\,yr$^{-1}$, in each grid. 
The errors of the angular offsets are about $1.7 \mathrm{mas~yr}^{-1}$, which are calculated from the 25, 50 and 75 percentiles of the angular offset distribution in the grid, 
using Eq.~(1) in \citet{2021MNRAS.501.5149Q}. 

\subsection{Proper motion measurement of stars}
\label{sec:pmmeasure}
\begin{figure*}
    \centering
    \includegraphics[width=0.95\hsize]{img/propermotionmaps.pdf}
    \caption{
    {\it Left panel}: The average ``net'' proper motion 
    of stars in each grid, which is estimated by averaging 
    the angular offsets of star positions in each grid between the HSC-SSP and HSC-PI. We convert it into the unit of proper motion, [mas~yr$^{-1}$], by dividing the angular offsets by the time baseline estimated for each object in Fig.~\ref{fig:visitdelta}.  
    These maps are constructed by 8720 nonzero grids in total, each of which has a size of $0.05 \times 0.05 ~\mathrm{deg}^2$. 
White-colour circles denote bright star masks which are not used for the following analyses. 
    {\it Middle panel}: Recalibration map of astrometric 
    solutions between the HSC-SSP and HSC-PI data, which 
    is estimated from the angular offsets of galaxies in each grid in units of [mas~yr$^{-1}$]. 
    Each grid contains about 25 galaxies on average. 
    {\it Right panel}: The proper motions of stars in each grid after subtracting the recalibration of astrometric solutions from the net motions (subtracting the middle from the left panels). 
    }
    \label{fig:pmmap}
\end{figure*}
As an initial estimate of proper motion for each of matched stars, we first estimate the angular offset of each 
matched star between the HSC-SSP and HSC-PI catalogues, 
and then divide it by the time difference between the two datasets. 
We adopt the separation time between the mean HSC-SSP and HSC-PI MJD dates for each matched object as the time baseline used for the proper motion measurements.
The time baseline ranges from 1.66 to 4.12 years for matched objects (stars and galaxies).
The average time baseline is about 2.66 years.

The left panel of Fig.~\ref{fig:pmmap} shows the map for the initial estimate of proper motions for matched stars. 
The map clearly shows spatial patterns reflecting the field-of-view  (FoV) of HSC.
This is because neither the HSC-SSP nor HSC-PI data corrected for the proper motions of stars before co-addition, 
and hence there are likely residual systematic errors in astrometric solutions in both catalogues.

In order to correct for the systematic errors of astrometric solutions, 
we use the catalogues of matched galaxies by making use of the fact that galaxies should not have any proper motions 
and can serve to define the reference coordinates in the two datasets: that is, the angular offset map introduced in Section~\ref{sec:recalibration}.
The middle panel of Fig.~\ref{fig:pmmap} shows the map of mean angular offsets for matched galaxies in each grid.
The spatial structure appears to display similar characteristics, 
such as the FoV of HSC pointings and depth variations to what is seen in the left panel.

To make a recalibration of the proper motion measurement for each of matched stars, we first find the closest grid of the galaxy offset map 
to each matched star, and then subtract the mean angular offset per year of galaxies from the measured proper motion of the star. 
This is our fiducial dataset of the proper motion measurements. 
The right panel of Fig.~\ref{fig:pmmap} shows the map for the proper motions of matched stars after the recalibration. 
It is now clear that the spatial structure of the map is gone. 
The proper motions display, except the Sextans region,
a spatially-homogeneous distribution, which reflects the proper motions of foreground halo stars in the MW  that should include both 
the intrinsic proper motion and the reflex motion due to our own motion of the Solar system.

\section{Results}
\label{sec:results}
In this section, we first estimate the structural parameters of Sextans using
three different density profile models.
We then estimate the systemic proper motion of Sextans. 

\subsection{Selection of Sextans member stars based on 
colour-magnitude diagram}
\label{sec:CMD}
\begin{figure}
    \centering
    \includegraphics[width=0.95\hsize]{img/CMD_membercut.png}
    %
    \caption{
    The colour-magnitude diagram (CMD)
    for stars in all data region (left panel)
    and the core region (right), where 
    the core region is defined by the region within 
    0.5~degree radius from the Sextan centre.
    The core region contains 12,377 stars in total, while our selection of Sextans member candidates, 
    inside the two solid lines, contains 7,951 stars.
    The points with error bars shown represent the width of star distribution in each magnitude.
    }
    \label{fig:member}
\end{figure}
To achieve a secure proper motion measurement of the Sextans galaxy, we make a separation of member star candidates of Sextans from the MW foreground stars
based on their colour-magnitude diagram (CMD) 
\citep[also see][for similar method]{10.1093/mnras/stw949}.
To define member star candidates in CMD, we look at the ``core'' region within $0.5^{\circ}$ radius from the centre of Sextans
that roughly corresponds to the half-light radius estimated in \cite{1995MNRAS.277.1354I}, i.e.  
about $27.6\arcmin$.
Through this section, we use the magnitude of each star where the Galactic dust extinction is corrected for.
\begin{figure}
    \centering
    \includegraphics[width=0.98\hsize]{img/spatialmap.png}
    %
    \caption{The spatial distribution of stars selected from the colour-magnitude cut given in the right panel 
    of Fig.~\ref{fig:member}. The figure clearly shows a concentration of stars due to member stars in Sextans.
    The blue point is the centre of Sextans. The red solid  circle shows 0.5~deg. radius from the centre, inside which stars are used to define the colour-magnitude cut in 
    Fig.~\ref{fig:member}. The red dashed outer circle denotes the radius inside which stars are used to estimate the structural parameters of Sextans. The orange solid elliptical line denotes twice the half-light radius $(R_h)$, where $R_h$ is estimated from the fitting to the Plummer profile (see text for details). We will use stars inside $2R_h$ 
     to measure the proper motion of Sextans. The orange dashed elliptical line denotes 
    the estimated tidal radius from the fitting to the King profile.
    }
    \label{fig:numspa}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=0.98\hsize]{img/maghist_corebg.pdf}
    %
    \caption{Blue histogram shows magnitude distribution of 7,747 stars in the inner region of Sextans, at radii $R_e<2R_h$ from the Sextans centre (see Fig.~\ref{fig:numspa}). This can be compared to the magnitude distribution of 3,581 stars outside Sextans, at separations $R_e>1.5$~deg. 
    \label{fig:maghist}
    }
\end{figure}
Fig.~\ref{fig:member} represents the $(g-i)$ colour and $i$-band magnitude diagram (CMD) of all stars in the core region of Sextans.
For comparison, we also show the CMD for stars in all the fields of our dataset.
The CMD reaches the main-sequence (MS) turn-off, and displays
a steep red giant branch (RGB) extending downwards from $i \sim 18$ to $i \sim 22.5$,
a blue and a red horizontal branch (BHB and RHB, respectively) at $i = 20 \sim 21$, 
a sub-giant branch (SGB), and blue stragglers (BSs) tracing down to $i\sim 23$.
As can be found from the left panel, there are many stars belonging to the colour-magnitude regions of
the Sextans RGB, SGB and MS stars, especially 
stars with colour $g-i\sim 0.3$-$0.5$ that are likely faint MS stars in the MW foregrounds \citep[also see][for the similar discussion]{10.1093/mnras/stx086}.
In contrast, stars belonging to RHB, BHB and BSs are relatively fewer.
Hence, we employ boundaries to preferentially select RGB, SGB and MS star candidates of the Sextans, as follows. 
We first use exclude horizontal branch stars, $i \sim 20$, and colour boundaries, $0 < \left(g-i\right) < 1$, to select stars inside the boundaries for the core region of Sextans. 
Secondly, we calculate the median of colour $(g-i)$ in each magnitude bin, stepped by $0.05$~mag, and 
fit the colour distribution for stars in each magnitude bin by a Gaussian function centred at the median of the colour.
Thirdly, we choose stars lying within $3\sigma$ of the Gaussian profile in each bin as member star candidates of the Sextans. 
More exactly, we perform a cubic interpolation of the $\pm 3\sigma$ width of the Gaussian profile as a function of the magnitude
and then choose stars inside the two bold solid lines in the right panel of Fig.~\ref{fig:member}.
This is our selection method of Sextans members.

Applying the selection of Sextans member candidates described above to the entire Sextans field including fields outside the core region, 
we have 13,792 stars.
The spatial distribution of the selected stars 
is shown in Fig.~\ref{fig:numspa}.
The figure clearly shows a significant concentration of stars around Sextans, meaning that our selection can properly pick up member stars. It is also clear that 
there are much fewer stars in the region outside Sextans. 
These remaining stars in the outside region are likely to be foreground halo stars in the MW, 
and we will use these stars to estimate contamination due to the foreground MW halo stars in our measurements, as 
will be described in detail. 

Fig.~\ref{fig:maghist} shows the magnitude distribution of stars in the core region of Sextans and outside Sextans, respectively, after the CMD selection. The figure clearly shows that the core region is dominated by faint stars 
at $i_{\rm PSF}>22.5$, reflecting member candidates of Sextans because Sextans is at a large heliocentric distance of about 86~kpc (its distance modulus 19.7) and therefore members stars tend to appear fainter than the foreground MW halo stars. 
On the other hand, stars outside the Sextans region (at separation greater than 1.5~deg.) are dominated by $i\sim 22$ stars, meaning that the stars tend to be brighter due to closer distances.
\subsection{Structural properties}
\label{sec:structure1d}

Before going to the proper motion measurement, 
we first estimate the structural parameters of Sextans following the maximum likelihood-based method in \citet{2008ApJ...684.1075M}
and \citet{2018ApJ...860...66M} \citep[also see][for the original work]{1998AJ....115.2359K,2006AJ....131..375W}.
The likelihood of estimating structural parameters, denoted as ${\bf p}$, from the observed angular positions of member stars, denoted
as $\{\boldsymbol{\theta}_i\}$,  
in the footprint is given as
\begin{align}
{\cal L}(\{\boldsymbol{\theta}_i\}|{\bf p})=\prod_{i}\ell_i(r_i |{\bf p})\, ,
\end{align}
where ${\bf p}$ is a set of parameters in the model profile of member stars (see below), $r_i$ is the angular separation of the $i$-th  star from the Sextans centre  (see below), $\ell_i$ is the likelihood for the $i$-th star,
and the product $\prod_i$ runs over all the stars used in the fitting. 
The likelihood $\ell_i$ for the $i$-th star is given as
\begin{align}
    \mathscr{l}_i(r_i|{\bf p})&=
    \frac{N_{\rm tot}-S_{\rm fit}\Sigma_b}{A}{\Sigma}_s(r_i,{\bf p})+\Sigma_b,
    \label{eq:likelihood_i}
\end{align}
where $N_{\rm tot}$ is the total number of stars used in the fitting, $S_{\rm fit}$ is the effective area ([deg$^2$]) used in the fitting, 
${\Sigma}_s(r_i,{\bf p})$ is the model of the normalised radial profile at the position $r_i$ from the centre, 
$\Sigma_b$ is the mean number density of foreground MW stars 
([deg$^{-2}$]). 
Note that we use dimension-less units for ${\Sigma}_s(r_i,{\bf p})$.
For the fitting, we use stars in the region within 1.5~deg. radius from the Sextans centre (see Fig.~\ref{fig:numspa}), where we have 10,591 stars in total after the colour-magnitude cut in Fig.~\ref{fig:member}.
The likelihood $\mathscr{l}_i$ for the $i$-th star takes into account both probabilities of the star being 
a member star of Sextans and a foreground MW star. 
For the model profile, we employ the ellipsoid profile throughout this paper, 
and in this case, the normalisation constant $A$ is defined by
\begin{align}
A\equiv (1-\epsilon) \int_0^{\infty}2\pi R\mathrm{d}R~ {\Sigma}_s(R,{\bf p}),
\end{align}
where $\epsilon$ is the ellipticity of the model profile. $A$ has a dimension of [deg$^2$]
and $R$ denote the radius from the centre in the elliptical coordinates that march the major and minor 
axes.
The elliptical radius $r_i$ for the $i$-th star is computed from its angular position as
\begin{eqnarray}
r_i^2  &=&  \frac{\left(X_i\cos\theta + Y_i\sin\theta\right)^2}{\left(1- \epsilon\right)^2} + \left(-X_i\sin\theta + Y_i\cos\theta\right)^2 \, 
\end{eqnarray}
where $X_i= (\alpha_i -\alpha_0)\cos\delta_0$ and $Y_i=\delta_i-\delta_0$, $\alpha_0$ and $\delta_0$ are the R.A. and Dec. of the assumed centre of Sextans, 
and $\theta$ is the position angle of the major axis of the ellipsoid profile, measured from the East-West direction in the sky coordinate. 

In this paper we consider the following three profiles:
a Plummer profile \citep{1911MNRAS..71..460P}, 
an exponential profile, and
a King profile \citep{1962AJ.....67..471K}, for $\Sigma_s$ in Eq.~(\ref{eq:likelihood_i}):
\begin{align}
&\Sigma_{\mathrm{exp}}(R;R_E) = \exp{\left(-\frac{R}{R_E}\right)},\nonumber\\
&\Sigma_{\mathrm{Plummer}}(R;R_p) = \left(1+\frac{R^2}{R_p^2}\right)^{-2}, \nonumber\\
&\Sigma_{\mathrm{King}}(R;R_c,R_t) = \left[
\left(1+\frac{R^2}{R_c^2}\right)^{-1/2}
-\left(1+\frac{R_t^2}{R_c^2}\right)^{-1/2}
\right]^2,
\end{align}
where $R_E$ and $R_p$ are the scale radius of the exponential and Plummer profiles, respectively, and $R_c$ and $R_t$ are the core and tidal radii for the King profile. 
Note that the King profile is defined only at $R<R_t$. 
The exponential and Plummer profiles are specified 
by 5 model parameters in total; 2 for the centre (R.A. and Dec.), the ellipticity ($\epsilon$), 
the position angle of the major axis ($\theta$) and the scale radius ($R_E$ or $R_p$). The King model 
is specified by 6 model parameters; in addition to the same first 4 parameters, $R_c$ and $R_t$. 
Note that $R_p$ is the same as the half-light radius $R_h$ (if the star distribution follows a Plummer profile), and we will often refer to $R_h$ throughout this paper. 

In this paper, we use stars in the region inside 
$1.5$ degrees in separation from the reference centre, (R.A., Dec.)$=(153.26,-1.614)$, to estimate the structural 
parameters from the model fitting, following \cite{2016MNRAS.460.4492D}.
To estimate
the effective area $S_{\rm fit}$,
we first divide the fitting region into grids, which each have 
$0.025\times 0.025~$deg$^2$ area, 
and then sum up all the grids containing galaxies inside to estimate $S_{\rm fit}$, based on the idea 
that we can efficiently remove the grids contained in bright star masks.
The reason to use galaxies to select the surviving grids is that there are much more galaxies than stars after 
the magnitude and colour cut selection of Sextans members (Section~\ref{sec:CMD} and the right panel of Fig.~\ref{fig:member}). Note that we used the smaller-size grid 
of $0.025\times 0.025$~deg$^2$ than the $0.05\times 0.05$~deg$^2$ grids used for the astrometry calibration (Section~\ref{sec:recalibration}), because 
the smaller-area grids allow us to more accurately estimate the effective area.
Similarly, we estimate
the foreground number density $\Sigma_b$ 
from the total number of stars outside the fitting region; 
$\Sigma_b\simeq 
218.41~\mathrm{deg^{-2}}$.

We perform 
the Markov chain Monte Carlo (MCMC) inference of the structural parameters using the publicly-available code
{\tt emcee} \citep{Foreman-Mackey_2013}.
The estimated parameters are summarised in Table~\ref{tab:structural1d}. 
We also show the best-fit Plummer and King models in the upper panel of Fig.~\ref{fig:radpro}, 
compared with the 1D radial profile of the projected star number density as a function of the radii from the centre. 
The figure shows that the best-fit model fairly well reproduces the measured radial profile, although we found the best-fit model from 
a 2D likelihood analysis and did not use the 1D radial profile. 
Our result can be compared to the previous works: 
$R_p=16.5\arcmin\pm0.10\arcmin$, $(R_c,R_t)=(20.1\arcmin\pm 0.5\arcmin,60.5\arcmin\pm 0.6\arcmin)$ \citep{2018ApJ...860...66M};
$R_p=26.60\arcmin\pm0.43$, $(R_c,R_t)=(17.91\arcmin\pm0.65\arcmin, 120.5\arcmin\pm7.7\arcmin)$ \citep{10.1093/mnras/stx086};
$R_p=23.0\arcmin\pm0.4\arcmin$, 
$(R_c,R_t)=(26.8\arcmin\pm1.2\arcmin, 83.2\arcmin\pm7.1\arcmin)$ \citep{10.1093/mnras/stw949}. 
Our result is consistent with \citet{10.1093/mnras/stw949}. 
Our analysis is advantageous in the sense that we use fainter stars down to $i_{\rm PSF}\simeq 24$ and 
use stars over a sufficiently wider area around Sextans, which allows us to estimate the contamination of foreground MW stars.  
The previous works are based on brighter stars or data limited to a narrower area around Sextans. 
We will use the estimated spatial distribution of Sextans members when estimating the systemic proper motion of Sextans.

\begin{table*}
    \centering
    \begin{tabular}{llllll} \hline\hline
    Model&RA$_c$ &DEC$_c$ & $\theta$[deg] & ellipticity: $\epsilon$ & $R_{\mathrm{scale}}$ \\ \hline
    Plummer&$153.280\pm0.003$&$-1.604\pm0.003$ & $60.566\pm1.296$ & $0.256\pm0.011$ & $R_p: 20\farcm142\pm0\farcm233$ \\
    Exponential&$153.280\pm0.003$&$-1.602\pm0.003$ & $61.548\pm1.336$& $0.251\pm 0.011$ & $R_E: 12\farcm352\pm0\farcm134$\\
    King&$153.281\pm0.003$&$-1.602\pm0.003$
    & $61.785\pm1.338$& $0.253\pm0.011$ & $R_c: 14\farcm718\pm0\farcm347$, 
     $R_t: 101\farcm600\pm1\farcm804$\\ \hline\hline
    \end{tabular}
    \caption{
    The top rows are the structural parameters of Sextans, obtained from the likelihood inference of the model parameters 
    from the observed spatial  distribution of member stars  (see Section~\ref{sec:structure1d}) for details.
    For the model profile, we consider the Plummer, exponential, and King models including the ellipticity. 
    The model parameters are the coordinates of the centre, RA$_c$ and DEC$_c$, the ellipticity ($\epsilon$), and the 
    scale radius that characterises each profile; $R_p$ for the Plummer profile, $R_E$ for the exponential profile, and 
    $R_c$ and $R_t$ for the core and tidal radii of the King profile, respectively. Here we give 
    the central value and the $\pm 1\sigma$ error of each quantity. 
    }
    \label{tab:structural1d}
\end{table*}
\subsection{An estimation of statistical errors in 
Sextans proper motion measurements}
\label{sec:error}

\begin{figure}
    \centering
    \includegraphics[width=0.95\hsize]{img/errorest.pdf}
    %
    \caption{Estimation of statistical errors in the proper motion measurements. 
    Circle symbols with error bars represent the dispersion proper motion in R.A. (upper panel) and Dec. (lower) directions, respectively, for stars 
    in each magnitude bin for the field {\it outside} Sextans by separation greater than 1.5~deg. (the region outside 
    the red dashed circle in Fig.~\ref{fig:numspa}). The error bars are estimated from the standard deviation of 100 bootstraps re-sampling in each bin. The Red dashed line denotes a best-fit model assuming the dispersion arises 
    from statistical errors of astrometry estimation (the centroid determination) of each star due to the photometric errors (see text for details). The Green dot-dashed line shows a best-fit model if the dispersion arises from the intrinsic, random velocity dispersion of MW halo stars (see text for details). The solid orange line is the sum of the two contributions.
    For comparison, the triangle symbols show the dispersion of apparent proper motions for the SDSS quasars in the Sextans field
    that are measured in the same way as that for stars. 
    }
    \label{fig:errorest}
\end{figure}
Before going to the measurement result of Sextans proper motion, we quantify statistical errors in the proper motion measurements. There are two contributions to the statistical errors.
One is a statistical error arising from measurement errors due to imperfect astrometry of individual stars.
The other is an error arising from the intrinsic random motions of stars that are especially important for foreground MW
halo stars. In this section, we show a method to discriminate the two contributions using the method described 
in Appendix~C of \citet{2021MNRAS.501.5149Q}.

Quasars, which are
distant objects from us and are point sources, should not have any apparent proper motion if we have perfect astrometry. In other words, apparent proper motions of quasars give an estimate of statistical errors in our proper motion measurements.
For this purpose, we use the quasar catalogue from the
SDSS-DR16Q \citep{2020ApJS..249....3A}. After matching each of the SDSS quasars to point sources of the HSC-SSP 
catalogue, we performed the apparent proper motion measurements using the same method used for the measurements 
of star proper motions. Note that the astrometry recalibration is done using the catalogue of galaxies in our method, 
and the quasar catalogue is not used in our method and is an independent catalogue.

Fig.~\ref{fig:errorest} shows the dispersion of proper motions for stars in the field outside the Sextans centre by more than 1.5~degree in radius, as a function of star magnitudes, where we estimate the dispersion 
from all stars sitting in each magnitude bin. The dispersion can be compared with the dispersion of quasars measured 
in the same way.
The figure shows that the dispersion of quasars is smaller than that of stars in magnitude bins of $i_{\rm PSF}\lesssim 22$, meaning that the dispersion of stars has an additional contribution to the statistical scatters indicated by the dispersion of quasars. The red dashed curve shows the best-fit model for the dispersion at $i_{\rm PSF}>22$, which is 
the error expected if the astrometry precision scales with 
flux as $1/{\rm flux}\propto 10^{0.4i_{\rm PSF}}$\footnote{Please see Section~6.2.3 in the document 
by Robert~Lupton, which describes how the statistical error or precision of centroid of each star depends on its flux or 
more exactly a signal-to-noise of the flux. Since our proper motion measurement depends primarily on the centroid 
determination of objects, the statistical errors of proper motions depend on the precision of centroid 
determination.}, where the normalisation is a model parameter. The red-dashed curve nicely reproduces the measured dispersion at the faint magnitudes, which is consistent with the statistical measurement errors. The dispersion of quasars is qualitatively consistent with the red dashed curve in bright magnitude bins, $i_{\rm PSF}<22$. On the other hand, the stars at $i_{\rm PSF}<22$ displays larger dispersion than
that of quasars or the red-dashed curve.

Hence we conclude that the dispersion of stars at $i_{\rm PSF}<22$ has a contribution from the intrinsic 
random proper motions of MW halo stars. MW halo stars are expected to have a constant random velocity dispersion of 
$80$--100~km$~$s$^{-1}$ \citep{2010ApJ...716....1B}. If the MW halo stars selected by our colour-magnitude cut (the 
right panel of Fig.~\ref{fig:member}) contain stars of a similar absolute magnitude in each $i_{\rm PSF}$ bin, 
the dispersion of proper motion is expected to scale with apparent magnitude $i_{\rm PSF}$ as 
$\sigma_\mu\propto \sigma_{v}/d\propto 10^{-0.2i_{\rm PSF}}$. The green dot-dashed curve in Fig.~\ref{fig:errorest} shows 
the best-fit model, where the normalisation parameter of $\sigma_\mu\propto 10^{-0.2i_{\rm PSF}}$ is estimated in that
the sum of the red dashed and green dot-dashed curves reproduce the measured dispersion. It is found that 
the best-fit model, shown by 
the solid orange curve, fairly well reproduces the measured dispersion. 
Thus HSC data would allow us to explore the random velocity 
dispersion of MW halo stars and the detailed study will be presented elsewhere.

\subsection{Measurement results for systemic proper motion of Sextans}
\label{sec:pm_measurements}

Now we show the results for measurements of a systemic proper motion of members stars 
in Sextans. To do this, we use two different methods; a method assuming 
the 1D radial profile of the spatial distribution of 
member stars and a method using the likelihood analysis of individual stars.

For the 1D fitting, we first measure the ``average'' proper motion of stars residing in each 
elliptical annulus from the Sextans centre, where we use the elliptical coordinates obtained from the best-fit Plummer 
model in Table~\ref{tab:structural1d}. We then estimate 
the proper motion of Sextans by minimising the following chi-square, for each of the R.A. or Dec. components:
\begin{align}
\chi^2=\sum_{i}\frac{\left[\mu_{\rm meas}(R_i)-\mu_{\rm model}(R_i)\right]^2}{\sigma(R_i)^2},
\label{eq:1dfitting}
\end{align}
where $R_i$ is the average of radial separations of stars in $i$-th radial bin, the summation runs over the annulus bins, 
$\mu_{\rm meas}(R_i)$ is the measured average proper motion of stars in the $i$-th bin, and $\mu_{\rm model}(R_i)$ is the model proper motion, given by
\begin{align}
\mu_{\rm model}(R_i)\equiv 
f(R_i) \mu_{\rm sxt}+\left[1-f(R_i)\right]\mu_{\rm MW}, 
\label{eq:mu_1D_model}
 \end{align}
with 
\begin{align}
f(R_i)\equiv 1-\frac{\Sigma_{b}S(R_i)}{N_{\rm obs}(R_i)}.
\end{align}
Here $\mu_{\rm sxt}$ is the systemic proper motion of Sextans, which is a model parameter, 
$\mu_{\rm MW}$ is the systemic proper motion of foreground MW stars, 
$f(R_i)$ is the fraction 
of Sextans member stars among the observed stars in the annulus of $R_i$,
$\Sigma_b$ is the number density of foreground MW halo stars, the same as that used in Eq.~(\ref{eq:likelihood_i}),
and $S(R_i)$ is the area of the $i$-th annulus.  
As we will show below, we employ radial binning so that each bin contains about 300 stars.
Throughout this paper, we use 7,446 stars at the annulus ranging from $R_0$ (inner radius of the first bin) to $2R_h$ (see Fig.~\ref{fig:compare}) 
for the fitting as our default choice. The lower limit is to prevent systematic errors because the first annulus where the region is so small that the galaxy correction is done by almost one grid. The upper limit is to include the region member stars dominate.
If we use different fitting regions to include stars outside the default region, the statistical uncertainties of proper motion measurement are not largely
improved, although the central values vary within the statistical errors. 
In the above chi-square, we properly consider contamination of foreground MW halo stars.
For $\mu_{\rm MW}$, we estimate the average coherent proper motions of MW halo stars in the regions outside 
Sextans, in separation greater than 1.5~deg from the centre: we estimated $(\mu^{\alpha}_{\rm MW},\mu^{\delta}_{\rm MW})=\left( -2.023, -3.585\right)$~[mas~yr$^{-1}$]. We accurately measure the proper motions of MW halo stars, so
we employ the fixed values in the following fitting. We checked that, even if we include the effect of the measurement errors of $\mu_{\rm MW}$,  the following results remain almost unchanged. 
For the annulus area $S(R_i)$ we estimate the effective area by summing the grids, each of which has an area of 
$0.025\times 0.025$~deg$^2$, contained in the annulus, excluding the grids due to bright star masks.

For the statistical uncertainties $\sigma(R_i)$ in the denominator of Eq.~(\ref{eq:1dfitting}), we model it as
\begin{align}
\sigma(R_i)^2=\sigma_{\mu, {\rm meas}}(R_i)^2+\frac{\left[\Sigma_bS(R_i)\right]^2}{N_{\rm obs}^3}
\left(\mu_{\rm sxt}-\mu_{\rm MW}\right)^2,
\end{align}
where $\sigma_{\mu,{\rm meas}}$ is the measurement error in each bin that is estimated from the distribution of 
proper motions of stars in the $R_i$ bin following the method in \citet{2021MNRAS.501.5149Q} and the second term 
is from the Poisson error of the number counts of stars in the $R_i$ bin; $\sigma(\delta N_{\rm obs})^2=N_{\rm obs}$. 

Fig.~\ref{fig:radpro} shows the systemic proper motion measurement of Sextans obtained from the 1D profile fitting 
(Eq.~\ref{eq:1dfitting}). At large separations $R\gtrsim 2$~deg., the measurements display a constant (separation-independent) proper motion, which means that the proper motions are dominated 
by a coherent motion of foreground MW halo stars, mainly due to the apparent reflex motion of stars by 
the proper motion of the Solar system. At radii $R_\lesssim 2$~deg., the measurements display scale-dependent 
profiles of star distribution and proper motions. The red-colour, shaded lines denote the best-fit model prediction 
and the $1\sigma$ uncertainty, which fairly well reproduce the measured proper motions. The obtained proper motion 
of Sextans is given by 
$\left(\mu^{\alpha}_\mathrm{sxt}, \mu^{\delta}_\mathrm{sxt}\right) = \left(-0.508^{+0.071}_{-0.071}, 0.083^{+0.071}_{-0.072}\right)$~mas~yr$^{-1}$.

\begin{figure}
    \centering
    \includegraphics[width=0.95\hsize]{img/pm_radial_profile.pdf}
    %
    \caption{{\it Upper panel}: The radial profile of number density of stars in each elliptical annulus of radius $R$ in the $x$-axis from the Sextans centre. Here we employ the elliptical coordinates (the axis ratio 
    and position angle) obtained from the fitting to the Plummer model in Table~\ref{tab:structural1d}. 
    The solid orange and blue lines denote the best-fit Plummer and King models, respectively.
    The vertical dashed line denotes twice the half-light radius, $2R_h$, inside which stars are used for 
    the proper motion measurements. The horizontal dashed line at outer radii denotes the density of MW halo stars that is estimated from the field outside Sextans in separations greater than 1.5~deg. from the centre of Sextans.
    Although we show the Poisson errors in each bin, the error bars are not visible because of small errors.  
    {\it Middle and lower panels}: The data points show the estimated systemic proper motion of Sextans, which is obtained by averaging the proper motions of stars in each annulus 
    in the R.A. and Dec. directions. The errors denote the statistical uncertainties that are estimated from the distribution of proper motions of stars in each bin. The shaded red region in each panel denotes the best-fit 
    model for the bulk proper motion of Sextans, where the width of the region denotes the $\pm 1\sigma$ range. 
    The horizontal dashed lines at outer radii denote the coherent proper motion of MW halo stars which is
    estimated from the averaged proper motion of stars in the region outside Sextans.
    }
    \label{fig:radpro}
\end{figure}
As an alternative method, we also perform a maximum likelihood analysis of Sextans proper motion measurement, 
inspired by  \citet{Pace_2019}. We define the likelihood function for the proper motions of all individual stars as 
\begin{align}
{\cal L}^{\rm PM}({\bf d}|
\boldsymbol{\mu}_{\rm model})\equiv \prod_i {\ell}_i^{\rm PM}({\bf d}_i|\boldsymbol{\mu}_{{\rm model}}),
\label{eq:likelihood_2d}
\end{align}
where $\ell_i^{\rm PM}$ is the likelihood for proper motion of the $i$-th star, given by 
\begin{align}
{\ell }^{\rm PM}_i&\equiv f(R_i){\cal N}(\boldsymbol{\mu}_{{\rm meas},i}|\boldsymbol{\mu}_{\rm sxt},\sigma_{\rm sxt})
\nonumber\\
&\hspace{4em}+\left[1-f(R_i)\right]{\cal N}(\boldsymbol{\mu}_{{\rm meas},i}|\boldsymbol{\mu}_{\rm MW},\sigma_{\rm MW}).
\end{align}
We again use 7,446 stars at $R<2R_h$ as in the 1D fitting. 
Here the first term denotes the likelihood for the star to be a member star of Sextans, and the second term 
is the likelihood for the star to be a foreground MW star,  
where $\boldsymbol{\mu}_{{\rm meas},i}$ is the measured proper motion for the $i$-th star 
and $\boldsymbol{\mu}_{\rm sxt}$ and $\boldsymbol{\mu}_{\rm MW}$ are the proper motions of Sextans and
MW halo stars, respectively. As we did for the 1D profile method, we use the measured, average proper motions of 
stars in the region sufficiently outside the Sextans, and we treat 
the Sextans proper motions, $\boldsymbol{\mu}_{\rm sxt}$, as fitting model parameters. We assume that the $i$-th 
star proper motion follows a Gaussian distribution, denoted as ${\cal N}(\mu_{\rm meas}|\mu_{\rm model}, \sigma)$:
\begin{align}
{\cal N}(\mu_{\rm meas}|\mu, \sigma)\equiv \frac{1}{\sqrt{2\pi}\sigma}
\exp\left[-\frac{(\mu_{\rm meas}-\mu_{\rm model})^2}{2\sigma^2}\right]. 
\end{align}
For a member star of Sextans, we use the dispersion denoted by the red dashed curve in Fig.~\ref{fig:errorest} 
as a function of $i_{\rm PSF}$ to model the error $\sigma_{\rm sxt}$. 
For a MW halo star, we use the dispersion denoted by the orange curve in Fig.~\ref{fig:errorest}
as a function of $i_{\rm PSF}$ to model the error $\sigma_{\rm MW}$. 
To model $f(R_i)$, which gives a probability of the star to be a member star of Sextans, we use the measured 1D profile in the upper panel of Fig.~\ref{fig:radpro}: after using the cubic interpolation of the measured 1D profile, 
we estimate $f(R_i)$ at an arbitrary separation $R_i$. 
Hence the fitting parameters in Eq.~(\ref{eq:likelihood_2d}) are only $\mu_{\rm sxt}^{\alpha}$ and $\mu_{\rm sxt}^{\delta}$, i.e. two model parameters. Note that the likelihoods for each of R.A. and Dec. components of the Sextans 
proper motions are separable, so we can perform the maximum likelihood analysis for the two components, separately. 
We find the best-fit parameters to be $(\mu_{\rm sxt}^{\alpha},\mu_{\rm sxt}^{\delta})= \left(-0.521^{+0.074}_{-0.071}, 0.047^{+0.070}_{-0.071}\right)$~mas~yr$^{-1}$. This 2D result is consistent with the result for the 1D profile fitting, although we use the same stars. The 2D likelihood analysis does not give an improved statistical uncertainty 
in the proper motions over the 1D profile method, meaning that including the 2D distribution of proper motions is not significant for the current data and statistical uncertainties.

\subsection{Discussion}
\label{sec:discussion}
\begin{figure}
    \centering
    \includegraphics[width=0.95\hsize]{img/pm_compare.pdf}
    %
    \caption{Comparison of our proper motion measurements of Sextans (star symbols) with other works. 
    Our results are taken from the 1D profile fitting and the 2D likelihood analysis presented 
    in Sections~\ref{sec:pm_measurements}. Other works are based on the $Gaia$ data of fewer, brighter stars in the 
    Sextans region than the stars used in this paper. 
    }
    \label{fig:compare}
\end{figure}
In Fig.\ref{fig:compare}, we compare the measured proper motion in this work with the results from other works, which are based on the $Gaia$ DR2 or EDR3 proper motions of bright stars in Sextans
\citep{2018AA...619A.103F, 2018AA...616A..12G, McConnachie_2020, Li_2021,  10.1093/mnras/stab1568}. 
Our results generally agree with the previous results to within the errors, although 
our results for the $\alpha$-component proper motion show a difference from the most recent results.  
The sample stars of $Gaia$ are limited to relatively bright BHB and RGB stars down to $G \sim 18$
while ours reach faint main sequence stars down to $i \sim 24$.

One may think the errors in our results are much bigger than $Gaia$'s. 
It is because of the difference in error definition.
Our errors are derived empirically from the proper motion dispersion, 
which accounts for the intrinsic proper motion dispersion of both the MW foreground and Sextans.
On the other hand, the previous results only include measurement errors.
Also, the error in this study is 
inversely proportional to the baseline time scale (about 2~yrs).
In other words, if we can combine new HSC data of Sextans with a longer time baseline, e.g. 5 years, combined with 
the data of this work, 
we can significantly reduce the statistical errors.

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

In this paper, we have used the two datasets of Subaru HSC, HSC-SSP, and HSC-PI called in this paper,  
for the Sextans field, which was taken at different epochs (about 2.66~years separation), to measure the proper motion of Sextans. The wide field of view and superb image quality 
(about 0\farcs6 seeing) brought an advantage to our study. First, the HSC data enabled us 
to make a clean separation of stars from galaxies, and the deep depth data enabled us to use stars down to 
$i\sim 24$ for the measurements. Then we used the positions of the matched galaxies in the two different-epoch HSC datasets to recalibrate the astrometry or equivalently define the reference frame for the proper motion measurements.

We selected secure candidates of Sextans member stars based on the colour-magnitude diagram (Fig.~\ref{fig:member}), which allowed us to select the candidates including main-sequence stars thanks to the HSC depth (the Sextans distance is about 86~kpc, and the distance modulus is about 19.8). We used the measured, spatial distribution of the member candidates to estimate the structural parameters of Sextans by performing a fitting of the assumed model profile. We found that the measured distribution of the member stars is well-fitted by 
Plummer, exponential, and King profiles. The model parameters are determined at a high precision thanks to a large number of member candidates, but the measured values show a slight difference from the results for some of the other works. Since the structure parameters reflect the distribution of dark matter for Sextans, our results would be important when combined with measurements of the kinematical structure \citep{2016MNRAS.461.2914H,2020ApJ...904...45H}.

We then measured the proper motion of Sextans and found that it is estimated as 
$(\mu_\alpha,\mu_\delta)= \left(-0.508^{+0.071}_{-0.071}, 0.083^{+0.071}_{-0.072}\right)$~mas~yr$^{-1}$. This result is basically consistent with the previous works within the error bars, but the precision of our measurement is limited by the time baseline (2.66~years). More importantly, our study gives a demonstration  
that ground-based data such as the Subaru HSC can be used to measure the proper motion of Sextans at a large distance (about 86~kpc), more generally proper motion of dwarf galaxies at such a large distance, if a recalibration of the astrometric solution is properly done. Our measurement can be improved by adding new data of HSC as it gives a longer time baseline. For example, if we have a factor of 2 longer time baseline, we can have a factor of 2 improvements in the statistical precision of the proper motion measurement. This is really promising because the upcoming LSST survey will give us a 10-year-long observation of the wide-solid angle sky (about 23,000~sq. degrees), which includes many dwarf galaxies or other clusters of stars such as globular clusters and stellar streams.

If we correct for the effect of reflection motion by the Solar system on our proper motion measurement, 
our results imply that Sextans has $v_t\simeq 141$~km~s$^{-1}$ for the tangential velocity assuming that its distance is at 86~kpc. Combining this with 
$v_r\simeq 224$~km~s$^{-1}$ for the radial velocity of Sextans in \citet{McConnachie_2020}, 
we find that Sextans has $v\simeq 265$~km~s$^{-1}$ for the three-dimensional velocity, where the radial motion towards 
the Galactic centre is dominated. This velocity value is consistent with the prediction for the MW halo model where the MW has a virial mass of $M_{\rm vir}\sim 10^{12}$~$M_\odot$ \citep{2018AA...619A.103F}.
We also find that Sextans has a relatively large impact parameter to the Galactic centre as $r_{\rm min}\sim 60$~kpc, implying that Sextans is not affected by the tidal disruption and preserves its own dark matter structure. Thus 
a proper motion measurement of dwarf galaxies is one of the key quantities to test the assembly history of the MW halo and also the nature of dark matter structure in the dwarf galaxies themselves and the MW halo. In fact, Sextans 
is one of the main dwarf targets with the upcoming Subaru wide-field-of-view multi-object spectrograph of Subaru, Prime Focus Spectrograph (PFS), for this purpose \citep{2014PASJ...66R...1T,2020ApJ...904...45H}, and our study gives the useful prior information on the dark matter analysis. 
We hope that our study gives useful guidance to the use of ground-based data for proper motion measurements.

\section*{Acknowledgements}
This work was supported in part by the World Premier International Research centre Initiative (WPI Initiative), MEXT, Japan, JSPS KAKENHI Grant Numbers JP18H05437, JP19H00677, JP20H00172, JP20H00181, JP20H05850, JP20H05855,  JP20H05856, JP21H05448, JP22J11959 and Basic Research Grant (Super AI) of Institute for AI and Beyond of the University of Tokyo. 

This paper is based on data collected at Subaru Telescope and obtained from the SMOKA, which is operated by the Astronomy Data Center, National Astronomical Observatory of Japan.We are honoured and grateful for the opportunity of observing the Universe from Maunakea, which has the cultural, historical, and natural significance in Hawaii.

The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University.  The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University.  Funding was contributed by the FIRST program from the Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency  (JST), the Toray Science  Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University.

This paper is based [in part] on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center (ADC) at NAOJ. Data analysis was in part carried out with the cooperation of Center for Computational Astrophysics (CfCA) at NAOJ.  We are honored and grateful for the opportunity of observing the Universe from Maunakea, which has the cultural, historical and natural significance in Hawaii.

This paper makes use of software developed for Vera C. Rubin Observatory. We thank the Rubin Observatory for making their code available as free software at \url{http://pipelines.lsst.io/}. 

The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg, and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

This work is based on data collected at Subaru Telescope and obtained from the SMOKA, which is operated by the Astronomy Data Center, National Astronomical Observatory of Japan.

\section*{Data Availability}

The data underlying this article are available in HSC-SSP at \url{https://doi.org/10.1093/pasj/psab122},
SDSS DR16 Quasar at \url{https://doi.org/10.3847/1538-4365/aba623}.
The data were derived from sources in the public domain: 
HSC-SSP at \url{https://hsc-release.mtk.nao.ac.jp/doc/}, 
Subaru Mitaka Okayama Kiso Archive at \url{https://smoka.nao.ac.jp/index.ja.jsp},
SDSS DR16 Quasar at \url{https://www.sdss.org/dr16/algorithms/qso_catalog}.


\bibliographystyle{mnras}
