11institutetext: Faculty of Computing and Informatics
Multimedia University, Malaysia

# Density-based Denoising of Point Cloud

Faisal Zaman  Ya Ping Wong  Boon Yian Ng

###### Abstract

Point cloud source data for surface reconstruction is usually contaminated
with noise and outliers. To overcome this deficiency, a density-based point
cloud denoising method is presented to remove outliers and noisy points.
First, particle-swam optimization technique is employed for automatically
approximating optimal bandwidth of multivariate kernel density estimation to
ensure the robust performance of density estimation. Then, mean-shift based
clustering technique is used to remove outliers through a thresholding scheme.
After removing outliers from the point cloud, bilateral mesh filtering is
applied to smooth the remaining points. The experimental results show that
this approach, comparably, is robust and efficient.

###### Keywords:

Point cloud, denoising, optimal bandwidth, particle swarm optimization,
bilateral filter

##  1 Introduction

In the past few decades, surface reconstruction has gained significant
attention due to the availability of commodity structured light sensors such
as Microsoft Kinect. Point cloud data acquired from these devices is usually
contaminated with noise and outliers [2], both of which could be caused by the
lighting or reflective nature of the surface or artifact in the scene. In
recent years, a number of point cloud denoising methods have been proposed [2,
1]. Most of the methods work only on certain type of noise or outliers.
Denoising algorithms such as data clustering [14, 15] are robust to outliers
but require prior knowledge about the input objects. Majority voting method
[16] can detect several types of outliers cluster, but high computation time
makes it infeasible for large datasets. Smoothing algorithms such as moving
least square [7] and mean curvature flow [5] are also used to remove noise and
outliers. These methods treat outliers as point with large noise and project
the noisy points to an estimated surface. However, these methods are sensitive
against large number of outliers and oversmooth the data points.

In this paper, a point cloud denoising method to remove noise and outliers
from the data is proposed. This method consists of the following steps. First,
the density of the input data point using kernel density estimation (KDE)
technique is evaluated. The performance of KDE is highly influenced by the
choice of smoothing parameter, also known as bandwidth. Most of the bandwidth
selection methods such as [11, 9, 17, 6] involved brute-force or exhaustive
search strategies have expensive computation. A particle swarm optimization
(PSO) (proposed by Kennedy and Eberhart, 1995) aided bandwidth selection
criterion using leave-one-out cross-validation (LOOCV) for multivariate kernel
density estimation is proposed. Because of the simple implementation and quick
convergence rates of PSO, it has been adapted in many practical applications
[3]. The optimal bandwidth selection method ensures the robust performance of
density estimation. After estimating the density, the mean-shift clustering
algorithm is applied to detect the local maxima of the constructed density
estimation. One of the main advantages of using mean-shift algorithm is its
ability to localize cluster modes automatically without any prior knowledge of
the number of clusters. Once clustering is done, outliers can be removed by
the proposed thresholding scheme. After removing the outliers, bilateral mesh
filtering is applied to smooth out the rest of the points. Experiments on
various point cloud datasets show that the method used in this research gives
better result than that of other approaches.

##  2 Our Method

In this section, an overview of the optimal bandwidth selection framework is
provided, starting from a review of the classical bandwidth selection methods.

###  2.1 Bandwidth in Kernel Density Estimation

Kernel density estimation (KDE), also known as Parzen window estimation, is a
non-parametric way of estimating the underlying distribution or probability
density function for a data set. Given $n$ data points $X_{i}$, $i=1,...,n$,
the kernel density estimate obtained with multivariate kernel function
$K_{H}(x)$ and bandwidth or smoothing parameter $H$, which is a symmetric
positive definite matrix, computed at the point $x$ is defined as:

$\hat{f}_{H}(x)=\frac{1}{n}\sum_{i=1}^{n}{K_{H}(x-X_{i})}$ (1)

where $x=(x_{1},x_{2},x_{3})^{T}$ is the $3$-dimension vector,
$K_{H}(x)=|H|^{-1/2}K(H^{-1/2}x)$. The resulting density through KDE depends
more heavily on the choice of bandwidth than that of the kernel. For the same
data, different bandwidths can produce different results. Optimal bandwidth
matrix — which minimizes the error between the estimated density
$\hat{f}_{H}(x)$ and the true density $f(x)$ is estimated. The performance of
$\hat{f}_{H}(x)$ can be determined by the risk function or integrated mean
squared error (IMSE) such as,

$R(\hat{f}_{H}(x),f(x))=E(L(\hat{f}_{H}(x),f(x)))$ (2)

where, $L$ is the loss function and
$L(\hat{f}_{H}(x),f(x))=\int{[\hat{f}_{H}(x)-f(x)]^{2}dx}$. The optimal
bandwidth can be obtained by minimizing the risk function, but since the true
density $f(x)$ is unknown leave-one-out cross-validation is used to estimate
the risk function. Therefore, loss function can be expressed as a function of
$H$.

$L(H)=\frac{1}{n}\sum_{i=1}^{n}log\,\hat{f}_{H,i}(x_{i})$ (3)

where, $\hat{f}_{H,i}(x_{i})$ is the leave-one-out estimator,
$\hat{f}_{-i}(x_{i})=\frac{1}{(n-1)}\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{n}{K_{H}(x_{i}-x_{j})}$ denotes the density estimation of
$x_{i}$ by using the other $(n-1)$ observed data points. Then the optimal
bandwidth matrix $H^{*}$ can be estimated by,

$H^{*}=arg\,\underset{H}{max}\,\frac{1}{n}L(H)$ (4)

Several classes of parameterizations of the bandwidth matrix are available in
the literature. Some popular choices include full bandwidth matrix, a diagonal
matrix with positive elements, or using a single bandwidth. In the
application, the diagonal bandwidth matrix,
$H=diag(h_{1}^{2},h_{2}^{2},h_{3}^{2})$, is employed since it is
computationally more efficient than the full bandwidth matrix.

Figure 1: The pipeline of the optimal bandwidth selection method, using PSO.

###  2.2 PSO for Optimal Bandwidth Selection

The PSO technique is proposed to carry out the bandwidth selection for better
approximation of kernel density estimation. PSO is a nature-inspired meta-
heuristic non-linear stochastic optimization technique. PSO is popular for its
ability to optimize complex non-linear functions and for its simple
implementation. PSO consists of a swarm of particles that fly through
hyperspace of the of objective function’s landscape. The pipeline of the
optimal bandwidth selection method shown in Fig. 1. The optimization problem
is solved from Eq. 4 using PSO, where $H=[h_{1},h_{2},h_{3}]^{T}$ is a
$3$-dimensional vector to be optimized and $L(\bullet)$ is the cost function.
Here are the steps of PSO algorithm to find the optimal bandwidth; where
$H_{i}^{t}$ is the potential solution of the $i^{th}$ particle at iteration
$t$.

  1. 1.

Swarm Initialization: Set the iteration index $t=0$. $S$ is the total number
of particles, where $\\{H_{i}^{t}\\}_{i=1}^{S}$ randomly generated particles
in solution space.

  2. 2.

Swarm Evaluation: Each particle denoted as $p_{i}^{t}$ remembers its best
position visited so far, which provides the cognitive information. Every
particle denoted as $p_{g}^{t}$ also knows the best position visited so far
among the entire swarm, which provides the social information. The cognitive
information $p_{i}^{t}$ and the social information $p_{g}^{t}$ are updated at
each iteration as follows:

if $L(H_{i}^{t})<L(p_{i}^{t-1})$ then

$p_{i}^{t}\leftarrow H_{i}^{t}$

end if

if $L(p_{i}^{t})<L(p_{g}^{t-1})$ then

$p_{g}^{t}\leftarrow p_{i}^{t}$

end if

  3. 3.

Swarm Update: Each particle $H_{i}^{t}$ has a velocity, denoted as
$v_{i}^{t}$. The velocity and position of the $i^{th}$ particle are updated in
each iteration according to:

$v_{i}^{t}=\omega{v_{i}^{t-1}}+C_{c}U[0,1](H_{i}^{t-1}-p_{i}^{t-1})+C_{s}U[0,1](H_{i}^{t-1}-p_{g}^{t-1})$
(5) $H_{i}^{t}=H_{i}^{t-1}+v_{i}^{t}$ (6)

Here $U[0,1]$ uniform random number between $0$ and $1$ and inertia weight
(proposed by Shi and Eberhart [13]) sets to $\omega=U[0,1]$. $C_{c}$ and
$C_{s}$ are constants and called cognitive learning rate and social learning
rate, respectively. It controls a particle’s learning behavior. Large value of
$C_{c}$ causes particles scatter around and slow convergence. On the other
hand, large value of $C_{s}$ causes fast convergence, and sometimes it can
cause the swarm to converge to local optima. These parameters are chosen
carefully so that global optima can be attained without compromising
convergence speed. In this application, the parameters $C_{c}$ reduced from,
$2.5$ to $0.5$ and $C_{s}$ varied between $0.5$ and $2.5$ according to
$C_{c}=2.5-2t/T$ and $C_{s}=0.5+2t/T$ [12].

The global best position $p_{g}$ of all particles during the previous three
steps is defined as $p_{g}=arg\,\underset{p_{i}}{max}[L(H_{i})]$, where
$p_{i}=H_{i}$. Once the maximum number of iterations reached algorithm stopped
and return the solution $p_{g}$.

Searching the solution depends highly on the number of selected particles. The
number of particles increases the chances of finding global optimal solution,
but it also increases the computation expenses significantly. To identify the
number of particles suitable for solving the current problem, a method also
known as population manager [10] is adopted. The idea behind the population
manager is to adjust the total number of particles in PSO according to the
solution-search status. The population manager determines the number of
suitable particles in the following ways: Assuming $k$ is population-manager
activating threshold. For $k$ consecutive generations, if there is no update
in $p_{g}$ and current population size does not exceed the maximum population
size, then a new particle will be added to the swarm. This new particle will
be added by combining the information from randomly selected two particles in
previous best solutions. If population size is equal to the maximum population
and there is no update in $p_{g}$, then remove a particle with poor
performance to make space for new potential particle. On the other hand, in
$k$ consecutive generations, if $p_{g}$ update multiple times, then remove a
particle with poor performance. Fig. 2a. demonstrates how population manager
dynamically adjust the particle in order to facilitate the searching process.

Figure 2: (a) Change in the number of particles using Population manager. (b)
Estimation performance of our optimal bandwidth selection method.

The initial population size in this application is set as one, and the maximum
population size is set as 45. To evaluate the PSO based bandwidth selection
method, a method to different types of point cloud models is applied (Fig.
2b). As the number of iterations increases, the mean square error (MSE)
gradually converges to zero and becomes steady once the optimal bandwidth is
found.

##  3 Noise and Outliers Removal

###  3.1 Algorithm for Outlier Removal

Mean-shift-based clustering algorithm [4] is a non-parametric unsupervised
clustering technique that does not require any prior information about the
number of clusters. It is applied to the kernel density estimate to identify
each local maxima or mode which represents one cluster. The mode of the
density function are located at the zeros of the gradient function
$\bigtriangledown{f}(x)=0.$ For each data point, say $x$, the mean-shift
procedure generates a sequence of points $\\{y_{j}\\}_{j=1,2,...}$ with

$y_{j+1}=\frac{\sum_{i=1}^{n}{x_{i}g\Big{(}\Big{\|}\frac{y_{j}-x_{i}}{h_{i}}\Big{\|}^{2}\Big{)}}}{\sum_{i=1}^{n}{g\Big{(}\Big{\|}\frac{y_{j}-x_{i}}{h_{i}}\Big{\|}^{2}\Big{)}}}\qquad
j=1,2,....$ (7)

and $y_{1}=x$. This sequence converges to a mode of density. The point $x$
belongs to the cluster corresponding to this mode.

Figure 3: Results of outliers removal on horse model. (a) The original point
cloud with outliers (b) The outliers detected by our method demonstrated in
red color, (c) The resultant point cloud after removing outliers (see the back
side of the horse in zoom-views)

After the clustering is done, the outliers need to be determined. For each
point $p$ in a cluster, where $x_{p}$ is the original location of $p$, the
average distance ($\bar{x}$) of all its neighboring points using k-nearest
neighbor(kNN)is calculated. Then the point $p$ is shifted iteratively to a new
position $\bar{x}$. After that, the distance between two positions,
$||\bar{x}$ \- $x_{p}||$, is calculated. If the difference between $||\bar{x}$
\- $x_{p}||$ and the average distance between its shifted neighbors below a
certain threshold, then point $p$ is considered outlier (Fig. 3b). The
outliers are discarded, leaving behind a new set of points, which have better
representation of the true shape of the cluster. Fig. 4 demonstrates the
results of the outliers removal algorithm in different models.

Figure 4: Results of the outliers removing algorithm: (a) Left: Bunny point
cloud with uniformly distributed noise, Right: Corresponding point cloud after
using the method used in this research, (b) Point cloud model of a car with
outliers data, Right: Car model after removing outliers.

###  3.2 Experimental Results and Analysis

In order to validate the performance and effectiveness of the proposed method,
the estimated performance in both synthetic and real world dataset is tested.
This program is written in Python, and the data obtained from the research is
from the tests running on Intel Core-i7 CPU 860 at 2.80 GHz plus 4GB RAM PC.
Fig. 4 demonstrates some of the results of this method(standford bunny and car
models). To evaluate the applicability of this method, both raw and synthetic
noise was tested. Table 1 summarizes timings and the parameters used to
generate the results. $T_{h}$ and $T_{f}$ is the computation time for the
optimal bandwidth selection using PSO and outliers removal using thresholding
scheme respectively.

Table 1: Computation time of the optimal bandwidth selection and outliers removal method for various models. Data | Input points | Filtered points | $T_{h}$ | $T_{f}$
---|---|---|---|---
Bunny | 330K | 293K | 29s | 30s
Car | 720K | 600K | 1m 45s | 1m 25s
Horse | 364K | 214K | 44.43s | 18.45s
Gargoyle | 2.1M | 796K | 3m 2s | 4m 44s
Head | 1.9M | 1.2M | 2m 5s | 3m 25s

Results of the denoising approach on scanned data show that it has good
performance on different types of scanned data. Experiments illustrate the
capability of this method in removing outliers.

###  3.3 Algorithm for Noise Removal

After removing outliers, triangular mesh for the rest of the points is
calculated. After that, the noise from the mesh by applying bilateral mesh
filter [8] is removed. Bilateral mesh filter is popular for its ability to
remove noise and preserve features at the same time. For denoising a point
$p$, other points around its neighborhood $\\{q_{i}\\}$ is identified; where,
$||p-q_{i}||<\rho=|2\sigma_{c}|$, here $\sigma_{c}$ is the radius of the
neighborhood of point $p$. To determine a desire value for $\sigma_{c}$, the
average distance of all adjacent triangles in an input mesh is used. The
bilateral filters in this method also have two parameters, $\sigma_{c}$ and
$\sigma_{s}$. The $\sigma_{s}$ is the standard deviation of the offsets in the
selected neighborhood. The result of bilateral mesh filtering operation on the
input model is illustrated in Fig. 5.

Figure 5: Point cloud denosing results: The first column is the point cloud
models after removing outliers; second column represents the triangular mesh
representation of the point cloud and third column is the models after
applying bilateral mesh denoising.

##  4 Conclusion

In this paper, a density-based approach for denoising point cloud data is
proposed. In particular, particle swarm optimization technique is used for
estimating the optimal bandwidth matrix of kernel density estimation, and
mean-shift clustering technique to remove outliers. To remove noise from
remaining points, bilateral mesh denoising method is applied. Our method is
very robust with highly noisy dataset, as can be seen from the examples shown.
While this research illustrates only a few applications of the method within
the domain of point cloud data, it is believed that its simplicity and
effectiveness will lead to its application in other domain.

## References

  * [1] Berger, M., Levine, J.A., Nonato, L.G., Taubin, G., Silva, C.T.: A benchmark for surface reconstruction. ACM Transactions on Graphics (TOG) 32(2), 20 (2013)
  * [2] Berger, M., Tagliasacchi, A., Seversky, L., Alliez, P., Levine, J., Sharf, A., Silva, C.: State of the art in surface reconstruction from point clouds. In: EUROGRAPHICS star reports. vol. 1, pp. 161–185 (2014)
  * [3] Blum, C., Groß, R.: Swarm intelligence in optimization and robotics. In: Springer Handbook of Computational Intelligence, pp. 1291–1309. Springer (2015)
  * [4] Comaniciu, D., Meer, P.: Mean shift: A robust approach toward feature space analysis. Pattern Analysis and Machine Intelligence, IEEE Transactions on 24(5), 603–619 (2002)
  * [5] Desbrun, M., Meyer, M., Schröder, P., Barr, A.H.: Anisotropic feature-preserving denoising of height fields and bivariate data. In: Graphics interface. vol. 11, pp. 145–152. Citeseer (2000)
  * [6] Duong, T., Hazelton, M.L.: Convergence rates for unconstrained bandwidth matrix selectors in multivariate kernel density estimation. Journal of Multivariate Analysis 93(2), 417–433 (2005)
  * [7] Fleishman, S., Cohen-Or, D., Silva, C.T.: Robust moving least-squares fitting with sharp features. In: ACM Transactions on Graphics (TOG). vol. 24, pp. 544–552. ACM (2005)
  * [8] Fleishman, S., Drori, I., Cohen-Or, D.: Bilateral mesh denoising. In: ACM Transactions on Graphics (TOG). vol. 22, pp. 950–953. ACM (2003)
  * [9] Guidoum, A.C.: Kernel estimator and bandwidth selection for density and its derivatives (2013)
  * [10] Hsieh, S.T., Sun, T.Y., Liu, C.C., Tsai, S.J.: Efficient population utilization strategy for particle swarm optimizer. Systems, Man, and Cybernetics, Part B: Cybernetics, IEEE Transactions on 39(2), 444–456 (2009)
  * [11] Hyndman, R.L., Zhang, X., King, M.L.: Bandwidth selection for multivariate kernel density estimation using mcmc. In: Econometric Society 2004 Australasian Meetings. No. 120, Econometric Society (2004)
  * [12] Ratnaweera, A., Halgamuge, S.K., Watson, H.C.: Self-organizing hierarchical particle swarm optimizer with time-varying acceleration coefficients. Evolutionary Computation, IEEE Transactions on 8(3), 240–255 (2004)
  * [13] Shi, Y., Eberhart, R.: A modified particle swarm optimizer. In: Evolutionary Computation Proceedings, 1998. IEEE World Congress on Computational Intelligence., The 1998 IEEE International Conference on. pp. 69–73. IEEE (1998)
  * [14] Song, Y.: Boundary fitting for 2d curve reconstruction. The Visual Computer 26(3), 187–204 (2010)
  * [15] Sotoodeh, S.: Outlier detection in laser scanner point clouds. International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences 36(5), 297–302 (2006)
  * [16] Wang, Y., Feng, H.Y.: Outlier detection for scanned point clouds using majority voting. Computer-Aided Design 62, 31–43 (2015)
  * [17] Zhang, X., King, M.L., Hyndman, R.J.: A bayesian approach to bandwidth selection for multivariate kernel density estimation. Computational Statistics & Data Analysis 50(11), 3009–3031 (2006)


# Multi-colour optical light curves of the companion star to the millisecond
pulsar PSR J2051$-$0827

V. S. Dhillon,1,2 M. R. Kennedy,3,4 R. P. Breton,4 C. J. Clark,4,5,6 D. Mata
Sánchez,2,4,7 G. Voisin,4,8 E. Breedt,9 A. J. Brown,1 M. J. Dyer,1 M. J.
Green,10 P. Kerry,1 S. P. Littlefair,1 T. R. Marsh,11 S. G. Parsons,1 I.
Pelisoli,11 D. I. Sahman,1 J. F. Wild,1 M. H. van Kerkwijk,12 B. W. Stappers4
1Department of Physics and Astronomy, University of Sheffield, Sheffield S3
7RH, UK
2Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain
3Department of Physics, University College Cork, Cork, Ireland
4Jodrell Bank Centre for Astrophysics, University of Manchester, Manchester
M13 9PL, UK
5Max Planck Institute for Gravitational Physics (Albert Einstein Institute),
D-30167 Hannover, Germany
6Leibniz Universität Hannover, D-30167 Hannover, Germany
7Departamento de Astrofísica, Universidad de La Laguna, E-38206 La Laguna,
Tenerife, Spain
8LUTHL, Observatoire de Paris, PSL Research University, 92195, Meudon, France
9Institute of Astronomy, University of Cambridge, Cambridge CB3 0HA, UK
10Department of Astrophysics, Tel Aviv University, Tel Aviv 6997801, Israel
11Department of Physics, University of Warwick, Coventry CV4 7AL, UK
12Department of Astronomy and Astrophysics, University of Toronto, Toronto, ON
M5S 3H4, Canada
E-mail<EMAIL_ADDRESS>
###### Abstract

We present simultaneous, multi-colour optical light curves of the companion
star to the black-widow pulsar PSR J2051$-$0827, obtained approximately 10
years apart using ULTRACAM and HiPERCAM, respectively. The ULTRACAM light
curves confirm the previously reported asymmetry in which the leading
hemisphere of the companion star appears to be brighter than the trailing
hemisphere. The HiPERCAM light curves, however, do not show this asymmetry,
demonstrating that whatever mechanism is responsible for it varies on
timescales of a decade or less. We fit the symmetrical HiPERCAM light curves
with a direct-heating model to derive the system parameters, finding an
orbital inclination of $55.9^{+4.8}_{-4.1}$ degrees, in good agreement with
radio-eclipse constraints. We find that approximately half of the pulsar’s
spin-down energy is converted to optical luminosity, resulting in temperatures
ranging from approximately $5150^{+190}_{-190}$ K on the day side to
$2750^{+130}_{-150}$ K on the night side of the companion star. The companion
star is close to filling its Roche lobe ($f_{\rm RL}=0.88^{+0.02}_{-0.02}$)
and has a mass of $0.039^{+0.010}_{-0.011}$ M⊙, giving a mean density of
$20.24^{+0.59}_{-0.44}$ g cm-3 and an apsidal motion constant in the range
$0.0036<k_{2}<0.0047$. The companion mass and mean density values are
consistent with those of brown dwarfs, but the apsidal motion constant implies
a significantly more centrally-condensed internal structure than is typical
for such objects.

###### keywords:

stars: neutron – pulsars: individual: PSR J2051$-$0827\.

††pubyear: 2022††pagerange: Multi-colour optical light curves of the companion
star to the millisecond pulsar PSR J2051$-$0827–Multi-colour optical light
curves of the companion star to the millisecond pulsar PSR J2051$-$0827

##  1 Introduction

Table 1: Journal of ULTRACAM and HiPERCAM observations of PSR J2051$-$0827\.
$n_{\rm exp}$ is the number of exposures and $t_{\rm exp}$ the exposure time
of each frame in seconds.

Instrument+telescope | Date start | UTC start | UTC end | $n_{\rm exp}$ | $t_{\rm exp}$ | Filters | Moon/transparency/$i$-band seeing
---|---|---|---|---|---|---|---
ULTRACAM+WHT | 2011/08/26 | 21:01 | 03:03 | 1085 | 20.0 | $u^{\prime}g^{\prime}i^{\prime}$ | Dark/non-photometric/0.9′′
ULTRACAM+WHT | 2011/08/27 | 20:50 | 03:27 | 1187 | 20.0 | $u^{\prime}g^{\prime}i^{\prime}$ | Dark/photometric/1.4′′
HiPERCAM+GTC | 2021/08/06 | 22:52 | 01:26 | 300∗ | 30.8∗ | $u_{\rm s}g_{\rm s}r_{\rm s}i_{\rm s}z_{\rm s}$ | Dark/photometric/0.8′′

∗Half the number of frames, each of double the exposure time, were obtained in
$u_{\rm s}$.

Table 2: Magnitudes of PSR J2051$-$0827 at light-curve maximum and minimum, measured with HiPERCAM – see Section 3.1 for details. Filter | Magnitude at
---|---
| max | min
$u_{\rm s}$ | $25.0\pm 0.1$ | $>25.9$
$g_{\rm s}$ | $23.47\pm 0.02$ | $>26.8$
$r_{\rm s}$ | $22.67\pm 0.01$ | $27.9\pm 0.4$
$i_{\rm s}$ | $22.30\pm 0.01$ | $25.6\pm 0.1$
$z_{\rm s}$ | $22.12\pm 0.02$ | $24.8\pm 0.2$

Pulsars are highly-magnetised, rotating neutron stars. Over 3000 are
known111http://www.atnf.csiro.au/research/pulsar/psrcat, most of which have
spin periods in the 0.1-1 s range that gradually increase with time due
primarily to the emission of magnetic dipole radiation. After tens of millions
of years, the spin slows to such an extent that the mechanism powering the
radio emission turns off and the pulsar ‘dies’. There exist a sub-set of more
than 500 known
pulsars222http://astro.phys.wvu.edu/GalacticMSPs/GalacticMSPs.txt (Manchester
et al., 2005), however, that have the fastest spin periods, of order
milliseconds, and are believed to be much older ($\sim 10^{9}$ yr) than
ordinary pulsars. These so-called ‘millisecond pulsars’ (MSPs) are believed to
be dead pulsars that have been spun up (or ‘recycled’) by the accretion of
mass from a companion star via Roche-lobe overflow (see Tauris & van den
Heuvel (2006) and references therein). During this accretion phase, the object
appears as an X-ray binary. Once accretion has stopped, the pulsar begins
emitting in the radio again. In some of the closer binaries, particle and
$\gamma$ radiation from the pulsar is then believed to ablate the companion
star, possibly evaporating it entirely to leave an isolated millisecond pulsar
(see Polzin et al. (2020) and references therein). Hence such systems are
sometimes referred to as ‘black-widow pulsars’. They consist of a millisecond
pulsar in a tight orbit ($P\lesssim 24$ h) with a very low mass companion star
($M_{2}\lesssim 0.05M_{\odot}$) and usually exhibit radio eclipses each orbit
due to the obscuration of the pulsar by the ablated material. Objects in each
of the evolutionary phases described above have been discovered, which lends
support to this general picture, although many uncertainties remain. For a
review of MSPs, see Lorimer (2008).

PSR J2051$-$0827 was the second black-widow pulsar to be discovered in the
Galactic disk (Stappers et al., 1996a), after PSR J1959+2048 (also known as
PSR B1957+20; Fruchter et al. 1988), and is a 4.5 ms pulsar in a 2.4 hr period
orbit with a low mass companion. The companion star was detected in the
optical by Stappers et al. (1996b), and orbital light curves were subsequently
obtained by Stappers et al. (1999) and Stappers et al. (2001). The nature of
the companion star in the PSR J2051$-$0827 system remains uncertain, due to
the fact that asymmetries were observed in the light curves of Stappers et al.
(2001) which resulted in fits that were unable to distinguish between
companion stars that were almost filling their Roche lobes and those that were
only half filling their Roche lobes. Not knowing the companion star radius
makes it difficult to determine if the companion star is a white dwarf, a
brown dwarf or a semi-degenerate helium star (Lazaridis et al., 2011). It also
makes it difficult to predict the future evolution of the system – if the
companion star is close to Roche-lobe filling, much less pulsar energy is
required to drive a wind from its surface, and mass may also be lost through
the inner-Lagrangian point, hastening the evaporation timescale to form an
isolated millisecond pulsar (e.g. Levinson & Eichler 1991). Knowing the Roche-
lobe filling factor also allows the apsidal motion constant of the companion
to be determined from the detection of its gravitational quadrupole moment by
Voisin et al. (2020), enabling the internal structure of the star to be
modelled.

In this paper, we present new simultaneous multi-colour light curves of PSR
J2051$-$0827 and use them to determine the system parameters and the nature of
the companion star.

##  2 Observations and data reduction

In 2011, we observed PSR J2051$-$0827 simultaneously in
$u^{\prime}g^{\prime}i^{\prime}$ using the high-speed, triple-beam camera
ULTRACAM (Dhillon et al., 2007) on the 4.2 m William Herschel Telescope (WHT)
on La Palma. In 2021, we observed PSR J2051$-$0827 again, this time
simultaneously in $u_{\rm s}g_{\rm s}r_{\rm s}i_{\rm s}z_{\rm s}$ using the
high-speed, quintuple-beam camera HiPERCAM (Dhillon et al., 2021) on the 10.4
m Gran Telescopio Canarias (GTC) on La Palma. Note that both filter sets use
the same cut points as the original Sloan Digital Sky Survey (SDSS) filters
(Fukugita et al., 1996), and the primed filters used by ULTRACAM also have
similar bandpass shapes to SDSS, but the HiPERCAM subscript-s filters are more
top-hat in shape and have significantly higher throughputs, particularly in
the $u_{\rm s}$ and $g_{\rm s}$ bands (see Dhillon et al. 2021). A journal of
observations is presented in Table 1. Both instruments were used in full-
frame, no-clear mode, giving a dead time between each frame in HiPERCAM and
ULTRACAM of 0.008 s and 0.024 s, respectively, where each HiPERCAM/ULTRACAM
frame is GPS time-stamped to a relative (i.e. frame-to-frame) accuracy of
0.1/50 $\mu$s and an absolute accuracy of 0.1/1 ms, respectively (Dhillon et
al. 2007, 2021).

The ULTRACAM and HiPERCAM data were reduced using their respective data
reduction pipelines (Dhillon et al. 2007, 2021). All frames were debiased and
then flat-fielded, the latter using the median of twilight-sky frames taken
with the telescope spiralling. The CCD fringing pattern was removed from the
$z_{\rm s}$ HiPERCAM frames using the median of night-sky frames taken with
the telescope spiralling. PSR J2051$-$0827 was invisible in $u^{\prime}$ in
the ULTRACAM data, so this band shall not be discussed further, and we skipped
the readout of every other $u_{\rm s}$ HiPERCAM frame using the NSKIP
parameter (see Dhillon et al. 2021) in order to double the exposure time (and
halve the number of frames) in this band.

We used software apertures that scaled in size with the seeing to extract the
counts from PSR J2051$-$0827 and a number of comparison stars in the same
field of view, the latter acting as the reference for the PSF fits,
transparency and extinction corrections. The comparison stars were also used
for flux calibration via their magnitudes given in the Pan-STARRS1 catalogue
(see Magnier et al. 2020 and references therein) and converted to SDSS
magnitudes (Finkbeiner et al., 2016). The aperture position of PSR
J2051$-$0827 relative to a bright comparison star was determined from a sum of
all the images, and this offset was then held fixed during the reduction so as
to avoid aperture centroiding problems during light-curve minimum. The effect
of atmospheric refraction on the relative aperture positions is negligible due
to the similarity in colour between the target and our chosen reference star,
and the fact that our observations on each night were approximately centred on
meridian transit. The sky level was determined from a clipped mean of the
counts in an annulus surrounding each star and subtracted from the object
counts.

##  3 Results

###  3.1 Light curves

Figure 1: Upper panel: HiPERCAM $u_{\rm s}g_{\rm s}r_{\rm s}i_{\rm s}z_{\rm
s}$ light curves (from bottom to top, the blue/green/orange/red/maroon points,
respectively) of PSR J2051$-$0827, phased according to the ephemeris given in
Eqn. LABEL:eqn:ephemeris. The folded light curve has been repeated three
times, the first cycle with error bars, the second without, and the third
cycle with every 6 data points binned, giving an effective exposure time of
approximately 3 min per point. The solid lines show the icarus fits to the
HiPERCAM data described in Section 3.2. The HiPERCAM data cover approximately
1.08 binary orbits – the ringed points around phase 1.6 show the overlapping
data from the second orbit, demonstrating the excellent agreement with the
data taken during the first orbit. Lower panel: Normalised residuals of the
icarus fits to the HiPERCAM data, offset vertically from each other by 0.6 and
with the $u_{\rm s}$ residuals divided by a factor of 10 for clarity.

PSR J2051$-$0827 was the subject of a 21 yr radio timing study by Shaifullah
et al. (2016). Since then, they have continued to monitor PSR J2051$-$0827 and
provided us with the following up-to-date binary ephemeris:

$\begin{array}[]{lrll}T_{\rm asc}=&{\rm BMJD}\ 59099.9673395&+\
0.09911025846&E\\\ &\pm\ 0.0000019&\pm\ 0.00000000014,&\\\ \end{array}$ (1)

where BMJD refers to the Modified Julian Date on the Barycentric Dynamical
Timescale (TDB) and $E$ is the cycle number. This ephemeris is stable and
precise enough to be applicable to both our HiPERCAM observations in 2021 and
our ULTRACAM observations in 2011. The HiPERCAM and ULTRACAM light curves we
obtained of PSR J2051$-$0827, folded on the above ephemeris, are shown in
Figs. 1 and 2, respectively. Note that $T_{\rm asc}$ in Eqn.
LABEL:eqn:ephemeris corresponds to the epoch of the pulsar’s ascending node.
In what follows, we have chosen to apply a phase offset of $-0.25$ so that
phase 0 corresponds to the superior conjunction of the pulsar.

The light curve of PSR J2051$-$0827 exhibits a single, broad hump that is
brightest at phase 0.5 and faintest at phase 0, indicating that the dominant
source of optical emission is the irradiated inner hemisphere of the companion
star; the neutron star itself is invisible at optical wavelengths in all (non-
transitional) MSPs333Optical pulsations that may originate directly from the
pulsar have been observed in the transitional system PSR J1023+0038 (Ambrosino
et al., 2017).. The $g_{\rm s}$ and $r_{\rm s}$ light curves show evidence for
a possible flaring event just prior to phase 0.5 – such flaring activity has
been seen before in black-widow systems, e.g. Romani (2012). The
$i^{\prime}$-band ULTRACAM light curve (Fig. 2) shows evidence for an
asymmetry, confirming the finding of Stappers et al. (2001), and suggesting
that the leading edge of the companion star is brighter than the trailing
edge, leading to a flux excess after the peak. At first glance, the
$g^{\prime}$-band ULTRACAM light curve does not appear to show this asymmetry,
but it is also present at a lower level (see Section 4.2). The HiPERCAM light
curves obtained a decade later (Fig. 1), on the other hand, are symmetrical.
For reference, the magnitudes at light-curve maximum/minimum of PSR
J2051$-$0827 are given in Table 2. These were determined from a clipped-mean
of the HiPERCAM fluxes in the phase ranges 0.45–0.55 at light-curve maximum
and 0.9–1.1 at minimum – note there was no significant flux detected at light-
curve minimum in $g_{\rm s}$ and $u_{\rm s}$ and so the measured 5$\sigma$
limiting magnitudes are quoted instead.

Figure 2: Upper panel: ULTRACAM $g^{\prime}i^{\prime}$ light curves (green
(bottom) and red (top) points, respectively) of PSR J2051$-$0827, phased
according to the ephemeris given in Eqn. LABEL:eqn:ephemeris. The folded light
curve has been repeated three times, the first cycle with error bars, the
second without, and the third cycle with every 45 data points binned, giving
an effective exposure time of approximately 15 min per point. The solid lines
show the result of fitting the HiPERCAM model presented in Fig. 1 to the
ULTRACAM data, allowing only the flux scaling factors to vary. Lower panel:
Normalised residuals of the icarus fits to the ULTRACAM data, offset
vertically from each other by 4.0.

###  3.2 Light-curve modelling

The multi-colour light curves of PSR J2051$-$0827 were fit using the icarus
modelling software (Breton et al., 2012). The companion star is assumed to be
tidally locked to the pulsar and its surface is modelled as a finite-element
grid, with the intensity of each element calculated from a BT-
Settl444https://phoenix.ens-lyon.fr/Grids/ phoenix model atmosphere (Allard,
2014) appropriate to the physical properties (temperature, gravity, velocity)
at that location on the stellar surface. The model atmospheres are folded
through the HiPERCAM and ULTRACAM filter profiles and the observed flux in
each photometric band is then obtained by integrating the specific intensity
emerging from the surface and visible to an observer located at a given
direction and distance. The best-fit model parameters and their errors are
determined using the MultiNest nested sampling algorithm (Feroz et al., 2013),
as implemented in the Python package PyMultiNest (Buchner et al., 2014).

The input parameters are: the ephemeris ($T_{\rm asc}$, $P$) given in Eqn.
LABEL:eqn:ephemeris; the projected semi-major axis of the pulsar orbit,
$x=a_{\rm psr}\sin i=0.04507837(398)$ lt-s, from the radio-timing data
described in Section 3.1; the gravity darkening exponent, $\beta$. For the
latter, we assumed a value of $\beta=0.08$, appropriate for stars with
convective envelopes, a reasonable assumption for the cool, low-mass companion
in PSR J2051$-$0827 (Stappers et al., 2001). As a check on the robustness of
our fit parameters to this assumption, we also modelled the light curves using
$\beta=0.25$, a value more appropriate for stars with radiative envelopes
(Lucy, 1967), and found negligible ($<1\sigma$) differences in the resulting
parameters.

The fit parameters are as follows:

  * •

$E(g-r)$ – the interstellar reddening. We adopted a Gaussian prior of
$E(g-r)=0.10\pm 0.02$, measured in the direction of PSR J2051$-$0827 from the
3D dust maps of Green et al. (2019)555http://argonaut.skymaps.info, which is
valid for $d>0.73$ kpc. The extinction in each band, $A$, is then calculated
from the extinction vectors, $R$, given by Green et al. (2019).

  * •

$d$ – the distance to the binary. No Gaia or radio timing parallax is
available for PSR J2051$-$0827, so we adopted a prior based on the dispersion
measure, $DM=20.7299$ pc cm-3 (Shaifullah et al., 2016), which corresponds to
a distance of 1.469 kpc using the Galactic free-electron density model of Yao
et al. (2017)666https://www.atnf.csiro.au/research/pulsar/ymw16. We adopted a
log-normal prior on this DM-derived distance, with a fractional error of 0.45
(Yao et al., 2017). Following Clark et al. (2021), we multiplied the DM-
distance prior by two additional priors. First, we adopted a prior based on
the Levin et al. (2013) model for the density of MSPs in the Galactic disk,
which has a Gaussian profile in Galactic radius with width $\sigma=4.5$ kpc,
and an exponential profile in height above the Galactic plane with scale
height $z=0.5$ kpc. Second, we adopted a prior based on the transverse
velocities of binary MSPs, which can be approximated by an exponential
distribution with a mean value of $93\pm 13$ km s-1 (Desvignes et al., 2016).
The radio timing of PSR J2051$-$0827 by Shaifullah et al. (2016) gives a total
proper motion of $\mu=6.1\pm 0.1$ mas yr-1, from which the transverse velocity
can be calculated: $v_{\rm T}=28.9\,d$ km s-1, where $d$ is in kpc. We
multiplied the distance prior by $e^{-v_{\rm T}}/93$ to take this into
account.

  * •

$T_{\rm base}$ – the base temperature of the surface of the companion star,
defined as the temperature at the pole of the star prior to irradiation. To
account for gravity darkening, the base temperature at each point on the
stellar surface is multiplied by the factor $(g/g_{\rm pole})^{\beta}$, where
$g$ is the acceleration due to gravity. $T_{\rm base}$ was constrained to lie
within the range $2100-5000$ K. The lower limit of the model atmospheres we
used is 2300 K; we verified that extrapolating them to 2100 K introduces no
significant systematic errors.

  * •

$T_{\rm irr}$ – the irradiating temperature, which accounts for the effect of
heating by the pulsar. $T_{\rm irr}$ is defined with respect to the centre of
mass of the companion star at a distance $a$ from the pulsar, which would
receive a flux of $\sigma T_{\rm irr}^{4}$, where $\sigma$ is the Stefan-
Boltzmann constant. Hence a location on the stellar surface at a distance $r$
from the pulsar, and which has a normal vector at an angle $\theta$ from the
vector pointing to the pulsar, receives heating power per unit area of $\sigma
T_{\rm irr}^{4}\cos\theta\,a^{2}/r^{2}$. Assuming that the pulsar’s
irradiating flux is immediately thermalised and re-radiated (Breton et al.,
2013), the surface temperature at this location is raised to $T=[T_{\rm
base}^{4}(g/g_{\rm pole})^{4\beta}+T_{\rm
irr}^{4}\cos\theta\,a^{2}/r^{2}]^{1/4}$. We also investigated the approach of
Romani & Sanchez (2016), who applied gravity darkening after rather than
before irradiation, and found that it did not alter the main conclusions of
this paper.

  * •

$f_{\rm RL}$ – the Roche-lobe filling factor, defined as the ratio between the
companion’s radius in the direction towards the pulsar and the distance
between the companion’s centre of mass and inner-Lagrangian point, $L_{\rm
1}$. A uniform prior was adopted, limiting the value between $0.2<f_{\rm
RL}<1$.

  * •

$i$ – the binary inclination angle. A prior that was uniform in $\cos i$ was
adopted to allow for the fact that high inclinations are more likely to be
observed if the binary orbits are randomly oriented.

  * •

$M_{\rm psr}$ – the mass of the pulsar. A uniform prior was adopted, limiting
the value between $1.0<M_{\rm psr}<2.5$ M⊙ (comfortably encompassing all
known, reliable neutron star masses777e.g. https://www3.mpifr-
bonn.mpg.de/staff/pfreire/NS_masses.html). The $M_{\rm psr}$ and $i$ fit
parameters, in conjunction with the $P$ and $x$ input parameters, were used to
derive the mass ratio ($q=M_{\rm psr}/M_{\rm c}$) from the binary mass
function:

$\frac{M_{\rm psr}\sin^{3}i}{(1+1/q)^{2}}=\frac{K_{\rm c}^{3}P}{2\pi
G}=\frac{q^{3}x^{3}4\pi^{2}}{GP^{2}},$ (2)

which in turn allowed the companion mass $M_{\rm c}$ and companion radial
velocity $K_{\rm c}$ to be derived using the above relations. The light curve
constrains $i$, and we chose to fit $M_{\rm psr}$ rather than $q$ or $K_{\rm
c}$ as we can place a more informative prior on $M_{\rm psr}$ based on our
knowledge of the observed neutron star mass distribution than we can on the
relatively unconstrained $q$ and $K_{\rm c}$.

As well as deriving $q$, $M_{\rm c}$ and $K_{\rm c}$ from the fit parameters,
we also used them to derive the volume-averaged radius, $R_{\rm c}$, and
density, $\rho_{\rm c}=M_{\rm c}/\frac{4}{3}\pi R_{\rm c}^{3}$, of the
companion star, and the heating efficiency $\varepsilon=L_{\rm irr}/\dot{E}$.
The latter compares the luminosity of the pulsar, $L_{\rm irr}=4\pi
a^{2}\sigma T_{\rm irr}^{4}$, where $a=x(1+q)/\sin{i}$ is the orbital
separation, to the pulsar’s spin-down power, $\dot{E}=4\pi^{2}I\dot{P}/P^{3}$,
assuming a canonical value for the neutron-star moment of inertia of $I={\rm
10}^{45}$ g cm2 (e.g. Abdo et al. 2013) and values for the pulsar spin period
and its derivative from the radio-timing data described in Section 3.1.

At each step in the sampling process, the input and fit parameters described
above are used to calculate the fluxes received on Earth from the modelled
companion star in each filter. To account for any systematic errors in the
flux calibration, extinction and atmosphere models, we allowed a flux scaling
factor (or equivalently, a magnitude offset) in each band, penalized by a
zero-mean Gaussian prior with a width of 0.1 mag. The latter value represents
the uncertainty in our flux calibration given that we did not account for the
colour terms relating to the differences between and within the HiPERCAM
Super-SDSS and ULTRACAM SDSS filter systems compared to SDSS (for a detailed
discussion of these differences, see Brown et al. 2022 and Wild et al. 2022,
respectively). Without such a prior, we would obtain unrealistically small
errors in the distance and reddening due to the degeneracy between these
parameters and the flux scaling factor. For reasons that will become apparent
below, we chose to fit the HiPERCAM $u_{\rm s}g_{\rm s}r_{\rm s}i_{\rm
s}z_{\rm s}$ and ULTRACAM $g^{\prime}i^{\prime}$ data separately, but for each
instrument we fit all of the filters simultaneously.

The icarus fits to the HiPERCAM $u_{\rm s}g_{\rm s}r_{\rm s}i_{\rm s}z_{\rm
s}$ light curves are shown as the solid curves in the upper panel of Fig. 1,
the resulting fit and derived parameters are given in Table 3, and the
posterior distributions of these parameters are shown in Fig. 3. The latter
figure shows that, of the three parameters with non-uniform priors, only the
prior on $E(g-r)$ has a noticeable effect on the posterior distributions
compared to a uniform prior. $E(g-r)$ is essentially unconstrained by the
data, and so the posterior distribution of this parameter closely follows the
prior, but $E(g-r)$ is slightly correlated with $T_{\rm irr}$ and $T_{\rm
base}$, and so the effect of the Gaussian prior on this parameter is to
slightly reduce the uncertainties on these temperatures. The posterior
distributions of $i$ and $d$, on the other hand, are much narrower than the
priors, indicating that their inferred values are dominated by the data and
not by the priors. Judging from the fit residuals shown in the lower panel of
Fig. 1, the reduced-$\chi^{2}$ value of 0.96, the near-unity flux scaling
factors given in Table 3 (which all lie well within the 0.1 mag Gaussian
priors we set), and the shape of the posterior distributions in Fig. 3, the
fit to the HiPERCAM data is acceptable, and the resulting parameter values
will be discussed in Section 4.1.

Table 3: Results of the icarus fits to the HiPERCAM $u_{\rm s}g_{\rm s}r_{\rm
s}i_{\rm s}z_{\rm s}$ light curves. The parameter values quoted are the median
of the marginalised posterior distributions shown in Fig. 3, with the 95%
confidence region given in sub- and superscript. $\chi^{2}$/degrees of freedom
= 1287.4/1341
---
$u_{\rm s}g_{\rm s}r_{\rm s}i_{\rm s}z_{\rm s}$ flux scaling factors =
1.06/0.98/0.96/0.98/0.94
Fit parameters
$E(g-r)$ | $0.096^{+0.038}_{-0.037}$
$d$ (kpc) | $2.48^{+0.39}_{-0.38}$
$T_{\rm base}$ (K) | $2750^{+130}_{-150}$
$T_{\rm irr}$ (K) | $5040^{+210}_{-200}$
$f_{\rm RL}$ | $0.88^{+0.02}_{-0.02}$
$i$ (∘) | $55.9^{+4.8}_{-4.1}$
$M_{\rm psr}$ (M⊙) | $1.77^{+0.69}_{-0.73}$
Derived parameters
$M_{\rm c}$ (M⊙) | $0.039^{+0.010}_{-0.011}$
$R_{\rm c}$ (R⊙) | $0.139^{+0.011}_{-0.015}$
$\rho_{\rm c}$ (g cm-3) | $20.24^{+0.59}_{-0.44}$
$K_{\rm c}$ (km s-1) | $454^{+66}_{-81}$
$q\equiv M_{\rm psr}/M_{\rm c}$ | $45.8^{+6.7}_{-8.2}$
$\varepsilon$ | $0.51^{+0.18}_{-0.17}$
$k_{2}$ | $0.0036-0.0047$
Figure 3: Corner plot showing the posterior distributions of the light-curve
fit parameters. The last six parameters ($M_{\rm c}$, $R_{\rm c}$, $\rho_{\rm
c}$, $K_{\rm c}$, $q$, $\varepsilon$) were derived from the values of the fit
parameters and the input parameters. The contours indicate the 1$\sigma$,
2$\sigma$ and 3$\sigma$ confidence regions. The dashed vertical lines on the
histograms show the mean and 95% confidence intervals of the parameter
distributions, and where non-uniform priors were used, these are shown as
curves on the histograms.

However, the same can not be said for the icarus fits to the ULTRACAM
$g^{\prime}i^{\prime}$ data. Like Stappers et al. (2001), we found that the
asymmetry in the light curves prevented a reliable model from being found.
This is unsurprising given that the simple, direct-heating model described
above can only produce symmetric light curves. So we experimented with
different modifications to the direct-heating model to fit the asymmetry:
allowing for a small orbital phase offset (e.g. Nieder et al. 2019); using
spherical harmonics or hot and cold and spots to account for brightness
asymmetries on the stellar surface (e.g. Clark et al. 2021); allowing for heat
redistribution due to convective flows on the stellar surface (e.g. Voisin et
al. 2020, Stringer et al. 2021). Although acceptable fits could be obtained by
adding free parameters in this way, the various methods did not point to a
consistent set of fit parameters: in particular, the best-fit inclinations and
Roche-lobe filling factors varied widely between models.

Given that we do not know which of the above methods of introducing light-
curve asymmetries in PSR J2051$-$0827 is physically correct, and that the
HiPERCAM data are of higher quality than the ULTRACAM data and do not show the
asymmetry, we decided to disregard the ULTRACAM fits. Instead, we assume that
the HiPERCAM light curve represents the underlying system (i.e. no asymmetric
heating or spots), and we adopt the HiPERCAM fit as correctly representing the
binary parameters. For reference, we took the best-fit model to the HiPERCAM
light curve and fit it to the ULTRACAM data, allowing only the flux scaling
factors to vary. The resulting fit to the ULTRACAM data ($\chi^{2}=5707$,
degrees of freedom $=4492$) is shown as the solid curves in Fig. 2 and
required flux scaling factors of $g^{\prime}=0.82$ and $i^{\prime}=0.93$, in
the sense that the best-fit HiPERCAM model is multiplied by these factors to
match the ULTRACAM fluxes. It can be seen that the HiPERCAM model provides a
reasonable match to the $g^{\prime}$ ULTRACAM light curve, but the fit to the
$i^{\prime}$ light curve is poor, particularly after phase 0.5 due to the
extra flux from the leading edge of the companion star. We shall discuss the
possible origins of this extra flux in Section 4.2.

##  4 Discussion

###  4.1 System parameters

With the exception of the Roche-lobe filling factor, which shall be discussed
further below, our light-curve fit parameters are in reasonable agreement with
those of Stappers et al. (2001). We are able to confirm that PSR J2051$-$0827
has only a moderate inclination of $55.9^{+4.8}_{-4.1}$ degrees, slightly
higher than the value of $i\sim 40^{\circ}$ found by Stappers et al. (2001)
and more consistent with the model of plasma-lensing during radio eclipse
derived by Lin et al. (2021). Like Stappers et al. (2001), we find that the
minimum night-side temperature ($\sim T_{\rm base}$) of the companion star is
cool ($2750^{+130}_{-150}$ K). This is close to the 2300 K minimum of our
model-atmosphere grid, which could be the cause of the slight over-prediction
of the $z_{\rm s}$ flux at minimum evident in Fig. 1. In contrast, the maximum
day-side temperature ($\sim[T_{\rm base}^{4}+T_{\rm irr}^{4}]^{1/4}$) is
$5150^{+190}_{-190}$ K. We find that approximately half of the pulsar’s spin-
down energy is converted to optical luminosity, which is somewhat higher than
the value of $\sim 20$% typically observed in black-widow systems (Breton et
al. 2013, Draghis et al. 2019).

PSR J2051$-$0827 is detected in $\gamma$-rays by the Fermi Large Area
Telescope (Wu et al., 2012), which measures an integrated energy flux above
100 MeV of $F_{\gamma}=(2.5\pm 0.3)\times 10^{-12}$ erg cm-2 s-1 (Fermi-LAT
collaboration et al., 2022). For our estimated $d=2.5\pm 0.2$ kpc (1$\sigma$
uncertainties), this corresponds to a luminosity of $L_{\gamma}=4\pi
F_{\gamma}d^{2}=(1.9\pm 0.4)\times 10^{33}$ erg s-1. This $\gamma$-ray
emission therefore represents $\sim 35$% of the pulsar’s spin-down power
budget, which is at the lower limit of the range of $\epsilon$ inferred from
our Icarus model. Beaming effects may result in the $\gamma$-ray flux being
stronger at the pulsar’s spin equator, which should be aligned with the
orbital plane since the pulsar is believed to have been spun-up by accretion,
than when viewed at $i\sim 56^{\circ}$, and so the pulsar’s $\gamma$-ray
emission may be sufficient to explain the irradiation in this system.

Without a measurement of the radial velocity of the companion star we are
unable to determine the pulsar mass, but by assuming it lies in the range
$1.0<M_{\rm psr}<2.5$ M⊙, the light-curve fit indicates that the companion
star has a mass in the range $0.039^{+0.010}_{-0.011}$ M⊙, similar to the
values derived for other black-widow systems (Roberts 2013, Draghis et al.
2019). The Roche-lobe filling factor lies in the range $0.88^{+0.02}_{-0.02}$,
resolving the ambiguity noted by Stappers et al. (2001) in favour of a
companion star that is close to filling its Roche lobe. The resulting volume-
averaged radius of the companion star ($0.139^{+0.011}_{-0.015}$ $R_{\odot}$)
implies a mean density of $20.24^{+0.59}_{-0.44}$ g cm-3. The mass and mean
density of the companion star in PSR J2051$-$0827 are consistent with those of
a hydrogen brown dwarf rather than a helium or carbon white dwarf (Tang et al.
2014, Hatzes & Rauer 2015, Kaplan et al. 2018).

With a reliable estimate of the Roche-lobe filling factor, we can now
determine the apsidal motion constant, $k_{2}$, from the orbital precession
measurement of PSR J2051$-$0827 by Voisin et al. (2020), constraining it to
the range $0.0036<k_{2}<0.0047$.888This value has been calculated using the
mass ratio given in Table 3, not by using Fig. 1 of Voisin et al. (2020). Note
also that the ordinate of the latter figure is incorrect due to an error of a
factor of 10 in the assumed mass ratio. The apsidal motion constant describes
how centrally condensed an object is, with higher central concentrations
corresponding to smaller values of $k_{2}$. Our measured value is two orders
of magnitude smaller than those typical of brown dwarfs and the gas giants of
the solar system (Heller et al., 2010). Detailed interpretation of our $k_{2}$
measurement requires integration of the equations of stellar structure and is
outside the scope of this paper.

###  4.2 Light-curve asymmetry

Asymmetric optical light curves have been seen in a number of other black-
widow systems, including PSR J1810+1744 (Romani et al., 2021), PSR
J1311$-$3430 (Romani et al., 2015), PSR J1653$-$0158 (Nieder et al., 2020),
PSR J0952$-$0607 (Nieder et al., 2019) and PSR J1959+2048 (Kandel & Romani,
2020). Models for the origin of the asymmetry include ducting of the particles
in the shock between the pulsar and companion star winds (the intra-binary
shock, IBS) onto the magnetic poles of the companion star (Sanchez & Romani,
2017), re-processed radiation from the IBS (Romani & Sanchez, 2016), hot and
cold spots on the companion star’s surface (van Staden & Antoniadis 2016,
Clark et al. 2021), and redistribution of energy on the companion star’s
surface by convection (Voisin et al. 2020, Kandel & Romani 2020).

Figure 4: HiPERCAM $g_{\rm s}i_{\rm s}$ (left) and ULTRACAM
$g^{\prime}i^{\prime}$ (right) light curves. To emphasize any asymmetries
present, the black points show the light curves mirrored around phase 0.5
(indicated by the vertical dashed line).

The $i^{\prime}$-band ULTRACAM light curve presented in Fig. 2 is clearly
asymmetric around phase 0.5. By comparing this light curve to the icarus fit
and the simultaneous $g^{\prime}$-band ULTRACAM light curve plotted in the
same figure, we can infer that the source of the asymmetry is predominantly
excess $i^{\prime}$-band light on the leading hemisphere of the companion
star. This asymmetry appears to be absent in the HiPERCAM light curves
obtained a decade later, as highlighted in Fig. 4. A close inspection of Fig.
4 shows that the asymmetry is actually present in the $g^{\prime}$-band
ULTRACAM light curve as well, but is weaker and covers a narrower phase range
than the $i^{\prime}$-band asymmetry.

Similar long-term changes in the light-curve symmetry of MSP binaries have
been observed before, e.g. in PSR J2339$-$0533 (Kandel et al., 2020) and PSR
J1723$-$2837 (van Staden & Antoniadis, 2016). These are both so-called ‘red-
back’ systems, which harbour more massive companion stars than black-widow
systems (Roberts, 2013), and the variations were accompanied by significant
changes in the overall luminosity of the system. In the case of PSR
J2051$-$0827, however, the flux scaling factors in the $g$ and $i$-bands are
consistent in the ULTRACAM and HiPERCAM data, to within the uncertainties in
the flux calibration (see Section 3.2), a fact that is also evident from the
direct comparison of the two light curves shown in Figure 4. This implies that
the luminosity at maximum remained approximately the same, implying a constant
irradiation power in the two sets of observations and ruling out any
significant variation in the pulsar emission or IBS, assuming that the latter
mediates the irradiation as proposed by Romani & Sanchez (2016). Luminosity
variations intrinsic to the companion star may also occur, for example due to
the Applegate mechanism (Applegate 1992, Applegate & Shaham 1994), but since
the luminosity at light-curve maximum is dominated by the irradiation power
this is unlikely to be detectable here. Besides, if intrinsic variations could
in principle be detected at minimum luminosity, none have been reported in
black-widow systems so far, to the best of our knowledge. Therefore, if one
assumes that the extrinsic (irradiation) and intrinsic power sources remain
stable, it seems plausible that the asymmetry observed in the ULTRACAM
observations results from a variation of the heat redistribution pattern on
the surface of the companion star, and the later HiPERCAM observations seem to
indicate that this redistribution pattern alters on timescales of a decade or
less. Repeat observations of the kind we present in this paper but with
significantly shorter intervals between them may help to determine if the
timescale on which the asymmetry appears is significantly shorter than this
and whether or not it is periodic.

Finally, our assumption that the symmetric HiPERCAM light curves represent the
underlying system provides us with a potential opportunity to compare and
constrain the various asymmetric heating models, free of the degeneracies
between the parameters that plagued the fits to the ULTRACAM data described at
the end of Section 3.2. To this end, the posterior distributions of the fit
parameters from the HiPERCAM photometry were modelled using a Gaussian mixture
model and then used as priors for the fit parameters when fitting the ULTRACAM
photometry (see Kennedy et al. 2022). We tried fitting both the spot model of
Clark et al. (2021) and the convection model of Voisin et al. (2020), but we
were unable to obtain satisfactory fits with either model, particularly of the
pre-maximum portion of the light curve. Either additional free parameters or a
different parameterisation of the spot and convection models will be required
to adequately model the ULTRACAM data.

##  5 Conclusions

Our light curves of PSR J2051$-$0827 have demonstrated that whatever mechanism
is responsible for the asymmetric heating observed on the companion stars in
black-widow systems is likely to be variable on timescales of a decade or less
and is most probably related to a change in the heat redistribution pattern on
the stellar surface rather than any change in the irradiation power. We find
that the companion star in PSR J2051$-$0827 is close to filling its Roche lobe
and has a mass and mean density consistent with a brown dwarf, but an apsidal
motion constant that implies a significantly more centrally-condensed internal
structure than is typical of such objects. We encourage continued monitoring
of PSR J2051$-$0827 to put tighter constraints on the heating-variability
timescale, detailed modelling of the internal structure of the companion star
to constrain the evolutionary history of the binary, and spectroscopy of the
companion star, perhaps with JWST or the coming generation of extremely large
telescopes, to measure the pulsar mass.

## Acknowledgements

We thank the anonmyous referee for their comments on the manuscript. The
design and construction of HiPERCAM was funded by the European Research
Council under the European Union’s Seventh Framework Programme (FP/2007-2013)
under ERC-2013-ADG Grant Agreement no. 340040 (HiPERCAM). VSD, ULTRACAM and
HiPERCAM operations are funded by the Science and Technology Facilities
Council (grant ST/V000853/1). RPB, CJC, DMS, MRK and GV acknowledge support
from the ERC under the European Union’s Horizon 2020 research and innovation
programme under Grant Agreement No. 715051 (Spiders). MRK acknowledges support
from the Irish Research Council in the form of a Government of Ireland
Postdoctoral Fellowship (GOIPD/2021/670: Invisible Monsters). DMS acknowledges
the Fondo Europeo de Desarrollo Regional (FEDER) and the Canary Islands
government for the financial support received in the form of grant number
PROID2020010104. SGP acknowledges the support of a STFC Ernest Rutherford
Fellowship. IP and TRM acknowledge support from the UK’s Science and
Technology Facilities Council (STFC), grant ST/T000406/1. The GTC is installed
at the Spanish Observatorio del Roque de los Muchachos (ORM) of the Instituto
de Astrofísica de Canarias (IAC), on the island of La Palma. The WHT is
operated on La Palma by the Isaac Newton Group in the Spanish ORM of the IAC.
For the purpose of open access, the author has applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

## Data availability

The data presented in this paper will be shared on reasonable request to the
corresponding author.

## References

  * Abdo et al. (2013) Abdo A. A., et al., 2013, ApJS, 208, 17
  * Allard (2014) Allard F., 2014, in Booth M., Matthews B. C., Graham J. R., eds, IAU Symposium Vol. 299, Exploring the Formation and Evolution of Planetary Systems. pp 271–272, doi:10.1017/S1743921313008545
  * Ambrosino et al. (2017) Ambrosino F., et al., 2017, Nature Astronomy, 1, 854
  * Applegate (1992) Applegate J. H., 1992, ApJ, 385, 621
  * Applegate & Shaham (1994) Applegate J. H., Shaham J., 1994, ApJ, 436, 312
  * Breton et al. (2012) Breton R. P., Rappaport S. A., van Kerkwijk M. H., Carter J. A., 2012, ApJ, 748, 115
  * Breton et al. (2013) Breton R. P., et al., 2013, ApJ, 769, 108
  * Brown et al. (2022) Brown A. J., et al., 2022, MNRAS, 513, 3050
  * Buchner et al. (2014) Buchner J., et al., 2014, A&A, 564, A125
  * Clark et al. (2021) Clark C. J., et al., 2021, MNRAS, 502, 915
  * Desvignes et al. (2016) Desvignes G., et al., 2016, MNRAS, 458, 3341
  * Dhillon et al. (2007) Dhillon V. S., et al., 2007, MNRAS, 378, 825
  * Dhillon et al. (2021) Dhillon V. S., et al., 2021, MNRAS, 507, 350
  * Draghis et al. (2019) Draghis P., Romani R. W., Filippenko A. V., Brink T. G., Zheng W., Halpern J. P., Camilo F., 2019, ApJ, 883, 108
  * Fermi-LAT collaboration et al. (2022) Fermi-LAT collaboration et al., 2022, arXiv e-prints, p. arXiv:2201.11184
  * Feroz et al. (2013) Feroz F., Hobson M. P., Cameron E., Pettitt A. N., 2013, Importance Nested Sampling and the MultiNest Algorithm (arXiv:1306.2144)
  * Finkbeiner et al. (2016) Finkbeiner D. P., et al., 2016, ApJ, 822, 66
  * Fruchter et al. (1988) Fruchter A. S., Stinebring D. R., Taylor J. H., 1988, Nature, 333, 237
  * Fukugita et al. (1996) Fukugita M., Ichikawa T., Gunn J. E., Doi M., Shimasaku K., Schneider D. P., 1996, AJ, 111, 1748
  * Green et al. (2019) Green G. M., Schlafly E. F., Zucker C., Speagle J. S., Finkbeiner D. P., 2019, A 3D Dust Map Based on Gaia, Pan-STARRS 1 and 2MASS (arXiv:1905.02734)
  * Hatzes & Rauer (2015) Hatzes A. P., Rauer H., 2015, ApJ, 810, L25
  * Heller et al. (2010) Heller R., Jackson B., Barnes R., Greenberg R., Homeier D., 2010, A&A, 514, A22
  * Kandel & Romani (2020) Kandel D., Romani R. W., 2020, ApJ, 892, 101
  * Kandel et al. (2020) Kandel D., Romani R. W., Filippenko A. V., Brink T. G., Zheng W., 2020, ApJ, 903, 39
  * Kaplan et al. (2018) Kaplan D. L., Stovall K., van Kerkwijk M. H., Fremling C., Istrate A. G., 2018, ApJ, 864, 15
  * Kennedy et al. (2022) Kennedy M. R., et al., 2022, MNRAS, 512, 3001
  * Lazaridis et al. (2011) Lazaridis K., et al., 2011, MNRAS, 414, 3134
  * Levin et al. (2013) Levin L., et al., 2013, MNRAS, 434, 1387
  * Levinson & Eichler (1991) Levinson A., Eichler D., 1991, ApJ, 379, 359
  * Lin et al. (2021) Lin F. X., Main R. A., Verbiest J. P. W., Kramer M., Shaifullah G., 2021, MNRAS, 506, 2824
  * Lorimer (2008) Lorimer D. R., 2008, Living Reviews in Relativity, 11, 8
  * Lucy (1967) Lucy L. B., 1967, Z. Astrophys., 65, 89
  * Magnier et al. (2020) Magnier E. A., et al., 2020, ApJS, 251, 6
  * Manchester et al. (2005) Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993
  * Nieder et al. (2019) Nieder L., et al., 2019, ApJ, 883, 42
  * Nieder et al. (2020) Nieder L., et al., 2020, ApJ, 902, L46
  * Polzin et al. (2020) Polzin E. J., Breton R. P., Bhattacharyya B., Scholte D., Sobey C., Stappers B. W., 2020, MNRAS, 494, 2948
  * Roberts (2013) Roberts M. S. E., 2013, in van Leeuwen J., ed., IAU Symposium Vol. 291, Neutron Stars and Pulsars: Challenges and Opportunities after 80 years. pp 127–132 (arXiv:1210.6903), doi:10.1017/S174392131202337X
  * Romani (2012) Romani R. W., 2012, ApJ, 754, L25
  * Romani & Sanchez (2016) Romani R. W., Sanchez N., 2016, ApJ, 828, 7
  * Romani et al. (2015) Romani R. W., Filippenko A. V., Cenko S. B., 2015, ApJ, 804, 115
  * Romani et al. (2021) Romani R. W., Kandel D., Filippenko A. V., Brink T. G., Zheng W., 2021, ApJ, 908, L46
  * Sanchez & Romani (2017) Sanchez N., Romani R. W., 2017, ApJ, 845, 42
  * Shaifullah et al. (2016) Shaifullah G., et al., 2016, MNRAS, 462, 1029
  * Stappers et al. (1996a) Stappers B. W., et al., 1996a, ApJ, 465, L119
  * Stappers et al. (1996b) Stappers B. W., Bessell M. S., Bailes M., 1996b, ApJ, 473, L119
  * Stappers et al. (1999) Stappers B. W., van Kerkwijk M. H., Lane B., Kulkarni S. R., 1999, ApJ, 510, L45
  * Stappers et al. (2001) Stappers B. W., van Kerkwijk M. H., Bell J. F., Kulkarni S. R., 2001, ApJ, 548, L183
  * Stringer et al. (2021) Stringer J. G., et al., 2021, MNRAS, 507, 2174
  * Tang et al. (2014) Tang S., et al., 2014, ApJ, 791, L5
  * Tauris & van den Heuvel (2006) Tauris T. M., van den Heuvel E. P. J., 2006, in Compact stellar X-ray sources. pp 623–665
  * Voisin et al. (2020) Voisin G., Clark C. J., Breton R. P., Dhillon V. S., Kennedy M. R., Mata-Sánchez D., 2020, MNRAS, 494, 4448
  * Wild et al. (2022) Wild J. F., et al., 2022, MNRAS, 509, 5086
  * Wu et al. (2012) Wu J. H. K., Kong A. K. H., Huang R. H. H., Takata J., Tam P. H. T., Wu E. M. H., Cheng K. S., 2012, ApJ, 748, 141
  * Yao et al. (2017) Yao J. M., Manchester R. N., Wang N., 2017, ApJ, 835, 29
  * van Staden & Antoniadis (2016) van Staden A. D., Antoniadis J., 2016, ApJ, 833, L12