# An Optical Study of the Black Widow Population

D. Kandel Department of Physics, Stanford University, Stanford, CA, 94305, USA
Roger W. Romani Department of Physics, Stanford University, Stanford, CA,
94305, USA

(Received 2022 June)

###### Abstract

The optical study of the heated substellar companions of ‘Black Widow’ (BW)
millisecond pulsars (MSP) provides unique information on the MSP particle and
radiation output and on the neutron star mass. Here we present analysis of
optical photometry and spectroscopy of a set of relatively bright BWs, many
newly discovered in association with Fermi $\gamma$-ray sources. Interpreting
the optical data requires sophisticated models of the companion heating. We
provide a uniform analysis, selecting the preferred heating model and
reporting on the companion masses and radii, the pulsar heating power and
neutron star mass. The substellar companions are substantially degenerate,
with average densities $15-30\times$ Solar, but are inflated above their zero
temperature radii. We find evidence that the most extreme recycled BW pulsars
have both large $>0.8M_{\odot}$ accreted mass and low $<10^{8}$G magnetic
fields. Examining a set of heavy BWs, we infer that neutron star masses larger
than $2.19M_{\odot}$ ($1\sigma$ confidence) or $2.08M_{\odot}$ ($3\sigma$
confidence) are required; these bounds exclude all but the stiffest equations
of state in standard tabulations.

pulsars: general — pulsars: individual (PSR J0023$+$0923, J0636$+$5128,
J0952$-$0607, J1301$+$0833, J1311$-$3430, J1653$-$0158, J1810$+$1744,
J1959$+$2048, J2052$+$1219)

##  1 Introduction

Black widows (BWs) with sub-stellar companions and redbacks (RBs) with low-
mass star companions, together known as ‘spiders’, are binary systems of MSP
in tight $<1$ d orbits, with the companion heated and evaporated by the pulsar
spindown power. Since the discovery of BW PSR J1959+2048’s companion(Fruchter
et al., 1988; Kulkarni et al., 1988), it has been known that the heating
pattern and the resulting light curve are important probes of the binary
geometry (Djorgovski & Evans 1988, Callanan et al. 1995) and the mass of the
MSP (Aldcroft et al. 1992, van Kerkwijk et al. 2011). These pulsars are
difficult to discover in the radio band, due to scattering and obscuration by
the evaporation wind. The pulsars’ peak photon output is in the penetrating
GeV gamma-rays and with the advent of the Fermi LAT sky survey and attendant
follow-up searches, the number of known ‘spiders’ has greatly increased. We
now have detailed optical light curves of many of these objects (Draghis et
al., 2019). Photometric data and its modeling provide an important probe of
the physics of these systems and can constrain the size and heating of the
companion. Most importantly, such modeling constrains the binary inclination
$i$, which together with spectroscopic measurements can help determine the
pulsar and companion masses.

Traditional modeling of these spider systems assumes direct pulsar irradiation
of the companion, which then spontaneously re-emits thermal radiation. For
such models, the optical light curve is symmetric with its maximum at pulsar
inferior conjunction. While some objects show optical modulation broadly
consistent with this picture, many light curves, especially those with higher
S/N, show significant peak asymmetries (Stappers et al., 2001) and phase-
shifts (Schroeder & Halpern, 2014). Past work has discussed the possible
physical origin of such asymmetries, including companion irradiation by an
asymmetric intrabinary shock (IBS; Romani & Sanchez 2016), IBS particles
ducting along magnetic field lines to companion poles (Sanchez & Romani,
2017), heat transfer from a global wind (Kandel & Romani 2020, hereafter KR20)
or general surface diffusion (Voisin et al., 2020). As described in Kandel &
Romani (2020), and Romani et al. (2021), properly accounting for these effects
is important to get an unbiased estimate of the NS mass. We caution that, in
many cases, masses based on simple direct heating modeling will not be
accurate.

In this paper, we discuss ten BW systems, presenting uniform LC modeling of
nine for which we have photometric data. Six out of these also have radial
velocity data, which gives us an opportunity to estimate the system masses. We
show that BW systems host heavy neutron stars, and by combining the inferred
neutron star masses, we can put a strong lower bound on the maximal neutron
star mass, appreciably above the values inferred from radio Shapiro delay
measurements. This should have significant ramifications for the equation of
state of dense matter.

Table 1: Summary of BW parameters and photometric observations Name |  $P_{s}$ (ms) |  $P_{b}$ (hr) |  $\dot{E}(10^{34}$erg/s) |  $x_{1}$(lt-s) |  $A_{V}$(mag) |  $B\,(10^{8}$G) | Bands | Ref.
---|---|---|---|---|---|---|---|---
J$0023+0923$ | 3.05 | 3.33 | 1.60 | 0.035 | 0.382 | 1.90 | $gri$ | Draghis et al. (2019)
J$0636+5128$ | 2.87 | 1.60 | 0.58 | 0.009 | 0.218 | 1.01 | $grizHK$ | Draghis & Romani (2018)
J$0952-0607$ | 1.41 | 6.42 | 6.65 | 0.063 | 0.137 | 0.82 | $ugriz$ | Romani et al. (2022)
J$1301+0833$ | 1.84 | 6.48 | 6.65 | 0.078 | 0.082 | 1.41 | $griz$ | Draghis et al. (2019)
J$1311-3430$ | 2.56 | 1.56 | 5.00 | 0.011 | 0.137 | 2.36 | $ugriz$ | Romani et al. (2015)
J$1653-0158$ | 1.97 | 1.25 | 1.20 | 0.011 | 0.710 | 0.68 | $u^{\prime}g^{\prime}r^{\prime}i^{\prime}$ | Nieder et al. (2020)
J$1810+1744$ | 1.66 | 3.56 | 3.97 | 0.095 | 0.390 | 0.88 | $ugriz$ | Romani et al. (2021)
J$1959+2048$ | 1.61 | 9.17 | 16.0 | 0.089 | 0.600 | 1.66 | $BVRIK$ | Reynolds et al. (2007)
J$2051-0827$ | 4.51 | 2.37 | 0.55 | 0.045 | 0.296 | 2.42 | $ugriz$ | Dhillon et al. (2022)
J$2052+1219$ | 1.99 | 2.75 | 3.37 | 0.061 | 0.328 | 1.17 | $gri$ | Draghis et al. (2019)

##  2 Observations

Table 1 summarises some basic properties of the ten BWs in our population
study. It also shows the photometric bands used; the acquisition and
processing of these data have been described in earlier papers. Photometry of
five systems – J0023+0923, J0636+5128, J1301+0833, J1959+2048 and J2052+1219 –
is described in Draghis et al. (2019) and references therein. Optical
photometry and spectroscopy of J0952$-$0607 is described in Romani et al.
(2022a). For J1311$-$3430, we use synthesized photometry from LRIS spectra, as
described in Romani et al. (2015). For J1653$-$0158, we use ULTRACAM data from
Nieder et al. (2020). Optical observations of J1810+1744 are described in
Romani et al. (2021). Photometry and fitting of PSR J2051$-$0827 are presented
in Dhillon et al. (2022).

##  3 Photometric Fitting

Table 2: Light Curve Fit Results Name | $i\,(\mathrm{deg})$ | $f$ | $T_{N}\,(\mathrm{K})$ | $L_{\mathrm{H}}\,(10^{33}\,\mathrm{erg/s})$ | $d\,({\rm kpc})$ | $\theta_{\rm hs}\,(\mathrm{deg})$ | $\phi_{\rm hs}\,(\mathrm{deg})$ | $\mathcal{A}_{\rm hs}$ | $\sigma_{\rm hs}\,(\mathrm{deg})$ | $\epsilon$ | $\chi^{2}/{\rm DoF}$
---|---|---|---|---|---|---|---|---|---|---|---
J$0023+0923$ | $74.1^{+8.8}_{-12}$ | $0.87^{+0.04}_{-0.05}$ | $2992^{+94}_{-92}$ | $1.22^{+0.30}_{-0.30}$ | $0.59^{+1.39}_{-0.39}$ | $353.0^{+2.7}_{-5.5}$ | $-22.6^{+75.8}_{-42.1}$ | $0.62^{+0.16}_{-0.15}$ | $11.4^{+3.8}_{-2.6}$ | $-$ | $105/53$
J$0636+5128$ | $34.2^{+2.9}_{-2.5}$ | $0.88\pm 0.02$ | $2284^{+100}_{-100}$ | $1.9^{+0.3}_{-0.3}$ | $0.63^{+1.44}_{-0.43}$ | $258.2^{+11.3}_{-12.4}$ | $66.1^{+6.9}_{-8.9}$ | $1.0^{+0.10}_{-0.13}$ | $11.0^{+3.0}_{-2.2}$ | - | $235/155$
J$0952-0607$ | $59.8^{+2.0}_{-1.9}$ | $0.79\pm 0.01$ | $3085^{+85}_{-80}$ | $3.81^{+0.46}_{-0.43}$ | $6.26^{+0.36}_{-0.40}$ | - | - | - | - | - | 286/287
J$1301+0833$ | $46.6^{+2.5}_{-2.2}$ | $0.61\pm 0.06$ | $2288^{+80}_{-84}$ | $4.7^{+0.5}_{-0.5}$ | $1.77^{+0.10}_{-0.11}$ | - | - | - | - | - | 157/121
J$1311-3430$ | $68.7^{+2.1}_{-2.0}$ | $0.99$ | $821^{+959}_{-402}$ | $228^{+56}_{-44}$ | $3.01\pm 0.15$ | $119.8^{+4.3}_{-4.3}$ | $27.5^{+15.7}_{-17.5}$ | $11.2^{+10.6}_{-5.9}$ | $16.8^{+3.4}_{-3.1}$ | $-0.072^{+0.001}_{-0.001}$ | $74/84$
J$1653-0158$ | $72.8^{+4.0}_{-4.0}$ | $0.84^{+0.02}_{-0.02}$ | $2253^{+499}_{-389}$ | $4.3^{+0.4}_{-0.6}$ | $0.87^{+0.07}_{-0.08}$ | $287.1^{+5.7}_{-5.2}$ | $-48.1^{+6.9}_{-5.2}$ | $1.23^{+0.47}_{-0.36}$ | $21.8^{+3.4}_{-3.2}$ | - | 1296/941
J$1810+1744$ | $66.3^{+0.5}_{-0.5}$ | 0.99 | $3101^{+77}_{-94}$ | $41.8^{+3.7}_{-4.8}$ | $2.86\pm 0.13$ | $86.9^{+1.1}_{-0.8}$ | $23.0^{+3.2}_{-2.6}$ | $0.86\pm 0.23$ | $23.4^{+1.2}_{-1.1}$ | $-0.19\pm 0.02$ | 229/168
J$1959+2048^{*}$ | $85.3^{+2.9}_{-1.2}$ | $0.74^{+0.03}_{-0.04}$ | $2710^{+25}_{-23}$ | $42.7^{+2.7}_{-2.2}$ | $2.25\pm 0.11$ | - | - | - | - | - | 128/88
J$2051-0837^{*}$ | $55.9^{+4.8}_{-4.1}$ | $0.88^{+-.02}_{-0.02}$ | $2750^{+120}_{-150}$ | $1.7^{+0.3}_{-0.3}$ | $2.48^{+0.39}_{-0.38}$ | - | - | - | - | - | 1288/1341
J$2052+1219$ | $59.3^{+2.3}_{-1.8}$ | $0.99$ | $2710^{+332}_{-1218}$ | $11.8^{+3.5}_{-3.4}$ | $5.57^{+0.47}_{-0.55}$ | - | - | - | - | $0.03\pm 0.01$ | 153/80

Our fits are performed with an outgrowth of the ICARUS light-curve model
(Breton et al., 2012) with additions described by Kandel et al. (2020). An
additional update replaces the simplified limb-darkening laws in the base code
with the more detailed limb-darkening coefficients computed by Claret &
Bloemen (2011) for two models, the ATLAS and PHOENIX atmospheres. We generally
find that the ATLAS coefficients perform better. Note that gravity- and limb-
darkening serve to rescale the local fluxes; to monitor the heating budget we
take care to integrate the emergent flux to determine the total companion
heating $L_{H}$, subtracting the thermal base emission (characterized by the
night-side temperature $T_{N}$). In all fits, we explore a simple direct
heating (DH) model, a wind model (WH), and a hotspot model (HS) and perform
model selection using Akaike Information Criterion (AIC, Akaike, 1974). In
general, asymmetry of the light-curve maximum and bluer colors to one side of
the peak indicate can the presence of a hotspot, whereas peak broadening and a
flux gradient across the maximum tend to indicate heat advection along
latitude lines from global winds.

In our fits, we choose the gravity darkening coefficient to be $\beta_{\rm
D}=$0.08 for companions with low front-side temperature ($T\lesssim 10,000$ K)
and 0.25 for front-side temperature $T\gtrsim 10,000$ K. To check whether
radiative conditions apply, for imtermediate front-side temperature
$5,000\lesssim T\lesssim 9,000$ K, we tried fitting the gravity darkening
coefficient and found that $\beta_{\rm D}$ of 0.08 was strongly preferred over
0.25; for such cases, $\beta_{\rm D}=0.08$ was fixed during model fitting.

The observed lightcurves are also affected by the intervening extinction
$A_{V}$. For many of our objects, we adopt the extinction estimates from the
three-dimensional dust model of Green et al. (2019) available as a query at
argonaut.skymaps.info, using the Bayestar2019 model and $A_{V}=2.73\,E(g-r)$
values at the estimated distance. For southern sources, we take the total
extinction column through the Galaxy from the NASA Extragalactic Database
(Schlafly & Finkbeiner, 2011) as an upper limit.

For the light curve modeling, the principal geometrical and physical fit
parameters are the inclination of the binary system $i$, the Roche lobe
filling factor $f$ (defined as the ratio of the companion nose radius to the
radius of $L_{1}$ point), the base temperature of the star $T_{N}$, the pulsar
irradiation power $L_{H}$, and the distance to the binary system $d$ in kpc.
For HS models, we also fit the Gaussian spot amplitude $\mathcal{A}_{\rm hs}$
(the multiplicative increase to the local temperature), radial size $r_{\rm
hs}$, and angular position $\theta_{\rm hs}$ (measured from the sub-pulsar
point on the companion), $\phi_{\rm hs}$ [measured from the orbital plane, 0
toward the dusk (East) terminator of the companion and $\pi/2$ toward the
North (spin) pole]. For WH models, we fit for the wind parameter
$\epsilon=\tau_{rad}/\tau_{\rm adv}$, the ratio of the radiation cooling and
advection timescales. For $i$, a uniform prior in $\cos\,i$ was applied; for
$L_{H}$, a uniform logarithmic prior was applied, whereas for all other
parameters a uniform prior was used. For all our fits, we explore models with
increasing model complexity: a direct heating (DH) model, DH+wind heating
(WH), DH+ hot spot(HS), and DH+WH+HS. To penalize increased model complexity,
the final model selection is based on Akaike Information Criterion.

Photometric fitting gives us constraints on heating parameters as well as
binary geometry. This can be combined with center-of-mass velocity obtained
from radial velocity (RV) fitting to estimate the NS and companion masses. In
general, the observed RV amplitude is different from the center-of-mass (CoM)
radial velocity amplitude $K_{\rm CoM}$ because strong pulsar irradiation
shifts the photometric center of the secondary star (“center of light” CoL)
away from its CoM. Moreover, different absorption lines have strengths varying
with the surface temperature; thus the CoL shift is different for each
spectral line. The temperature sensitivity of the line profile can be
characterized by monitoring its equivalent-width $EW(T)$ variation across the
companion surface. Therefore, we model the CoL RVs for a given $K_{\rm CoM}$,
by averaging up the radial velocities over the companion surface elements,
weighted by flux and temperature-dependent equivalent-width $EW(T)$ of the
lines that dominate the radial velocity template. We fit these predicted CoL
RVs to the cross-correlation velocities, estimating $K_{\rm CoM}$ and hence
the binary masses. Although we do not perform a simultaneous photometry fit,
we do marginalize the spectroscopic fit over the geometrical parameters from
the end of the photometric Markov Chain Monte Carlo chains, sampling $\sim
2\sigma$ uncertainties. Thus, the mass errors do include all uncertainties in
the model fitting, spectroscopic and photometric.

Below, we describe photometric and RV fitting of individual BW objects. Our
population study includes J0952-0607 and J1810+1744, but since these objects
were already fitted with the most updated prescription in Romani et al. (2022)
and Romani et al. (2021), we skip the modeling details and present only the
fit results. The values for J2051$-$0827 are from Dhillon et al. (2022).

### PSR J0023+0923

Lacking radial velocity measurements, we do not have enough constraints for a
detailed NS mass estimate, so we fix at $M_{\rm NS}=1.5\,M_{\odot}$. We find
that a HS model fits the data the best, with $\chi^{2}/$DoF $=107/57$ and a
poorly constrained $i(^{\circ})=74.5^{+9.7}_{-15.5}$. With this inclination,
the companion mass is rather small $M_{c}\sim 0.018M_{\odot}\,(M_{\rm
NS}/1.5\,M_{\odot})^{2/3}$, typical of other BWs. The maximum is relatively
well sampled, thus the fill-factor is modestly constrained at $0.87\pm 0.05$,
implying a volume-averaged companion radius of $0.10\,R_{\odot}$. As noted in
Draghis & Romani (2018), BW companions appear to be inflated by the heating
and hence we expect companion radius to be larger than the radius $R_{\rm
cold}$ of a cold (degenerate) remnant of a stellar core. This heating-induced
inflation is discussed in §4.2. Due to limited photometry near optical
minimum, our inclination and hot-spot position estimates have substantial
uncertainty. While the hotspot’s nose angle $\theta_{\rm hs}$ is near the
nose, its phase angle is highly unconstrained. Additional photometry covering
the minimum, especially in the $i$ and $z$ bands, would help constrain $T_{N}$
and $i$, and $u$ band data would help constrain hot-spot parameters.

Figure 1: Four BW LCs. $\phi_{\rm B}=0$ is defined as pulsar TASC (time of the
pulsar ascending node). Two cycles are shown. For J0952 and J1301 direct
heating provides the best-fit model and is shown in both cycles. For the
others the first cycle shows the best-fit direct heating model, and the second
is the best-fit hotspot or wind model. Lower panels show residuals to the
model (for a red band and a blue band) of the corresponding cycle – the model
differences are often subtle, best seen in the residuals, but are
statistically significant (see text). Colors denote the various photometric
bands (blue=$u$, green=$g$ red=$r$, magenta=$i$, black=$z$, yellow=$H$,
cyan=$K$.

Figure 2: As for figure 1, except for J1959+2048, where blue=$B$, green=$V$,
red=$R$, magenta=$I$.

### PSR J0636+5128

Fig 1 shows the $grizHK$ LC of this object. The LC is relatively shallow for
all colors, indicating a low binary inclination. Since radial velocity
measurements are unavailable, we cannot fit the NS mass, so we fix at $M_{\rm
NS}=1.5\,M_{\odot}$. The best model is achieved with a $11.0^{\circ}$ Gaussian
hotspot located close to the pole of the ‘northern’ hemisphere, with a
substantial temperature excess of 64%. The direct heating power is $1.9\times
10^{33}\,$erg/s representing about $36\%$ of the spin-down power. With the HS
model, the inclination is $i\approx 34.0^{\circ}$, slightly larger than the DH
estimate in Draghis et al. (2019). With this inclination, the companion mass
is rather small $M_{c}\sim 0.01M_{\odot}(M_{\rm NS}/1.5\,M_{\odot})^{2/3}$, so
for a modest $M_{\rm NS}>1.5\,M_{\odot}$, J0636’s has a BW-type companion. The
timing parallax (Stovall et al., 2014) of J0636 gave a very small distance,
but this has been superseded by more extended timing, giving
$d=1.1^{+0.6}_{-0.3}$ kpc (Arzoumanian et al., 2018), and our photometric
distance fit of $0.63^{+1.44}_{-0.43}$ is in good agreement with this value
although the constraint is rather loose. At this distance, the fill factor of
$0.88$ gives a companion radius of $0.095\,R_{\odot}$.

### PSR J1301+0833

Romani et al. (2016) found that for J1301, the combination of faint magnitude
and $\sim 4500$ K effective temperature companion requires a distance of $\sim
6$ kpc for direct full-surface heating of a Roche-lobe filling star. This is
dramatically larger than the $d=1.23$ kpc implied by the pulsar’s
$DM=13.2\,{\rm pc\,cm}^{-3}$ (Ray et al., 2012) in the YMW17 model. We have
fit with DH, HS and WH models and find that all prefer a companion that
substantially under-fills its Roche lobe at $f\sim 0.61$. This is the smallest
fill factor found for any of our BWs, but with this value the fit distance is
a plausible 1.77 kpc. While the HS and WH models do slightly decrease
$\chi^{2}$ (by 10 and 4, respectively), with the extra parameters AIC prefers
the basic DH model. All models give a similar $i\approx 46.6^{\circ}$. This
small inclination is quite consistent with the relatively small RV amplitude
of Romani et al. (2016). At our fit distance the Shklovskii-corrected spin-
down power is $\dot{E}_{c}=\dot{E}I_{45}\left[1-(0.31\pm 0.07)d_{\rm
kpc}\right]=3.0\times 10^{34}\,$erg/s for the moment of inertia
$10^{45}I_{45}$ g/cm2. Thus, our best-fit $L_{H}=4.7\times 10^{33}$ erg/s is a
modest $16\%$ of ${\dot{E}}$.

### PSR J1311-3430

J1311 is particularly interesting since Romani et al. (2012) found that the NS
might be very massive. Romani et al. (2015) used early versions of several HS
models to show that the fit $i$ could range between $64^{\circ}-85^{\circ}$,
with the NS mass as low as $1.8\,M_{\odot}$ and as large as $2.7\,M_{\odot}$.
With our improved modeling we can reduce these ranges.

Figure 2 shows the LCs of this object from simultaneous photometry synthesized
from Keck LRIS spectra. Simultaneous light curves are essential for this
object, since J1311 is known to have intense optical flaring, with detectable
flares occurring nearly every orbit and large flares appearing every dozen
orbits. With simultaneous colors and multiple orbits we can prune these flares
leaving the underlying quiescent thermal flux to be described by the heating
model. J1311’s LCs show three main features: i) the maximum is particularly
flat, indicating strong gravity darkening, ii) the maximum is asymmetrical,
especially in bluer colors, iii) a slight gradient is present across maximum.
All fits find that the companion is very close to Roche-lobe filling, so we
set the fill factor at $f=0.99$. The large fill factor (low surface gravity
near the $L_{1}$ point) and the high average front-side temperature $>10,000$
K should require a large $\beta_{\rm D}$. Indeed, if freed in the fits,
$\beta_{\rm D}=0.30\pm 0.06$, consistent with the $\beta_{\rm D}=0.25$
expected for a radiative atmosphere.

In our fitting, the best model invokes both HS + WH, with a sub- Alfvénic wind
with $\epsilon=-0.07\pm 0.01$, and a hotspot just past the day-night
terminator with a large $\approx 12\times$ excess over the local nightside
temperature. The resulting binary inclination is quite well constrained at
$68.7^{\circ}\pm 2.1^{\circ}$. We also find that the base temperature $T_{N}$
is very low, with samples converging below the lowest temperature limit of
1000K of our spectral library. Such extrapolation to low temperature will
likely be imprecise as unmodelled effects such as dust settling and clouds are
important for low-temperature atmosphere spectra (Husser et al., 2013). Thus,
we consider $T_{N}$ low, but poorly constrained. Note that the inclination is
hardly affected by this uncertainty as the model light curves are dominated by
hot surface elements resulting from extreme pulsar irradiation.

Using the DM$=37.84$ pc cm-3 (Ray et al., 2013), Antoniadis (2021) estimated
$2.43\pm 0.48$ kpc. Our distance estimate of $3.01\pm 0.15$ kpc is consistent
with this DM estimate. The best-fit $L_{H}=2.28\times 10^{35}$ erg/s is
$\approx 4$ larger than the nominal spin-down power. For J1311, a substantial
portion of the flux seems to be from IBS particle precipitation (a strong hot
spot) and some leakage from this particle-mediate heating may also contribute
to the large $L_{H}$.

### PSR J1653-0158

The lightcurve shows shallow modulation implying either a low inclination or a
strong veiling flux. However, with the companion RV amplitude $K_{CoL}=669\pm
7.5\,{\rm km\,s^{-1}}$ we find a large mass function $f(M)=1.60\pm
0.05M_{\odot}$ (Romani et al., 2014), so small inclination would lead to an
unphysically large NS mass. Noting that the minimum is flat and that the
modulation decreases for bluer colors, we infer that a strong blue veiling
flux dominates at orbital minimum. This is also visible in the phase-resolved
spectra (Romani et al., 2014). Although this veiling flux is likely associated
with the IBS, we model it here as a simple power-law with form
$f_{\nu}=f_{A}(\nu/10^{14}{\rm Hz})^{-p}$, with $f_{A}\sim 101\pm 10\,\mu$Jy
and $p=0.50\pm 0.03$. This flux is fairly constant through the orbit, although
there are hints of sharp phase structures in the light curve, e.g. in $r$ and
$i$ at $\phi_{\rm B}=0.72$. Any model without such a hard-spectrum component
is completely unacceptable. These fits prefer an $A_{V}\sim 1.0$ mag slightly
higher than, but consistent with the maximum in this direction (which is found
for all $d>300$pc).

We find that the HS model best fits the data. However, with the fine structure
near maximum, the model is not yet fully acceptable ($\chi^{2}/DoF\sim 1.38$).
More detailed models, including modulated veiling flux from the IBS, may be
needed to fully model the light curves. Such modeling would be greatly helped
by light curves over an even broader spectral range, with IBS effects
increasingly dominant in the UV and low-temperature companion emission better
constrained in the IR. With many cycles, we could also gauge the reality (and
stability) of the apparent fine structure and test for hot spot motion.

### PSR J1810+1744

Our fit follows the assumptions of Romani et al. (2021), except that we allow
for possible errors in the photometric zero points. This results in small
changes in the fit parameters and a substantial increase in the distance and
$L_{H}$ uncertainties. The neutron star mass is decreased from that paper’s
value by $\sim 0.5\sigma$.

### PSR J1959+2048

This, the original BW pulsar, is of special interest, since van Kerkwijk et
al. (2011), with an approximate treatment of DH effects, infer that it may
have a pulsar mass as large as $2.4\pm 0.12M_{\odot}$. Recently, evidence of
an eclipse of the pulsed $\gamma-$ray emission from J1959 has been presented
in (Clark et al., 2021). With a small $\sim 0.1R_{\odot}$ companion, this
requires a binary inclination $i>83^{\circ}$, far from the result of optical
LC modeling in van Kerkwijk et al. (2011). Some support for this edge-on view
also comes from evidence for an X-ray eclipse in Kandel et al. (2021).

With our improved treatment of gravity darkening, a simple DH model gives
$i=65.7^{\circ}\pm 2.4^{\circ}$ and $\chi^{2}/$DoF = $138/89$. However,
looking at the light curves, one can see a small phase shift in the maximum,
with the peak somewhat broadened and excess flux at phases $\phi\sim
0.75–0.9$, especially in the bluer bands. In Kandel & Romani (2020), we showed
how this excess can be explained by both a Gaussian hotspot and a banded wind.
Because the LC maxima are relatively flat, the best-fit $i$ is $\sim
64^{\circ}$. Strong gravity darkening might also produce such a flat maximum,
but the low front-side temperature of J1959 ($\lesssim 5000\,$K) together with
a relatively low fill-factor of $\sim 0.8$ precludes this possibility.

For our analysis, we impose a prior on $i$ of $(83^{\circ},90^{\circ})$, where
the bounds are the result of geometrical constraints imposed by the
observation of $\gamma-$ray eclipse. At such a high inclination, the
lightcurve becomes relatively narrow at any reasonable value of the fill
factor. Thus, to match the observational data, we require a spread of heat
away from the maximum. While a banded wind can achieve this, a model with heat
diffusion fits the data the best, resulting in $i=85.2^{+2.9}_{-1.2}$, with
$\chi^{2}/$DoF=$128/88$. The angular scale of diffusion is $17.1^{\circ}\pm
4.2^{\circ}$. While this model fits data reasonably well, it is statistically
worse than a model with an inclination as a free parameter. Fig 2 shows the
light curves for the best-fit model with inclination bound by $\gamma$-ray
eclipse (first cycle) and one with free inclination (second cycle). $V$ and
$R$ fit residuals for the two models are shown in the lower panels. This
substantial tension between the gamma-ray eclipse and the optical light curve
modeling is worrisome, and shows that this binary should be re-measured with
modern multi-band photometry, to see if systematics affect the limited
existing data or if additional physical effects are required to obtain a good
fit.

### PSR J2051$-$0827

Recently Dhillon et al. (2022) have described simultaneous multi-color
observations and ICARUS fitting of this BW. We do not attempt to re-fit here,
but for comparison report the parameters of that study. They find significance
evidence for a hot spot before 2011 (parameters not given), but that this spot
is no longer prominent in 2021.

### PSR J2052+1219

All model fits require a fill factor $\sim 1$, so we fix $f=0.99$. Lacking RV
data, we fix $M_{\rm NS}=1.5\,M_{\odot}$. A DH model gives $\chi^{2}/$DoF of
$165/81$, with best-fit $i=59.8^{\circ}\pm 2.2^{\circ}$. Adding a hotspot
reduces the $\chi^{2}$/DoF to $145/76$, with $i=60.6^{\circ}\pm 2.0^{\circ}$.
However, the hotspot is rather large with a radius $\sim 39^{\circ}$ and a
temperature-excess of $\sim 200\%$ of the base temperature. A WH model gives
$\chi^{2}$/DoF to $153/80$, with a very similar inclination estimate of
$i=59.3^{\circ}\pm 2.0^{\circ}$. After penalizing for the extra DoF in the HS
model, the WH model is preferred at a marginal level of 60% over the HS model.
Physically, very shallow gradients (large effective spot radius) seem more
natural for wind flow than a magnetic pole, so we prefer WH on these grounds.
With very similar $i$, either model appears to capture the overall heating
pattern. With good overall light curve matches (Fig 2), we infer that the
large $\chi^{2}/$ DoF is the result of low-level stochastic flaring or
possibly under-estimation of photometry errors. With this inclination, the
companion mass is $M_{c}\sim 0.041(M_{\rm NS}/1.5\,M_{\odot})^{2/3}$, typical
for a BW.

Table 3: Mass estimates of some BWs and RBs for which RV measurements are available. Name | $i\,(\mathrm{deg})$ | $K_{\rm CoM}\,(\mathrm{km/s})$ | $M_{\rm PSR}\,(M_{\odot})$ | $M_{\rm C}\,(M_{\odot})$
---|---|---|---|---
J$0952-0607$ | $59.8^{+2.0}_{-1.9}$ | $376.1\pm 5.1$ | $2.35\pm 0.17$ | $0.032\pm 0.002$
J$1301+0833$ | $46.6^{+2.5}_{-2.2}$ | $274.8\pm 8.1$ | $1.60_{-0.25}^{+0.22}$ | $0.036\pm 0.006$
J$1311-3430$ | $68.7^{+2.1}_{-2.0}$ | $641.2\pm 3.6$ | $2.22\pm 0.10$ | $0.012\pm 0.0006$
J$1653-0158$ | $72.8^{+4.0}_{-4.0}$ | $700.2\pm 8.0$ | $2.15^{+0.16}_{-0.16}$ | $0.014\pm 0.001$
J$1810+1744$ | $66.3^{+0.5}_{-0.5}$ | $462.9\pm 2.2$ | $2.11\pm 0.04$ | $0.064\pm 0.001$
J$1959+2048$ | $85.3^{+2.3}_{-1.2}$ | $334.6\pm 3.6$ | $1.55_{-0.05}^{+0.06}$ | $0.024\pm 0.001$

##  4 Discussion and Summary

###  4.1 Thermal emission from the companion

We find that a basic direct (photon) heating model is insufficient to
represent the light curves of many BWs and RBs. There is statistically
significant improvement in most model fits when adding a localized hot spot
and/or global winds, and improved treatments of gravity and limb darkening. We
compute multiple models for each object, including HS, WH, and diffusion,
before deciding on the best, using the AIC. In the cases where the AIC did not
select a preferred model and the competing models fit the light curve well,
both give very similar binary kinematic parameters. In other words a well-
matched light curve gives a reliable mass, even when the fit model is not
unique.

However, variability can be an important factor in this modeling. The hot spot
phases and brightness can vary, especially for RB, and may even show secular
trends (e.g van Staden & Antoniadis, 2016). When correctly modeled, the
underlying heating pattern, and the fit binary kinematic parameters should be
stable. Other variability issues arise from the optical/X-ray flaring activity
seen in strongly heated spiders. Such flares are best isolated in simultaneous
multi-color photometric observations. It is important to excise these events
before forming the ‘quiescent’ light curve needed for the steady heating model
and fit of the orbital parameters. An important, but often subtle effect is
the presence of a non-thermal veiling flux. If associated with the IBS this
may be phase dependent and may also exhibit secular variations. When bright
enough to dominate at binary minimum (as for J1653), this too can be an
important complication in fitting the quiescent heating pattern.

We find that $\sim$ half of the BW sources in our study prefer models with a
hotspot. In Table 5 we separate the companions’ integrated thermal surface
emission into day side $L_{\rm comp}$ and hotspot $L_{\rm hs}$ components. The
hotpots often appear far from the direct heating maximum at the sub-pulsar
point, and so can represent substantial changes in the local surface flux (and
corresponding light curve features) with only modest energy input. The largest
$L_{\rm hs}/L_{\rm comp}$ ratio is $\sim 0.25$ for J1653. Such hotspots are
even more common in the RBs where the heavier ($\geq 0.08\,M_{\odot}$)
companions should have core fusion and strongly convective envelopes. Since
tidal locking ensures rapid rotation, the RBs should support strong
$\alpha-\Omega$ dynamos resulting in large (but constantly refreshing) dipole
magnetic fields. Such fields can direct IBS-generated particles to the
magnetic poles (Sanchez & Romani, 2017), resulting in substantial hotspots,
whose positions may vary with epoch. Although the BW companions lack core
fusion, the very large front-back temperature gradients may well drive
internal convection, resulting in similar dynamo-supported global fields.
Indeed the objects preferring hotspots tend to have stronger heating. A recent
two-epoch photometric study of BW J2051$-$0827 (Dhillon et al., 2022) suggests
that hot spot variability may occur in BWs, as well.

We conclude that direct heating dominates energetically. Since pulsar SEDs are
dominated by the GeV $\gamma$-rays we naturally assume that $\gamma$-rays will
typically dominate the companion heating. But it is important to recall that
the Fermi-observed $\gamma$-ray flux is from a single slice through the
$\gamma$-ray beam along the Earth line-of-sight (for pulsars with aligned spin
and orbital angular momenta this will be at co-latitude $i$). This cut may not
represent the $\gamma$-ray flux intercepted near the orbital plane by the
companion; for spin aligned pulsars this is the spin equator. Most modern
outer magnetosphere models have the $\gamma$-ray flux concentrated to this
equatorial plane, so in general we expect the companion-intercepted heating
will be larger than that seen by Fermi, although the correction factors are
highly model dependent (Draghis et al., 2019).

For most MSP, the observed photon (i.e. $\gamma$-ray) flux represents only a
modest fraction of the spin-down power; the remainder is assumed to be carried
off as the $e^{\pm}/B$ pulsar wind. For the BWs in our sample if we (naively,
incorrectly) assume that the $\gamma$-ray flux is isotropic, we find that it
represents $<0.25I\Omega{\dot{\Omega}}={\dot{E}}/4$. Table 5 lists this
isotropic $\gamma$-ray efficiency $f_{iso}$ for a standard moment of inertia
$I=10^{45}{\rm g\,cm^{2}}$ (i.e. $I_{45}=1$). The very low $\gamma$-ray flux
of J0636 indicates that our small $i$ view is outside of the main $\gamma$-ray
beam. Three sources have a large isotropic $\gamma$-ray efficiency: For J2052,
our photometric $d$ is substantially larger than the Table 4 $DM$ distance
values and may be an overestimate; at the 2.4 kpc DM distance it would have a
more typical $f_{iso}=0.13$. From our fitting below, J1810 has a large pulsar
mass (Table 3), so that we expect $I_{45}=2-3$, reducing $f_{iso}$ to
0.19-0.28. For J1311, the very large $f_{iso}=1.4$ may well be produced by a
combination of the two factors; the inferred NS mass is large and $d$ exceeds
the $DM$ estimates. Here an independent (eg. optical) parallax distance would
be quite valuable. Of course all these large values should be mitigated by
$\gamma$-ray beaming preferentially to the Earth line-of-sight.

Since in modern models the pulsar wind power (as well as the $\gamma$-ray
pulsed emission) is concentrated to the spin equator, we approximate this in
our direct heating models by distributing the heating power as $\propto{\rm
sin}^{2}\theta$, with $\theta$ measured from the spin axis. With this
distribution we can write the model-fit heating flux’s value for the Earth
line-of-sight $f_{H}$. It now becomes interesting to examine
$\eta_{\gamma}=f_{H}/f_{\gamma}$. If we assume that the direct heating is
$\gamma$-ray pulsed flux $\eta_{\gamma}>1$ tells us that there is more
$\gamma$-emission heating the companion near the spin equator than at the
Earth sightline, and vice-versa. Of course, we can also compare the integral
heating flux required with ${\dot{E}}$, via $\eta_{\dot{E}}=L_{H}/{\dot{E}}$.
At $\eta_{\dot{E}}=1$, ${\rm sin}^{2}\theta$ heating would require the full
$I_{45}=1$ spindown power. For $\eta_{\dot{E}}>1$ we infer some combination of
spindown power more equatorially focused, large $I_{45}$ and smaller $d$.
J1311 certainly requires some or all of these effects to reconcile the
observed heating with the pulsar energy loss rate.

Table 4: BW distances and Shklovskii-corrected magnetic fields Name |  $d$ (kpc) | CL02(kpc) | Y17(kpc) | Parallax(kpc) |  $B_{\rm int}(10^{8}$G)
---|---|---|---|---|---
J$0023+0923$ | $0.59^{+1.39}_{-0.39}$ | 0.69 | 1.25 | $-$ | $1.82^{+0.05}_{-0.17}$
J$0636+5128$ | $0.63^{+1.44}_{-0.43}$ | 0.49 | 0.21 |  $1.1^{+0.6}_{-0.3}$aaTiming parallax Arzoumanian et al. (2018); bGAIA DR3; cVLBI parallax Romani et al. (2022b) | $1.0^{+0.01}_{-0.06}$
J$0952-0607$ | $6.26^{+0.36}_{-0.40}$ | 0.97 | 1.74 | $-$ | $0.61\pm 0.02$
J$1301+0833$ | $1.77^{+0.10}_{-0.11}$ | 0.67 | 1.23 | $-$ | $0.96\pm 0.03$
J$1311-3430$ | $3.01\pm 0.15$ | 1.4 | 2.43 | $-$ | $2.29\pm 0.03$
J$1653-0158$ | $0.87^{+0.07}_{-0.08}$ |  |  |  $0.57^{+0.44}_{-0.18}$bbBolometric luminosity of the hotspot and the companion star without hotspot. Uncertainties generated through MCMC chain of the photometric fitting. | $0.36\pm 0.04$
J$1810+1744$ | $2.86\pm 0.13$ | 2.00 | 2.36 |  $1.54^{+7.61}_{-0.70}$bbBolometric luminosity of the hotspot and the companion star without hotspot. Uncertainties generated through MCMC chain of the photometric fitting. | $0.77\pm 0.01$
J$1959+2048$ | $2.25\pm 0.11$ | 2.49 | 1.73 |  $2.57^{+1.84}_{-0.77}$ccfootnotemark:  | $1.19\pm 0.01$
J$2052+1219$ | $5.57^{+0.47}_{-0.55}$ | 2.4 | 3.92 | $-$ | $0.26\pm 0.22$
Table 5: BW Derived and Observed Heating Name | $f_{\gamma}$ | $\eta_{iso}$ |  $f_{\rm H}$aa$f_{H}=3{\rm sin^{2}}i\,L_{h}/(8\pi d^{2})$, i.e. assuming $\sin^{2}\theta$ $\gamma$-ray beaming, orbital and spin momenta aligned. | $\eta_{\dot{E}}$ | $\eta_{\gamma}$ |  $L_{\rm comp}$bbBolometric luminosity of the hotspot and the companion star without hotspot. Uncertainties generated through MCMC chain of the photometric fitting. |  $L_{\rm hs}$bbBolometric luminosity of the hotspot and the companion star without hotspot. Uncertainties generated through MCMC chain of the photometric fitting.
---|---|---|---|---|---|---|---
|  $10^{-12}$erg/s/cm2 |  |  $10^{-12}$erg/s/cm2 |  |  |  $10^{30}$erg/s |  $10^{30}$erg/s
J$0023+0923$ | $7.3\pm 1.4$ | 0.019 | $40.6^{+430}_{-38.3}$ | 0.076 | 5.6 | $8.0\pm 1.0$ | $1.5^{+0.7}_{-0.5}$
J$0636+5128$ | $<0.5$ | $<0.004$ | $18.9^{+230}_{-17.7}$ | 0.33 |  $>$38 | $4.1^{+0.8}_{-0.7}$ | $0.1^{+0.04}_{-0.02}$
J$0952-0607$ | $2.2\pm 0.5$ | 0.154 | $0.91^{+0.23}_{-0.22}$ | 0.057 | 0.41 | $108^{+10}_{-11}$ | $-$
J$1301+0833$ | $10.6\pm 1.5$ | 0.060 | $10.8^{+3.9}_{-2.84}$ | 0.071 | 1.0 | $16.6^{+1.7}_{-1.8}$ | $-$
J$1311-3430$ | $64.7\pm 1.9$ | 1.40 | $273^{+106}_{-85.6}$ | 4.56 | 4.2 | $321^{+80}_{-62}$ | $16^{+19}_{-7}$
J$1653-0158$ | $33.7\pm 1.8$ | 0.26 | $65.0^{+26.3}_{-18.1}$ | 0.11 | 1.9 | $7.3\pm 0.9$ | $1.9^{+0.9}_{-0.6}$
J$1810+1744$ | $22.4\pm 1.3$ | 0.56 | $47.8^{+5.9}_{-5.6}$ | 1.05 | 2.1 | $187^{+17}_{-19}$ | $19\pm 2$
J$1959+2048$ | $17.9\pm 1.9$ | 0.068 | $105\pm 9.5$ | 0.27 | 5.9 | $110^{+7}_{-8}$ | $-$
J$2052+1219$ | $6.5\pm 1.0$ | 0.720 | $3.5^{+2.3}_{-1.5}$ | 0.35 | 0.54 | $72.1^{+13.8}_{-12.0}$ | $-$

###  4.2 Mass distributions

Figure 3: The companion mass versus NS mass for all BWs (black points) and RBs
(red points) for which a reliable mass measurement exists. Masses for objects
with circled points come from this paper; cross points are RBs from Strader et
al. (2019) and BW J1555 from Kennedy et al. (2022).

Figure 3 shows the distribution of NS and companion masses for both BWs and
RBs, with most reliable measurements of BWs coming from this paper. The NS
masses of RBs have a wide spread and are on average lower than those of BWs.
For companion masses, the BW/RB separation remains clear, with BWs having a
low companion mass $\sim 0.01-0.06\,M_{\odot}$, and RBs having a higher
companion mass $\sim 0.2-0.5\,M_{\odot}$. However, the detection of relatively
heavy BW companions in J1810 and J1555 suggests that the mass gap between
these two populations is a factor $\sim 3\times$.

In Chen et al. (2013), it was argued that RBs and BWs are distinct
evolutionary populations, with irradiation efficiency determining the
formation channel; in this picture RBs do not evolve into BWs. On the other
hand, numerical simulations in Benvenuto et al. (2014) suggests that $P_{B}<1$
d RBs will indeed evolve into BWs. Our sample does not distinguish these two
scenarios, although a decrease in the BW-RB companion mass gap could be seen
as supporting an evolutionary connection.

Table 6: BW Companion Mean Densities Name |  $P_{b}$(h) |  $\rho$ (g/cm3) | $R/R_{\rm Y_{e}=Y_{\odot}}$
---|---|---|---
J$0023+0923$ | 3.33 | 22.8$\pm 0.8$ | 0.92
J$0636+5128$ | 1.60 | 21.3$\pm 0.3$ | 1.77aaUCXRB progenitor; cold radius for $Y_{e}=0.5$.
J$0952-0607$ | 6.42 | 40.1$\pm 2.5$ | 1.12
J$1301+0833$ | 6.48 | 36.1$\pm 7.5$ | 1.25
J$1311-3430$ | 1.56 | 28.6$\pm 1.4$ | 1.52aaUCXRB progenitor; cold radius for $Y_{e}=0.5$.
J$1653-0158$ | 1.25 | 33.3$\pm 2.5$ | 1.61aaUCXRB progenitor; cold radius for $Y_{e}=0.5$.
J$1810+1744$ | 3.56 | 30.8$\pm 0.5$ | 1.94
J$1959+2048$ | 9.17 | 28.6$\pm 1.9$ | 1.04
J$2051-0827$ | 2.37 | 20.2$\pm 0.6$ | 1.60bbDhillon et al. (2022)
J$2052+1219$ | 2.75 | 33.3 | 1.64

In Table 6, we list the inferred mean density of our companions. With our fit
parameters, these are surprisingly similar, varying by less than $2\times$. We
also list the ratio of their volume equivalent radius to that of a cold
degenerate object of the companion mass. The $T=0$ radius for a $\Gamma=5/3$
degenerate object with $Y_{e}$ electrons per nucleon is $R_{\rm
cold}=0.0126\left(\frac{M_{C}}{M_{\odot}}\right)^{-1/3}(2Y_{e})^{5/3}R_{\odot}$.
In most BW evolutionary scenarios, mass transfer and subsequent companion
evaporation are initiated well before core exhaustion. In this case we might
expect a $\sim$ solar metalicity for the companion with $Y_{e}\approx 0.833$.
We see that most companions are inflated above this value (up to nearly
$2\times$ for the strongly heated J1810). A few objects have small
$R/R_{Y_{\odot}}$, (e.g. J0023); for these we might assume some core
enrichment before evaporation commences, giving smaller cold $R$ and stronger
inflation. On the other hand the shortest period BWs are believed to be the
descendants of ultra-compact X-ray binaries, with hydrogen depleted
secondaries, so for the three objects with $P_{b}<2$ h we can assume
$Y_{e}=0.5$. This is supported by spectroscopy of J1311 and J1653, where the
absence of H features shows that even the photosphere is H depleted (Romani et
al., 2014).

Figure 4: Shklovskii-corrected surface dipole magnetic field vs. accreted mass
for some known BWs (blue) and RBs (red), assuming photometric distances except
when more accurate parallax distances are available. Note that for three PSRs
(J1311, J2039 and J2215) we assume that they were born as higher mass neutron
stars. While the RBs show no trend the BWs seem to have lower magnetic fields
at higher accreted (and final) mass.

###  4.3 Relationship between NS mass and B field

Binary pulsars have lower inferred dipole fields than isolated neutron stars,
and it has long been proposed (e.g., Bisnovatyi-Kogan, 1974; Romani, 1990)
that past accretion has decreased young pulsar magnetic fields to millisecond-
pulsar values. It is unclear whether the amount of mass accreted (as in simple
burial models) or the duration of the accretion phase (as in heating-driven
ohmic decay) is most important (see Mukherjee, 2017, for a review). Note that
small MSP $P_{s}$ are not, of themselves, precise probes of accretion history.
Since the angular momentum needed for spin-up to breakup periods is only $\sim
0.1\,M_{\odot}$ and since birth masses seem to vary by at least this amount it
is unlikely that even high precision mass measurements can determine the mass
gain with sufficient accuracy to probe angular moment addition (although
statistical studies might prove useful). Instead $P_{s}$ tracks the neutron
star dipole field in equilibrium spin-up models, thus telling us something
about the accretion rate at the end of spin-up, and subsequent spindown. The
best hope for a probe of dipole field reduction occurs if large amounts of
mass are needed to decrease the field from initial Terra-Gauss values to $\sim
10^{8}$G levels required for small $P_{s}$. Even then the situation is
complicated, since MSP population studies suggest that the initial masses
$M_{i}$ of these pulsars are bimodal with a dominant component of $M_{i}\leq
1.4M_{\odot}$ and a 20% sub-population with $M_{i}\sim 1.8M_{\odot}$
(Antoniadis et al., 2016, and references therein).

We make a first attempt at such comparison in Figure 4 where the Shklovskii-
corrected intrinsic dipole surface field $B_{\rm int}$ for an assumed moment
of inertia $I_{45}=1$ is plotted against inferred mass increase. Table 4
$B_{\rm int}$ values are computed assuming the nominal photometric distance
uncertainty; for the plot we add an 10% additional uncertainty in quadrature
to acknowledge $I_{45}$ variation with the uncertain mass. Three BWs that have
particularly high accreted mass (J0952 $\Delta M\approx 1.0M_{\odot}$,
J1653$-$0158 $\Delta M\approx 0.8M_{\odot}$, and J1810+1740 $\Delta M\approx
0.76M_{\odot}$, assuming a start at typical $1.35\,M_{\odot}$) also have some
of the lowest $B_{\rm int}$ known. However, J1723 has a relatively low surface
field and has apparently accreted very little mass (although its mass
measurement has large error bars, and this object has not been subject to the
uniform fits of our study). Conversely the well-measured J1311, J2039 and
J2215 are heavy neutron stars. Their $B_{\rm int}$ fields also seem
substantial at $>1.5\times 10^{8}$ G, although with the possible photometric
distance errors above, the Shklovskii corrections remain quite uncertain. A
possible solution is that these systems may be from the 20% subset of MSP
starting the accretion phase from $\sim 1.8\,M_{\odot}$. This is assumed in
Figure 4. At this point we can only conclude that the largest $M$, smallest
$B_{\rm int}$ and shortest $P_{s}$ appear associated with substantial mass
increase – it is likely that other factors are important in determining
$B_{\rm int}$ and we will need more, and better ($B_{\rm int}$, $\Delta M$)
pairs to draw detailed conclusions.

Figure 5: Normalized probability distributions of mass estimates for MSPs with
$M_{\rm NS}>2.0\,M_{\odot}$. The dashed curves show the mass estimates for the
two most massive radio-selected pulsars with white dwarf (WD) companions,
measured from pulse timing (supplemented by WD atmosphere modeling for J0348).
The solid curves show spider binary mass estimates, relying on companion
spectrophotometry. The bottom panel shows the cumulative probability
distributions for $M<M_{\rm max}$, for the radio objects, the spiders, and all
six pulsars. Lower bounds on $M_{\rm max}$ ($1\sigma$, $2\sigma$, and
$3\sigma$) for these distributions are listed in the legend at top.

###  4.4 Mass implications on NS EoS

Shapiro delay measurements of radio pulsar PSR J0348+0432 give $2.01\pm
0.04\,M_{\odot}$ (Antoniadis et al., 2013) and PSR J0740+6620 give $2.08\pm
0.07\,M_{\odot}$ (Fonseca et al., 2021). However, the quoted errors are at the
$1\sigma$ level, and so these data do not require $M_{\rm max}>2\,M_{\odot}$
at high statistical confidence. Since the impact on EoS modeling is driven by
the minimum required $M_{\max}$, higher masses, with tighter estimates should
provide even greater physics impact. The spiders pulsars are good candidates
for such measurements and, indeed, J1810 is the first individual object for
which a $\sim 3\sigma$ lower bound on the mass exceeds $2\,M_{\odot}$ (Romani
et al., 2021). We have argued above that our improved, uniform treatment of BW
LC and spectroscopy data results in more stable and better-determined mass
measurements, relatively immune to modeling uncertainty. Our task now is to
combine these into a global constrain on $M_{\max}$. Figure 5 shows the mass
uncertainty ranges for the objects in this paper with central values
$>2M_{\odot}$.

Several approaches can be used to estimate $M_{\rm max}$. One option is to
model the full distribution of (binary) NS masses and see if an upper cutoff
is required; Alsing et al. (2018), for example, determine that $M_{\rm max}$
is within a $1\sigma$ range of 2.0–2.2 $M_{\odot}$. Here we only attempt to
determine a lower bound to $M_{\rm max}$. Assuming that our heavy NS sample is
drawn from a uniform population of lower bound $M_{1}\equiv 1.8\,M_{\odot}$
and upper bound $M_{\rm max}$, $U(M_{1},M_{\rm max})$, the marginal
distribution of individual observations $\\{M_{i},\sigma_{i}\\}$ is:

$\displaystyle p(M_{i})$ $\displaystyle=\int_{0}^{\infty}{\rm
d}\mu_{i}\,N(M_{i}|\mu_{i};\sigma_{i})\,U(\mu_{i}|M_{1},M_{\rm max})$ (1)
$\displaystyle=\frac{{\rm
Erf}\left[\frac{M_{i}-M_{1}}{\sqrt{2}\sigma_{i}}\right]-{\rm
Erf}\left[\frac{M_{i}-M_{\rm max}}{\sqrt{2}\sigma_{i}}\right]}{2(M_{\rm
max}-M_{1})}$ (2)

For $n$ observations, the log(likelihood) is

$\displaystyle\log\mathcal{L}$ $\displaystyle=-n\log(M_{\rm max}-M_{1})$
$\displaystyle+\sum_{i=1}^{n}\log\left({\rm
Erf}\left[\frac{M_{i}-M_{1}}{\sqrt{2}\sigma_{i}}\right]-{\rm
Erf}\left[\frac{M_{i}-M_{\rm max}}{\sqrt{2}\sigma_{i}}\right]\right)$ (3)

The log(likelihood) curves for different sample sets (Radio and Spiders) are
plotted in the lower panel of Fig. 5 and the median value, as well as the
$1-,2-$ and $3-\sigma$ lower bounds on $M_{\rm max}$ are listed for each
sample. At the $1\sigma$ level, $M_{\rm max}$ inferred from spiders is higher
than that inferred from radio measurements by $0.18\,M_{\odot}$, suggesting
that NSs in the spider systems are significantly heavier that those in NS-
white dwarf binaries. In our uniform re-fit of these spiders we have
marginalized over the photometric parameters. We have also shown that with
good fits to high quality LC data, we can select between various physically-
motivated heating models. When different models’ fit-quality is similar, we
find that the orbital parameters determining the fit mass are similar too.
Thus we have reduced systematic uncertainties, and it becomes worth discussing
statistical bounds at higher confidence levels. We can now say with high
($3\sigma$ one-sided lower bound) statistical confidence that $M_{\rm
max}>2.08\,M_{\odot}$, and that at $\sim 1\sigma$ significance $M_{\rm
max}>2.19\,M_{\odot}$ is preferred.

Figure 6: Mass-radius curves for a range of EoS. Figure adapted from Özel &
Freire (2016); see that paper for the EoS labels. Note that any viable EOS
must rise above the $M_{\rm max}$ lower bound. As this increases many
otherwise viable EoS are excluded at high confidence level.

For each EoS, integration of the Oppenheimer-Volkoff equations gives a
different $M-R$ curve. Curves for various literature EoS are shown in Figure
6. Although some of these are already deprecated by radius and tidal
deformation measurements of common $\sim 1.4\,M_{\odot}$ NS, if one is
interested in the physics of extreme conditions, the most important aspect of
these curves is probed by the roll-off at their maximum. This is because only
high mass objects achieve the core densities to probe such physics. ‘Soft’,
compressible EoS are already strongly excluded by the $\sim 2M_{\odot}$ radio
MSP. But the compressibility at $\rho>4\rho_{\rm sat}$ is best probed at
higher masses. In general EOS with non-nucleonic phases, such as hyperons or
condensates cannot reach these masses. Conversely, to reach masses far above
$2M_{\odot}$ a phase transition from nucleonic to even less compressible
matter might be required. In any event, the bounds from figure 5 are the
tightest constraints in this region. We have transferred the 1$\sigma$ and
$3\sigma$ limits from all spiders to Figure 6. Note that only five of the
listed models (two with implausibly large radii at $\sim 1.4M_{\odot}$)
survive the $1\sigma$ bound. Extending to $3\sigma$, only two more survive. It
is clear that astrophysics community’s efforts to pin down neutron stars in
the $M-R$ plane, substantially enhanced by our spider measurements, are
dramatically narrowing the options for the dense matter EoS.

We thank Alex Filippenko and colleagues for collaboration on many of the
observations re-analyzed in this paper. We also thank Rene Breton for useful
discussions on spider light curve fitting. DK and RWR were supported in part
by NASA grants 80NSSC17K0024390 and 80NSSC17K0502.

## References

  * Akaike (1974) Akaike, H. 1974, IEEE Transactions on Automatic Control, 19, 716
  * Aldcroft et al. (1992) Aldcroft, T. L., Romani, R. W., & Cordes, J. M. 1992, The Astrophysical Journal, 400, 638
  * Alsing et al. (2018) Alsing, J., Silva, H. O., & Berti, E. 2018, MNRAS, 478, 1377
  * Antoniadis (2021) Antoniadis, J. 2021, Monthly Notices of the Royal Astronomical Society, 501, 1116
  * Antoniadis et al. (2016) Antoniadis, J., Tauris, T. M., Ozel, F., et al. 2016, arXiv e-prints, arXiv:1605.01665.  https://arxiv.org/abs/1605.01665
  * Antoniadis et al. (2013) Antoniadis, J., Freire, P. C., Wex, N., et al. 2013, Science, 340
  * Arzoumanian et al. (2018) Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2018, The Astrophysical Journal Supplement Series, 235, 37
  * Benvenuto et al. (2014) Benvenuto, O. G., De Vito, M. A., & Horvath, J. E. 2014, The Astrophysical Journal Letters, 786, L7
  * Bisnovatyi-Kogan (1974) Bisnovatyi-Kogan, G. 1974, Soviet Astronomy, 18, 261
  * Breton et al. (2012) Breton, R. P., Rappaport, S. A., van Kerkwijk, M. H., & Carter, J. A. 2012, The Astrophysical Journal, 748, 115
  * Callanan et al. (1995) Callanan, P. J., Van Paradijs, J., & Rengelink, R. 1995, The Astrophysical Journal, 439, 928
  * Chen et al. (2013) Chen, H.-L., Chen, X., Tauris, T. M., & Han, Z. 2013, The Astrophysical Journal, 775, 27
  * Claret & Bloemen (2011) Claret, A., & Bloemen, S. 2011, Astronomy & Astrophysics, 529, A75
  * Clark et al. (2021) Clark, C., M., K., & Breton, R. P. 2021, in 9th International Fermi Symposium
  * Dhillon et al. (2022) Dhillon, V. S., Kennedy, M. R., Breton, R. P., et al. 2022, arXiv e-prints, arXiv:2208.09249.  https://arxiv.org/abs/2208.09249
  * Djorgovski & Evans (1988) Djorgovski, S., & Evans, C. R. 1988, The Astrophysical Journal, 335, L61
  * Draghis & Romani (2018) Draghis, P., & Romani, R. W. 2018, The Astrophysical Journal Letters, 862, L6
  * Draghis et al. (2019) Draghis, P., Romani, R. W., Filippenko, A. V., et al. 2019, The Astrophysical Journal, 883, 108
  * Fonseca et al. (2021) Fonseca, E., Cromartie, H., Pennucci, T. T., et al. 2021, The Astrophysical Journal Letters, 915, L12
  * Fruchter et al. (1988) Fruchter, A. S., Stinebring, D. R., & Taylor, J. H. 1988, Nature, 333, 237, doi: 10.1038/333237a0
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, The Astrophysical Journal, 887, 93
  * Husser et al. (2013) Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, Astronomy & Astrophysics, 553, A6
  * Kandel & Romani (2020) Kandel, D., & Romani, R. W. 2020, The Astrophysical Journal, 892, 101
  * Kandel et al. (2021) Kandel, D., Romani, R. W., & An, H. 2021, The Astrophysical Journal Letters, 917, L13
  * Kandel et al. (2020) Kandel, D., Romani, R. W., Filippenko, A. V., Brink, T. G., & Zheng, W. 2020, ApJ, 903, 39, doi: 10.3847/1538-4357/abb6fd
  * Kennedy et al. (2022) Kennedy, M. R., Breton, R. P., Clark, C. J., et al. 2022, MNRAS, 512, 3001, doi: 10.1093/mnras/stac379
  * Kulkarni et al. (1988) Kulkarni, S. R., Djorgovski, S., & Fruchter, A. S. 1988, Nature, 334, 504, doi: 10.1038/334504a0
  * Mukherjee (2017) Mukherjee, D. 2017, Journal of Astrophysics and Astronomy, 38, 1
  * Nieder et al. (2020) Nieder, L., Clark, C., Kandel, D., et al. 2020, The Astrophysical journal letters, 902, L46
  * Özel & Freire (2016) Özel, F., & Freire, P. 2016, Annual Review of Astronomy and Astrophysics, 54, 401
  * Ray et al. (2012) Ray, P., Abdo, A., Parent, D., et al. 2012, arXiv preprint arXiv:1205.3089
  * Ray et al. (2013) Ray, P., Ransom, S., Cheung, C., et al. 2013, The Astrophysical journal letters, 763, L13
  * Reynolds et al. (2007) Reynolds, M. T., Callanan, P. J., Fruchter, A. S., et al. 2007, Monthly Notices of the Royal Astronomical Society, 379, 1117
  * Romani (1990) Romani, R. W. 1990, Nature, 347, 741
  * Romani et al. (2014) Romani, R. W., Filippenko, A. V., & Cenko, S. B. 2014, The Astrophysical Journal Letters, 793, L20
  * Romani et al. (2015) —. 2015, The Astrophysical Journal, 804, 115
  * Romani et al. (2012) Romani, R. W., Filippenko, A. V., Silverman, J. M., et al. 2012, The Astrophysical Journal Letters, 760, L36
  * Romani et al. (2016) Romani, R. W., Graham, M. L., Filippenko, A. V., & Zheng, W. 2016, The Astrophysical Journal, 833, 138
  * Romani et al. (2021) Romani, R. W., Kandel, D., Filippenko, A. V., Brink, T. G., & Zheng, W. 2021, The Astrophysical Journal Letters, 908, L46
  * Romani et al. (2022) —. 2022, The Astrophysical Journal Letters, 934, L18
  * Romani et al. (2022a) Romani, R. W., Kandel, D., Filippenko, A. V., Brink, T. G., & Zheng, W. 2022a, ApJ, 934, L17, doi: 10.3847/2041-8213/ac8007
  * Romani & Sanchez (2016) Romani, R. W., & Sanchez, N. 2016, ApJ, 828, 7, doi: 10.3847/0004-637X/828/1/7
  * Romani et al. (2022b) Romani, R. W., Deller, A., Guillemot, L., et al. 2022b, ApJ, 930, 101, doi: 10.3847/1538-4357/ac6263
  * Sanchez & Romani (2017) Sanchez, N., & Romani, R. W. 2017, ApJ, 845, 42, doi: 10.3847/1538-4357/aa7a02
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, The Astrophysical Journal, 737, 103
  * Schroeder & Halpern (2014) Schroeder, J., & Halpern, J. 2014, The Astrophysical Journal, 793, 78
  * Stappers et al. (2001) Stappers, B., van Kerkwijk, M., Bell, J., & Kulkarni, S. 2001, The Astrophysical Journal Letters, 548, L183
  * Stovall et al. (2014) Stovall, K., Lynch, R., Ransom, S., et al. 2014, The Astrophysical Journal, 791, 67
  * Strader et al. (2019) Strader, J., Swihart, S., Chomiuk, L., et al. 2019, The Astrophysical Journal, 872, 42
  * van Kerkwijk et al. (2011) van Kerkwijk, M., Breton, R., & Kulkarni, S. 2011, The Astrophysical Journal, 728, 95
  * van Staden & Antoniadis (2016) van Staden, A. D., & Antoniadis, J. 2016, The Astrophysical Journal Letters, 833, L12
  * Voisin et al. (2020) Voisin, G., Kennedy, M., Breton, R., Clark, C., & Mata-Sánchez, D. 2020, MNRAS, 499, 1758


# On the spontaneous magnetization of two-dimensional ferromagnets.

D. Pescia  Laboratory for Solid State Physics, ETH Zurich, 8093 Zurich,
Switzerland<EMAIL_ADDRESS>A. Vindigni  Laboratory for Solid State
Physics, ETH Zurich, 8093 Zurich, Switzerland.  Department of Evolutionary
Biology and Environmental Studies, University of Zurich, 8057 Zurich,
Switzerland.

###### Abstract

Ferromagnetism is typically discussed in terms of the exchange interaction and
magnetic anisotropies. Yet real samples are inevitably affected by the
magnetostatic dipole-dipole interaction. Because of this interaction, a
theorem (R.B. Griffiths, Free Energy of interacting magnetic dipoles, Phys.
Rev. 176, 655 (1968)) forbids a spontaneous magnetization in, nota bene,
three-dimensional bodies. Here we discuss perpendicularly and in-plane
magnetized ferromagnetic bodies in the shape of a slab of finite thickness. In
perpendicularly magnetized slabs, magnetic domains are energetically favored
when the lateral size is sufficiently large, i.e. there is no spontaneous
magnetization. For in-plane magnetization, instead, spontaneous magnetization
is possible below a critical thickness which, in real thin films, could be as
small as few monolayers. At this critical thickness, we predict a genuine
phase transition to a multi-domain state. These results have implications for
two-dimensional ferromagnetism.

A situation familiar to ferromagnetism foresees that, below the Curie
temperature, the graph of the the free-energy as a function of the
magnetization has a flat portionIsing  (graph ”a” in Fig.1). This flatness
defines a situation in which the magnetization can acquire a value between
zero and a so called ”spontaneous magnetization” $\pm{M}_{0}$ without any
change in the free energy. This ”flatness” is a property of the thermodynamic
limit, i.e. of infinite bodies. In finite bodies, the free energy assumes a
shape that resembles the graph ”b” in Fig.1, $\pm M_{0}$ being the values at
which the free energy has minima. Systems with spontaneous magnetization are,
for example, the Ising and classical Heisenberg ferromagnets in three
dimensions (3D)Ising ; Jurg , the 2D-Ising model or the 2D-planar and
classical 2D-Heisenberg models with symmetry breaking single ion
interactionsJose . In real ferromagnetic bodies, however, the inevitable
dipole-dipole interaction, originating within Maxwell equation of
magnetostatics, must be considered alongside the main exchange interaction (of
purely quantum mechanical origin) and the single ion magnetic anisotropies
(produced by the spin-orbit interaction). The dipole-dipole interaction is,
typically, much weaker than the exchange interaction (by about two orders of
magnitude). Yet, it is long-ranged, as it decays only with the third power of
the distance between two magnetic moments. Because of the dipole-dipole
interaction, a theorem, proved by GriffithsGriff  for bodies with linear
dimension $L$ approaching infinity along all three spatial dimensions, implies
that any non-zero magnetization produces an increase of the free energy, i.e.
the graph of the free energy as a function of $M$ has a minimum at $M=0$ at
any temperature (graph ”c” in Fig.1). Accordingly, ferromagnetic order can
only be local and, globally, the spontaneous magnetization is exactly
vanishing.
This no-spontaneous magnetization rule is somewhat similar to the absence of
long-range order foreseen for the isotropic 2D-planar and 2D-Heisenberg
ferromagnetsMermin  but it refers, remarkably, to a 3D-body. In fact, it
appears that an important assumption underlying Griffith’s theorem is the size
of the body approaching infinity along all three spatial dimensions. In this
paper, we discuss ferromagnetism in the presence of exchange, magnetic
anisotropy and dipole-dipole interaction but in a slab-geometry, where only
two spatial dimensions are allowed to increase and the third is assigned a
finite thickness. Our results should be relevant for discussing ferromagnetism
in the new class of monolayer thin materials obtained by mechanical
exfoliationGong ; Huang ; Burch ; Li ; Novo ; Dai . They are known to be
perfectly flat over large distances and have been shown to be vertically
engineerableWua . As experiments are often analyzed in term of abstract
models, exact theorems such as Griffiths’s oneGriff  or scaling arguments such
as those presented here should allow experiments to distinguish those features
that are general and universal from those ones that originates from less known
details of a sample (such as defects).
Perpendicular Magnetization. We first analyze the situation of perpendicular
magnetization. In the state of spontaneous perpendicular magnetization, all
magnetic moments in the slab point along one of the two $z$-directions
perpendicular to the slab (Fig.2a), e.g. the $+z$-direction. In Fig.2a, this
state is rendered with a white color and the magnetization vector with
absolute value $M_{0}$ is given by a black arrow. $M_{0}$ represents the
magnetic moment per volume of the unit cell, i.e.
$M_{0}\doteq\frac{g\cdot\mu_{B}\cdot S}{a^{3}}$ (with $S$ being the spin in
units of $\hbar$ and $a$ the lattice constant). A possible, elementary state
of vanishing spontaneous magnetization is shown in Fig.2b. One half of the
slab is still filled with magnetic moments pointing upwards ”$\uparrow$” but
the other half (gray in Fig.2b) contains magnetic moments pointing downwards
”$\downarrow$” (indicated as state $\uparrow\downarrow$ henceforth). Assuming
that Griffith’s theorem is valid for the perpendicular magnetization in the
slab geometry as well, the $\uparrow\downarrow$ state should have a lower
total energy than the state of spontaneous uniform magnetization (indicated as
state $\uparrow\uparrow$ henceforth). We, therefore, proceed to compare the
energy of the two states in a situation where the size $L$ of the slab is much
larger than its thickness $d$.
Magnetostatic energy $E_{M}$. The magnetostatic energy $E_{M}$ for the
perpendicular magnetization configuration is most appropriately computed as
the energy of the interacting Ampèrian effective currect densities
$\vec{\nabla}\times\vec{M}(\vec{r})$ resulting from the magnetization vector
$\vec{M}(\vec{r})$ (see Section I of the Supplementary Material (SM)). The
current density vectors arising in the $\uparrow\uparrow$ and
$\uparrow\downarrow$ configurations are summarized by red arrows in Fig.2. We
find that the formation of the two neighboring domains with opposite
magnetization lowers the total magnetostatic energy and is thus the driving
force behind the suppression of the spontaneous magnetization, i.e.
$E_{M}(\uparrow\downarrow)-E_{M}(\uparrow\uparrow)$ is negative. The self-
energies of the Ampèrian currents circulating along the perimeter of the slab
cancel out exactly from $E_{M}(\uparrow\downarrow)-E_{M}(\uparrow\uparrow)$.
Their mutual interaction provides terms of the order $L\cdot d^{2}$. The
leading contribution to $E_{M}(\uparrow\downarrow)-E_{M}(\uparrow\uparrow)$ is
the self-energy of the current flowing along the wall that separates the
domains. Assuming that the domain wall has a finite thickness $w$, the leading
contribution writes (Section I, SM), in the limit $d\\!<\\!<w\\!<\\!<\\!L$

$E_{M}(\uparrow\downarrow)\\!-\\!E_{M}(\uparrow\uparrow)\\!\approx\\!\frac{L\\!\cdot\\!d}{a^{2}}\\!\cdot\\!\left[-\frac{2}{\pi}\cdot(\Omega\\!\cdot\\!\frac{d}{a})\\!\cdot\\!\ln\frac{L}{w}\right]+\\!{\cal
O}(\frac{d^{2}\\!\cdot\\!L}{a^{3}})$ (1)

In Eq.1, $L\\!\cdot\\!d$ is the surface of the domain wall.
$\Omega\doteq\frac{\mu_{0}}{2}\cdot M_{0}^{2}\cdot a^{3}$ is a parameter used
for expressing the strength of the magnetostatic energy per unit cell. As an
example, $\Omega$ for metallic Fe amounts to $\approx 0.28$ meV (see Section I
of SM)omega . The slab model (see Section I of SM) shows explicitly that the
relevant coupling constant entering
$E_{M}(\uparrow\downarrow)\\!-\\!E_{M}(\uparrow\uparrow)$ is not $\Omega$
itself but $\Omega\cdot\\!\frac{d}{a}$, i.e. the characteristic magnetostatic
energy per unit surface cell. The logarithmic term in Eq. 1 provides,
formally, a divergence of
$E_{M}(\uparrow\downarrow)\\!-\\!E_{M}(\uparrow\uparrow)$ with the size $L$.
It is universal in the sense that it does not depend on the exact geometry of
the wall separating the domains: both the shape of this wall and the exact
shape of slab contribute only terms of the order
$\Omega\\!\cdot\\!L\cdot\\!d^{2}$.
Energy of the wall between $\uparrow$-and $\downarrow$-domains. The formation
of a domain wall in the $\uparrow\downarrow$-state increases the total energy
of the ferromagnetic slab and therefore promotes the state of spontaneous
magnetization. Within the wall, the magnetic moments rotate away from the
$z$-direction. For simplicity, we assume the wall to run parallel to the
$y$-direction and the rotation to take place along the $x$-direction. Let the
rotation be characterized by an angle $\theta$, which increases from $0$ to
$\pi$ when moving along $x$-within the wall. The misalignement is associated,
in the first place, with an increase of the single ion magnetic anisotropy
energy that favors the perpendicular magnetization, introduced conceptually by
NéelNeel  and computed for the first time from first principles by Gay and
Richter Gay  for the monolayer of Fe. This term originates from the breaking
of translational symmetry perpendicular to the slab plane. In ultrathin slabs,
it is only weakly dependent on $d$GR  as it arises from the two surfaces
bounding the slab. Using the convention of Ref.Gay , we write this term as
$-\lambda\cdot\cos^{2}\theta(x)$, the parameter $\lambda$ to be intended, as
in Ref.Gay , as an energy per surface unit cell. For $\lambda$, a value of
$\approx 0.4$ meV is reportedGay ; Shi  for the one monolayer of Fe.
The Néel-anisotropy is not the only contribution to the magnetic anisotropy
affecting the magnetic moment rotation in the wall. The magnetostatic energy
itself favors the magnetic moment to lie within the slab plane and contributes
a term $+\Omega\\!\cdot\\!\frac{d}{a}\cdot\\!\cos^{2}\theta(x)$ (Section II of
SM). The coefficient of this contribution scales with $d$, see Section II of
SM and Ref. GR , in contrast to $\lambda$.
Finally, for the building of the wall, one must also consider that the
misalignment of two neighboring magnetic moments at the sites $x$ and $x\pm a$
increases the energy bySmall  $J\cdot S^{2}\cdot\cos\left(\theta(x\pm
a)-\theta(x)\right)$, $J$ being the exchange coupling energy per spin couple.
For bulk Fe, Ref.Small  estimates $J\approx 48$ meV. The total energy of a
wall is proportional to its surface $L\cdot d$ multiplied by the geometric
mean of the relevant coupling constants (see Section III of SM for a short
reminding on how the term $+\Omega\\!\cdot\\!d\\!\cdot\\!\cos^{2}\theta(x)$ is
embedded into the wall energy):

$\propto\frac{L\cdot
d}{a^{2}}\\!\cdot\\!\left[2\cdot\\!\sqrt{\left(\lambda-\Omega\\!\cdot\\!\frac{d}{a}\right)\\!\cdot\\!2\\!\cdot\\!J\cdot
S^{2}}\right]$ (2)

Absence of spontaneous magnetization and crossover length. Comparing the
domain wall energy cost to the magnetostatic energy gain, we recognize that
the logarithmic term always favors the building of domains for sufficiently
large $L$. Accordingly, Griffith’s theorem about the absence of spontaneous
magnetization in the thermodynamic limit holds true in the slab geometry with
perpendicular magnetization.
One interesting outcome of our argument is the estimate of the cross-over
length $L_{c}$ at which a slab will transit from a monodomain state to a
multi-domain state, resulting from equating the magnetostatic energy gain to
the wall energy:

$L_{c}\approx w\cdot
e^{\frac{\pi\cdot\sqrt{\left(\lambda-\Omega\\!\cdot\\!\frac{d}{a}\right)\cdot
2\cdot J\cdot S^{2}}}{\Omega\\!\cdot\\!\frac{d}{a}}}$ (3)

There are three aspects of this result that deserve amplification. First, if
we insert the values for $J$, $\Omega$ and $\lambda$ given previously, one
recognizes that $L_{c}$ assumes astronomically large values for the monolayer
limit (In Section III of SM we find that the ”monolayer limit” corresponds, in
the slab model, to $d\\!\approx\\!a$; in this limit, the expression $L_{c}$ is
consistent with known resultsPighin ).
Second, one recognizes a threshold thickness
$\frac{d_{R}}{a}=\frac{\lambda}{\Omega}$ at which, formally, the argument of
the exponential function vanishes. $d_{R}$ is in the subnanometer rangeBader ;
Shi . We therefore propose that $L_{c}$ decreases exponentially with $d$ and,
toward $d_{R}$, it assumes values of few tens of micrometers. These are the
lateral lenghts over which exfoliated two-dimensional magnets are believed to
be almost perfectly flatGong ; Huang ; Burch ; Li ; Novo ; Dai . Accordingly,
a sequence of exfoliated samples with suitable thickness and with increasing
lateral size $L$ should allow an insight into the yet unexplored mechanism of
penetration of magnetic domains in two-dimensional ferromagnetic elements as a
function of their size $L$. Some preliminary results in this direction were
reported in Ref.Oliver  on epitaxially grown ultrathin films.
Third, the spin wave excitations produce a renormalization of the various
coupling constants $J$, $\Omega$ and $\lambda$ as a function of temperature
Politi , so that $d_{R}$ is itself a function of the temperature. One finds
that $d_{R}(T)$ defines a line of phase transitions at which the perpendicular
magnetization turns into the plane of the slabBader ; Shi ; Politi .
Accordingly, the exponential decay of $L_{c}$ should become observable when
the temperature $T$ is increased and $d_{R}(T)$ approaches the actual
thickness $d$ of the slab.
A final comment is dedicated to the stability against a perpendicular magnetic
field of the domain phase that should appear at sufficiently large $L$. In
Section IV of the SM, we analyze this problem by considering the energetics of
one stripe of reversed magnetization $-M_{0}$, embedded into a background with
magnetization $+M_{0}$, subject to a perpendicularly applied field with
strength $+B_{0}$. We find that the state of uniform magnetization becomes the
energetically favored one when $B_{0}$ exceeds a threshold strength
$B_{t}\propto\left(\mu_{0}\cdot M_{0}\cdot\frac{d}{L_{c}}\right)$. Far away
from the $d_{R}(T)$-transition line, the threshold field might be as small as
few $nT$. Close to the transition line we expect this field to be in the
$mT$-rangeNiculin ; NiculinPRB .
In-plane magnetization. We now analyze the slab geometry with in-plane
magnetization. The two configurations considered are one of uniform
magnetization along e.g. the $+x$-direction (Fig.3a), and one where one half
of the slab has magnetization along the $+x$-direction and the other half has
magnetization along $-x$ (Fig.3b). In this situation, the magnetostatic energy
is most appropriately computed as the Coulomb interaction between effective
”charges” produced by $\vec{\nabla}\cdot\vec{M}(\vec{r})$. The charges
resulting from the two spin configurations are indicated in Fig.3 in red. The
change of magnetostatic energy produced by the building of in-plane domains is
negative, i.e. the magnetostatic energy favor a state of vanishing spontaneous
magnetization. However, the logarithmic term produced by the self-energies
cancel out exactly when the energies of the two states are subtracted,
provided the wall is parallel to the magnetization, i.e. the wall is not
charged. The remaining contributions provide terms proportional to
$-(\Omega\cdot d)\cdot L\cdot d$ (Section V, SM). Again, there is a wall
between the two domains, in which spins rotate away from the $x$-direction. We
assume, for simplicity, an in-plane uniaxial anisotropy with the strength
$\Lambda$ such as the one encountered in ultrathin Fe films on W(110)Grad .
This uniaxial anisotropy provides and energy barrier against rotations away
from the $x$ in-plane direction. A typical value for $\Lambda$, desumed from
Ref.Grad , is $0.04$ meV per unit surface cellGrad2 . Notice that $\Lambda$
originates from the rectangular nature of the surface unit cell and is about
two-orders of magnitude smaller than the Néel magnetic anisotropy constant
forthcoming in the perpendicular magnetization configuration. $\Lambda$ is
rather of the same order of magnitude as the quartic in-plane magnetic
anisotropy constant computed e.g. in Ref.Gay . Given $\Lambda$, the energy of
the wall per wall surface unit cell is then proportional to the geometric
average of the exchange coupling $J$ and $\Lambda$, by virtue of the same
arguments exposed in Section 3 of SM, i.e. $\approx L\cdot d\cdot
2\\!\cdot\sqrt{\Lambda\\!\cdot\\!2\cdot J\cdot S^{2}}$.
Equating the total energy change due to the formation of a domain wall to zero
provides an estimate of the critical thickness $d_{c}$ below which a state of
spontaneous magnetization is favored:

$d_{c}\propto a\cdot\frac{\sqrt{\Lambda\cdot J}}{\Omega}$ (4)

As $L$ cancels out, we argue that it should be possible to find a rigorous
proof of spontaneous magnetization even in the thermodynamic limit
$L\rightarrow\infty$ in a truly 2D system with in-plane magnetization. As both
states above and below $d_{c}$ are stable for $L\rightarrow\infty$, the
transition at $d_{c}$ from a single-domain state with spontaneous
magnetization to a multi-domain stripe state should be a genuine phase
transition. Using the values for $\Lambda,J,\Omega$ introduced in this Letter,
we obtain $d_{c}\approx$ five lattice constants or less. This small number
means that a sample must be fabricated with uniform thickness over large
lateral distances, in order for this transition to be observed. The new class
of two dimensional magnetsGong ; Huang ; Burch ; Li ; Novo ; Dai  might
provide this kind of precision. A final remark: provided that $\Omega$, $J$,
and $\Lambda$ renormalize slightly differently with temperature, we might
expect a $d_{c}(T)$ line of phase transitions which can also be crossed at a
fixed thickness $d$ by varying the temperature. This situation would represent
the analogon to a similar phase transition observed in perpendicularly
magnetized filmsNiculin .

Supplementary Material. Details of the computations used to obtain Eq.1,2,3
and 4 are given in the ”Supplemental Material”.
Author contribution: D.P. wrote the paper and A. V. contributed to the
scientific content and to the writing.
Competing Interests: The authors declare that they have no competing financial
interests.
Data availability statement: The informations that support the findings of
this study are available within the article and its supplementary material.

## References

  * (1) R.B.Griffiths, Peierls Proof of Spontaneous Magnetization in a Two-Dimensional Ising Ferromagnet, Phys. Rev. 136, A437-A439 (1964)
  * (2) J. Fröhlich and B. Simon, Infrared Bounds, Phase Transitions and Continuous Symmetry Breaking, Commun. Math. Phys. 50, 79-85 (1976)
  * (3) Jorge V. Jose, Leo P. Kadanoff, S. Kirkpatrick, David R. Nelson, Renormalization, vortices, and symmetry-breaking perturbations in the two-dimensional planar model, Phys. Rev. B. 16, 1217-1241 (1977)
  * (4) R.B. Griffiths (1968), Free Energy of interacting magnetic dipoles, Phys. Rev. 176, 655-659, (1968)
  * (5) N. D. Mermin, H. Wagner, Absence of Ferromagnetism or Antiferromagnetism in One- or Two-Dimensional Isotropic Heisenberg Models, Phys. Rev. Lett. 17, 1133-1136 (1966)
  * (6) C. Gong et al., Discovery of intrinsic ferromagnetism in two-dimensional van der Waals crystals, Nature 546, 265-269 (2017)
  * (7) B. Huang et al., Layer-dependent ferromagnetism in a van der Waals crystal down to the monolayer limit. Nature 546, 270-273 (2017).
  * (8) K.S. Burch, D. Mandrus and Je-Geun Park, Magnetism in two-dimensional van der Waals materials, Nature 563, 47-52 (2018)
  * (9) Yang Li, Baishun Yang, Shengnan Xu, Bing Huang, and Wenhui Duan, Emergent Phenomena in Magnetic Two-Dimensional Materials and van der Waals Heterostructures, ACS Appl. Electron. Mater. 4 3278-3302 (2022)
  * (10) M. Gibertini, M. Koperski, A. F. Morpurgo, and K. S. Novoselov, Magnetic 2D materials and heterostructures, Nature Nanotechnology 408, 408-419 (2019)
  * (11) Chuying Dai et al, Sci. Research progress of two-dimensional magnetic materials, China Mater. 66 859-876 (2023)
  * (12) Yecun Wua et al., Interlayer engineering of Fe3GeTe2: From 3D superlattice to 2D monolayer, Proc.Nat.Ac.Sci. 121, (4) e2314454121 (2024)
  * (13) In bulk Fe ($a\approx 0.28$nm) there are two magnetic moments per unit cell, each carrying a magnetic moment of about $2.2$ $\mu_{B}$ ($g\approx 2$, $S\approx 1.1$Small ), leading to $\Omega\approx 0.28$ meV (see also Section I of SM). For the monolayer of Fe, Gay and RichterGay  give $\Omega\approx 0.3$ meV.
  * (14) L. M. Small and V Heine, A couple method for calculating interatomic interactions in itinerant electron magnetic systems, J. Phys. F: Met. Phys. 14, 3041-3052 (1984)
  * (15) J. G. Gay and Roy Richter, Spin Anisotropy of Ferromagnetic Films, Phys. Rev. Lett. 56, 2728-2731 (1986)
  * (16) J. G. Gay and Roy Richter, Spin anisotropy of ferromagnetic slabs and overlayers, J. Appl. Phys. 61 3362-3365 (1987)
  * (17) M.L. Néel, Anisotropie magnetique superficielle et surstructures d’orientation, J. Phys. Radium, 15, 225-239 (1954)
  * (18) S. A. Pighin, O. V. Billoni, and S. A. Cannas, Finite-temperature phase diagram of ultrathin magnetic films without external fields, Phys. Rev. E. 86, 051119 (2012)
  * (19) Long-Pei Shi, Perpendicular magnetic anisotropy in ultrathin fcc iron films and surfaces at finite temperature,Physics Letters A 189 409-414 (1994)
  * (20) Z.Q. Qiu, J. Pearson and S,.D. Bader, Asymmetry of the Spin Reorientation Transition in Ultrathin Fe Films and Wedges Grown on Ag(100), Phys. Rev. Lett. 70 1006-1009 (1993)
  * (21) O. Portmann t al., Micromagnetism in the ultrathin limit, Thin Solid Films 505, 2-9 (2006)
  * (22) P. Politi, A. Rettori, M.G. Pini and D. Pescia, Magnetic Phase Diagram of a Thin Film with a Reorientation Transition, Europhysics Letters 28 71 (1994); DOI 10.1209/0295-5075/28/1/013
  * (23) N. Saratz et al, Critical exponents and scaling invariance in the absence of a critical point, Nat. Commun. 7, 13611 (2016);
  * (24) N. Saratz et al, Irreversibility, reversibility, and thermal equilibrium in domain patterns of Fe films with perpendicular magnetization, Phys. Rev. B 82, 184416 (2010)
  * (25) U. Gradmann, J. Korecki, and G. Waller, In-Plane Magnetic Surface Anisotropies in Fe(110), Appl. Phys. A 39, 101-108 (1986)
  * (26) In Ref.Grad , the anisotropy parameter is given as $0.105\cdot\frac{\text{erg}}{\text{cm}^{2}}$. To convert this energy in meV per unit cell we use $1\text{erg}=6.242\cdot 10^{11}eV$. In the [110]-plane of bcc Fe we have $\frac{\sqrt{2}}{2.83^{2}}\cdot\frac{\text{Atom}}{\text{Angstrom}^{2}}$ so that the magnitude of the two-fold in-plane magnetic anisotropy is $\approx 4\cdot 10^{-2}\cdot\text{meV}$.

Figure 1: The free energy as a function of $M$. a: the graph for an infinite
ferromagnetic body with spontaneous magnetization has a flat portion between
$\pm M_{0}$ below the Curie temperature, see Ref.Ising . b: for a finite
ferromagnetic body, the free energy has minima at $\pm M_{0}$. c: A theorem by
Griffiths Griff  implies that the free energy has a minimum at $M=0$ at any
temperature.  Figure 2: a: the state of uniform perpendicular magnetization
(represented in white) in a slab. The magnetization vector is represented by
the vertical black arrow. Red arrows represent the effective current density
vectors flowing along the perimeter of the slab. b: The slab is filled by two
domains with magnetization vector parallel (white domain) and antiparallel
(gray domain) to the vertical $z$-axis. Red arrows represent the effective
current density vectors.  Figure 3: a: the state of uniform in-plane
magnetization (represented in white) in a slab. The magnetization vector is
represented by the horizontal black arrow. Effective charge densities (their
sign being given in red) appear along the perimeter. b: The slab is filled by
two domains with magnetization vector parallel (white domain) and antiparallel
(gray domain) to the horizontal $x$-axis. The sign of the effective charge
densities is given in red.