#  JWST Detects Neon Line Variability in a Protoplanetary Disk

C. C. Espaillat Institute for Astrophysical Research, Department of Astronomy,
Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA  T.
Thanathibodee Institute for Astrophysical Research, Department of Astronomy,
Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA  C. V.
Pittman Institute for Astrophysical Research, Department of Astronomy, Boston
University, 725 Commonwealth Avenue, Boston, MA 02215, USA  J. A. Sturm Leiden
Observatory, Leiden University, PO Box 9513, NL–2300 RA Leiden, The
Netherlands  M. K. McClure Leiden Observatory, Leiden University, PO Box 9513,
NL–2300 RA Leiden, The Netherlands  N. Calvet Department of Astronomy,
University of Michigan, 1085 South University Avenue, Ann Arbor, MI 48109, USA
F. M. Walter Department of Physics and Astronomy, Stony Brook University,
Stony Brook, NY 11794-3800, USA  R. Franco-Hernández Instituto de Astronomía y
Meteorología, Universidad de Guadalajara, Avenida Vallarta No. 2602, Col.
Arcos Vallarta Sur, CP 44130, Guadalajara, Jalisco, Mexico  J. Muzerolle Page
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218,
USA

###### Abstract

We report the first detection of variability in the mid-infrared neon line
emission of a protoplanetary disk by comparing a JWST MIRI MRS spectrum of SZ
Cha taken in 2023 with a Spitzer IRS SH spectrum of this object from 2008. We
measure the [Ne III]-to-[Ne II] line flux ratio, which is a diagnostic of the
high-energy radiation field, to distinguish between the dominance of EUV- or
X-ray–driven disk photoevaporation. We find that the [Ne III]-to-[Ne II] line
flux ratio changes significantly from $\sim 1.4$ in 2008 to $\sim 0.2$ in
2023. This points to a switch from EUV-dominated to X-ray–dominated
photoevaporation of the disk. We present contemporaneous ground-based optical
spectra of the H$\alpha$ emission line that show the presence of a strong wind
in 2023. We propose that this strong wind prevents EUV radiation from reaching
the disk surface while the X-rays permeate the wind and irradiate the disk. We
speculate that at the time of the Spitzer observations, the wind was
suppressed and EUV radiation reached the disk. These observations confirm that
the MIR neon emission lines are sensitive to changes in high-energy radiation
reaching the disk surface. This highlights the [Ne III]-to-[Ne II] line flux
ratio as a tool to gauge the efficiency of disk photoevaporation in order to
provide constraints on the planet-formation timescale. However,
multiwavelength observations are crucial to interpret the observations and
properly consider the star-disk connection.

accretion disks, stars: circumstellar matter, planetary systems:
protoplanetary disks, stars: formation, stars: pre-main sequence

##  1 Introduction

The protoplanetary disk lifetime is influenced by the rate at which gas is
eroded by photoevaporative winds created by high-energy radiation from the
central star (see review by Pascucci et al., 2022). However, photoevaporation
models predict significantly different mass-loss rates depending on whether
the photoevaporation is driven by the extreme ultraviolet (EUV) or the X-ray
radiation field. These different mass-loss rates lead to different timescales
for disk clearing and hence the timescales by which giant planets must form.
Models of EUV photoevaporation predict low mass-loss rates of $\sim$10-10
$M_{\sun}$ yr-1 (Font et al., 2004; Alexander et al., 2006) while X-ray
photoevaporation achieves high mass-loss rates up to $\sim$10-8 $M_{\sun}$
yr-1 (Gorti & Hollenbach, 2009; Owen et al., 2011). Note that the mass-loss
rates can change based on model assumptions such as chemical coolants or the
adopted ultraviolet and X-ray spectra (Pascucci et al., 2022). Initially,
viscous evolution dominates disk evolution, leading to a decrease in the mass
accretion rate with time as disk material is depleted onto the star (i.e.,
less gas in the disk leads to a slower accretion rate onto the star). Once the
rate of mass accretion inward through the disk equals the outward
photoevaporative mass-loss rate, the inner disk cannot be replenished and is
accreted onto the star in $<10^{5}$ yrs (Clarke et al., 2001).

The onset of rapid disk clearing is reached later for the EUV photoevaporation
model since it takes longer for the disk accretion rate to decline to $\sim
10^{-10}$ $M_{\sun}$ yr-1 (Hartmann et al., 1998). In the case of X-ray–driven
photoevaporation, disks should quickly dissipate after reaching a mass
accretion rate of $\sim 10^{-8}$ $M_{\sun}$ yr-1. The average mass accretion
rate of young ($\sim 1$ Myr old) low-mass pre-main sequence stars (i.e., T
Tauri stars; TTS) with disks is $\sim 10^{-8}$ $M_{\sun}$ yr-1 (Hartmann et
al., 1998), suggesting that X-ray–driven photoevaporation often dominates disk
clearing. However, the existence of young stars surrounded by substantial
disks with mass accretion rates $<10^{-10}$ $M_{\sun}$ yr-1 indicates that
X-ray–driven photoevaporation may not always dominate (Ingleby et al., 2011;
Thanathibodee et al., 2023). The observed disk frequency can also be explained
by low photoevaporation rates (Manzo-Martínez et al., 2020). The main obstacle
in distinguishing between photoevaporation models is the lack of observational
constraints on the high-energy emission of stars, particularly since our
knowledge of the strength of EUV emission from TTS is limited due to
absorption of EUV photons by neutral hydrogen along the line of sight.

Infrared (IR) fine structure lines can determine the high-energy radiation
fields impinging on the disk surface, which can be used to distinguish between
photoevaporation models. Many theoretical works have examined the importance
of the neon line luminosities and their connection to high-energy radiation
(Glassgold et al., 2007; Meijerink et al., 2008; Ercolano et al., 2008; Gorti
& Hollenbach, 2008; Hollenbach & Gorti, 2009; Schisano et al., 2010; Ercolano
& Owen, 2010). Due to the high first and second ionization potentials of neon
(21.56 eV and 41.0 eV, respectively), the [Ne II] and [Ne III] lines lie
$>1,000$ K above the ground state, making them ideal to probe high-energy
radiation fields.

Both EUV and X-ray radiation have been proposed to be responsible for neon
forbidden lines in disks around TTS (Glassgold et al., 2007; Hollenbach &
Gorti, 2009), while far-ultraviolet photons have insufficient energy to ionize
neon (Glassgold et al., 2007). [Ne II] was detected with the Infrared
Spectrograph (IRS) on Spitzer in more than 50 TTS (Espaillat et al., 2007;
Pascucci et al., 2007; Lahuis et al., 2007; Flaccomio et al., 2009; Güdel et
al., 2010; Baldovin-Saavedra et al., 2011; Szulágyi et al., 2012). However,
comparing observations and theories has not led to clear results. High-
resolution ground-based observations of the [Ne II] line are fit equally well
by different photoevaporation models (e.g., Pascucci et al., 2011).

One way to distinguish between EUV and X-ray creation of neon fine-structure
emission is by measuring the [Ne III]-to-[Ne II] line flux ratio. A hard EUV
spectrum ($L_{\nu}$ $\propto$ $\nu^{-1}$) can lead to [Ne III]-to-[Ne II] line
flux ratios $>1$ (see Figure 1 in Hollenbach & Gorti, 2009). Neon ion
production by X-rays and/or a soft EUV spectrum can lead to [Ne III]-to-[Ne
II] line flux ratios of the order $\sim 0.1$ (Glassgold et al., 2007;
Meijerink et al., 2008). Of the more than 50 disks with [Ne II] line
detections, only five have [Ne III] line detections; of these five [Ne III]
line detections, four have [Ne III]/[Ne II]$<0.1$ (Lahuis et al., 2007; Najita
et al., 2010; Espaillat et al., 2013). One of the five disks (SZ Cha) has [Ne
III]/[Ne II] $>1$, showing that EUV radiation can dominate the creation of
neon forbidden lines (Espaillat et al., 2013).

SZ Cha is a TTS with a spectral type of K0–K2 (Rydgren, 1980; Manara et al.,
2014), a stellar accretion rate of $\sim 0.2$–$1.6\times 10^{-8}$ $M_{\sun}$
yr-1 (Espaillat et al., 2011; Manara et al., 2014; Espaillat et al., 2019b),
and an X-ray luminosity of $\sim 1\times 10^{30}$ erg s-1 (Ingleby et al.,
2011; Espaillat et al., 2019b). The IR spectral energy distribution (SED)
indicates that SZ Cha is surrounded by a pre-transitional disk (i.e., an inner
disk separated from an outer disk by a large gap $\sim$ tens of astronomical
units wide; Kim et al., 2009; Espaillat et al., 2014), and SED modeling points
to a gap size of about 20 au (Espaillat et al., 2011; Ribas et al., 2016).
SPHERE scattered-light near-infrared images detect disk rings at 40, 70, and
110 au from the star (Asensio-Torres et al., 2021). Lastly, SZ Cha displays
“seesaw” variability in the mid-infrared (MIR) as seen with Spitzer (Espaillat
et al., 2011). Seesaw variability is characterized by a decrease in emission
at shorter IR wavelengths with a corresponding increase at longer IR
wavelengths and visa versa (Muzerolle et al., 2009; Espaillat et al., 2011;
Flaherty et al., 2012); this behavior has been attributed to variable
shadowing of the outer disk by the inner disk (Espaillat et al., 2011).

Here we analyze a JWST spectrum of SZ Cha as well as a CHIRON optical spectrum
taken within 15 hours of the JWST observations in 2023. We compare the JWST
data to a Spitzer spectrum of SZ Cha taken in 2008. In Section 2, we present
the new MIR and optical spectra as well as review an archival Spitzer spectrum
and several archival optical spectra of SZ Cha. In Section 3, we measure the
MIR atomic emission lines and model the optical H$\alpha$ emission line. In
Section 4, we compare the 2023 and 2008 data sets and interpret the results in
the context of theoretical models. We end with a summary and conclusions in
Section 5.

Figure 1:  Left: JWST MIRI MRS spectrum from 2023 (black) and Spitzer IRS SH
spectrum from 2008 (purple) of SZ Cha. We label the [Ar II], H2, [Ne II], and
[Ne III] emission lines and a polycyclic aromatic hydrocarbon feature around
11 µm. Directly to the left of the [Ar II] line is an H2 line. To the right of
the [Ar II] line is a feature where [He II], [H I], and [Ni II] emission lines
may overlap. Right: Comparisons between the [Ne II] and [Ne III] lines within
the same observation.

##  2 Observations and Data Reduction

Here we present a new JWST Mid-InfraRed Instrument Medium Resolution
Spectrometer (MIRI MRS) spectrum of SZ Cha and a new high-resolution CHIRON
optical spectrum taken contemporaneously (within 15 hrs) of the JWST
observations in 2023. In addition, we revisit a previously published Spitzer
IRS spectrum of SZ Cha from 2008 and archival optical spectra taken at
Magellan and ESO between 2007 to 2014.

Figure 2:  Line fits (orange) to MIR spectra of [Ne II] 12.81 µm, [Ne III]
15.55 µm, and [Ar II] 6.99 µm (blue lines). The top row corresponds to _JWST_
and the bottom row to _Spitzer_.

###  2.1 MIR Spectra

SZ Cha was observed with JWST on 2023 April 21 starting at 10:12 UT using MIRI
(Rieke et al., 2015; Wright et al., 2023) MRS (Wells et al., 2015) as part of
GO program 1676 (PI: C. Espaillat) in Cycle 1. Dedicated background
observations were also obtained. The target and background observations had an
exposure time of 336 s and used the default four-point dither.

We use the JWST Science Calibration Pipeline v.1.11.4 (Bushouse et al., 2023)
to reduce the uncalibrated raw MRS data using the calibration reference file
version jwst_1121.pmap. In Stage 1 of the pipeline the raw data or ramps are
converted into uncalibrated slope images using the default Detector1 class. In
Stage 2 the dedicated background is subtracted from the slope images, which
are subsequently calibrated using the Spec2 class with additional residual
fringe correction. Finally, the composite data cube from all dithered
exposures is created from the calibrated slope images in Stage 3.

The final reduced spectrum is presented in Figure 1 (left) along with a
Spitzer IRS Short-High (SH) spectrum from GO Program 40247, which was
published in Espaillat et al. (2013). There are significant differences in the
continuum level between the _JWST_ spectrum and the Spitzer IRS SH spectrum.
These differences are especially large in channel 4 (17.7–27.9 µm). We do not
show data beyond about 24 µm, at which point the spectrum appears to be
affected by the reported sensitivity decline in MIRI MRS. As noted earlier,
the MIR continuum of SZ Cha is known to be variable, but further analysis of
channel 4 is necessary to disentangle intrinsic variability from flux
calibration issues and we will address this in future work. In Figure 1
(right), we highlight the change in the Neon line ratio between the Spitzer
(top) and JWST observations (bottom).

###  2.2 Optical Spectra

A high-resolution ($R\sim 25,000$) optical spectrum of SZ Cha was obtained on
2023 April 22 starting at 01:23 UT with the CHIRON bench-mounted, fiber-fed,
cross-dispersed echelle spectrograph (Schwab et al., 2010) on a 1.5-m
telescope that is part of the Small and Moderate Aperture Research Telescope
System (Tokovinin et al., 2013) at Cerro Tololo Inter-American Observatory.
The exposure time was 3600 s. The CHIRON observations were taken within $\sim
15$ hrs after the _JWST_ observations. The data were obtained in fiber mode
with $4\times 4$ on-chip binning, leading to $R\sim 27,800$. The data were
reduced following the standard steps using an IDL
pipeline.111http://www.astro.sunysb.edu/fwalter/SMARTS/CHIRON/ch_reduce.pdf

We also searched for archival optical spectra of SZ Cha. High-resolution
($R\sim 35,000$) optical spectra of SZ Cha were taken with the Magellan
Inamori Kyocera Echelle (MIKE; Bernstein et al., 2003) on the Magellan-Clay
telescope at Las Campanas Observatory in Chile on the following UT dates: 2007
February 10, February 11; 2008 February 15, February 18; 2009 January 18; 2011
January 5; 2012 February 17. The Spitzer SH spectrum was taken on 2008 August
16, so we have MIKE spectra taken about six months before (2008 February 15,
February 18) the Spitzer spectrum and another taken about five months after
(2009 January 18). The data were reduced using the Image Reduction and
Analysis Facility (IRAF; Tody, 1993) tasks CDPROC, APFLATTEN, and DOECSLIT. We
also retrieved spectra from the ESO Science Archive
Facility.222http://archive.eso.org The ESO data were taken with X-shooter
(Vernet et al., 2011) on 2010 January 19, UVES (Dekker et al., 2000) on 2014
April 17, and FEROS (Kaufer et al., 1999) on 2014 May 31. We present these
optical spectra in Section 3.2.

Table 1: MIR Atomic Line Fluxes and Ratios for SZ Cha  Date | Instrument | [Ar II] 6.99 µm | [Ne II] 12.81 µm | [Ne III] 15.55 µm | [Ne III]-to-[Ne II] | [Ne II]-to-[Ar II]
---|---|---|---|---|---|---
(UT) |  | (10-15 erg cm-2 s-1) | (10-15 erg cm-2 s-1) | ( 10-15erg cm-2 s-1) | Ratio | Ratio
2008 Aug 16 |  Spitzer IRS SH | … | $15.5\pm 2.9$ | $22.1\pm 3.4$ | $1.4\pm 0.4$ | …
2023 Apr 21 |  JWST MIRI MRS | $3.2\pm 0.6$ | $6.0\pm 1.6$ | $1.2\pm 0.2$ | $0.20\pm 0.06$ | $1.9\pm 0.6$
Figure 3:  H$\alpha$ profiles of SZ Cha (blue) and the accretion flow modeling
fits (black/grey). We show the average accretion flow model fit (black line)
from the top 100 best-fitting models (gray lines). The horizontal dotted line
denotes the continuum, and the vertical dotted line marks the line center.
Horizontal green lines mark regions included in the fit.

##  3 Analysis and Results

In Section 3.1, we identify and measure MIR atomic emission lines. No evidence
of any strong molecular emission besides H2 in SZ Cha has emerged. In Section
3.2, we model the H$\alpha$ emission lines in the CHIRON optical spectra to
obtain mass accretion rates.

###  3.1 MIR Atomic Lines

As depicted in Figure 1, we detect the following atomic lines in the _JWST_
spectrum: [Ar II] 6.99 µm, [Ne II] 12.81 µm, and [Ne III] 15.55 µm. We note
that [S III] 18.71 µm and [Ar III] 21.83 µm appear to be present as well, but
these are located in channel 4, and due to the issues noted in Section 2.1, we
defer analysis of these lines to future work.

We fit each line and its underlying continuum with a composite model of a
Gaussian and a linear function within $\pm 700$ km s-1 from the line center
using SciPy’s curve_fit. The fits to the [Ar II], [Ne II], and [Ne III] lines
in the _JWST_ data are shown in the top panels of Figure 2. Line fluxes,
calculated from the width and the amplitude of the Gaussians, and line flux
ratios are reported in Table 1. The uncertainty in the fluxes is calculated by
propagating the error of the fitted Gaussian parameters returned from
curve_fit. The uncertainties of the fluxes and ratios are calculated following
the standard error propagation.

We revisit a _Spitzer_ SH spectrum taken on 2008 August 16 and measure the [Ne
II] and [Ne III] lines (bottom panels of Figure 2, Table 1). Note that the
wavelength range of _Spitzer_ SH does not cover the [Ar II] line. We measure
neon line fluxes (Table 1) that are consistent with previous work (Espaillat
et al., 2013).

###  3.2 H$\alpha$ Emission Lines

In Figure 3 we show H$\alpha$ emission line profiles of SZ Cha obtained
between 2007 to 2023. There is significant variability in the H$\alpha$
emission line strength and in many epochs there is strong blueshifted
absorption, likely from a strong, moderate-velocity wind.

Here we model all of the H$\alpha$ profiles using the magnetospheric accretion
flow code from Hartmann et al. (1994); Muzerolle et al. (1998, 2001). These
models assume alignment of the magnetic, stellar rotation and disk rotation
axes. An axisymmetric accretion flow arises from the gas disk that deposits
material onto the star. The flow geometry follows that of a dipolar magnetic
field and is characterized by an inner radius (Ri) and the width of the flow
(Wr) at the disk plane. Additional parameters include the maximum temperature
in the flow (Tmax) and the given viewing inclination ($i$).

We created a grid of 133,848 models by varying $\dot{M}$ ($1-7\times 10^{-9}$
$M_{\sun}$ yr-1), Ri ($1.6-6.0$ R∗), Wr ($0.2-2.2$ R∗), Tmax ($8-11\times
10^{3}$ K), and $i$ ($20\arcdeg-70\arcdeg$). We convolved the model profiles
with a Gaussian instrumental profile of each instrument’s resolution. We
follow the procedure of Thanathibodee et al. (2023) to determine the best fits
and the average model parameters for each observation. Specifically, we
calculate the $\chi^{2}$ for each combination of the model and observed
profile. For each observed profile, we selected the models where the
normalized likelihood is $\geq 0.5$ and calculated the weighted mean of
$\dot{M}$, Ri, Wr, Tmax, and $i$. We show the best-fitting model to the
H$\alpha$ profiles in Figure 3 and the best fit parameters in Table 2. Winds
are not included in the model; the regions displaying blueshifted absorption
and low-velocity redshifted absorption in 2014 were excluded from the fit.

Table 2: Results of the Magnetospheric Accretion Model  Obs. Date | Instrument | $\dot{M}$ | Ri | Wr | Tmax | $i$
---|---|---|---|---|---|---
(UT) |  | ($10^{-9}$ $M_{\sun}$ yr-1) | (R⋆) | (R⋆) | ($10^{3}\,$K) | (deg)
2007-02-10 | MIKE | 3.8$\pm$1.8 | 2.3$\pm$0.5 | 1.3$\pm$0.5 | 10.6$\pm$0.3 | 65$\pm$ 4
2007-02-11 | MIKE | 4.7$\pm$1.5 | 2.4$\pm$0.5 | 0.7$\pm$0.4 | 9.8$\pm$0.8 | 53$\pm$ 8
2008-02-15 | MIKEaaObservation taken about 6 months before the Spitzer observation. | 3.3$\pm$1.8 | 2.1$\pm$1.3 | 0.5$\pm$0.5 | 9.1$\pm$0.8 | 54$\pm$16
2008-02-18 | MIKEaaObservation taken about 6 months before the Spitzer observation. | 3.2$\pm$1.8 | 2.1$\pm$1.2 | 0.5$\pm$0.5 | 9.3$\pm$0.8 | 51$\pm$16
2009-01-18 | MIKEbbObservation taken about 5 months after the Spitzer observation. | 4.4$\pm$1.8 | 2.0$\pm$0.4 | 0.8$\pm$0.4 | 10.1$\pm$0.7 | 58$\pm$10
2010-01-19 | X-Shooter | 4.8$\pm$1.6 | 2.4$\pm$0.5 | 0.9$\pm$0.6 | 10.1$\pm$0.7 | 57$\pm$11
2011-01-05 | MIKE | 3.5$\pm$1.8 | 3.1$\pm$1.5 | 1.2$\pm$0.5 | 9.6$\pm$0.8 | 39$\pm$17
2012-02-17 | MIKE | 4.3$\pm$1.8 | 1.8$\pm$0.6 | 0.4$\pm$0.3 | 9.3$\pm$0.8 | 58$\pm$11
2014-04-17 | UVES | 4.0$\pm$1.4 | 2.2$\pm$0.3 | 0.6$\pm$0.4 | 9.7$\pm$0.8 | 44$\pm$10
2014-05-31 | FEROS | 2.7$\pm$1.4 | 1.6$\pm$0.1 | 0.2$\pm$0.0 | 10.0$\pm$0.6 | 61$\pm$ 9
2023-04-22 | CHIRONccObservation taken about 15 hrs after the JWST observation. | 5.2$\pm$1.3 | 2.1$\pm$0.2 | 1.6$\pm$0.4 | 10.6$\pm$0.4 | 65$\pm$ 3

##  4 Discussion

We do not detect any significant molecular emission besides H2 from the pre-
transitional disk of SZ Cha. This is consistent with previous MIR molecular
line observations of transitional disks (i.e., a disk with a central cavity
devoid of dust and $\sim$ tens of astronomical units wide; Espaillat et al.,
2014) that find no significant molecular emission (Najita et al., 2010;
Banzatti et al., 2020); this is in contrast to observations of disks without
large cavities (i.e., “full” disks) that are rich in molecular emission (Carr
& Najita, 2008; Pontoppidan et al., 2010; Banzatti et al., 2023).

Pre-transitional and transitional disks are ideal candidates to search for [Ne
III] line emission. Only five protoplanetary disks have [Ne III] line
detections to date (Lahuis et al., 2007; Najita et al., 2010; Espaillat et
al., 2013). Of these, four were from pre-transitional and transitional disks.
As noted above, little molecular emission is found in pre-transitional and
transitional disks, which helps with the detection of the [Ne III] line since
there are H2O lines in the [Ne III] wavelength region. In addition, there is
less MIR continuum emission in pre-transitional and transitional disks, which
facilitates line detection. Also, pre-transitional and transitional disks are
not known to have strong outflows or jets as seen in full disks that can
contribute to neon emission (Güdel et al., 2010). High-resolution ground-based
MIR spectroscopy shows that the [Ne II] line emission comes from a region
close ($\sim$ tens of astronomical units) to the star for pre-transitional and
transitional disks (Sacco et al., 2012).

Aside from the Spitzer SH observation of SZ Cha, no evidence has been found to
support the possibility that EUV radiation dominates the creation of MIR
atomic emission lines. Szulágyi et al. (2012) have studied a sample of $\sim
60$ pre-transitional and transitional disks with Spitzer and find [Ne
III]-to-[Ne II] line flux ratios $<0.1$, which point to a soft EUV or X-ray
spectrum. Indirect measurements of the EUV luminosity in TTS using radio
emission suggest that the EUV emission is not sufficient to reproduce the
observed [Ne II] 12.81 µm luminosities (Pascucci et al., 2014, 2020).

To explain the large difference between the [Ne III]-to-[Ne II] line flux
ratio measured by Spitzer and JWST, we consider the possibility that the
Spitzer measurement of the [Ne III] line flux of SZ Cha is contaminated by
unresolved line emission. As a test, we convolved our _JWST_ spectrum down to
the resolution of Spitzer IRS SH and see that the [Ne III]-to-[Ne II] line
flux ratio does not change. Therefore, it is unlikely that line contamination
played a significant role in the [Ne III] line flux ratio measurement of SZ
Cha obtained with Spitzer. We then consider if at the time of the Spitzer
observations strong emission existed from the 13${}_{6~{}8}$ to 12${}_{3~{}9}$
transition of H2O, which is located at $\sim 15.57$ µm and may be weakly
present in the JWST spectrum (see top middle panel of Figure 2). However, this
is unlikely given that no water emission was seen in the Spitzer SH spectrum
at any other wavelengths.

Next we consider that the significant change in the [Ne III]-to-[Ne II] line
flux ratio of SZ Cha between 2008 and 2023 is due to a jet that was present at
the time of the Spitzer observations and that is no longer active at the time
of the JWST observations. Pascucci et al. (2014) did not find evidence of a
jet when studying centimeter emission of SZ Cha in 2012 October. The
variability of jet emission is not currently well understood (Espaillat et
al., 2019a, 2022; Curone et al., 2023), and while we cannot rule out this
scenario, we consider it to be unlikely.

To explain the variability in the [Ne III]-to-[Ne II] line flux ratio of SZ
Cha between 2008 and 2023, we propose a variable wind. A signature of a wind
is blueshifted absorption in the H$\alpha$ profile, and in Figure 3 it is
clear that SZ Cha’s wind is variable on short timescales (i.e., days). The
presence of strong blueshifted absorption is evident in the H$\alpha$ profile
of SZ Cha that was obtained about 15 hrs after the JWST observations (bottom
row, right panel, Figure 3). To the best of our knowledge, no H$\alpha$
observations of SZ Cha have been taken contemporaneously with Spitzer. There
are MIKE spectra taken about six months before (top row, two right panels in
Figure 3) and about five months after (second row, leftmost panel in Figure 3)
the Spitzer observations. We propose that at the time of the Spitzer
observations, the wind was suppressed, and that during the JWST observations,
a strong wind existed. A strong wind would absorb the EUV from the star and
prevent EUV emission from reaching the disk, but X-rays would permeate the
wind and reach the disk and generate the neon emission.

We can consider the [Ne II]-to-[Ar II] line flux ratio to further investigate
if a strong wind was present at the time of the JWST observations. Hollenbach
& Gorti (2009) find that [Ne II]/[Ar II]$\sim$1 in the case that EUV or soft
X-ray photons generate these emission lines and that this ratio is about 2.5
in the case of hard X-rays. The measured [Ne II]-to-[Ar II] line flux ratio of
1.9$\pm$0.6 is roughly compatible with either scenario within the
uncertainties. More JWST data are necessary to achieve a better signal-to-
noise ratio. Additional contemporaneous multiwavelength observations of SZ Cha
can test how the [Ne III]-to-[Ne II] line flux ratio varies with the wind.

##  5 Summary and Conclusions

We present a 2023 JWST MIRI MRS spectrum of the pre-transitional disk of SZ
Cha. We find no significant molecular emission, in line with previous Spitzer
IRS observations of this object specifically and disks with large gaps in
general. We detect atomic argon and neon emission, particularly the [Ar II]
6.99 µm, [Ne II] 12.81 µm, and [Ne III] 15.55 µm fine-structure emission
lines. The JWST [Ne III]-to-[Ne II] line flux ratio point to X-ray–dominated
photoevaporation of the disk.

We compare the neon line flux ratios from JWST in 2023 to those measured with
Spitzer in 2008 and find significant variability. With Spitzer, the [Ne
III]-to-[Ne II] line flux ratio pointed to EUV-dominated photoevaporation of
the disk. The change in the [Ne III]-to-[Ne II] line flux ratio between the
Spitzer and JWST observations indicates a switch between EUV and X-ray disk
photoevaporation. To explore the reason behind this switch, we analyze a
ground-based optical spectrum of the H$\alpha$ emission line taken within 15
hrs of the JWST observations. We find strong blueshifted absorption, evidence
of a wind. We propose that a strong wind was present at the time of the JWST
observations and that the EUV radiation was prevented from reaching the disk
while X-rays were able to permeate the wind and irradiate the disk. By
extension, we propose that a weak wind was present at the time of the Spitzer
observations and therefore EUV photons reached the disk and dominated the disk
photoevaporation.

This study shows that disk photoevaporation can be dominated by either EUV or
X-ray radiation and that this is impacted by the strength of winds that play a
role in dictating how much EUV reaches the disk. More observations of MIR
atomic lines in a larger sample of protoplanetary disks are necessary to
better understand how MIR atomic emission can be used to study disk
photoevaporation. In addition, more expansive multiepoch multiwavelength
contemporaneous data sets with JWST are critical to probe the star-disk
connection.

We thank Todd Henry, Wei-Chun Jao, and the CHIRON team for their prompt
scheduling of the CHIRON observations. We acknowledge support from JWST grant
GO-01676. M.K.M. acknowledges support from the Dutch Research Council (NWO;
grant VI.Veni.192.241). The JWST data presented in this paper were obtained
from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope
Science Institute. The specific observations analyzed can be accessed via
https://doi.org/10.17909/092y-8624 (catalog 10.17909/092y-8624).

## References

  * Alexander et al. (2006) Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006, MNRAS, 369, 229, doi: 10.1111/j.1365-2966.2006.10294.x
  * Asensio-Torres et al. (2021) Asensio-Torres, R., Henning, T., Cantalloube, F., et al. 2021, A&A, 652, A101, doi: 10.1051/0004-6361/202140325
  * Baldovin-Saavedra et al. (2011) Baldovin-Saavedra, C., Audard, M., Güdel, M., et al. 2011, A&A, 528, A22, doi: 10.1051/0004-6361/201015622
  * Banzatti et al. (2020) Banzatti, A., Pascucci, I., Bosman, A. D., et al. 2020, ApJ, 903, 124, doi: 10.3847/1538-4357/abbc1a
  * Banzatti et al. (2023) Banzatti, A., Pontoppidan, K. M., Carr, J., et al. 2023, arXiv e-prints, arXiv:2307.03846, doi: 10.48550/arXiv.2307.03846
  * Bernstein et al. (2003) Bernstein, R., Shectman, S. A., Gunnels, S. M., Mochnacki, S., & Athey, A. E. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4841, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, ed. M. Iye & A. F. M. Moorwood, 1694–1704, doi: 10.1117/12.461502
  * Bushouse et al. (2023) Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2023, JWST Calibration Pipeline, 1.11.4, Zenodo, doi: 10.5281/zenodo.8247246
  * Carr & Najita (2008) Carr, J. S., & Najita, J. R. 2008, Science, 319, 1504, doi: 10.1126/science.1153807
  * Clarke et al. (2001) Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485, doi: 10.1046/j.1365-8711.2001.04891.x
  * Curone et al. (2023) Curone, P., Testi, L., Macias, E., et al. 2023, arXiv e-prints, arXiv:2307.10798, doi: 10.48550/arXiv.2307.10798
  * Dekker et al. (2000) Dekker, H., D’Odorico, S., Kaufer, A., Delabre, B., & Kotzlowski, H. 2000, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4008, Optical and IR Telescope Instrumentation and Detectors, ed. M. Iye & A. F. Moorwood, 534–545, doi: 10.1117/12.395512
  * Ercolano et al. (2008) Ercolano, B., Drake, J. J., Raymond, J. C., & Clarke, C. C. 2008, ApJ, 688, 398, doi: 10.1086/590490
  * Ercolano & Owen (2010) Ercolano, B., & Owen, J. E. 2010, MNRAS, 406, 1553, doi: 10.1111/j.1365-2966.2010.16798.x
  * Espaillat et al. (2011) Espaillat, C., Furlan, E., D’Alessio, P., et al. 2011, ApJ, 728, 49, doi: 10.1088/0004-637X/728/1/49
  * Espaillat et al. (2007) Espaillat, C., Calvet, N., D’Alessio, P., et al. 2007, ApJ, 664, L111, doi: 10.1086/520879
  * Espaillat et al. (2013) Espaillat, C., Ingleby, L., Furlan, E., et al. 2013, ApJ, 762, 62, doi: 10.1088/0004-637X/762/1/62
  * Espaillat et al. (2014) Espaillat, C., Muzerolle, J., Najita, J., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 497, doi: 10.2458/azu_uapress_9780816531240-ch022
  * Espaillat et al. (2019a) Espaillat, C. C., Macías, E., Hernández, J., & Robinson, C. 2019a, ApJ, 877, L34, doi: 10.3847/2041-8213/ab2193
  * Espaillat et al. (2022) Espaillat, C. C., Macías, E., Wendeborn, J., et al. 2022, ApJ, 924, 104, doi: 10.3847/1538-4357/ac365a
  * Espaillat et al. (2019b) Espaillat, C. C., Robinson, C., Grant, S., & Reynolds, M. 2019b, ApJ, 876, 121, doi: 10.3847/1538-4357/ab16e6
  * Flaccomio et al. (2009) Flaccomio, E., Stelzer, B., Sciortino, S., et al. 2009, A&A, 505, 695, doi: 10.1051/0004-6361/200810972
  * Flaherty et al. (2012) Flaherty, K. M., Muzerolle, J., Rieke, G., et al. 2012, ApJ, 748, 71, doi: 10.1088/0004-637X/748/1/71
  * Font et al. (2004) Font, A. S., McCarthy, I. G., Johnstone, D., & Ballantyne, D. R. 2004, ApJ, 607, 890, doi: 10.1086/383518
  * Glassgold et al. (2007) Glassgold, A. E., Najita, J. R., & Igea, J. 2007, ApJ, 656, 515, doi: 10.1086/510013
  * Gorti & Hollenbach (2008) Gorti, U., & Hollenbach, D. 2008, ApJ, 683, 287, doi: 10.1086/589616
  * Gorti & Hollenbach (2009) —. 2009, ApJ, 690, 1539, doi: 10.1088/0004-637X/690/2/1539
  * Güdel et al. (2010) Güdel, M., Lahuis, F., Briggs, K. R., et al. 2010, A&A, 519, A113, doi: 10.1051/0004-6361/200913971
  * Hartmann et al. (1998) Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385, doi: 10.1086/305277
  * Hartmann et al. (1994) Hartmann, L., Hewett, R., & Calvet, N. 1994, ApJ, 426, 669, doi: 10.1086/174104
  * Hollenbach & Gorti (2009) Hollenbach, D., & Gorti, U. 2009, ApJ, 703, 1203, doi: 10.1088/0004-637X/703/2/1203
  * Ingleby et al. (2011) Ingleby, L., Calvet, N., Hernández, J., et al. 2011, AJ, 141, 127, doi: 10.1088/0004-6256/141/4/127
  * Kaufer et al. (1999) Kaufer, A., Stahl, O., Tubbesing, S., et al. 1999, The Messenger, 95, 8
  * Kim et al. (2009) Kim, K. H., Watson, D. M., Manoj, P., et al. 2009, ApJ, 700, 1017, doi: 10.1088/0004-637X/700/2/1017
  * Lahuis et al. (2007) Lahuis, F., van Dishoeck, E. F., Blake, G. A., et al. 2007, ApJ, 665, 492, doi: 10.1086/518931
  * Manara et al. (2014) Manara, C. F., Testi, L., Natta, A., et al. 2014, A&A, 568, A18, doi: 10.1051/0004-6361/201323318
  * Manzo-Martínez et al. (2020) Manzo-Martínez, E., Calvet, N., Hernández, J., et al. 2020, ApJ, 893, 56, doi: 10.3847/1538-4357/ab7ead
  * Meijerink et al. (2008) Meijerink, R., Glassgold, A. E., & Najita, J. R. 2008, ApJ, 676, 518, doi: 10.1086/527411
  * Muzerolle et al. (1998) Muzerolle, J., Calvet, N., & Hartmann, L. 1998, ApJ, 492, 743, doi: 10.1086/305069
  * Muzerolle et al. (2001) —. 2001, ApJ, 550, 944, doi: 10.1086/319779
  * Muzerolle et al. (2009) Muzerolle, J., Flaherty, K., Balog, Z., et al. 2009, ApJ, 704, L15, doi: 10.1088/0004-637X/704/1/L15
  * Najita et al. (2010) Najita, J. R., Carr, J. S., Strom, S. E., et al. 2010, ApJ, 712, 274, doi: 10.1088/0004-637X/712/1/274
  * Owen et al. (2011) Owen, J. E., Ercolano, B., & Clarke, C. J. 2011, MNRAS, 412, 13, doi: 10.1111/j.1365-2966.2010.17818.x
  * Pascucci et al. (2022) Pascucci, I., Cabrit, S., Edwards, S., et al. 2022, arXiv e-prints, arXiv:2203.10068, doi: 10.48550/arXiv.2203.10068
  * Pascucci et al. (2014) Pascucci, I., Ricci, L., Gorti, U., et al. 2014, ApJ, 795, 1, doi: 10.1088/0004-637X/795/1/1
  * Pascucci et al. (2007) Pascucci, I., Hollenbach, D., Najita, J., et al. 2007, ApJ, 663, 383, doi: 10.1086/518535
  * Pascucci et al. (2011) Pascucci, I., Sterzik, M., Alexander, R. D., et al. 2011, ApJ, 736, 13, doi: 10.1088/0004-637X/736/1/13
  * Pascucci et al. (2020) Pascucci, I., Banzatti, A., Gorti, U., et al. 2020, ApJ, 903, 78, doi: 10.3847/1538-4357/abba3c
  * Pontoppidan et al. (2010) Pontoppidan, K. M., Salyk, C., Blake, G. A., et al. 2010, ApJ, 720, 887, doi: 10.1088/0004-637X/720/1/887
  * Ribas et al. (2016) Ribas, Á., Bouy, H., Merín, B., et al. 2016, MNRAS, 458, 1029, doi: 10.1093/mnras/stw333
  * Rieke et al. (2015) Rieke, G. H., Wright, G. S., Böker, T., et al. 2015, PASP, 127, 584, doi: 10.1086/682252
  * Rydgren (1980) Rydgren, A. E. 1980, AJ, 85, 444, doi: 10.1086/112694
  * Sacco et al. (2012) Sacco, G. G., Flaccomio, E., Pascucci, I., et al. 2012, ApJ, 747, 142, doi: 10.1088/0004-637X/747/2/142
  * Schisano et al. (2010) Schisano, E., Ercolano, B., & Güdel, M. 2010, MNRAS, 401, 1636, doi: 10.1111/j.1365-2966.2009.15799.x
  * Schwab et al. (2010) Schwab, C., Spronck, J. F. P., Tokovinin, A., & Fischer, D. A. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7735, Proc. SPIE, 77354G, doi: 10.1117/12.856709
  * Szulágyi et al. (2012) Szulágyi, J., Pascucci, I., Ábrahám, P., et al. 2012, ApJ, 759, 47, doi: 10.1088/0004-637X/759/1/47
  * Thanathibodee et al. (2023) Thanathibodee, T., Molina, B., Serna, J., et al. 2023, ApJ, 944, 90, doi: 10.3847/1538-4357/acac84
  * Tody (1993) Tody, D. 1993, in Astronomical Society of the Pacific Conference Series, Vol. 52, Astronomical Data Analysis Software and Systems II, ed. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 173
  * Tokovinin et al. (2013) Tokovinin, A., Fischer, D. A., Bonati, M., et al. 2013, PASP, 125, 1336, doi: 10.1086/674012
  * Vernet et al. (2011) Vernet, J., Dekker, H., D’Odorico, S., et al. 2011, A&A, 536, A105, doi: 10.1051/0004-6361/201117752
  * Wells et al. (2015) Wells, M., Pel, J. W., Glasse, A., et al. 2015, PASP, 127, 646, doi: 10.1086/682281
  * Wright et al. (2023) Wright, G. S., Rieke, G. H., Glasse, A., et al. 2023, PASP, 135, 048003, doi: 10.1088/1538-3873/acbe66


# A Pedagogical Review of the Vacuum Retarded Dipole Model of Pulsar Spin Down

J. C. Satherley  School of Physical and Chemical Sciences, University of
Canterbury, Christchurch, New Zealand  [ C. Gordon  School of Physical and
Chemical Sciences, University of Canterbury, Christchurch, New Zealand

(dd Mmm YYYY; dd Mmm YYYY; dd Mmm YYYY)

###### Abstract

Pulsars are rapidly spinning highly magnetised neutron stars. Their spin
period is observed to decrease with time. An early analytical model for this
process was the vacuum retarded dipole (VRD) by Deutsch (1955). This model
assumes an idealised star and it finds that the rotational energy is radiated
away by the electromagnetic fields. This model has been superseded by more
realistic numerical simulations that account for the non-vacuum like
surroundings of the neutron star. However, the VRD still provides a reasonable
approximation and is a useful limiting case that can provide some qualitative
understanding. We provide detailed derivations of the spin down and related
electromagnetic field equations of the VRD solution. We also correct
typographical errors in the general field equations and boundary conditions
used by Deutsch (1955).

###### doi:

%

###### keywords:

stars: pulsars – stars: magnetic field – stars: neutron

J. C<EMAIL_ADDRESS>19 March 2021

##  1 INTRODUCTION

Pulsars are highly magnetised rapidly rotating neutron stars that emit
electromagnetic radiation almost across the entire spectrum, from radio to
$\gamma$-rays (Travelle, 2011). Because of the extreme magnetic and electric
fields, the environment around a pulsar proves to be interesting and is
believed to be the main cause of a pulsar’s energy emission. This energy
emission results in a decreasing spin-rate for the neutron star. An analytical
approximation for this process is given by the Vacuum Retarded Dipole (VRD)
model which was originally derived by Deutsch (1955) (hereafter referred to as
D55).

The derivations of the results in D55 are very concise, with many details
missing. The purpose of this review is to rederive the main results of D55 and
in doing so provide details on how this derivation takes place mathematically
and physically. We also uncover three typographical errors in D55, the most
important of which is in our Eqn. (31), the polar component of the electric
field. It is worth noting the first papers we found that mentioned the error
in the general field equations are Melatos (1997) and Michel & Li (1999).
However, they do not show the derivations for the general field equations and
how the typographical error is found. A reasonably detailed derivation of the
VRD solution, without the typos, is shown in the textbook by Michel (1991).
However, the method used by Michel is different from D55. While D55 uses
methods outlined by Stratton (1941), Michel just considers multipole
expansions at an inclination and does not derive the normalisations. Also the
following working carried out by us goes into greater detail for the
derivations. In this review we use methods outlined in Stratton (1941)
(hereafter referred to as S41). In particular, we refer frequently to Chapters
VII and VIII, which consider spherical waves and radiation respectively. We
use Griffiths (2017) (hereafter referred to as G17) as a reference for
standard electrodynamics theory.

In Sec. 2, we begin by outlying the definition of an idealised star, along
with fundamental equations that apply due to this assumption, which include
Maxwell’s equations. We also note down the Alfvén theorem and a brief
description of a vector field. In Sec. 3, we define the _magnetic field
symmetry axis_ and its associated coordinate system. The magnetic field is
assumed to take the form of a dipole about this axis. The transformation
between the two coordinate systems is shown and used to derive the
electromagnetic field boundary conditions in the rotation symmetry axis
coordinates. In Sec. 4, we show the derivations for the time-dependent and
static components of the electromagnetic fields using series expansions and
matching coefficients to the boundary conditions. In Sec. 5, we simplify the
form of the electromagnetic fields for a large radius, then the power and spin
down equations are derived.

##  2 The Internal Fields of An Idealised Star

For the VRD solution, we consider the electric and magnetic fields of an
idealised star. An idealised star is a star that is a sharply bounded sphere,
is isotropic, is perfectly conducting, and has permittivity and permeability
of a vacuum everywhere inside. The boundary is not continuous but rather an
immediate change between the star’s medium and the surrounding vacuum. Which,
in the case of a real star is not true. However, it greatly simplifies the
situation.

###  2.1 Maxwell’s Equation

The Maxwell equations in SI units are,

$\displaystyle\nabla\cdot\bm{E}=\frac{\rho}{\epsilon_{0}},$ (1)
$\displaystyle\nabla\cdot\bm{H}=0,$ (2)
$\displaystyle\nabla\times\bm{E}=-\mu_{0}\frac{\partial\bm{H}}{\partial t},$
(3)
$\displaystyle\nabla\times\bm{H}=\bm{J}+\epsilon_{0}\frac{\partial\bm{E}}{\partial
t},$ (4)

where $\bm{E}$ is the electric field, $\bm{H}$ is the magnetic field (related
to $\bm{B}$ by $\bm{B}=\mu_{0}\bm{H}$), $\epsilon_{0}$ is the permittivity of
free space, and $\mu_{0}$ the permeability of free space. $\bm{H}$ will
typically be referred to as the $H$ field. In D55, there is a missing minus
sign in his equivalent of equation (3). This is clearly a typo as he takes
into account the minus sign in later parts of his article.

###  2.2 The H field

Firstly, D55 describes the typical spherical polar coordinates around the
star’s axis of rotation, $(r,\theta,\varphi)$, with rotation axis
$\bm{\omega}$ aligned with the $z$-axis and origin at the centre of the star.
Where $r$ is the radius, $\theta$ is the angle from the pole and $\varphi$ is
the azimuth angle measured counterclockwise from the $x$-axis.

Let $\bm{H}(r,\theta,\varphi,t)$ be an arbitrary solenoidal function within
the star. A solenoidal function is one for which the divergence is zero at all
points in the field ($\nabla\cdot\bm{H}=0$). D55 states that a necessary and
sufficient condition that $\bm{H}$ be “frozen into” the star is that

$\bm{H}(r,\theta,\varphi,t)=r_{0}H_{r}(r,\theta,\lambda)+\theta_{0}H_{\theta}(r,\theta,\lambda)+\varphi_{0}H_{\varphi}(r,\theta,\lambda),$

(5)

where $\lambda$ is an azimuthal coordinate measured from a meridian fixed in
the star, and $r_{0}$, $\theta_{0}$ and $\phi_{0}$ are the unit vectors in the
respective coordinate directions. That is the time and azimuthal angle
$\varphi$ can be replaced by $\lambda$ as the $H$ field and the meridian move
with the star. The variable $\lambda$ is defined more clearly in eqn. (6) and
Fig. 2(b). $H_{r}$, $H_{\theta}$ and $H_{\varphi}$ give the magnitude of the
field in each basis vector’s direction. The concept of the magnetic field
being “frozen into” a medium comes from Alfvén’s Theorem (Roberts, 2007).

####  2.2.1 Vector Fields

As $\bm{H}$ is a vector field, at a particular point $(r,\theta,\varphi)$ in
space, the field can be represented as a vector with some magnitude in each
unit vector’s direction. In D55, the directions are given by the subscript $0$
terms, and the magnitudes are given by the functions proceeding the unit
vector terms. In polar coordinates, $r_{0}$ is the radial component, pointing
away from the centre of the sphere. The polar component is $\theta_{0}$, its
direction is tangent to the polar angle. Lastly, $\varphi_{0}$ is the
azimuthal component which has a direction tangent to the azimuthal angle. This
notation continues throughout D55. However, we will here on forth elect to use
$\,\mathbf{\hat{r}}$, $\,\bm{\hat{\theta}}$ and $\,\bm{\hat{\varphi}}$ to
represent these components of the vector field. A visual representation of the
unit vectors of a field in polar coordinates is shown in Fig. 1(a). Hence, (5)
in the notation used for the rest of the report becomes,

$\bm{H}(r,\theta,\varphi,t)=H_{r}(r,\theta,\lambda)\,\mathbf{\hat{r}}+H_{\theta}(r,\theta,\lambda)\,\bm{\hat{\theta}}+H_{\varphi}(r,\theta,\lambda)\,\bm{\hat{\varphi}}.$

$z,\bm{\omega}$$\,\mathbf{\hat{r}}$$\,\bm{\hat{\theta}}$$\,\bm{\hat{\varphi}}$
(a)

$z,\bm{\omega}$$\bm{e}$$\,\mathbf{\hat{r}}$$\,\bm{\hat{\psi}}$$\,\bm{\hat{\nu}}$
(b)

Figure 1: (a) The unit vectors of a vector field at a point on a sphere in
$(r,\theta,\varphi)$. (b) The unit vectors of a vector field at a point on a
sphere in $(r,\psi,\nu)$.

##  3 The External Field in Vacuum

It is convenient to define a new set of polar coordinates $(r,\psi,\nu)$ with
axis $\bm{e}$ such that the magnetic field is symmetrical about $\bm{e}$
(known as the _magnetic field symmetry axis_). The axis $\bm{e}$ will be in
fixed rotation with the star due to Alfvén’s theorem. This new coordinate
system is shown in Fig. 2(a). The magnetic field symmetry here refers to, for
a fixed $r$ and $\psi$, the magnitude of the field does not change as $\nu$ is
changed and the direction does not change relative to the local unit vectors.
The angle between the rotation axis $\bm{\omega}$ and $\bm{e}$ is $\chi$. This
relationship is shown in Fig. 2(b). The axis $\bm{e}$ has polar angle $\chi$
and azimuthal angle $\omega t$ in $(r,\theta,\varphi)$, $\bm{e}=(r,\chi,\omega
t)$, at some time $t$ during the rotation of the star. The azimuthal angle
$\lambda$ mentioned earlier is measured from $\bm{e}$ along $\varphi$, so that

$\lambda=\varphi-\omega t.$ (6)

$x$$y$$z$$\bm{e}$$\psi$$\nu$ (a)

$x$$y$$z$$\bm{\omega}$$\bm{e}$$\chi$$\omega t$$\lambda$ (b)

Figure 2: The relation between the two coordinate systems $(r,\theta,\varphi)$
and $(r,\psi,\nu)$. $x$, $y$, and $z$ are the Cartesian coordinates. (a) The
coordinate system $(r,\psi,\nu)$ portrayed on a sphere. (b) The star’s
rotation axis $\bm{\omega}$ and the magnetic field’s symmetry axis $\bm{e}$.
$\chi$ is the angle of inclination between $\bm{\omega}$ and $\bm{e}$, and
$\lambda$ is measured from the projection of $\bm{e}$ onto the xy-plane.

We shall now consider a magnetic field which within the star has the form,

$\bm{H}=R_{1}(r)S_{1}(\psi)\,\mathbf{\hat{r}}+R_{2}(r)S_{2}(\psi)\,\bm{\hat{\psi}}+R_{3}(r)S_{3}(\psi)\,\bm{\hat{\nu}},$
(7)

where $R_{i}$ and $S_{i}$ are arbitrary functions of $r$ and $\psi$
respectively, with $i=1,2$, or $3$. Note that $\bm{H}$ has no dependence on
$\nu$ for the magnitude of the vectors and hence the magnitude is symmetrical
about axis $\bm{e}$ as stated earlier. In this situation D55 assumes the
interior magnetic field to be a _dipole_ about $\bm{e}$ by taking
$S_{1}(\psi)=\cos\psi$.

###  3.1 Matrix Transformation between Vector Field Basis Vectors

A vector transformation between the two coordinate systems mentioned is
required whilst working out the boundary conditions for the fields in both
coordinate systems. It is convenient to figure out the transformation now,
then approach the boundary conditions. Define the transformation matrix
$R:(\,\mathbf{\hat{r}},\,\bm{\hat{\theta}},\,\bm{\hat{\varphi}})\to(\,\mathbf{\hat{r}},\,\bm{\hat{\psi}},\,\bm{\hat{\nu}})$
which rotates between the two unit vector systems at a point in space. Figs.
1(a) and 1(b) show the two sets of basis vectors at the same point. It is
clear that $\,\mathbf{\hat{r}}$ points radially out, normal to the sphere’s
surface, in both figures. Hence, at the same point in space
$\,\mathbf{\hat{r}}$ is the same in both coordinate systems. However, the
angular components undergo a rotation as
$\,\bm{\hat{\theta}},\,\bm{\hat{\varphi}},\,\bm{\hat{\psi}}$ and
$\,\bm{\hat{\nu}}$ all lie in a plane tangent to the sphere. The rotation
between these unit vectors will be a two dimensional rotation about
$\,\mathbf{\hat{r}}$. Let the angle of this rotation be $\xi$, the
transformation matrix $R$ is then

$R=\begin{bmatrix}1&0&0\\\ 0&\cos\xi&-\sin\xi\\\
0&\sin\xi&\cos\xi\end{bmatrix}.$ (8)

This matrix transformation takes a basis of vectors at some point in space
from $(\,\mathbf{\hat{r}},\,\bm{\hat{\theta}},\,\bm{\hat{\varphi}})$ and
rotates them anticlockwise by an angle $\xi$ with $\,\mathbf{\hat{r}}$ as axis
to the unit vectors $(\,\mathbf{\hat{r}},\,\bm{\hat{\psi}},\,\bm{\hat{\nu}})$.
To reverse this transformation the angle of rotation is replaced by $-\xi$.

The angle $\xi$ can be found in terms of already defined coordinate parameters
and values. This is done by considering triangles on the surface of a sphere.
The spherical triangle for this situation can be found in Fig. 4(b) in F. Then
from solutions of spherical triangles (53) and (54) the following
relationships can be found,

$\displaystyle\sin\xi$ $\displaystyle=\frac{\sin\chi\sin\lambda}{\sin\psi},$
(9) $\displaystyle\cos\xi$
$\displaystyle=\frac{\cos\chi-\cos\psi\cos\theta}{\sin\psi\sin\theta}.$

An interesting relationship that is obtained from the standard spherical
triangle relation (52) is

$\cos\psi=\cos\chi\cos\theta+\sin\chi\sin\theta\cos\lambda,$ (10)

which will be required momentarily.

###  3.2 Boundary Conditions of the Fields

As the surface of the star is a rigid boundary between the interior of the
star and free space, the necessary components of the external electromagnetic
fields must satisfy the boundary conditions at $r=a$ (see G17). These boundary
conditions are set by the internal electromagnetic fields. From equations
(7.61) and (7.62) of G17, the magnetic fields orthogonal to the surface must
be equal and the electric fields tangent to the surface must be equal. Hence,
using (7), the first condition at the surface of the star that must be
satisfied is,

$H_{r}^{\text{ext}}=H_{r}^{\text{int}}=R_{1}(a)S_{1}(\psi).$

Next we look at the internal electric field. To do this we use Ohm’s law from
G17, Eqn. 7.2, that relates the electric field and current density, $\bm{J}$
and is given by,

$\bm{E}+\mu_{0}\bm{V}\times\bm{H}=\frac{\bm{J}}{\sigma}$ (11)

where $\sigma$ is the conductivity and $\bm{V}$ is the local velocity of
rotation of the star. Assuming $\sigma\rightarrow\infty$ gives

$\bm{E}=-\mu_{0}\bm{V}\times\bm{H}.$ (12)

The velocity in
$(\,\mathbf{\hat{r}},\,\bm{\hat{\theta}},\,\bm{\hat{\varphi}})$ is

$\bm{V}=\bm{\omega}\times\bm{r}=r\omega\sin\theta\,\bm{\hat{\varphi}}.$ (13)

However, the vector field (7) must be converted from
$(\,\mathbf{\hat{r}},\,\bm{\hat{\psi}},\,\bm{\hat{\nu}})$ to
$(\,\mathbf{\hat{r}},\,\bm{\hat{\theta}},\,\bm{\hat{\varphi}})$ so that
$\bm{V}$ and $\bm{H}$ are defined on the same basis. This is done by rotating
(7) using (8). However, as the reverse transformation is needed, the signs of
$R_{2,3}$ and $R_{3,2}$ elements are switched.

$\displaystyle R(-\xi)\bm{H}=$ $\displaystyle
R_{1}(r)S_{1}(\psi)\,\mathbf{\hat{r}}$
$\displaystyle+\left[R_{2}(r)S_{2}(\psi)\cos\xi+R_{3}(r)S_{3}(\psi)\sin\xi\right]\,\bm{\hat{\theta}}$
$\displaystyle+\left[-R_{2}(r)S_{2}(\psi)\sin\xi+R_{3}(r)S_{3}(\psi)\cos\xi\right]\,\bm{\hat{\varphi}}.$

Applying the previous result to (12) with (13) and using Fig. 4(b) of F we can
find the $E$ field inside the star,

$\displaystyle\bm{E}=$ $\displaystyle-\mu_{0}\bm{V}\times(R(-\xi)\bm{H})$
$\displaystyle=$
$\displaystyle\mu_{0}r\omega\sin\theta\left(R_{2}(r)S_{2}(\psi)\cos\xi+R_{3}(r)S_{3}(\psi)\sin\xi\right)\,\mathbf{\hat{r}}$
$\displaystyle-\omega\mu_{0}rR_{1}(r)S_{1}(\psi)\sin\theta\,\bm{\hat{\theta}}$
$\displaystyle+0\,\bm{\hat{\varphi}},$

and evaluating this at $r=a$ gives the tangential components as

$\displaystyle E_{\theta}$
$\displaystyle=-\omega\mu_{0}aR_{1}(a)S_{1}(\psi)\sin\theta,$ $\displaystyle
E_{\varphi}$ $\displaystyle=0.$

All but the third boundary condition match what D55 obtained. There is a
typographical error in D55’s boundary conditions (Deutsch’s Eqn. (10)), the
third condition should be written as
$E_{\varphi}^{\text{ext}}=E_{\varphi}^{\text{int}}=0$ rather than
$E_{r}^{\text{ext}}=E_{\varphi}^{\text{int}}=0$.

A more useful form of the boundary conditions can be found. In particular, the
field equations are solved for when $S_{1}=\cos\psi$. So by recalling (10),
the boundary conditions then results in,

$\begin{aligned}
H_{r}^{\text{ext}}=H_{r}^{\text{int}}&=R_{1}(a)(\cos\chi\cos\theta+\sin\chi\sin\theta\cos\lambda),\\\
E_{\theta}^{\text{ext}}=E_{\theta}^{\text{int}}&=-\omega\mu_{0}aR_{1}(a)(\cos\chi\cos\theta+\sin\chi\sin\theta\cos\lambda)\sin\theta\\\
&=-\frac{1}{2}\omega\mu_{0}aR_{1}(a)(\sin{2\theta}\cos\chi+(1-\cos{2\theta})\sin\chi\cos\lambda),\\\
E_{\varphi}^{\text{ext}}=E_{\varphi}^{\text{int}}&=0.\end{aligned}$

(14)

From the boundary conditions in the $(r,\theta,\phi)$ coordinate system, there
is a static (time-independent) component and a time-dependent component for
the fields. The terms in (14) which relate to the time-dependent component
have a dependency on $\lambda$, which depends on time from (6).

##  4 Derivation of General Field Equations

The general field equations for the VRD solution external to the star are
provided in the Appendix of D55. However, with minimal explanation given in
D55 about the derivation of the field equations, we wish to show how they are
found. We use results found in Chapters VII and VIII of S41 for this section.
Chapter VII proves beneficial in supplying general solutions to the vector
wave equation which, must be satisfied by the time evolving electromagnetic
fields around the neutron star. We apply these solutions to the neutron star
situation by firstly considering the propagation factor of the waves outside
the star. In S41, on page 392, the propagation factor (complex wavenumber) is
given as,

$k^{2}=\epsilon\mu\omega^{2}+i\sigma\mu\omega,$ (15)

and as we are now considering the fields external of the star, in a vacuum,
$\epsilon=\epsilon_{0}$, $\mu=\mu_{0}$, $\sigma=0$ and
$c=1/\sqrt{\mu_{0}\epsilon_{0}}$. Hence,

$k=\omega/c.$ (16)

Then for convenience define,

$\begin{gathered}\rho=kr=(\omega/c)r,\\\ \alpha=ka=(\omega/c)a,\end{gathered}$
(17)

where $\omega$ is the rotational velocity of the star, $c$ the speed of light,
and $a$ is the radius of the star. Here $\rho$ and $\alpha$ are unitless
numbers. A trivial relationship that will be required is

$\frac{\partial}{\partial\rho}=\frac{\omega}{c}\frac{\partial}{\partial r}.$

Next consider the EM fields surrounding the star in two parts, a static time-
independent component and a time-dependent component, where the complete field
will be the linear combination of the two components (G17). The time-
independent component will be axially symmetric about the axis $\bm{\omega}$
and hence have no $\varphi$ dependence ($m=0$, $m$ will be defined clearer
shortly). Whilst the time-dependent component will be non-axially symmetric
and hence have a dependence on $\varphi$ ($m\neq 0$). This can be expressed as

$\begin{gathered}\bm{H}=\bm{H}_{\text{static}}(r,\theta)+\bm{H}_{\text{time-
dependent}}(r,\theta,\lambda),\\\
\bm{E}=\bm{E}_{\text{static}}(r,\theta)+\bm{E}_{\text{time-
dependent}}(r,\theta,\lambda).\end{gathered}$

We first consider the time-dependent component in the following section on the
wave solution. Then use a multipole expansion to provide the static component.

###  4.1 Wave Solution

Due to the rotating nature of the magnetic field symmetry axis around the
star’s axis of rotation, the vector fields $\bm{E}$ and $\bm{H}$ satisfy the
same vector differential wave equation from S41, Sec. 7.1, Eqn. (1), which is
given as

$\nabla^{2}\bm{C}-\mu\epsilon\frac{\partial^{2}\bm{C}}{\partial
t^{2}}-\mu\sigma\frac{\partial\bm{C}}{\partial t}=0,$ (18)

where $\bm{C}$ represents either vector $\bm{E}$ or $\bm{H}$. We can express
$\bm{E}$ and $\bm{H}$ as a series of solutions $\bm{\mathrm{M}}_{n}$ and
$\bm{\mathrm{N}}_{n}$ which satisfy the vector differential wave equation (18)
for $\bm{C}$. From S41, Sec. 7.1, Eqn. (12) this representation is,

$\begin{gathered}\bm{E}=-\sum_{n}\left(a_{n}\bm{\mathrm{M}}_{n}+b_{n}\bm{\mathrm{N}}_{n}\right),\\\
\bm{H}=-\frac{k}{i\omega\mu_{0}}\sum_{n}\left(a_{n}\bm{\mathrm{N}}_{n}+b_{n}\bm{\mathrm{M}}_{n}\right)\end{gathered}$
(19)

where $\bm{\mathrm{M}}_{n}$ and $\bm{\mathrm{N}}_{n}$ have more recently been
referred to as the Magnetic and Electric harmonics respectfully. In
particular, the functions are expressed as
$\bm{\mathrm{M}}_{n}=\nabla\times\bm{a}\psi$ and
$\bm{\mathrm{N}}_{n}=(1/k)\nabla\times\bm{\mathrm{M}}_{n}$, where $\bm{a}$ is
a constant vector and $\psi$ is a scalar function satisfying
$\nabla^{2}\psi+k^{2}\psi=0$, which is known as the Helmholtz equation. The
derivation for the general solutions of $\bm{\mathrm{M}}_{n}$ and
$\bm{\mathrm{N}}_{n}$ is provided in S41 and the results are quoted below.
These solutions will also contain a time dependence. Due to the linearity of
the wave equation, the time dependency can be split off such that
$\bm{\mathrm{M}}=\bm{\mathrm{m}}e^{-i\omega t}$ and
$\bm{\mathrm{N}}=\bm{\mathrm{n}}e^{-i\omega t}$ without any loss of
generality. Then from S41, Sec. 7.11, Eqns. (11) and (12), and rewriting them
in a slightly more useful form with $\rho$ instead of $kr$,

$\displaystyle\bm{\mathrm{m}}_{\genfrac{}{}{0.0pt}{}{e}{o}{m\,n}}=$
$\displaystyle\mp\frac{m}{\sin\theta}h^{(1)}_{n}(\rho)P_{n}^{m}(\cos\theta)\genfrac{}{}{0.0pt}{}{\sin}{\cos}(m\varphi)\,\bm{\hat{\theta}}$
(20) $\displaystyle-h^{(1)}_{n}(\rho)\frac{\partial
P_{n}^{m}}{\partial\theta}\genfrac{}{}{0.0pt}{}{\cos}{\sin}(m\varphi)\,\bm{\hat{\varphi}},$
$\displaystyle\bm{\mathrm{n}}_{\genfrac{}{}{0.0pt}{}{e}{o}{m\,n}}=$
$\displaystyle\frac{n(n+1)}{\rho}h^{(1)}_{n}(\rho)P_{n}^{m}(\cos\theta)\genfrac{}{}{0.0pt}{}{\cos}{\sin}(m\varphi)\,\mathbf{\hat{r}}$
(21)
$\displaystyle+\frac{1}{\rho}\frac{\partial}{\partial\rho}\left[\rho\,h^{(1)}_{n}(\rho)\right]\frac{\partial}{\partial\theta}P_{n}^{m}(\cos\theta)\genfrac{}{}{0.0pt}{}{\cos}{\sin}(m\varphi)\,\bm{\hat{\theta}}$
$\displaystyle\mp\frac{m}{\rho\sin\theta}\frac{\partial}{\partial\rho}\left[\rho\,h^{(1)}_{n}(\rho)\right]P_{n}^{m}(\cos\theta)\genfrac{}{}{0.0pt}{}{\sin}{\cos}(m\varphi)\,\bm{\hat{\varphi}},$

where $h^{(1)}_{n}$ is the spherical Bessel function of the third kind (refer
to C), $P_{n}^{m}(\cos\theta)$ are the associated Legendre polynomials (refer
to D) , and $m$ and $n$ are associated to a particular solution to the
Helmholtz equation. In particular, $m$ defines the azimuthal dependence of the
$E$ and $H$ fields. Whilst $n$ controls the polar dependence. These functions
(20) and (21) are now known as vector spherical harmonics. The subscript $e$
or $o$ define the even and odd solution, this chooses the sign of some terms
and whether the trigonometric function for $\varphi$ is $\sin$ or $\cos$. From
here on $h^{(1)}_{n}$ will be referred to as $h_{n}$. However, D55 does not
use the even-odd notation from S41 but rather D55 multiplies the function by a
complex exponential then takes the real part of the result to give the wave
solutions. This can be expressed as an operator that first multiplies a
function $f(x)$ by $e^{im\varphi}$, then takes the real part of the result,

$\displaystyle f(x)\genfrac{}{}{0.0pt}{}{\cos}{\sin}(m\varphi)$
$\displaystyle=\Re(f(x)e^{im\varphi}),$ (22) $\displaystyle\mp
f(x)\genfrac{}{}{0.0pt}{}{\sin}{\cos}(m\varphi)$
$\displaystyle=-\Re(f(x)ie^{im\varphi}),$

where $f(x)$ may be complex-valued function and $\Re$ denotes the operator
which returns the real part of an expression. The operation (22) essentially
selects $\cos(m\varphi)$ or $\sin(m\varphi)$ depending on the complex-valued
nature of $f(x)$.

We wish to solve for the field equations in polar coordinates around rotation
axis $\bm{\omega}$. Due to the dependence of $\varphi$ in our electromagnetic
boundary conditions (14) only being $\cos(\varphi)$ and $\sin(\varphi)$, we
are limited to $m=1$ for our equations (20) and (21). Otherwise terms which
contain $\cos(m\varphi)$ and $\sin(m\varphi)$, where $m\neq 1$, will appear in
the solutions of (19) and cause disagreement of the external and internal
electromagnetic fields at the surface of the star. By also comparing the
series expression (19) and the boundary conditions (14) we see that if a
trigonometric function of $\theta$ does not appear in the boundary conditions,
its weighting coefficient will be zero (as will be noticed in the derivations
shortly). Hence, it is apparent from inspection of the boundary conditions and
the Associated Legendre Polynomials that only $n=1$ and $n=2$ are required. So
evaluating (20) and (21) for $m=1$, $n=1$ and $n=2$ and including the new
notation from (22) gives,

$\displaystyle\bm{\mathrm{m}}_{1,1}=$ $\displaystyle
ih_{1}(\rho)e^{i\varphi}\,\bm{\hat{\theta}}-h_{1}(\rho)\cos\theta
e^{i\varphi}\,\bm{\hat{\varphi}},$ (23) $\displaystyle\bm{\mathrm{m}}_{1,2}=$
$\displaystyle 3ih_{2}(\rho)\cos\theta
e^{i\varphi}\,\bm{\hat{\theta}}-3h_{2}(\rho)\cos{2\theta}e^{i\varphi}\,\bm{\hat{\varphi}},$
$\displaystyle\bm{\mathrm{n}}_{1,1}=$ $\displaystyle
2\frac{h_{1}(\rho)}{\rho}\sin\theta
e^{i\varphi}\,\mathbf{\hat{r}}+\frac{1}{\rho}(\rho
h_{1}^{\prime}(\rho)+h_{1}(\rho))\cos\theta e^{i\varphi}\,\bm{\hat{\theta}}$
$\displaystyle+\frac{i}{\rho}(\rho
h_{1}^{\prime}(\rho)+h_{1}(\rho))e^{i\varphi}\,\bm{\hat{\varphi}},$
$\displaystyle\bm{\mathrm{n}}_{1,2}=$ $\displaystyle
9\frac{h_{2}(\rho)}{\rho}\sin{2\theta}e^{i\varphi}\,\mathbf{\hat{r}}+\frac{3}{\rho}(\rho
h_{2}^{\prime}(\rho)+h_{2}(\rho))\cos{2\theta}e^{i\varphi}\,\bm{\hat{\theta}}$
$\displaystyle+\frac{3i}{\rho}(\rho
h_{2}^{\prime}(\rho)+h_{2}(\rho))\cos\theta e^{i\varphi}\,\bm{\hat{\varphi}},$

where we have neglected the operator that extracts the real part of the
result. This will be added during the calculation for the coefficients in
(19).

###  4.2 Multipole Expansion

A static field symmetric about the rotation axis can be expressed using a
multipole expansion. Upon inspection of the static terms in the boundary
conditions (14) and the multipole expansions given in G17 it can be seen, with
some working, that only a magnetic dipole and electric quadrupole are present
in this situation based on the conditions that multipole terms are linearly
independent of each other and that the boundary conditions will uniquely
determine the fields (electromagnetism uniqueness theorem, see G17).

####  4.2.1 Magnetic Dipole

We can begin by taking a magnetic dipole and check whether it satisfies the
$H$ field boundary condition from (14). From equation (5.88) of G17, a
magnetic dipole is given as,

$\bm{H}_{\text{dip}}(\bm{r})=\frac{m}{4\pi
r^{3}}\left(2\cos\theta\,\mathbf{\hat{r}}+\sin\theta\,\bm{\hat{\theta}}\right),$
(24)

where $m$ is the magnetic dipole moment of the field. Now compare the radial
component of (24) to the static term of the radial $H$ field boundary
condition from (14).

$\frac{m}{4\pi a^{3}}2\cos\theta=R_{1}(a)\cos\chi\cos\theta.$

We see that the function depending on $\theta$ in the boundary condition is
satisfied fully by the dipole term. Hence, a magnetic dipole is the only
static term present due to the boundary conditions. This implies the ratio
$m/(4\pi)$ will satisfy,

$\frac{m}{4\pi}=\frac{1}{2}a^{3}R_{1}(a)\cos\chi.$

Therefore, the static component of the magnetic field will be

$\bm{H}_{\text{static}}(\bm{r})=\frac{1}{2}\frac{a^{3}}{r^{3}}R_{1}(a)\cos\chi\left(2\cos\theta\,\mathbf{\hat{r}}+\sin\theta\,\bm{\hat{\theta}}\right).$
(25)

####  4.2.2 Electric Quadrupole

We begin by assuming the solution for the static electric field to be a
quadrupole, after which we show that it is the only component required to
satisfy the boundary conditions due to the dependency on $\theta$. From G17,
equation (3.65), the quadrupole potential term is

$V_{\text{quad}}(r)=\frac{B_{2}}{r^{3}}\frac{1}{2}\left(3\cos^{2}\theta-1\right).$

where $B_{2}$ is the weighting coefficient from the series expansion in G17.

From equation (2.23) of G17, the electric field is the negative gradient of
the potential. Hence, the polar component of the $E$ field is

$E_{\text{quad},\,\theta}=-\frac{1}{r}\frac{\partial
V}{\partial\theta}=\frac{B_{2}}{r^{4}}3\cos\theta\sin\theta.$

The functions depending on the variable $\theta$ match those of the boundary
condition (14). Therefore, the boundary condition only requires the static
component of the electric field to be a quadrupole. Comparing the above result
to (14) we can uniquely determine $B_{2}$ for the situation,

$\begin{gathered}\frac{B_{2}}{a^{4}}3\cos\theta\sin\theta=-\omega\mu_{0}aR_{1}(a)\cos\chi\cos\theta\sin\theta,\\\
\implies B_{2}=-\frac{1}{3}\omega\mu_{0}a^{5}R_{1}(a)\cos\chi.\end{gathered}$

and on applying the gradient function to the potential in spherical
coordinates (see E),

$\displaystyle E_{\text{quad},\,r}$ $\displaystyle=-\frac{\partial V}{\partial
r}=-\frac{1}{2}\omega\mu_{0}a\frac{a^{4}}{r^{4}}R_{1}(a)\cos\chi\left(3\cos^{2}\theta-1\right),$
$\displaystyle E_{\text{quad},\,\theta}$
$\displaystyle=-\frac{1}{r}\frac{\partial
V}{\partial\theta}=-\omega\mu_{0}aR_{1}(a)\frac{a^{4}}{r^{4}}\cos\chi\cos\theta\sin\theta,$
$\displaystyle E_{\text{quad},\,\varphi}$
$\displaystyle=-\frac{1}{r\sin\theta}\frac{\partial V}{\partial\varphi}=0.$

The trigonometric identities
$3\cos^{2}\theta-1=\frac{1}{2}\left(3\cos{2\theta}+1\right)$ and
$\cos\theta\sin\theta=\frac{1}{2}\sin 2\theta,$ can be used to write the
quadrupole terms in a similar form as D55,

$\displaystyle E_{\text{quad},\,r}$
$\displaystyle=-\frac{1}{4}\omega\mu_{0}a\frac{a^{4}}{r^{4}}R_{1}(a)\cos\chi\left(3\cos
2\theta+1\right),$ (26) $\displaystyle E_{\text{quad},\,\theta}$
$\displaystyle=-\frac{1}{2}\omega\mu_{0}aR_{1}(a)\frac{a^{4}}{r^{4}}\cos\chi\sin
2\theta,$ $\displaystyle E_{\text{quad},\,\varphi}$ $\displaystyle=0.$

###  4.3 Derivation of General Field Equations

We are now in a position to derive the equations in the Appendix of D55. The
general field equations will be the sum of the multipoles derivations (25) and
(26), and the series expressions (19),

$\begin{gathered}\bm{H}=\bm{H}_{\text{dip}}-\frac{1}{ic\mu_{0}}(a_{1}\bm{\mathrm{n}}_{1,1}+a_{2}\bm{\mathrm{n}}_{1,2}+b_{1}\bm{\mathrm{m}}_{1,1}+b_{2}\bm{\mathrm{m}}_{1,2})e^{-i\omega
t},\\\
\bm{E}=\bm{E}_{\text{quad}}-(a_{1}\bm{\mathrm{m}}_{1,1}+a_{2}\bm{\mathrm{m}}_{1,2}+b_{1}\bm{\mathrm{n}}_{1,1}+b_{2}\bm{\mathrm{n}}_{1,2})e^{-i\omega
t},\\\ \end{gathered}$

(27)

where $a_{1}$, $a_{2}$, $b_{1}$ and $b_{2}$ are constants to be determined by
the boundary conditions (14). These constants may not necessarily be real
valued.

####  4.3.1 H Field Radial Component

By substituting the radial components from (23) and (25) into (27) we have,

$\displaystyle H_{r}=$
$\displaystyle\frac{a^{3}}{r^{3}}R_{1}(a)\cos\chi\cos\theta$
$\displaystyle-\frac{1}{ic\mu_{0}}\left(2a_{1}\frac{h_{1}(\rho)}{\rho}\sin\theta
e^{i\varphi}+9a_{2}\frac{h_{2}(\rho)}{\rho}\sin{2\theta}e^{i\varphi}\right)e^{-i\omega
t},$

which, when $r=a$ ($\rho=\alpha$), will be equal to the boundary condition on
$H_{r}$ from (14). We can use $e^{i\varphi}e^{-i\omega t}=e^{i(\varphi-\omega
t)}=e^{i\lambda}$ to write the exponent in terms of $\lambda$ and that
$e^{i\lambda}=\cos\lambda+i\sin\lambda$. The constants can then be determined
by matching trigonometric functions. For $\sin\theta$,

$\Re\left(-\frac{1}{ic\mu_{0}}a_{1}\frac{h_{1}(\alpha)}{\alpha}e^{i\lambda}\right)=R_{1}(a)\sin\chi\cos\lambda,$

where the $\Re$ operator has come from (22). Hence, $a_{1}$ must be purely
imaginary to cancel out $i$ in the bottom of the fraction so that
$\cos\lambda$ is obtained instead of $\sin\lambda$.

$\displaystyle-\frac{2}{ic\mu_{0}}a_{1}\frac{h_{1}(\alpha)}{\alpha}\cos\lambda=R_{1}(a)\sin\chi\cos\lambda,$
$\displaystyle\implies
a_{1}=-\frac{R_{1}(a)ic\mu_{0}\sin\chi}{2}\frac{\alpha}{h_{1}(\alpha)}.$ (28)

As there is no $\sin{2\theta}$ term in the boundary condition, we will have
$a_{2}=0$. The gives the radial component of the $H$ field as the real part of

$H_{r}=R_{1}(a)\left(\frac{a^{3}}{r^{3}}\cos\chi\cos\theta+\frac{\alpha}{h_{1}(\alpha)}\frac{h_{1}(\rho)}{\rho}\sin\chi\sin\theta
e^{i\lambda}\right).$

####  4.3.2 E Field Polar Component

Now considering the $E_{\theta}$ component which can be found by combing the
$\,\bm{\hat{\theta}}$ components of results (23) and (26) into (27),

$\displaystyle E_{\theta}=$
$\displaystyle-\frac{1}{2}\omega\mu_{0}aR_{1}(a)\frac{a^{4}}{r^{4}}\cos\chi\sin
2\theta-a_{1}ih_{1}(\rho)e^{i\varphi}e^{-i\omega t}$ $\displaystyle-
b_{1}\frac{1}{\rho}(\rho h_{1}^{\prime}(\rho)+h_{1}(\rho))\cos\theta
e^{i\varphi}e^{-i\omega t}$ $\displaystyle-b_{2}\frac{3}{\rho}(\rho
h_{2}^{\prime}(\rho)+h_{2}(\rho))\cos{2\theta}e^{i\varphi}e^{-i\omega t},$

and when $r=a$ this is equal to the $E_{\theta}$ component from (14). Once
again, matching functions of $\theta$ and using $e^{i\varphi}e^{-i\omega
t}=e^{i\lambda}$, we find that for the term excluding any function of
$\theta$,

$\begin{gathered}\Re\left(-\frac{c\mu_{0}R_{1}(a)}{2}\frac{\alpha}{h_{1}(\alpha)}h_{1}(\alpha)\sin\chi
e^{i\lambda}\right)\\\
=-\frac{1}{2}\omega\mu_{0}aR_{1}(a)\sin\chi\cos\lambda,\end{gathered}$

then it is convenient to change $\alpha$ to $(\omega/c)a$,

$-\frac{\omega\mu_{0}aR_{1}(a)}{2}\sin\chi\cos\lambda=-\frac{1}{2}\omega\mu_{0}aR_{1}(a)\sin\chi\cos\lambda,$

the LHS and RHS match, which validates the result for $a_{1}$ found in (28).
As the boundary condition does not contain any $\cos\theta$ terms,

$\begin{gathered}-b_{1}\frac{1}{\rho}(\rho
h_{1}^{\prime}(\rho)+h_{1}(\rho))\cos\theta e^{i\lambda}=0,\\\ \implies
b_{1}=0.\end{gathered}$ (29)

Lastly, the $\cos{2\theta}$ terms,

$\displaystyle\Re\left(-b_{2}\frac{3}{\alpha}(\alpha
h_{2}^{\prime}(\alpha)+h_{2}(\alpha))\cos{2\theta}e^{i\lambda}\right)$
$\displaystyle=-\frac{1}{2}\omega\mu_{0}aR_{1}(a)\cos{2\theta}\sin\chi\cos\lambda,$

again we want $\cos\lambda$ instead of $\sin\lambda$ so $b_{2}$ must be purely
real. Hence,

$\displaystyle-b_{2}\frac{3}{\alpha}(\alpha
h_{2}^{\prime}(\alpha)+h_{2}(\alpha))\cos{2\theta}\cos\lambda$
$\displaystyle=\frac{1}{2}\omega\mu_{0}aR_{1}(a)\cos{2\theta}\sin\chi\cos\lambda,$
$\displaystyle\implies
b_{2}=-\frac{1}{2}\omega\mu_{0}aR_{1}(a)\frac{1}{3}\frac{\alpha}{\alpha
h_{2}^{\prime}(\alpha)+h_{2}(\alpha)}\sin\chi.$ (30)

All the results for the constants thus far lead to,

$\begin{aligned}
E_{\theta}=&\frac{1}{2}\omega\mu_{0}aR_{1}(a)\left[-\frac{a^{4}}{r^{4}}\cos\chi\sin{2\theta}\right.\\\
&+\left.\left(\frac{\alpha}{\alpha
h_{2}^{\prime}(\alpha)+h_{2}(\alpha)}\frac{\rho
h_{2}^{\prime}(\rho)+h_{2}(\rho)}{\rho}\cos{2\theta}-\frac{h_{1}(\rho)}{h_{1}(\alpha)}\right)\sin\chi
e^{i\lambda}\right].\end{aligned}$

(31)

D55 has a typographical error in his Appendix for eqn. (31). The
$\frac{\alpha}{\alpha h_{2}^{\prime}(\alpha)+h_{2}(\alpha)}\frac{\rho
h_{2}^{\prime}(\rho)+h_{2}(\rho)}{\rho}$ term is instead written as
$\frac{\alpha h_{2}^{\prime}(\alpha)+h_{2}(\alpha)}{\alpha}\frac{\rho}{\rho
h_{2}^{\prime}(\rho)+h_{2}(\rho)}$.

####  4.3.3 E Field Azimuthal Component

The final boundary condition is on the $E_{\varphi}$ component. We combine the
results from above with azimuthal component of (27), which yields

$\begin{aligned}
E_{\varphi}=&-\left(-a_{1}h_{1}(\rho)\cos\theta+b_{2}\frac{3i}{\rho}(\rho
h_{2}^{\prime}(\rho)+h_{2}(\rho))\cos\theta\right)e^{i\lambda}\\\
=&\frac{1}{2}\omega\mu_{0}aR_{1}(a)\left(\frac{\alpha}{\alpha
h_{2}^{\prime}(\alpha)+h_{2}(\alpha)}\frac{\rho
h_{2}^{\prime}(\rho)+h_{2}(\rho)}{\rho}-\frac{h_{1}(\rho)}{h_{1}(\alpha)}\right)i\sin\chi\cos\theta
e^{i\lambda},\end{aligned}$

and clearly at $r=a$ ($\rho=\alpha$) this reduces to zero and matches the
required boundary condition $E_{\varphi}(a,\theta,\varphi)=0$, implying that
the values for the constants are consist with all the boundary conditions in
this situation.

####  4.3.4 Remaining Vector Field Components

Noting that the expressions for the constants $a_{1}$ and $b_{2}$ can be
rewritten using $\alpha=(\omega/c)a$, the field components are given by the
real part of the following,

$\begin{aligned}
H_{\theta}&=\frac{1}{2}\frac{a^{3}}{r^{3}}R_{1}(a)\cos\chi\sin\theta-\frac{1}{ic\mu_{0}}\left(a_{1}\frac{\rho
h_{1}^{\prime}(\rho)+h_{1}(\rho)}{\rho}\cos\theta+b_{2}(3ih_{2}(\rho)\cos\theta)\right)e^{i\lambda}\\\
&=\frac{1}{2}R_{1}(a)\left[\frac{a^{3}}{r^{3}}\cos\chi\sin\theta+\left(\frac{\alpha}{h_{1}(\alpha)}\frac{\rho
h_{1}^{\prime}(\rho)+h_{1}(\rho)}{\rho}\right.\right.\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\left.\left.+\frac{1}{c}\omega
a\frac{\alpha}{\alpha
h_{2}^{\prime}(\alpha)+h_{2}(\alpha)}h_{2}(\rho)\right)\sin\chi\cos\theta
e^{i\lambda}\right]\\\
&=\frac{1}{2}R_{1}(a)\left[\frac{a^{3}}{r^{3}}\cos\chi\sin\theta+\left(\frac{\alpha}{h_{1}(\alpha)}\frac{\rho
h_{1}^{\prime}(\rho)+h_{1}(\rho)}{\rho}+\frac{\alpha^{2}h_{2}(\rho)}{\alpha
h_{2}^{\prime}(\alpha)+h_{2}(\alpha)}\right)\sin\chi\cos\theta
e^{i\lambda}\right],\\\ \\\
H_{\varphi}&=-\frac{1}{ic\mu_{0}}\left(a_{1}\frac{i}{\rho}(\rho
h_{1}^{\prime}(\rho)+h_{1}(\rho))-3b_{2}h_{2}(\rho)\cos{2\theta}\right)e^{i\lambda}\\\
&=\frac{1}{2}R_{1}(a)\left(i\frac{\alpha}{h_{1}(\alpha)}\frac{\rho
h_{1}^{\prime}(\rho)+h_{1}(\rho)}{\rho}\sin\chi-\frac{1}{i}\frac{\alpha^{2}}{\alpha
h_{2}^{\prime}(\alpha)+h_{2}(\alpha)}h_{2}(\rho)\sin\chi\cos{2\theta}\right)e^{i\lambda}\\\
&=\frac{1}{2}R_{1}(a)\left(\frac{\alpha}{h_{1}(\alpha)}\frac{\rho
h_{1}^{\prime}(\rho)+h_{1}(\rho)}{\rho}+\frac{\alpha^{2}}{\alpha
h_{2}^{\prime}(\alpha)+h_{2}(\alpha)}h_{2}(\rho)\cos{2\theta}\right)i\sin\chi
e^{i\lambda},\\\ \\\
E_{r}&=-\frac{1}{4}\omega\mu_{0}a\frac{a^{4}}{r^{4}}R_{1}(a)\cos\chi\left(3\cos
2\theta+1\right)-9b_{2}\frac{h_{2}(\rho)}{\rho}\sin{2\theta}e^{i\lambda}\\\
&=\frac{1}{2}\omega\mu_{0}aR_{1}(a)\left(-\frac{1}{2}\frac{a^{4}}{r^{4}}\cos\chi(3\cos{2\theta}+1)+3\frac{\alpha}{\alpha
h_{2}^{\prime}(\alpha)+h_{2}(\alpha)}\frac{h_{2}(\rho)}{\rho}\sin\chi\sin{2\theta}e^{i\lambda}\right).\end{aligned}$

The components of the general field equations are gathered in A for easier
reference. A plot111For animations of the solution see
https://github.com/jsa113/General_Field_Equation_gifs of the magnetic field is
shown in Fig. 3.

Figure 3: The magnetic field from (41) using typical star parameters is
plotted as streamlines showing the direction and the log magnitude of the
magnetic field strength represented by the colour. Directly vertical is the
star’s rotation axis, and the magnetic field axis is set as $\chi=\pi/4$. Note
the dashed red circle is the light cylinder, where $\rho=1$ or $r=c/\omega$.

##  5 The Radiation Fields of the Star

Far from a star, the electromagnetic fields take on their radiation form,
allowing us to find the fields which ‘leak’ out of the magnetosphere and
propagate away from the electromagnetic source. To find the radiation form of
the electromagnetic fields we consider the limit in which $r\gg c/\omega$
(Griffiths, 2017). Hence, from (17) if $r\gg c/\omega$ then $\rho\gg 1$. Then
the limit of the spherical Bessel functions when $\rho$ is large (47) can be
substituted into (41). However, we find that $\alpha=(\omega/c)a\ll 1$, which
means $h_{1}(\alpha)$ can be approximated using (46). Using
$c\simeq$3\text{\times}{10}^{8}\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$
and the typical star values given by D55,

$\begin{gathered}a\simeq${10}^{9}\text{\,}\mathrm{m}$,\quad\omega\simeq${10}^{-5}\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{/}\mathrm{s}$,\end{gathered}$

we get $\alpha\simeq 3\times 10^{-5}\ll 1$. While a typical neutron star has

$\begin{gathered}a\simeq${10}^{4}\text{\,}\mathrm{m}$,\quad\omega\simeq
2\pi\,$\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{/}\mathrm{s}$,\end{gathered}$

which gives $\alpha\simeq 2\times 10^{-4}\ll 1$. We can also express this as a
scaled equation,

$\alpha=2\times
10^{-4}\left(\frac{a}{${10}^{4}\text{\,}\mathrm{m}$}\right)\left(\frac{\omega}{2\pi\,$\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{/}\mathrm{s}$}\right)$

The fastest recorded neutron star has
$a\simeq$1.6\text{\times}{10}^{4}\text{\,}\mathrm{m}$$ and
$\omega\simeq$4.5\text{\times}{10}^{3}\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{/}\mathrm{s}$$
(Hessels et al., 2006). We see from these values that $\alpha\simeq 0.2<1$.
Though this limit is not as strong as slower neutron stars, it shows an
appreciable upper limit that is still less than one. Suggesting that taking
the approximation $\alpha\ll 1$ in the following derivations is valid for
slower stars but can still apply to faster stars, but with some error.

###  5.1 Radiation Form Derivation

To derive the radiation form we take the limit in which $r\gg c/\omega$. In
taking this limit, we only consider terms that will be of the largest
magnitude for each electromagnetic field component, neglecting terms
containing high powers of $r$ in their denominator. Recall $\Re$ denotes the
real part of an expression from (22) and $\lambda=\varphi-\omega t$ from (6).
Working firstly on the $H$ field’s radial component, after neglecting smaller
magnitude terms and substituting the limit forms of $h_{1}(\rho)$ and
$h_{1}(\alpha)$ in, $H_{r}$ becomes

$\displaystyle H_{r}$
$\displaystyle\simeq\Re\left({R_{1}(a)\frac{e^{i\rho}}{\rho^{2}}}\frac{\alpha^{3}}{i}\sin\chi\sin\theta
e^{i\lambda}\right)$
$\displaystyle=\Re\left(R_{1}(a)\frac{\omega}{c}\frac{a^{3}}{r^{2}}\frac{1}{i}\sin\chi\sin\theta
e^{i(\lambda+\rho)}\right)$
$\displaystyle=\Re\left(R_{1}(a)\frac{\omega}{c}\frac{a^{3}}{r^{2}}\frac{1}{i}\sin\chi\sin\theta\left[\cos(\lambda+\rho)+i\sin(\lambda+\rho)\right]\right)$
$\displaystyle=\frac{\omega}{c}a^{3}R_{1}(a)\frac{1}{r^{2}}\sin\chi\sin\theta\sin\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right],$

where $\lambda+\rho=\varphi-\omega t+\omega r/c=\omega(r/c-t)+\varphi$.
Repeating the same limit and substituting the approximations in for the polar
component yields,

$\begin{aligned}
H_{\theta}&\simeq\Re\left(\frac{1}{2}R_{1}(a)\left[\left(\frac{\alpha^{2}}{\alpha\frac{9i}{\alpha^{4}}-\frac{3i}{\alpha^{3}}}\right)\frac{ie^{i\rho}}{\rho}+\frac{\alpha^{3}}{-i}\left(\frac{-ie^{i\rho}}{\rho}-\frac{e^{i\rho}}{\rho^{2}}\right)\right]\sin\chi\cos\theta
e^{i\lambda}\right)\\\
&=\Re\left(\frac{1}{2}R_{1}(a)\left[\frac{\alpha^{5}}{6\rho}+\frac{\alpha^{3}}{\rho}+\frac{\alpha^{3}}{i\rho^{2}}\right]\sin\chi\cos\theta
e^{i(\lambda+\rho)}\right).\end{aligned}$

The term with the largest magnitude in the line above is the $\alpha^{3}/\rho$
term as $\alpha\ll 1$ and $\rho\gg 1$. So,

$\displaystyle H_{\theta}$
$\displaystyle\simeq\Re\left(\frac{1}{2}R_{1}(a)\frac{\alpha^{3}}{\rho}\sin\chi\cos\theta
e^{i(\lambda+\rho)}\right)$
$\displaystyle=\Re\left(\frac{1}{2}R_{1}(a)\frac{\alpha^{3}}{\rho}\sin\chi\cos\theta\left[\cos(\lambda+\rho)+i\sin(\lambda+\rho)\right]\right)$
$\displaystyle=\frac{1}{2}\frac{\omega^{2}}{c^{2}}a^{3}R_{1}(a)\frac{1}{r}\sin\chi\cos\theta\cos\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right].$

The azimuthal component has similar intermediate steps as above. Hence
skipping intermediate working gives,

$\begin{aligned}
H_{\varphi}&\simeq\Re\left(\frac{1}{2}R_{1}(a)\left[\left(\frac{\alpha^{2}}{\frac{9i}{\alpha^{4}}-\frac{3i}{\alpha^{3}}}\right)\frac{ie^{i\rho}}{\rho}\cos{2\theta}-\frac{\alpha^{3}}{i}\left(\frac{-ie^{i\rho}}{\rho}+\frac{-e^{i\rho}}{\rho^{2}}\right)\right]i\sin\chi
e^{i\lambda}\right)\\\
&=\Re\left(\frac{1}{2}R_{1}(a)\left[\left(\frac{\alpha^{6}}{9-3\alpha}\right)\frac{1}{\rho}\cos{2\theta}+\frac{\alpha^{3}}{\rho}+\frac{\alpha^{3}}{i\rho^{2}}\right]i\sin\chi
e^{i(\lambda+\rho)}\right).\end{aligned}$

The largest term in the line above is the $\alpha^{3}/\rho$ term. Removing all
smaller magnitude terms yields,

$\displaystyle H_{\varphi}$
$\displaystyle\simeq\Re\left(\frac{1}{2}R_{1}(a)\frac{\alpha^{3}}{\rho}i\sin\chi
e^{i(\lambda+\rho)}\right)$
$\displaystyle=\frac{1}{2}\frac{\omega^{2}}{c^{2}}a^{3}R_{1}(a)\frac{1}{r}\sin\chi\sin\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right].$

Now consider the $E$ field. Using similar intermediate steps and the same
principle as before where the largest order term is the only term considered
and the limit forms of $h_{1}(\rho)$ and $h_{1}(\alpha)$ are used. The E
field’s radial component will be,

$\displaystyle E_{r}$
$\displaystyle\simeq\Re\left(\frac{1}{2}\omega\mu_{0}aR_{1}(a)3\left(\frac{\alpha}{\alpha\frac{9i}{\alpha^{4}}-\frac{3i}{\alpha^{3}}}\right)\frac{ie^{i\rho}}{\rho^{2}}\sin\chi\sin{2\theta}e^{i\lambda}\right)$
$\displaystyle=\Re\left(\frac{3}{2}\omega\mu_{0}aR_{1}(a)\frac{\alpha^{4}}{6}\frac{1}{\rho^{2}}\sin\chi\sin{2\theta}e^{i(\lambda+\rho)}\right)$

however, due to the $\alpha^{4}/\rho^{2}$ term, the radial component is
equivalently equal to zero when compared to the other $E$ field components,
this implies that $E_{r}\approx 0$. The polar component is,

$\begin{aligned}
E_{\theta}&\simeq\Re\left(\frac{1}{2}\omega\mu_{0}aR_{1}(a)\left[\frac{\alpha^{4}}{6i}\left(\frac{-e^{i\rho}+\frac{ie^{i\rho}}{\rho}}{\rho}\right)\cos{2\theta}-\frac{\alpha^{2}}{i\rho}e^{i\rho}\right]\sin\chi
e^{i\lambda}\right)\\\
&=\Re\left(\frac{1}{2}\omega\mu_{0}aR_{1}(a)\left[\frac{\alpha^{4}}{6i}\left(\frac{-1}{\rho}+\frac{i}{\rho^{2}}\right)\cos{2\theta}-\frac{\alpha^{2}}{i\rho}\right]\sin\chi
e^{i(\lambda+\rho)}\right)\\\
&\simeq-\frac{1}{2}\frac{\omega^{2}\mu_{0}}{c}a^{3}R_{1}(a)\frac{1}{r}\sin\chi\sin\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right].\end{aligned}$

By using some of the previous intermediate steps from above, the azimuthal
component is,

$\begin{aligned}
E_{\varphi}&\simeq\Re\left(\frac{1}{2}\omega\mu_{0}aR_{1}(a)\left[\frac{\alpha^{4}}{6i}\left(\frac{-1}{\rho}+\frac{i}{\rho^{2}}\right)-\frac{\alpha^{2}}{i\rho}\right]i\sin\chi\cos\theta
e^{i(\lambda+\rho)}\right)\\\
&\simeq-\frac{1}{2}\frac{\omega^{2}\mu_{0}}{c}a^{3}R_{1}(a)\frac{1}{r}\sin\chi\cos\theta\cos\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right].\end{aligned}$

The results of the radiation field equations are gathered in B as equations
(42) and (43). The results obtained here match the results obtained by D55.

These results are interesting as we can see the _retarded time_ affect,
$(\frac{r}{c}-t)$, on the magnetic field in the radiation zone. That is we
would not expect an instantaneous change to the electromagnetic fields when it
is far from the surface from the star as electromagnetic “news” travels at the
speed of light (G17). Instead we have a delayed affect which depends on the
distance from the star.

###  5.2 Power Equation - D55 Equation (15)

The power $P$ passing through a surface is the integral of the Poynting vector
($\bm{S}=\bm{E}\times\bm{H}$) over said surface. From G17, Eqn. (11.1),

$P(r,t)=\oint{\bm{S}\cdot d\bm{a}}.$ (32)

This is the power passing through the surface, so a negative must be
introduced to give the rate at which energy is radiated away from the system,

$-\frac{dW}{dt}=\oint{\bm{S}\cdot d\bm{a}}.$ (33)

For the $E$ and $H$ fields G17 advocates taking the fields when $r$ is very
large, when they are in radiation form. Therefore, by taking the field
equations (42) and (43) (which apply when $r\gg c/\omega$) leads to the
Poynting vector,

$\bm{S}=\begin{pmatrix}0\\\ E_{\theta}\\\
E_{\varphi}\end{pmatrix}\times\begin{pmatrix}H_{r}\\\ H_{\theta}\\\
H_{\varphi}\end{pmatrix}=\begin{pmatrix}E_{\theta}H_{\varphi}-E_{\varphi}H_{\theta}\\\
E_{\varphi}H_{r}\\\ E_{\theta}H_{r}\end{pmatrix}$

and by also taking the surface to be a sphere of radius $r$, this implies that
$S_{r}\parallel\bm{a}$, $S_{\theta}\perp\bm{a}$ and $S_{\varphi}\perp\bm{a}$,
as $\bm{a}$ is the normal of the surface (see Fig. 1(a)). So (33) becomes,

$-\frac{dW}{dt}=\int_{0}^{2\pi}\int_{0}^{\pi}(E_{\theta}H_{\varphi}-E_{\varphi}H_{\theta})r^{2}\sin(\theta)\,d\theta\,d\varphi.$

Upon applying expressions for the $H$ and $E$ fields from (42) and (43) we
get,

$\begin{aligned}
-\frac{dW}{dt}=\int_{0}^{2\pi}\int_{0}^{\pi}&\frac{1}{4}\frac{\omega^{4}\mu_{0}}{c^{3}}a^{6}R_{1}^{2}(a)\sin^{2}(\chi)\left\\{\sin^{2}\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right]\right.\\\
&\left.+\cos^{2}(\theta)\cos^{2}\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right]\right\\}\sin(\theta)\,d\theta\,d\varphi.\end{aligned}$

(34)

Now we evaluate the integrals over the functions that depend on $\theta$ and
$\varphi$. The double integral can firstly be expressed as,

$\displaystyle\begin{aligned}
\int_{0}^{2\pi}\int_{0}^{\pi}&\sin(\theta)\sin^{2}\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right]\\\
&+\sin(\theta)\cos^{2}(\theta)\cos^{2}\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right]\,d\theta\,d\varphi,\end{aligned}$
(35) and using the integration identities,
$\displaystyle\int_{0}^{2\pi}\sin(\theta)\,d\theta=2,$
$\displaystyle\int_{0}^{2\pi}\sin(\theta)\cos^{2}(\theta)\,d\theta=\frac{2}{3},$
$\displaystyle\begin{gathered}\int_{0}^{\pi}\sin^{2}\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right]\,d\varphi\\\
=\int_{0}^{\pi}\cos^{2}\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right]\,d\varphi=\pi,\end{gathered}$
(38)

which leads to _the integral_ (35) yielding $8\pi/3$. On combining this result
with (34) we get

$-\frac{dW}{dt}=\frac{2\pi\omega^{4}\mu_{0}}{3c^{3}}a^{6}R_{1}^{2}(a)\sin^{2}(\chi).$
(39)

This is the same result obtained by D55 in his equation (15). The RHS is
always positive. Hence, $dW/dt$ must always be negative. That is, energy is
always being radiated away from the star. The energy that is radiated comes
from the rotational kinetic energy causing the star’s spin down.

####  5.2.1 Spin Down Rate - D55 Equation (16)

To derive the spin down rate of a star in the VRD solution, we require two
standard relationships (see for example Moebs et al. (2016), Eqns. (10.31) and
(11.8)),

$\bm{\tau}=\frac{d\bm{L}}{dt},\qquad\frac{dW}{dt}=P=\bm{\tau}\cdot\bm{\omega},\quad\implies\quad\frac{dW}{dt}=\frac{d\bm{L}}{dt}\cdot\bm{\omega},$

where $\bm{L}$ is the angular momentum and $\bm{\tau}$ is the torque. As the
angular velocity and the angular momentum are in the same direction, we can
rearrange the equation above to become,

$\frac{d\bm{L}}{dt}=\frac{\bm{\omega}}{\left\lVert\bm{\omega}\right\rVert^{2}}\frac{dW}{dt},$

which is the result obtained by D55.

However, another typo is present in D55. It is not obviously shown that his
$\omega$ is actually a vector and that the $\omega$ in the bottom of the
fraction is the norm of $\omega$. Explaining why he has $\omega/\omega^{2}$
rather than just $1/\omega$. Lastly, D55 displays it with a negative so it can
be compared easier with (39),

$-\frac{d\bm{L}}{dt}=\frac{\bm{\omega}}{\left\lVert\bm{\omega}\right\rVert^{2}}\left(-\frac{dW}{dt}\right).$
(40)

It can be seen from this equation that the angular momentum also decreases
with time. The star will be brought to rest by its own radiation as D55
states.

##  6 Conclusion

By considering an idealised star and assuming a dipolar form for the internal
magnetic field, the general field equations of the electromagnetic fields
external of the star were found analytically. These general field equations
were then able to be reduced to their radiation limit ($r\gg c/\omega$). From
the radiation form, the relationship for the power radiated by the rotating
magnetic field was found.

It was found that the power equation matches the equation found in D55 and the
derivation for which is shown in this report for scrutiny. However, a typo
appears in D55’s general field equations. In particular, the $E_{\theta}$ (31)
component appears to evaluate the wrong fraction at $\alpha$. This error in
D55 causes trouble when reducing the general field equations to their
radiation form. When using the $E_{\theta}$ component found in this report the
reduction to radiation form matches D55.

A plot of the magnetic field using the general field equations was also
produced, showing the summation of the magnetic dipole and the wave solution.
This generalises a similar plot given in D55 which was zoomed in around the
$r\ll c/\omega$ region.

A 1969 paper by Goldreich and Julian argues that a pulsar cannot be surrounded
by a vacuum but rather a magnetosphere containing plasma (Goldreich & Julian,
1969). A related approach can be taken when estimating the power equation. The
Force-Free (FF) magnetosphere, which is the solution to Maxwell’s Equations
for a rotating star with a dipole field. It ignores the plasma pressure and
requires that the electric field parallel to the magnetic field is zero
everywhere (Contopoulos et al., 1999). This solution requires approximations
just like the VRD solution. The most realistic way is to use numerical
simulations. The results from Kalapotharakos et al. (2018) find their power
equations ranging from the near FF solution to near the VRD solution given in
equation (39). The VRD solution proves to be an interesting limiting case of
numerical simulations.

The VRD solution discussed in this review has been extended in a number of
ways. For example, general relativistic corrections to the Deustch solution
can be added (Rezzolla et al., 2001; Pétri, 2013). Also, considering an
electromagnetic field that consists of more multipole terms than just a dipole
adds extra correction terms that are potentially required to match reality
(Pétri, 2015; Bonazzola et al., 2020).

## References

  * Bonazzola et al. (2020) Bonazzola, S., Mottez, F., & Heyvaerts, J. 2020, arXiv e-prints, arXiv:2007.02539
  * Contopoulos et al. (1999) Contopoulos, I., Kazanas, D., & Fendt, C. 1999, The Astrophysical Journal, 511, 351, arXiv: astro-ph/9903049
  * Deutsch (1955) Deutsch, A. J. 1955, Annales d’Astrophysique, 18
  * Goldreich & Julian (1969) Goldreich, P., & Julian, W. H. 1969, ApJ, 157, 869
  * Griffiths (2017) Griffiths, D. J. 2017, Introduction to Electrodynamics, 4th edn. (Cambridge University Press), doi:10.1017/9781108333511
  * Hessels et al. (2006) Hessels, J. W. T., Ransom, S. M., Stairs, I. H., et al. 2006, Science, 311, 1901
  * Kalapotharakos et al. (2018) Kalapotharakos, C., Brambilla, G., Timokhin, A., Harding, A. K., & Kazanas, D. 2018, The Astrophysical Journal, 857, 44
  * Melatos (1997) Melatos, A. 1997, Monthly Notices of the Royal Astronomical Society, 288, 1049
  * Michel & Li (1999) Michel, F., & Li, H. 1999, Physics Reports, 318, 227
  * Michel (1991) Michel, F. C. 1991, Theory of neutron star magnetospheres, Theoretical astrophysics (Chicago: University of Chicago Press)
  * Moebs et al. (2016) Moebs, W., Ling, S. J., & Sanny, J. 2016, University Physics Volume 1 (Houston, Texas: OpenStax)
  * Pétri (2013) Pétri, J. 2013, MNRAS, 433, 986
  * Pétri (2015) Pétri, J. 2015, Monthly Notices of the Royal Astronomical Society, 450, 714
  * Rezzolla et al. (2001) Rezzolla, L., Ahmedov, B. J., & Miller, J. C. 2001, Monthly Notices of the Royal Astronomical Society, 322, 723
  * Roberts (2007) Roberts, P. H. 2007, in Encyclopedia of Geomagnetism and Paleomagnetism, ed. D. Gubbins & E. Herrero-Bervera (Dordrecht: Springer Netherlands), 7–11
  * Stratton (1941) Stratton, J. A. 1941, Electromagnetic Theory, International series in physics (McGraw-Hill book Company, Incorporated)
  * Travelle (2011) Travelle, P. A., ed. 2011, Pulsars: discoveries, functions, and formation (Hauppauge, N.Y: Nova Science Publishers)
  * Van Brummelen (2012) Van Brummelen, G. 2012, Heavenly Mathematics (Princeton University Press)

##  Appendix A General Field Equations

The results of the general field equations found in Sec. 4 are gathered here
for reference. The fields external a star for the VRD solution is the _real_
part of the following,

$\displaystyle H_{r}$
$\displaystyle=R_{1}(a)\left(\frac{a^{3}}{r^{3}}\cos\chi\cos\theta+\frac{\alpha}{h_{1}(\alpha)}\frac{h_{1}(\rho)}{\rho}\sin\chi\sin\theta
e^{i\lambda}\right),$ (41) $\displaystyle H_{\theta}$
$\displaystyle=\frac{1}{2}R_{1}(a)\left[\frac{a^{3}}{r^{3}}\cos\chi\sin\theta+\left(\frac{\alpha}{h_{1}(\alpha)}\frac{\rho
h_{1}^{\prime}(\rho)+h_{1}(\rho)}{\rho}+\frac{\alpha^{2}}{\alpha
h_{2}^{\prime}(\alpha)+h_{2}(\alpha)}h_{2}(\rho)\right)\sin\chi\cos\theta
e^{i\lambda}\right],$ $\displaystyle H_{\varphi}$
$\displaystyle=\frac{1}{2}R_{1}(a)\left(\frac{\alpha}{h_{1}(\alpha)}\frac{\rho
h_{1}^{\prime}(\rho)+h_{1}(\rho)}{\rho}+\frac{\alpha^{2}}{\alpha
h_{2}^{\prime}(\alpha)+h_{2}(\alpha)}h_{2}(\rho)\cos{2\theta}\right)i\sin\chi
e^{i\lambda},$ $\displaystyle E_{r}$
$\displaystyle=\frac{1}{2}\omega\mu_{0}aR_{1}(a)\left(-\frac{1}{2}\frac{a^{4}}{r^{4}}\cos\chi(3\cos{2\theta}+1)+3\frac{\alpha}{\alpha
h_{2}^{\prime}(\alpha)+h_{2}(\alpha)}\frac{h_{2}(\rho)}{\rho}\sin\chi\sin{2\theta}e^{i\lambda}\right),$
$\displaystyle E_{\theta}$
$\displaystyle=\frac{1}{2}\omega\mu_{0}aR_{1}(a)\left[-\frac{a^{4}}{r^{4}}\cos\chi\sin{2\theta}+\left(\frac{\alpha}{\alpha
h_{2}^{\prime}(\alpha)+h_{2}(\alpha)}\frac{\rho
h_{2}^{\prime}(\rho)+h_{2}(\rho)}{\rho}\cos{2\theta}-\frac{h_{1}(\rho)}{h_{1}(\alpha)}\right)\sin\chi
e^{i\lambda}\right],$ $\displaystyle E_{\varphi}$
$\displaystyle=\frac{1}{2}\omega\mu_{0}aR_{1}(a)\left(\frac{\alpha}{\alpha
h_{2}^{\prime}(\alpha)+h_{2}(\alpha)}\frac{\rho
h_{2}^{\prime}(\rho)+h_{2}(\rho)}{\rho}-\frac{h_{1}(\rho)}{h_{1}(\alpha)}\right)i\sin\chi\cos\theta
e^{i\lambda}.$

##  Appendix B Radiation Field Equations

The results of the radiation field equations found in Sec. 5 are gathered here
for reference. The magnetic field when $r\gg c/\omega$ is,

$\displaystyle H_{r}$
$\displaystyle=\frac{\omega}{c}a^{3}R_{1}(a)\frac{1}{r^{2}}\sin\chi\sin\theta\sin\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right],$
(42) $\displaystyle H_{\theta}$
$\displaystyle=\frac{1}{2}\frac{\omega^{2}}{c^{2}}a^{3}R_{1}(a)\frac{1}{r}\sin\chi\cos\theta\cos\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right],$
$\displaystyle H_{\varphi}$
$\displaystyle=\frac{1}{2}\frac{\omega^{2}}{c^{2}}a^{3}R_{1}(a)\frac{1}{r}\sin\chi\sin\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right].$

The electric field when $r\gg c/\omega$ is,

$\displaystyle E_{r}$ $\displaystyle=0,$ (43) $\displaystyle E_{\theta}$
$\displaystyle=-\frac{1}{2}\frac{\omega^{2}\mu_{0}}{c}a^{3}R_{1}(a)\frac{1}{r}\sin\chi\sin\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right],$
$\displaystyle E_{\varphi}$
$\displaystyle=-\frac{1}{2}\frac{\omega^{2}\mu_{0}}{c}a^{3}R_{1}(a)\frac{1}{r}\sin\chi\cos\theta\cos\left[\omega\left(\frac{r}{c}-t\right)+\varphi\right].$

##  Appendix C Analysis of Spherical Bessel Functions

The general solution for the fields includes spherical Bessel functions of the
third kind (first kind Hankel functions). This appendix highlights what these
functions look like and what they are approximately equal to when considered
at the limits $\rho\ll 1$ and $\rho\gg 1$. From Ref. Stratton (1941) page 405,

$\begin{gathered}h_{1}^{(1)}(\rho)=-e^{i\rho}\left(\frac{\rho+i}{\rho^{2}}\right),\\\
h_{2}^{(1)}(\rho)=ie^{i\rho}\left(\frac{\rho^{2}+3i\rho-3}{\rho^{3}}\right).\end{gathered}$
(44)

We will omit the superscript denoting the functions as the third kind for
convenience. Differentiating these functions with respect to $\rho$ yields,

$\begin{gathered}h_{1}^{\prime}=-e^{i\rho}\left(\frac{i\rho^{2}-2\rho-2i}{\rho^{3}}\right),\\\
h_{2}^{\prime}=-e^{i\rho}\left(\frac{\rho^{3}+4i\rho^{2}-9\rho-9i}{\rho^{4}}\right).\end{gathered}$
(45)

In the limit $\rho\ll 1$ eqns. (44) and (45) become,

$\displaystyle h_{1}\approx$ $\displaystyle\frac{-i}{\rho^{2}},\qquad$
$\displaystyle h_{2}\approx$ $\displaystyle\frac{-3i}{\rho^{3}},$ (46)
$\displaystyle h_{1}^{\prime}\approx$ $\displaystyle\frac{2i}{\rho^{3}},$
$\displaystyle h_{2}^{\prime}\approx$ $\displaystyle\frac{9i}{\rho^{4}}.$

In the limit $\rho\gg 1$ eqns. (44) and (45) become,

$\displaystyle h_{1}\approx$ $\displaystyle\frac{-e^{i\rho}}{\rho},\qquad$
$\displaystyle h_{2}\approx$ $\displaystyle\frac{ie^{i\rho}}{\rho},$ (47)
$\displaystyle h_{1}^{\prime}\approx$ $\displaystyle\frac{-ie^{i\rho}}{\rho},$
$\displaystyle h_{2}^{\prime}\approx$ $\displaystyle\frac{-e^{i\rho}}{\rho}.$

The form these Bessel functions take at these limits are relatively simple and
are required in the derivation of a eqns. (13), (14), (18) and (19) from D55
(equations (42) and (43) in this report).

##  Appendix D Associated Legendre Functions

In Appendix IV of Stratton (1941), the associated Legendre functions are given
as,

$\displaystyle P_{0}(\cos\theta)$ $\displaystyle=1,$ (48) $\displaystyle
P_{1}(\cos\theta)$ $\displaystyle=\cos\theta,$ $\displaystyle
P_{1}^{1}(\cos\theta)$ $\displaystyle=\sin\theta,$ $\displaystyle
P_{2}(\cos\theta)$
$\displaystyle=\frac{1}{2}(3\cos^{2}\theta-1)=\frac{1}{4}(3\cos{2\theta}+1),$
$\displaystyle P_{2}^{1}(\cos\theta)$
$\displaystyle=3\cos\theta\sin\theta=\frac{3}{2}\sin{2\theta},$

where $\theta\in[0,\pi]$, the subscript denotes the degree of the polynomial,
and the superscript denotes the order of the derivative.

##  Appendix E Vector Calculus in Spherical Coordinates

Let $f(r,\theta,\varphi)$ be a scalar function and
$\bm{A}(r,\theta,\varphi)=(A_{r},A_{\theta},A_{\varphi})$ be a vector-valued
function in spherical coordinates $(r,\theta,\varphi)$. From appendix A of
G17, the gradient of a scalar function is,

$\nabla f=\frac{\partial f}{\partial
r}\,\mathbf{\hat{r}}+\frac{1}{r}\frac{\partial
f}{\partial\theta}\,\bm{\hat{\theta}}+\frac{1}{r\sin\theta}\frac{\partial
f}{\partial\varphi}\,\bm{\hat{\varphi}}.$ (49)

The divergence of a vector field is,

$\nabla\cdot\bm{A}=\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}A_{r}\right)+\frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}\left(A_{\theta}\sin\theta\right)+\frac{1}{r\sin\theta}\frac{\partial
A_{\varphi}}{\partial\varphi}.$

(50)

The curl of a vector field is,

$\displaystyle\nabla\times\bm{A}=$
$\displaystyle\frac{1}{r\sin\theta}\left(\frac{\partial}{\partial\theta}\left(A_{\varphi}\sin\theta\right)-\frac{\partial
A_{\theta}}{\partial\varphi}\right)\,\mathbf{\hat{r}}$ (51)
$\displaystyle+\frac{1}{r}\left(\frac{1}{\sin\theta}\frac{\partial
A_{r}}{\partial\varphi}-\frac{\partial}{\partial
r}\left(rA_{\varphi}\right)\right)\,\bm{\hat{\theta}}$
$\displaystyle+\frac{1}{r}\left(\frac{\partial}{\partial
r}\left(rA_{\theta}\right)-\frac{\partial
A_{r}}{\partial\theta}\right)\,\bm{\hat{\varphi}}.$

##  Appendix F Solutions of Spherical Triangles

From Van Brummelen (2012) on page 98 the spherical Law of Cosines is given as,

$\cos c=\cos a\cos b+\sin a\sin b\cos C,$ (52)

and on page 63 the spherical Law of Sines is given as,

$\frac{\sin a}{\sin A}=\frac{\sin b}{\sin B}.$ (53)

Upon rearranging the spherical Law of Cosines for $\cos C$ we get,

$\cos C=\frac{\cos c-\cos a\cos b}{\sin a\sin b}.$ (54)

$a$$b$$A$$B$$C$$c$ (a)

$z,\bm{\omega}$$\chi$$\theta$$\tiny-\nu$$\tiny\xi$$\small\lambda$$\,\bm{\hat{\theta}}$$\bm{e}$$\psi$$\,\bm{\hat{\psi}}$$\xi$
(b)

Figure 4: Spherical triangles with general notation and the labels according
to parameters defined in the report. (a) Spherical Triangle with side lengths
and angles labelled according to the notation in this appendix. Note that
sides and angles labelled with the same letter are opposite each other. (b)
Description according to coordinate parameters and inclination of $\bm{e}$
axis. Note that $\nu$ is an angle which is $2\pi$ periodic so $-\nu\equiv
2\pi-\nu$ (note angles increase though anticlockwise).