# Quantitative coarse graining of laminar fluid flow penetration in rough
boundaries

Akankshya Majhi  Physical Chemistry and Soft Matter,
Wageningen University and Research, Wageningen, The Netherlands  Lars Kool
Physical Chemistry and Soft Matter,
Wageningen University and Research, Wageningen, The Netherlands  Laboratoire
de Physique et Mécanique des Milieux Hétèrogènes,
École supérieure de physique et de chimie industrielles de la Ville de Paris,
Paris, France  Jasper van der Gucht  Physical Chemistry and Soft Matter,
Wageningen University and Research, Wageningen, The Netherlands  Joshua A.
Dijksman****Corresponding author<EMAIL_ADDRESS>Physical Chemistry
and Soft Matter,
Wageningen University and Research, Wageningen, The Netherlands

###### Abstract

The interaction between a fluid and a wall is described with a certain
boundary condition for the fluid velocity at the wall. To understand how
fluids behave near a rough wall, the fluid velocity at every point of the
rough surface may be provided. This approach requires detailed knowledge of,
and likely depends strongly on the roughness. Another approach of modeling the
boundary conditions of a rough wall is to coarse grain and extract a
penetration depth over which on average the fluid penetrates into the
roughness. In this work we show that for a broad range of periodic roughness
patterns and relative flow velocities, a universal penetration depth function
can be obtained. We obtain these results with experiments and complementary
numerical simulations. Our results show that wall roughness boundary
conditions can be captured with an average “slip length” and so indicate that
surface patterning yields extensive control over wall slip.

##  I Introduction

The interaction of a fluid with a wall can be very complex, especially when
the wall possesses small-scale features [1]. Typically, complexities arise in
the form of slip layers [2, 3] and other complicated boundary conditions [4].
These wall-fluid interactions can, however, be tuned to control the fluid flow
behaviour further away from the wall. The walls of the flow geometry can be
modified by introducing roughness in the form of riblets or grooves on their
surface.

The concept of a slip length goes back to Navier, but the study on explicitly
patterned ribbed surfaces came to prominence around 40 years ago, with
inspiration from nature [1, 5, 6, 7]. Sharks have denticles on their skin
which reduces drag and allows for easy swimming. To understand the effect of
drag reduction in animals, work has been done on simplified model systems with
various designs for ridges such as L-shaped, V-shaped and U-shaped [8, 9, 10].
These studies primarily deal with variations in ridge dimensions such as ridge
depth, ridge spacing, ridge width and angle of ridge orientation and in
particular, their effect on the performance of the ridged surfaces [8, 9, 11,
12]. Since then, several studies have been carried out on engineered patterned
ribbed geometries in parallel plate and concentric cylindrical systems. These
patterned geometries have been studied in the context of reducing viscous
(friction) drag [13, 14, 15, 16, 17, 18, 19, 7, 20, 9, 21, 22] or suppressing
wall slip [23, 2, 24]. Walsh [14] indicated that microsurface geometry
variations change the near-wall structure of the flow turbulent boundary layer
and thus, are effective in reducing (viscous) drag. Bechert et al. [17, 5]
carried out an extensive parametric study on surfaces with longitudinal ribs,
where they argued that the velocity profile in between the ridges penetrates
to a distance below the ridge tips which they refer as “protrusion height”.
This protrusion height depended on the ridge dimensions and was indicative of
the drag reduction. Despite showing interesting results, their geometries
posed a problem in terms of manufacturing and durability. Davies et al. [18]
numerically studied the effect of patterned channel walls with alternating
microribs and cavities. They investigated the influence of the vapour cavity
depth in the entire laminar flow regime and showed significant reduction in
the frictional resistance in laminar fluid flow. Their experimental studies
were documented in [25]. Experimental work by Maynes et al. [26] investigated
the laminar flow in a parallel plate microchannel with ultrahydrophobic top
and bottom walls and indicated dramatic decrease in the overall flow
resistance. Djenidi et al. [10, 27] showed that riblets cause a reduction in
frictional drag in laminar flows. These drag reducing techniques find
applications in aerospace industry for saving fuel costs [13, 14, 15, 16] and
also serve as a non-additive drag reduction technology in pipe flow.

Even very recently, McKinley and co-workers [21] used different surface
microtextures to study the drag reducing effect on the fluid flow over a range
of flow speeds. They examined the interaction of Taylor vortices with the
riblets in case of Newtonian fluids. Nickerson and Kornfield [23] have shown
that cleated surfaces on parallel plate geometry can be used to suppress wall
slip. Although most of the previous researchers have established that rough
surfaces can often be used to prevent wall slip, for instance in rheological
measurements, so far they have not considered that this wall roughness also
modifies the flow of the fluid near the wall, either by introduction of
secondary flows or by penetration of the fluid flow into the gaps between the
ridged wall. Additionally, riblets used are typically very complex in
structure and the flow fields in such geometries are not easy to study
systematically. Such complexities make interpretation of rheological
measurements on the effect of rough walls difficult.

In the present paper, we show that the role of the wall roughness on the
boundary layer of a laminar fluid flow can be effectively quantified by an
effective penetration depth for a range of different roughness conditions. We
use simple 3D printed ridged concentric cylinder geometries with ridges of
different depth, spacing and orientation. This way of creating rough surfaces
gives the ability to systematically investigate the role of a ridge in fluid
flow, to demonstrate how the flow penetrates between the ridges and to examine
whether there are effects of secondary flows, orthogonal to the primary flow
direction. Fig. 1(a) shows examples of patterned geometries. $D$, $S$ and $w$
are the ridge depth, ridge spacing, ridge width, respectively. The ridge angle
$\theta$ can be varied with respect to the vertical axis. A helical wall
patterning is shown to induce secondary flows at all flow rates, similar to
Taylor vortices which are usually observed at high rotation rates as indicated
by the Taylor number [28].

Figure 1: (a) 3D printed geometries with patterned surfaces: vertical ridges
(left), angled ridges (centre) and horizontal ridges (right). (b) Schematic
representation of a concentric cylinder (CC) geometry with angled ridges.
$R_{i}$ and $R_{o}$ are the radii of the inner and outer cylinder of the
ridged geometry, respectively. Here the inner radius $R_{i}$ corresponds to
the distance from the center of the coordinate system to the farthest point on
the ridges. $D$ is the ridge depth, $S$ is spacing between the ridges, $w$ is
the ridge width.

In this work, we investigate the influence of 3D printed ridged geometries on
the flow of a Newtonian fluid at low Reynolds number using standard concentric
cylinder rheology. In Section II, we first describe the materials, and
experimental and numerical methods used for the determination of protrusion
height. Later in Section III, we discuss both numerical and experimental
results of this flow profile penetration. We first present the simulated flow
profiles in different ridged geometries. Then we explore the flow profile
penetration, both experimentally and numerically, for two limiting cases of
horizontal and vertical ridges and then discuss the angled ridge case.

##  II Materials and Methods

###  II.1 Design and fabrication of ridged geometries

Custom 3D printed ridged geometries were fabricated on a Stratasys Objet30
Scholar and a Formlabs 2 SLA (Formlabs, Inc.) printer to fit the Anton Paar
rheometer shaft for disposable geometries.

This way of patterning with 3D printing allows to create well-defined
geometries, compared to other manufacturing methods, with a good control to
produce fine ($<$1 mm) features with a dimensional accuracy of $\sim$100
$\mu$m. The additional benefit of using the 3D printed ridged geometries is
that they can be easily incorporated into a flow imaging technique such as MRI
[29, 30], which allows to fully quantify the flow profiles and extract the
wall stresses present in the system.

The geometries with ridges perpendicular to the rotational axis were designed
by making cutouts of various spacings (0.3$\\--$3 mm) and depths (0.5$\\--$3
mm) from the standard CC17 geometry with constant $R_{i}=8.5$ mm and $R_{o}=9$
mm. They were printed on the Objet30 Scholar, using the Vero Black
photopolymer (Stratasys) with water-soluble supports. The geometries were
printed vertically with a layer thickness of 16 $\mu$m.

The geometries with ridges at an angle with respect to the rotational axis
were designed with the angle ranging from 10∘ to 60∘ with 10∘ intervals.
Geometries with ridge angles between 60∘ to 90∘ could not be printed due to
limitations of the printing technique. The geometries were designed such that
the shortest distance (spacing) between the ridges was 2 mm, with a ridge
depth of 1 mm and ridge thickness of 0.3 mm. This was achieved by changing the
width and the number of cutouts. The width and the number of cutouts were
calculated using a Python script with the length and the diameter of the
ridged geometry, ridge angle and shortest distance between the ridges as
input. These geometries with angled ridges were fabricated on a Formlabs 2 SLA
printer, using the glass reinforced Rigid resin (Formlabs, Inc.) with a layer
height of 50 $\mu$m and xy-resolution of 140 $\mu$m. The support needed for
the printing was only attached to the shaft and the top chamber to prevent the
introduction of artifacts coming from small residuals of the support material.

After printing, the geometries were washed twice for 10-15 minutes in 95$\%$
isopropanol. After drying, the geometries were post-cured under a 366 nm lamp
for 1.5-2 hours. The geometries were rotated 180∘ halfway through to cure
evenly. The support material was removed manually using flush side-clippers
after curing. The shaft of all geometries was printed with a diameter of 8 mm
and milled down on a lathe to 6.95 mm to fit the Anton Paar shaft for
disposable geometries. We determined that this lathe post-processing ensured
the best concentricity of the geometry with the Anton Paar rheometer shaft and
cup.

###  II.2 Rheological experiments

Rheology experiments were performed on Anton Paar MCR 300 and MCR 501
rheometers, using the custom geometries in standard cups and cup holders
(Anton Paar CC17). The experiments were performed at 20.5 ∘C, using either a
Peltier heat exchange element and waterbath at 20 ∘C as heat sink (for MCR
501) or a high flowrate waterbath at 20.5 ∘C (for MCR 300). All the printed
geometries were tested on a Newtonian fluid. We used castor oil as a Newtonian
fluid, due to its high viscosity ($\sim$1 Pa s at 20 ∘C compared to 50 mPa s
for most oils) and its stability over time. Unlike other viscous fluids like
glycerol, castor oil is non-hygroscopic and it does not evaporate at room
temperature, does not degrade over time, and does not swell nor degrade the 3D
printed geometries (contrary to most organic solvents).

Before each measurement, the geometry was placed in a container with the test
fluid to eliminate air bubbles between the ridges. The geometries were then
visually inspected for the presence of air bubbles, which, if present, were
removed by rubbing the geometry against the cup while it was submerged in the
fluid.

We performed steady shear measurements with average shear rates ranging from
$\dot{\gamma}$ = 10-3 to 102 s-1, resulting in a range of Re from 2.5 $\times$
10-5 to 2.5 for castor oil (only the laminar regime was probed). Each
measurement consisted of 26 datapoints, distributed logarithmically, and each
datapoint started with a stepwise increase of $\dot{\gamma}$, followed by an
equilibration period of maximum 2 min, after which the data was averaged for 5
s. For each measurement, the samples were presheared at $\dot{\gamma}$ = 1000
s-1 for 60 s followed by a recovery period of 60 s, to ensure that the stress
from loading the sample in the measurement geometry did not influence the
results and also to establish reproducible initial conditions.

To test the viability of using 3D printed geometries in rheology measurements,
a solid CC geometry was designed with an outer diameter as close as possible
to that of a commercially available smooth standard stainless steel CC
geometry. These geometries were checked by comparing the flow curves of the
test fluid. Both the 3D printed and the standard geometries yielded identical
flow curves indicating that the accuracy of the 3D printer was high enough to
produce complex rheology geometries and no deviations were expected in the
results due to the intrinsic roughness of the 3D printed geometries.

###  II.3 Simulations

Finite element simulations were done using COMSOL Multiphysics 5.6 to
complement the experiments. To simplify the calculations, the ridged
geometries were approximated as infinitely large parallel plate geometries. We
performed simulations for varying $S$, $D$, and $\theta$, while keeping the
gap size $R_{\textnormal{gap,0}}$ (i.e. the distance between the outer wall
and the tip of the ridges) and ridge width $w$ constant and equal to the
experimentally used values. We then impose a constant velocity difference
$\Delta v$ between the inner and outer wall, and solve the Stokes equation
(hence, ignoring inertial effects) numerically using quadratic quadrilateral
elements for the velocity components and linear elements for the pressure.
Since the flat plate limit does not account for the finite curvature of the
actual CC geometry, we performed additional simulations for curved geometries
with horizontal and vertical ridges to confirm that the effects of curvature
were small.

###  II.4 Determination of penetration depth

When a Newtonian fluid is sheared between two surfaces at a certain shear rate
$\dot{\gamma}$, a shear stress $\sigma$ is generated and it can be defined as
$\sigma=\eta\dot{\gamma}$, with $\eta$ being the viscosity of the fluid. The
fluid flow behaviour can be modified using different boundary conditions,
imposed by employing different geometries in a rheometer. Concentric cylinder
(CC) is one of the most commonly used geometries, where the fluid is sheared
in the gap between the inner and outer cylinder of radius $R_{i}$ and $R_{o}$,
respectively. The inner cylinder is made to rotate at a rotational velocity
$\Omega_{i}$ in revolutions per second (rps). The average shear rate
$\dot{\gamma}$, in s-1, is defined as the surface velocity of the moving inner
cylinder divided by the gap size $R_{\textnormal{gap}}$ (Eq. (1)). It should
be noted that the shear rate is not constant within the gap.

$\dot{\gamma}=\frac{2\pi R_{i}\Omega_{i}}{R_{\textnormal{gap}}}\hskip
28.45274pt\text{where}\hskip 28.45274ptR_{\textnormal{gap}}=R_{o}-R_{i}$ (1)

The torque $M$ acting across the cylindrical surface is given by Eq. (2),
where $\omega_{i}$ is the angular velocity in rad s-1 [31, 32].

$M=\frac{4\pi\eta\omega_{i}L(R_{o}R_{i})^{2}}{R_{o}^{2}-R_{i}^{2}}\hskip
14.22636pt\text{where}\hskip 14.22636pt\omega_{i}=2\pi\Omega_{i}$ (2)

The shear force in the azimuthal direction $F_{\textnormal{shear}}$ induced by
the fluid can be expressed as the ratio of $M$ and the radius of the inner
cylinder $R_{i}$. Thus, the shear stress in the azimuthal direction $\sigma$
can be expressed as $F_{\textnormal{shear}}$ divided by the surface area of
the inner cylinder $A_{i}$ (Eq. (3)).

$\sigma=\frac{F_{\textnormal{shear}}}{A_{i}}=\hskip 2.84544pt\frac{M}{2\pi
R_{i}^{2}L}\hskip 2.84544pt\text{where}\hskip
5.69046ptF_{\textnormal{shear}}=\frac{M}{R_{i}}\hskip 2.84544pt,\hskip
8.5359ptA_{i}=2\pi R_{i}L$ (3)

When roughness is introduced into the CC geometry in the form of ridges (as
shown in Fig. 1(b)), the flow profiles penetrate in the space between the
ridges. We now make the step to quantify this flow profile penetration not via
the flow field, but via the effect it has on the effective stress on the wall.
We denote the extent of this penetration by a parameter called the penetration
depth $\delta$, similar to the protrusion height as given in [5, 33, 17, 24,
34]. Typically, the penetration depth is determined by inducing fluid flow
over a flat surface with ridges on it. In their theoretical investigation,
Bechert and coworkers [5] simulated velocity distributions for various ridge
configurations and found that the apparent origin of the velocity profile in
between the ridges lies in the gap between the inner flat surface and the tip
of the ridge. Lee and Lee [34] experimentally investigated the flow structures
inside the semi-circular ridges. They referred to the distance between this
origin and the ridge tip as the protrusion height. Another equivalent way to
determine penetration depth is to create ridges on concentric cylinder
geometries and introduce continuous flow by rotating the ridged geometry with
respect to the fluid. We use the second method to generate different flow
patterns. Due to the introduction of the ridged geometries, the effective gap
size ${R_{\textnormal{gap}}}$ changes as a function of the ridge spacing $S$
and the ridge depth $D$, which is denoted as

$R_{\textnormal{gap}}(S,D)=R_{o}-(R_{i}-\delta)$ (4)

where $R_{i}$ is assumed to be identical to the radius of the inner cylinder
without the cutouts (as in Eq. (1)). Due to the penetration of the flow in
between the ridges, the effective gap size increases. The ridged geometry can
thus be considered similar to the standard solid CC geometry with an effective
inner radius $R_{i,\textnormal{eff}}=R_{i}-\delta$.

For every geometry with a particular ridge spacing $S$ and ridge depth $D$,
the penetration depth $\delta$ can be measured from standard rheology
experiments. Rewriting the torque-rotation rate conversion equation (Eq. (2)),
we get

$M=\frac{4\pi\eta\omega_{i}L(R_{o}R_{i,\textnormal{eff}})^{2}}{R_{o}^{2}-R_{i,\textnormal{eff}}^{2}}\hskip
14.22636pt\text{where}\hskip 14.22636ptR_{i,\textnormal{eff}}=R_{i}-\delta$
(5)

In this equation, $R_{i}$, $R_{o}$, $L$ are known from the geometry and
$\Omega_{i}$, $M$ were measured from the rheology experiments. $\eta$ was
determined by fitting the flow curves in a non-ridged (solid) CC geometry with
$\sigma=\eta\dot{\gamma}$. The penetration depth for every rotation rate can
be obtained by solving Eq. (5). For Newtonian fluids, the penetration depth
does not depend on the rotation rate $\Omega_{i}$, therefore, we averaged the
penetration depths corresponding to ten higher rotation rates to give a single
penetration depth value for a geometry with a specific ridge spacing and ridge
depth.

From the simulations in the flat plate limit, the penetration depth was
calculated from the average shear stress on the outer wall using
$\sigma=\eta\dot{\gamma}_{\textnormal{eff}}=\eta\Delta
v/(R_{\textnormal{gap,0}}+\delta)$.

##  III Penetration of the flow profiles between the ridges

We investigated the effect of different ridged geometries on the flow
behaviour of a Newtonian fluid. We numerically predicted the flow profiles for
the three cases of ridged geometries: horizontal, vertical and angled ridges
using finite element simulations. In the horizontal ridge case, the imposed
flow direction is parallel to the ridges and the flow is purely in the
azimuthal directon while in the vertical ridge case, the azimuthal flow is
perturbed by the ridges, and the flow includes a radial velocity component.
Finally, in the angled ridge case, the ridges exert an axial force on the
fluid, therefore, there is also an axial velocity component (Fig. 2(c)).

From the azimuthal velocity profiles (Fig. 2(a,b)), we observed that the fluid
flow profile penetrates into the gap between the ridges up to a certain
extent. The extent of this penetration is denoted by a parameter that we call
the penetration depth ($\delta$). This flow profile penetration will influence
the torque measured on the rheometer, thus, by measuring the torque, we can
quantify the penetration depth as described in Subsection II.4. We then
characterised the penetration depth for different ridge spacing, ridge depth
and ridge angles.

Figure 2: Simulated flow profiles of a Newtonian fluid in different ridged
geometries: (a) horizontal (b)-(c) angled ridges. Colorbar shows velocity in
the respective directions. Here the azimuthal and axial velocities are of the
order of $\sim$ 10-3 and 10-4 in m s-1, respectively. For these calculations,
the curved geometries were approximated as flat (infinite radius) parallel
plate geometries.

We first evaluated the penetration depth for the two limiting cases of
horizontal and vertical ridges, i.e. with ridges that are at 90∘ and 0∘ angles
with respect to the vertical axis, respectively.

###  III.1 Limiting case I: Horizontal ridges

We measured the flow curves from the rotational tests on multiple horizontal
ridged geometries. On a log-log scale, the torque $M$ increased linearly with
the rotation rate $\Omega$ for over four orders of magnitude in rotation rate
(Fig. 3(a)), indicating that the presence of ridges did not change the
Newtonian shear rate dependence of the flow. To investigate the variations due
to the different roughness parameters, we looked more closely on a linear
scale (Fig. 3(b)) and found that the measured torque decreases with increasing
ridge depths $D$ for a given ridge spacing $S$. This is due to the influence
of the flow profile penetration on the measured torque, i.e., the penetration
depth ($\delta$) increases with ridge depth, leading to an increase in the
effective gap size and a decrease in the measured torque. It should be noted
that in all these experiments, the rotation rates were low enough to avoid
Taylor vortices.

Fig. 3(c) shows the penetration depth, obtained from the measured torques
using Eq. (5), as a function of $D$ for various ridge spacings $S$. We see
that $\delta$ first increases with $D$, but then levels off to a plateau that
depends on the value of $S$. Hence, the flow can penetrate between the ridges
only up to a certain maximum depth, $\delta_{\infty}$, which increases
approximately linearly with the ridge spacing $S$ (Fig. 3(d)).

Figure 3: Experimental flow curves in (a) double logarithmic scale and (b)
linear scale of castor oil measured in different horizontal ridged geometries.
For all geometries, $S$ = 3 mm and $D$ is indicated in the legend. (c)
Variation of $\delta$ in a horizontal ridged geometry with ridge depth for
different ridge spacing (legends) Experiments: symbols+lines, simulations:
lines. (d) Plateau value of $\delta$ as a function of ridge spacing. A linear
fit through the origin give a of slope of 0.18. (e) Master curve showing the
collapse of the curves onto one exponential function when plotted as
normalised $\delta$ versus normalised ridge depth. Experiments: symbols,
simulations: lines.

To complement the experiments and to gain more insight into the flow profile
penetration, finite element simulations were performed. The geometries were
recreated as described in Subsection II.3, and the curvature of the geometry
was neglected. For this geometry, the only non-zero velocity component is in
the azimuthal ($\phi$) direction, while all gradients in the $\phi$-direction
are zero. The pressure was also constant in this case.

The simulated penetration depths as a function of different ridge depths and
ridge spacings is plotted in Fig. 3(c) as dashed lines. The experimental and
simulation results match quite well, especially at larger $D$. Similar to the
experiments, the $\delta$ values plateau after a certain ridge depth; as in
the experiments, the plateau value $\delta_{\infty}$ is found to increase
linearly with $S$, see Fig. 3(d).

For a solid cylinder case ($S=0$), $\delta_{\infty}$ must be zero as the fluid
cannot penetrate into a solid cylinder, and therefore, we would normally
expect $\delta_{\infty}$ to scale linearly with the ridge spacing and the
linear plot (dashed line in Fig. 3(d)) to extrapolate to the origin. However,
this is not observed and $\delta_{\infty}$ seems to deviate from the linear
dependence on $S$ at small ridge spacings below 0.5 mm. This is because in
this small $S$ limit, the ridge spacing becomes comparable to the ridge width,
which also starts to play a role in determining the penetration depth.
Strikingly, the normalised penetration depth $\delta/\delta_{\infty}$ for
different ratios of $D/S$ collapses onto a master curve (Fig. 3(e)). This
master curve can be described by a simple exponential function:

$\frac{\delta}{\delta_{\infty}}=1-e^{kD/S}\hskip 14.22636pt\text{with}\hskip
14.22636ptk\approx-6.65$ (6)

The initial part of this master curve, corresponding to small $D$ or large
$S$, can be linearized to give $\delta/\delta_{\infty}\approx kD/S$. It is
expected that $\delta\rightarrow D$ in this regime, so that we expect
$k=1/k_{2}$ with $k_{2}$ the slope of the curve of $\delta_{\infty}$ versus
$S$. We find $1/k_{2}\approx 5.7$ from Fig. 3(d), which is indeed close to
(but slightly below) the value found for $k$. This can also be understood as
follows: a static boundary layer is created due to the presence of the small
ridges. The thickness of this static boundary layer is $D-\delta$ and it goes
to zero as the spacing between the ridges becomes very large. The prefactor
$k$ takes into account the other geometric effects, most likely the effect
from the other length scale, i.e., the ridge width $w$. At a constant ridge
spacing $S$, the thicker the ridges, the lower will be the flow penetration,
hence, $\delta$ will be smaller. We expect that as the ridge width increases,
$\delta$ will be small and hence, $\delta_{\infty}$ will also drop. If the
normalised penetration depth $\delta/\delta_{\infty}$ decreases, then $k$
increases. But if both $\delta$ and $\delta_{\infty}$ decrease by the same
factor, then $\delta/\delta_{\infty}$ will not change, hence, $k$ will not
change. Therefore, $k$ can be written as a function of normalized $w$ as:
$k=f(w/S)$.

For smaller $D$, experimental results deviate from that of simulations and
this is a consistent deviation for all ridge spacings. There are some possible
sources of error. It could be that the small inaccuracies during the
fabrication of the geometries have a bigger influence at smaller cutouts.
Moreover, the extraction of penetration depths from experiments is affected by
several sources of error caused while measuring the actual inner and outer
radii and the length of the cylindrical geometry and also due to end effects
from the top and bottom conical parts of the geometry. For every geometry, the
actual dimensions of the 3D printed geometries were measured. The contribution
of end effects to the torque was calibrated on a solid cylindrical geometry.
The torques for the ridged geometries were corrected for these end effects.
Additionally, boundary conditions like wall slip and end effects are
uncontrolled in the experiments especially for large $S$ and small $D$, which
may contribute to small errors. Another possibility is that small air bubbles
may be trapped between the ridges, especially for small ridge spacings. Due to
these reasons, we took the plateau value of delta in experimental results as
the average of the last three $\delta$ values. By linking the experiments and
the simulations, it is possible to determine how changing the $D$ and $S$ of
the ridges changes the penetration depth ($\delta$).

###  III.2 Limiting case II: Vertical ridges

The second limiting case for patterned walls is a cylindrical geometry with
vertical ridges. Here, we show only the numerical results in the limit where
the curvature can be neglected (i.e. for large cylinder radii), see Subsection
II.3. In this case, the fluid velocity has components both in the $\phi$ and
the $r$-directions, and gradients in the $z$-direction are zero.

The penetration depth for vertical ridges with different ridge depths and
ridge spacings was calculated in the same way as before and plotted in Fig.
4(a). As for the horizontal ridges, the flow profiles penetrate more deeply in
between the ridges with increasing ridge depth and spacing, therefore,
$\delta$ increases. The penetration depth values are smaller than for the
horizontal ridge case because the main flow direction is perpendicular to the
ridges while for the horizontal ridges, the flow is along the ridges. Also,
contrary to the horizontal ridge case, we see two (quasi-)plateau values of
$\delta$ as a function of $D$: the first plateau value is at around $D$ = 0.5
mm, where the ridge depth is close to the width of the ridge and most likely,
the plateau value is caused by the effect of the ridge width. The second
plateau value is at large ridge depths similar to what was observed in case of
horizontal ridges. We consider the second plateau value and denote this as
$\delta_{\infty}$.

As for the horizontal ridges, $\delta_{\infty}$ increases linearly with the
ridge spacing $S$, but the slope is about 2.5 times smaller (Fig. 4(b)). As
observed in the previous case, the linear plot does not extrapolate exactly to
the origin and there seems to be non-linear dependence at small ridge spacings
below 0.5 mm.

The normalised penetration depth $\delta/\delta_{\infty}$ for different ratios
of $D/S$ again collapses quite well onto a master curve (Fig. 4(c)). However,
in this case, a single exponential fit does not work so well compared to what
was seen in the horizontal case and there seems to be a two step function fit.
This is mostly due to the presence of two plateaus in the $\delta$ vs D
variation (Fig. 4(a)). For $S\gg D$ i.e. for small values of $D/S$, we again
expect $\delta\rightarrow D$; the plot is linear in this region and the slope
of this plot is inversely related to the slope of the variation of
$\delta_{\infty}$ with $S$.

Lastly, the relation between $\delta_{\infty}$ and $S$ might depend on the
width of the ridges $w$ and on the gap size ($R_{o}-R_{i}$). There are a few
length scales to explore, which may give different regimes (for example: two
plateau values in $\delta$ for vertical ridges). The current plots shown are
for $S\gg w$. This is a very general rescaling that to our knowledge has not
been reported in literature, and it is helpful to get a better understanding
of the angled ridge case, as discussed below.

Figure 4: (a) Variation of simulated $\delta$ in a vertical ridged geometry
(considered as flat plate) with ridge depth and different ridge spacing
(legends). Here we see a first plateau at the dashed line, followed by a
second plateau at large ridge depths. (b) Plateau value of $\delta$ as a
function of ridge spacing (slope = 0.067). (c) Normalised $\delta$ as a
function of normalised ridge depth.

###  III.3 Angled ridge case

We can now look at the influence of geometries with ridges at finite angles on
the flow behaviour of a Newtonian fluid. For different angles, the torque
required to sustain a certain $\Omega$ again shows a linear dependence on the
rotation rate even at higher rotation rates, indicating that the angled ridges
do not affect the Newtonian shear rate dependence of the flow (Fig. 5(a)). All
the angled ridges have the same ridge spacing and ridge depth. However, the
slope of the torque versus rotation rate changes as a function of ridge angle
(Fig. 5(b)), indicating that the extent of flow profile penetration into the
gap between the ridges changes with the angle. To examine the effect of
different angles, we calculated the penetration depths from the $M$ vs
$\Omega$ relationship as previously done for the limiting cases.

Interestingly, $\delta$ shows a monotonic dependence on the ridge angle (Fig.
5(c)) and increases as $\theta$ increases, with the largest $\delta$ for the
horizontal ridges ($\theta=90^{\circ}$). There are a few outliers at smaller
ridge angles. One of the plausible reasons is $\delta$ is derived from the
assumption that there is a shear stress only in the azimuthal direction. In
reality, there is a normal force (and hence, a shear stress in the axial
direction) and it is not clear whether the normal force is completely
decoupled from the azimuthal shear stress.

We again compare our experimental results with numerical results. In this
case, the velocity has non-zero components in the $\phi$, $r$, and
$z$-directions, and also the gradients in all three directions are non-zero,
so, a full 3D-calculation is necessary. As we restrict ourselves to a flat
geometry (corresponding to large cylinder radii), we can make use of
symmetries to limit the computational domain. In particular, we model a domain
with five ridges and impose a zero-outflow condition at the bottom and top of
the domain, and periodic boundary conditons along the azimuthal directions. To
reduce the influence of end effects, we exclude the top and bottom regions
when calculating the average shear stress on the wall.

As shown in Fig. 5(c)), the simulated $\delta$ values for angled ridges in the
flat plate regime follow a similar trend as the experimental values, with a
$\delta$ that increases monotonically when the angle increases from 0∘
(vertical) to 90∘ (horizontal). However, the experimental values are
systematically above the simulated ones, with an offset of $\sim$ 0.1 mm. It
is not clear what determines this offset; potential factors include partial
wall slip (which would reduce the measured torque), curvature effects, or end
effects.

Figure 5: Experimental flow curves in double logarithmic scale (a) and linear
scale (b) of castor oil measured in angled ridged geometries with ridges at
various $\theta$ with respect to the vertical axis. For all geometries, ridge
spacing ($S$) = 2 mm and ridge depth ($D$) = 1 mm. (c) Variation of $\delta$
as a function of the ridge angle. The simulation results are for the two
limiting cases with the same ridge spacing ($S$) = 2 mm and ridge depth ($D$)
= 1 mm.

##  IV Conclusions

We have introduced textured surfaces in the form of ridged geometries that can
be readily fabricated using 3D printing technique. This way of patterning
gives more freedom to explore an entire range of wall roughness features that
can be easily characterised. These ridged geometries introduce an anisotropic
wall roughness and we have shown that this has an influence on the flow
behaviour in case of simple Newtonian fluids even at low Reynolds number. Due
to the presence of ridges, the flow profiles penetrate in between the ridges
and the extent of this penetration was defined as penetration depth
($\delta$). We have measured how this penetration depth depends on the various
roughness parameters ($D$, $S$, $\theta$) from standard rheological
measurements and we found a general scaling onto an exponential function. The
characteristic flow behaviour of the fluid did not change due to the ridges
but their magnitudes were affected. We also performed complementary simulation
studies to support our experiments.

## Acknowledgments

The authors thank the industrial partners and collaborators for their useful
comments. We thank Bob Mulder for introducing 3D printing in rheology to us.
This work is a part of the Industrial Partnership Programme Controlling
Multiphase Flow (Project No. WP-30-01) that is carried out under an agreement
between Shell, Unilever Research and Development B.V., Evodos and the
Netherlands Organisation for Scientific Research (NWO). This project is co-
funded by the Dutch Research Council (NWO) and TKI-E$\&$I with the
supplementary grant ‘TKI-Toeslag’ for Top Consortia for Knowledge and
Innovation (TKI’s) of the Ministry of Economic Affairs and Climate Policy.
This work is within the framework of the Institute of Sustainable Process
Technology.

## References

  * Bottaro [2019] A. Bottaro, Flow over natural or engineered surfaces: an adjoint homogenization perspective, Journal of Fluid Mechanics 877, P1 (2019).
  * Chen  _et al._ [2012] Y. Chen, C. Zhang, M. Shi, and G. P. Peterson, Slip boundary for fluid flow at rough solid surfaces, Applied Physics Letters 100, 074102 (2012).
  * Gao  _et al._ [2018] W. Gao, X. Zhang, X. Han, and C. Shen, Role of solid wall properties in the interface slip of liquid in nanochannels, Micromachines 9 (2018).
  * Achdou  _et al._ [1998] Y. Achdou, O. Pironneau, and F. Valentin, Effective boundary conditions for laminar flows over periodic rough boundaries, Journal of Computational Physics 147, 187 (1998).
  * Bechert and Bartenwerfer [1989] D. W. Bechert and M. Bartenwerfer, The viscous flow on surfaces with longitudinal ribs, Journal of Fluid Mechanics 206, 105–129 (1989).
  * Bechert  _et al._ [2000] D. W. Bechert, M. Bruse, and W. Hage, Fluid mechanics of biological surfaces and their technological application, Naturwissenschaften 87, 157–171 (2000).
  * Bixler and Bhushan [2013] G. D. Bixler and B. Bhushan, Fluid drag reduction with shark-skin riblet inspired microstructured surfaces, Advanced Functional Materials 23, 4507 (2013).
  * Abdulbari  _et al._ [2013] H. A. Abdulbari, R. Yunus, N. Abdurahman, and A. Charles, Going against the flow—a review of non-additive means of drag reduction, Journal of Industrial and Engineering Chemistry 19, 27 (2013).
  * Raayai-Ardakani and McKinley [2019] S. Raayai-Ardakani and G. H. McKinley, Geometric optimization of riblet-textured surfaces for drag reduction in laminar boundary layer flows, Physics of Fluids 31, 053601 (2019).
  * Djenidi  _et al._ [1989] L. Djenidi, F. Anselmet, J. Liandrat, and L. Fulachier, Numerical and experimental investigation of the laminar boundary layer over riblets, Applied Scientific Research 46, 263 (1989).
  * Soleimani and Eckels [2021] S. Soleimani and S. Eckels, A review of drag reduction and heat transfer enhancement by riblet surfaces in closed and open channel flow, International Journal of Thermofluids 9, 100053 (2021).
  * Cafiero and Iuso [2022] G. Cafiero and G. Iuso, Drag reduction in a turbulent boundary layer with sinusoidal riblets, Experimental Thermal and Fluid Science 139, 110723 (2022).
  * Walsh [1982] M. J. Walsh, Turbulent boundary layer drag reduction using riblets, in  _20th Aerospace Sciences Meeting_ (1982) p. 169.
  * Walsh [1983] M. J. Walsh, Riblets as a viscous drag reduction technique, AIAA Journal 21, 485 (1983).
  * Walsh and Lindemann [1984] M. J. Walsh and A. M. Lindemann, Optimization and application of riblets for turbulent drag reduction, in  _22nd Aerospace Sciences Meeting_ (1984) p. 347.
  * Walsh [1990] M. J. Walsh, Riblets, in  _Viscous Drag Reduction in Boundary Layers_, Vol. 123, edited by D. M. Bushnell and J. N. Hefner (AIAA, 1990) pp. 203–261.
  * Bechert  _et al._ [1997] D. W. Bechert, M. Bruse, W. Hage, J. G. T. Van Der Hoeven, and G. Hoppe, Experiments on drag-reducing surfaces and their optimization with an adjustable geometry, Journal of Fluid Mechanics 338, 59–87 (1997).
  * Davies  _et al._ [2006] J. Davies, D. Maynes, B. W. Webb, and B. Woolford, Laminar flow in a microchannel with superhydrophobic walls exhibiting transverse ribs, Physics of Fluids 18, 087110 (2006).
  * Grüneberger and Hage [2011] R. Grüneberger and W. Hage, Drag characteristics of longitudinal and transverse riblets at low dimensionless spacings, Experiments in Fluids 50, 363 (2011).
  * Raayai-Ardakani and McKinley [2017] S. Raayai-Ardakani and G. H. McKinley, Drag reduction using wrinkled surfaces in high reynolds number laminar boundary layer flows, Physics of Fluids 29, 093605 (2017).
  * Raayai-Ardakani and McKinley [2020] S. Raayai-Ardakani and G. H. McKinley, Geometry mediated friction reduction in taylor-couette flow, Phys. Rev. Fluids 5, 124102 (2020).
  * Owens  _et al._ [2020] C. E. Owens, A. J. Hart, and G. H. McKinley, Improved rheometry of yield stress fluids using bespoke fractal 3d printed vanes, Journal of Rheology 64, 643 (2020).
  * Nickerson and Kornfield [2005] C. S. Nickerson and J. A. Kornfield, A “cleat” geometry for suppressing wall slip, Journal of Rheology 49, 865 (2005).
  * Luchini  _et al._ [1991] P. Luchini, F. Manzo, and A. Pozzi, Resistance of a grooved surface to parallel flow and cross-flow, Journal of Fluid Mechanics 228, 87–109 (1991).
  * Woolford  _et al._ [2009] B. Woolford, D. Maynes, and B. W. Webb, Liquid flow through microchannels with grooved walls under wetting and superhydrophobic conditions, Microfluidics and Nanofluidics 7, 121–135 (2009).
  * Maynes  _et al._ [2007] D. Maynes, K. Jeffs, B. Woolford, and B. W. Webb, Laminar flow in a microchannel with hydrophobic surface patterned microribs oriented parallel to the flow direction, Physics of Fluids 19, 093603 (2007).
  * Djenidi  _et al._ [1994] L. Djenidi, F. Anselmet, J. Liandrat, and L. Fulachier, Laminar boundary layer over riblets, Physics of Fluids 6, 2993 (1994).
  * Faber [1995] T. E. Faber,  _Fluid Dynamics for Physicists_ (Cambridge University Press, 1995).
  * Callaghan [1999] P. T. Callaghan, Rheo-NMR: nuclear magnetic resonance and the rheology of complex fluids, Reports on Progress in Physics 62, 599 (1999).
  * Coussot [2020] P. Coussot, Progress in rheology and hydrodynamics allowed by nmr or mri techniques, Experiments in Fluids 61, 207 (2020).
  * Landau and Lifshitz [1987] L. D. Landau and E. M. Lifshitz, Chapter ii - viscous fluids, in  _Fluid Mechanics (Second Edition)_ (Pergamon, 1987) pp. 44–94.
  * Macosko [1994] C. W. Macosko,  _Rheology: Principles, Measurements, and Applications_ (Wiley-VCH, New York, 1994).
  * Dean and Bhushan [2010] B. Dean and B. Bhushan, Shark-skin surfaces for fluid-drag reduction in turbulent flow: a review, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 368, 4775 (2010).
  * Lee and Lee [2001] S. J. Lee and S. H. Lee, Flow field analysis of a turbulent boundary layer over a riblet surface, Experiments in Fluids 30, 153–166 (2001).


11institutetext: Institute of Astronomy, University of Cambridge, Madingley
Road, Cambridge CB3 0HA, UK
email<EMAIL_ADDRESS>
# Gamma-ray spectra of the Crab, Vela and Geminga pulsars
fitted with SED of the emission from their current sheet

Houshang Ardavan

We show that the spectral energy distribution (SED) of the tightly focused
radiation generated by the superluminally moving current sheet in the
magnetosphere of a non-aligned neutron star fits the gamma-ray spectra of the
Crab, Vela and Geminga pulsars over the entire range of photon energies so far
detected by Fermi-LAT, MAGIC and H.E.S.S. from them: over $10^{2}$ MeV to $20$
TeV. While emblematic of any emission that entails caustics, the SED
introduced here radically differs from those of the disparate emission
mechanisms currently invoked in the literature to fit the data in different
sections of these spectra. We specify, moreover, the connection between the
values of the fit parameters for the analysed spectra and the physical
characteristics of the central neutron stars of the Crab, Vela and Geminga
pulsars and their magnetospheres.

###### Key Words.:

pulsars: individual: J0534+2200, J0835-4510, J0633+1746 – gamma-rays: stars –
stars: neutron – methods: data analysis – radiation mechanisms: non-thermal

##  1 Introduction

It is now firmly established that the magnetosphere of a non-aligned neutron
star contains a rigidly rotating current sheet outside its light cylinder (see
the review article Philippov & Kramer 2022 and the references therein). That
this current sheet moves faster than light in vacuum is not incompatible with
the requirements of special relativity because its superluminally moving
distribution pattern is created by the correlated motion of subluminally
moving charged particles (Bolotovskii & Bykov, 1990). The field generated by a
constituent volume element of this current sheet embraces a synergy between
the superluminal version of the field of synchrotron radiation and the vacuum
version of the field of Čerenkov radiation (Ardavan, 2023, Sect. 2): the waves
emanating from it possess a two-sheeted cusped envelope that spirals outward
around the light cylinder all the way into the far zone (Ardavan, 2023, Figs.
3a and 4).

Only a single wave front propagates past an observer outside the envelope at
any given time. But an observer inside the envelope receives three
intersecting wave fronts, emanating from three distinct retarded positions of
the source element simultaneously (Ardavan, 2023, Figs. 3b and 4). On each
sheet of the envelope, two of these wave fronts coalesce and interfere
constructively to generate an infinitely large radiation field in the time
domain (Ardavan, 2021, Sect. 3.4). On the cusp locus of the envelope, where
all three of these wave fronts (together with the three corresponding images
of the source element) coalesce, the resulting radiation field has even a
higher-order singularity (Ardavan, 2023, Figs. 5 and 6). Divergence of the
field of a source point on the envelope and its cusp, which stems from the
relativistic restrictions inherent in electrodynamics, reflects the fact that
no superluminal source can be point-like. Superposition of the contributions
from the constituent volume elements of the current sheet results in a field
that, while sharply peaked, is finite (Ardavan, 2021, Sect. 4).

The field of the entire volume of the source receives its main contributions
from those volume elements that approach the observer along the radiation
direction with the speed of light and zero acceleration at the retarded time
(Ardavan, 2021, Sect. 4.2). Phases of the waves emitted by such volume
elements of the source possess two nearby stationary points that draw closer
to each other the farther the waves are from their source (Ardavan, 2023, Fig.
7) and eventually coalesce at infinity. As a result, the waves in question
constitute progressively more focused pulses, which by virtue of having time-
domain profiles whose peaks get narrower the farther they propagate away from
their source, embody the high-frequency component of the radiation from the
current sheet (Ardavan, 2021, Sect. 5). The lower-frequency emission in the
radio loud gamma-ray pulsars arises from the isolated stationary points of the
phases of these waves (Ardavan, 2021, Sect. 5.4).

A detailed analysis of the radiation that is generated by the superluminally
moving current sheet in the magnetosphere of a non-aligned neutron star can be
found in Ardavan (2021). A heuristic account of the mathematical results of
that analysis in more transparent physical terms is presented in Ardavan
(2022) and Ardavan (2023, Sect. 2).

We begin here by deriving the spectral energy distribution (SED) of the most
tightly focused component of the radiation that is emitted by the
magnetospheric current sheet (Sect. 2). We then specify the values of the free
parameters in the derived expression for which this SED best fits the data on
the observed gamma-ray spectra of the Crab, Vela and Geminga pulsars (Sect.
3). The specified values of the fit parameters will be used, in conjunction
with the results of the analysis in Ardavan (2021), to determine certain
attributes of the central neutron stars of these pulsars and their
magnetospheres in Sect. 4. The radical departure of the single explanation
given here for the entire breadths of the analysed spectra from those normally
given in terms of disjointed sets of ad hoc spectral distribution functions
(such as simple or broken power-law functions with exponential cutoffs; see,
e.g. Zanin 2017) and disparate emission mechanisms (such as synchro-curvature
processes, inverse Compton scattering, or magnetic reconnection; see, e.g. The
H.E.S.S. Collaboration et al. 2023) is briefly discussed in the final section
of the paper (Sect. 5).

##  2 SED of the caustics generated by the superluminally moving current sheet

The frequency spectrum of the radiation that is generated by the
superluminally moving current sheet in the magnetosphere of a non-aligned
neutron star was presented, in its general form, in Ardavan (2021, Sect. 5.3).
Here we derive the SED of the most tightly-focused component of this radiation
from the general expression given in Eq. (177) of that paper.

In a case where the magnitudes of the vectors denoted by $\boldsymbol{\cal
P}_{l}$ and $\boldsymbol{\cal Q}_{l}$ in Eq. (177) of Ardavan (2021) are
appreciably larger than those of their counterparts, ${\bar{\boldsymbol{\cal
P}}}_{l}$ and ${\bar{\boldsymbol{\cal Q}}}_{l}$, and the dominant contribution
towards the Poynting flux $S_{\nu}$ of the radiation at the observation point
is made by the term corresponding to $l=2$, Eq. (177) of Ardavan (2021) can be
written as

$S_{\nu}=\kappa_{0}\,k^{-2/3}\left|\boldsymbol{\cal P}_{2}\,{\rm
Ai}(-k^{2/3}\sigma_{21}^{2})-{\rm i}k^{-1/3}\boldsymbol{\cal Q}_{2}\,{\rm
Ai}^{\prime}(-k^{2/3}\sigma_{21}^{2})\right|^{2},$ (1)

where Ai and ${\rm Ai}^{\prime}$ are the Airy function and the derivative of
the Airy function with respect to its argument, respectively,
$k=2\pi\nu/\omega$ is the frequency $\nu$ of the radiation in units of the
rotation frequency $\omega/2\pi$ of the central neutron star, and $\kappa_{0}$
and $\sigma_{21}$ are two positive scalars. In the high-frequency regime $k\gg
1$, the coefficients of the Airy functions in the above expression stand for
the complex vectors $\boldsymbol{\cal P}_{2}=k^{-1/2}\boldsymbol{\cal
P}_{2}^{(2)}$ and $\boldsymbol{\cal Q}_{2}=k^{-1/2}\boldsymbol{\cal
Q}_{2}^{(2)}$ in which $\boldsymbol{\cal P}_{2}^{(2)}$ and $\boldsymbol{\cal
Q}_{2}^{(2)}$ are defined by Eqs (138)-(146) of Ardavan (2021). The variable
$\sigma_{21}$ has a vanishingly small value at those privileged colatitudes
(relative to the star’s spin axis) where the high-frequency radiation is
observable (Ardavan, 2021, Sect. 4.5). For the purposes of the analysis in
this paper, we may therefore replace $\boldsymbol{\cal P}_{2}^{(2)}$ and
$\boldsymbol{\cal Q}_{2}^{(2)}$ by their limiting values for $k\gg 1$ and
$0\leq\sigma_{21}\ll 1$ and treat them as constant parameters.

Multiplying Eq. (1) by the radiation frequency $\nu=\omega k/2\pi$ and
expressing $\boldsymbol{\cal P}_{2}$ and $\boldsymbol{\cal Q}_{2}$ in terms of
$\boldsymbol{\cal P}_{2}^{(2)}$ and $\boldsymbol{\cal Q}_{2}^{(2)}$, we find
that the SED of this emission is given by

$\displaystyle\nu S_{\nu}$ $\displaystyle=$
$\displaystyle\frac{\omega\kappa_{0}}{2\pi}k^{-2/3}\left|\boldsymbol{\cal
P}^{(2)}_{2}\,{\rm Ai}(-k^{2/3}\sigma_{21}^{2})\right.$ (2)
$\displaystyle\left.-{\rm i}k^{-1/3}\boldsymbol{\cal Q}^{(2)}_{2}\,{\rm
Ai}^{\prime}(-k^{2/3}\sigma_{21}^{2})\right|^{2}.$

Evaluation of the right-hand side of Eq. (2) results in

$\displaystyle\nu S_{\nu}$ $\displaystyle=$
$\displaystyle\kappa_{1}\,k^{-2/3}\left[{\rm
Ai}^{2}(-k^{2/3}\sigma_{21}^{2})+\zeta_{1}^{2}k^{-2/3}{{\rm
Ai}^{\prime}}^{2}(-k^{2/3}\sigma_{21}^{2})\right.$ (3)
$\displaystyle\left.+2\zeta_{1}\cos\beta\,k^{-1/3}{\rm
Ai}(-k^{2/3}\sigma_{21}^{2}){\rm
Ai}^{\prime}(-k^{2/3}\sigma_{21}^{2})\right],$

where

$\kappa_{1}=\frac{\omega\kappa_{0}}{2\pi}\left|\boldsymbol{\cal
P}^{(2)}_{2}\right|^{2},\,\,\,\zeta_{1}=\frac{\left|{\boldsymbol{\cal
Q}^{(2)}_{2}}\right|}{\left|{\boldsymbol{\cal
P}^{(2)}_{2}}\right|},\,\,\,\cos\beta=\frac{\Im\left(\boldsymbol{\cal
Q}^{(2)}_{2}\cdot\boldsymbol{\cal P}^{(2)*}_{2}\right)}{\left|\boldsymbol{\cal
Q}^{(2)}_{2}\right|\left|\boldsymbol{\cal P}^{(2)}_{2}\right|},$ (4)

and $\Im{}$ and $*$ denote an imaginary part and the complex conjugate,
respectively. The above spectrum is emblematic of any radiation that entails
caustics (see Stamnes, 1986).

To take account of the fact that the parameter $\sigma_{21}$ assumes a non-
zero range of values across the (non-zero) latitudinal width of the detected
radiation beam (Ardavan, 2021, Sect. 4.5), we must integrate $\nu S_{\nu}$
with respect to $\sigma_{21}$ over a finite interval
$\rho\sigma_{0}\leq\sigma_{21}\leq\sigma_{0}$ with $\sigma_{0}\ll 1$ and
$0\leq\rho<1$. Performing the integration of the Airy functions in Eq. (3)
with respect to $\sigma_{21}$ by means of Mathematica, we thus obtain

$\displaystyle{\cal S}_{\nu}$ $\displaystyle=$
$\displaystyle\nu\int_{\rho\sigma_{0}}^{\sigma_{0}}S_{\nu}\,{\rm
d}\sigma_{21}$ (5) $\displaystyle=$
$\displaystyle\kappa\left[f_{1}(\chi,\rho)+\frac{\zeta^{2}}{4\sqrt{3}}f_{2}(\chi,\rho)-\frac{\zeta\cos\beta}{2\sqrt{3}}f_{3}(\chi,\rho)\right],$

where

$\kappa=\left(\frac{\sigma_{0}}{\sqrt{3\pi}}\right)^{3}\kappa_{1},\qquad\zeta=\sigma_{0}\zeta_{1},\qquad\chi=\frac{2}{3}{\sigma_{0}}^{3}k,$
(6) $\displaystyle f_{1}$ $\displaystyle=$
$\displaystyle\bigg{[}3\Gamma\left(\frac{7}{6}\right)\eta\chi^{-2/3}{}_{2}F_{3}\left(\begin{matrix}1/6&1/6&{}\\\
1/3&2/3&7/6\end{matrix};-\eta^{6}\chi^{2}\right)$
$\displaystyle+\pi^{1/2}\eta^{3}{}_{2}F_{3}\left(\begin{matrix}1/2&1/2&{}\\\
2/3&4/3&3/2\end{matrix}\,;-\eta^{6}\chi^{2}\right)$
$\displaystyle+\frac{9}{20}\Gamma\left(\frac{5}{6}\right)\eta^{5}\chi^{2/3}{}_{2}F_{3}\left(\begin{matrix}5/6&5/6&{}\\\
4/3&5/3&11/6\end{matrix}\,;-\eta^{6}\chi^{2}\right)\bigg{]}_{\eta=\rho}^{\eta=1},$
$f_{2}=\eta\chi^{-4/3}\,{}_{24}G^{31}\left(-\eta^{2}\chi^{2/3},\frac{1}{3}\left|\,\begin{matrix}5/6&7/6&{}&{}\\\
0&2/3&4/3&-1/6\end{matrix}\right)\right.\bigg{|}_{\eta=\rho}^{\eta=1},$ (8)
$f_{3}=\eta\chi^{-1}\,{}_{24}G^{31}\left(-\eta^{2}\chi^{2/3},\frac{1}{3}\left|\,\begin{matrix}5/6&1/2&{}&{}\\\
0&1/3&2/3&-1/6\end{matrix}\right)\right.\bigg{|}_{\eta=\rho}^{\eta=1},$ (9)

and ${}_{2}F_{3}$ and ${}_{24}G^{31}$ are respectively the generalised
hypergeometric function (see Olver et al., 2010) and the generalised Meijer
G-Function 111 https://mathworld.wolfram.com/MeijerG-Function.html. The
variable $\chi$ that appears in the above expressions is related to the
frequency $\nu$ of the radiation via $\chi=4\pi{\sigma_{0}}^{3}\nu/(3\omega)$.

The scale and shape of the SED given in Eq. (5) depend on the five parameters
$\kappa$, $\zeta$, $\beta$, $\sigma_{0}$ and $\rho$: parameters whose values
are dictated by the characteristics of the magnetospheric current sheet (see
Sect. 4). The parameters $\zeta$, $\beta$ and $\rho$ determine the shape of
the energy distribution while the parameters $\kappa$ and $\sigma_{0}$
determine the position of this distribution along the energy-flux and photon-
energy axes, respectively.

Figure 1: Phase-averaged gamma-ray spectrum of the Crab pulsar. The FERMI-LAT
data points below $25$ GeV (coloured black) are those reported in Abdo et al.
(2010). The MAGIC data points between $25$ and $100$ GeV (coloured blue) are
those reported in Aleksić et al. (2011). The MAGIC data points beyond $100$
Gev (coloured green) are extracted, by the procedure described in Ardavan
(2023, Sect. 4), from those reported by Ansoldi et al. (2016). The curve
(coloured red) is a plot of the SED described by Eq. (5) for the parameters
given in Eq. (10).

##  3 Fits to the gamma-ray data on the SEDs of the Crab, Vela and Geminga
pulsars

In this section we use Mathematica’s ‘NonlinearModelFit’ procedure222
https://reference.wolfram.com/language/ref/NonlinearModelFit.html and the
statistical information that it provides to determine the values of the fit
parameters and their standard errors. Where, owing to the complexity of the
expression in Eq. (5), this procedure breaks down and so the fits to the data
are obtained by trial and error, only the values of these parameters are
specified. The fit residuals in the case of each pulsar turn out to be smaller
than the corresponding observational errors for almost all data points (see
Figs. 1–3). Since there are no two values of any of the fit parameters for
which the present SED has the same shape and position, the specified values of
the fit parameters are moreover unique.

###  3.1 Spectrum of the Crab pulsar (PSR J0534+2200)

Figure 1 shows the data on the energy spectrum of the phase-averaged gamma-ray
emission from the Crab pulsar (Abdo et al., 2010; Aleksić et al., 2011;
Ansoldi et al., 2016) and a plot of the SED described by Eq. (5) that fits
them best. The parameters for which the function ${\cal S}_{\nu}(\chi)$
(depicted by the red curve in this figure) is plotted have the following
values and standard errors:

$\displaystyle\kappa$ $\displaystyle=$ $\displaystyle(8.15\pm 0.75)\times
10^{-11}\quad{\rm erg}\,{\rm cm}^{-2}\,{\rm s}^{-1},$ $\displaystyle\chi$
$\displaystyle=$ $\displaystyle(0.576\pm 0.095)\times(h\nu)_{\rm GeV},$
$\displaystyle\sigma_{0}$ $\displaystyle=$ $\displaystyle(4.76\pm 2.61)\times
10^{-8},\,\,\,\zeta=0.394\pm 0.076,$ $\displaystyle\beta$ $\displaystyle=$
$\displaystyle 0.258\pm 0.249,\quad\rho=0,$ (10)

where $h$ stands for the Planck constant and $2\pi/\omega$ has been set equal
to the period of the Crab pulsar ($0.033$ s). Figure 1 was first reported in
Ardavan (2023) but with an incorrect label on its vertical axis.

As the photon energy $h\nu$ increases beyond $1.2$ TeV, the slope of the curve
in Fig. 1 continues to decrease at the relatively slow rate in which it
decreases past $50$ GeV. Thus the upper limits given by Ansoldi et al. (2016)
for the flux density of gamma-rays whose energies exceed $1.2$ TeV all lie
above the continuation of the curve shown in Fig. 1.

Figure 2: Gamma-ray spectrum of the Vela pulsar. The FERMI-LAT data points
below $10^{2}$ GeV (coloured black) are those reported in Abdo et al. (2013).
The H.E.S.S. data above $260$ GeV (shown as the magenta area between $1\sigma$
uncertainty bands and the magenta upper limits) are those reported in The
H.E.S.S. Collaboration et al. (2023). The curve (coloured red) is a plot of
the SED described by Eq. (5) for the parameters given in Eqs. (11) and (12).

###  3.2 Spectrum of the Vela pulsar (PSR J0835-4510)

The data on the Vela pulsar both in the GeV band detected by Fermi-LAT (Abdo
et al., 2013) and in the TeV band detected by H.E.S.S. (The H.E.S.S.
Collaboration et al., 2023) are shown in Fig. 2. The SED that best fits these
data (represented by the red curve in Fig. 2) is described by the expression
in Eq. (5) for the following values of its free parameters:

$\displaystyle\kappa$ $\displaystyle=$ $\displaystyle(3.31\pm 0.047)\times
10^{-10}\quad{\rm erg}\,{\rm cm}^{-2}\,{\rm s}^{-1},$ $\displaystyle\chi$
$\displaystyle=$ $\displaystyle(0.436\pm 0.071)\times(h\nu)_{\rm GeV},$
$\displaystyle\sigma_{0}$ $\displaystyle=$ $\displaystyle(3.12\pm 1.70)\times
10^{-8},\quad\zeta=0.518\pm 0.065,$ $\displaystyle\beta$ $\displaystyle=$
$\displaystyle 0.244\pm 0.073,\quad\rho=0.409\pm 0.198,$ (11)

over the range $0.1\leq(h\nu)_{\rm GeV}\leq 10$ and

$\displaystyle\kappa$ $\displaystyle=$ $\displaystyle 7.94\times
10^{-13}\quad{\rm erg}\,{\rm cm}^{-2}\,{\rm s}^{-1},$ $\displaystyle\chi$
$\displaystyle=$ $\displaystyle 2.51\times 10^{-6}\times(h\nu)_{\rm GeV},$
$\displaystyle\sigma_{0}$ $\displaystyle=$ $\displaystyle 5.59\times
10^{-10},\,\,\zeta=0.203,\,\,\beta=0.105,\,\,\rho=0.99,$ (12)

over the range $10\leq(h\nu)_{\rm GeV}\leq 2\times 10^{4}$ of photon energies.

The difference between the values of the free parameters in the above two
ranges of photon energies stems from the fact that the higher-frequency
radiation is appreciably more focused than its lower-frequency counterpart in
this case. This is reflected not only in the much lower value of $\sigma_{0}$
in Eq. (12) compared to that in Eq. (11), but also in the higher value of
$\rho$ (i.e. the shorter range of values of $\sigma_{21}$): the smaller the
value of $\sigma_{21}$, the closer to each other are the stationary points of
the phases of the emitted waves and so the more focused is the radiation
(Ardavan, 2021, Sect. 4.5).

Figure 3: Gamma-ray spectrum of the Geminga pulsar. The FERMI-LAT data points
below $35$ GeV (coloured black) are those reported in Abdo et al. (2013). The
MAGIC data points above $35$ GeV (coloured blue) are those reported in The
MAGIC Collaboration et al. (2020). The curve (coloured red) is a plot of the
SED described by Eq. (5) for the parameters given in Eqs. (13) and (14).

###  3.3 Spectrum of the Geminga pulsar (PSR J0633+1746)

The Fermi-LAT (Abdo et al., 2013) and MAGIC (The MAGIC Collaboration et al.,
2020) data on the gamma-ray emission from the Geminga pulsar are shown in Fig.
3. The fit to these data (the red curve in Fig 3) is described by the
expression in Eq. (5) for

$\displaystyle\kappa$ $\displaystyle=$ $\displaystyle(1.47\pm 0.23)\times
10^{-10}\quad{\rm TeV}\,{\rm cm}^{-2}\,{\rm s}^{-1},$ $\displaystyle\chi$
$\displaystyle=$ $\displaystyle(0.392\pm 0.038)\times(h\nu)_{\rm GeV},$
$\displaystyle\sigma_{0}$ $\displaystyle=$ $\displaystyle(2.17\pm 0.99)\times
10^{-8},\quad\zeta=0.646\pm 0.028$ $\displaystyle\beta$ $\displaystyle=$
$\displaystyle 0.194\pm 0.017,\quad\rho=0.568\pm 0.131,$ (13)

over the range $0.1\leq(h\nu)_{\rm GeV}\leq 7.94$ and

$\displaystyle\kappa$ $\displaystyle=$ $\displaystyle 1.01\times
10^{-13}\quad{\rm TeV}\,{\rm cm}^{-2}\,{\rm s}^{-1},$ $\displaystyle\chi$
$\displaystyle=$ $\displaystyle 3.98\times 10^{-6}\times(h\nu)_{\rm GeV},$
$\displaystyle\sigma_{0}$ $\displaystyle=$ $\displaystyle 4.71\times
10^{-10},\,\,\zeta=0.101,\,\,\beta=0.071,\,\,\rho=0.511,$ (14)

over the range $7.94\leq(h\nu)_{\rm GeV}\leq 63.1$ of photon energies. As in
the case of the spectrum of the Vela pulsar, the change in the values of the
fit parameters across $(h\nu)_{\rm GeV}=7.94$ reflects the fact that the
higher-frequency radiation is appreciably more focused than its lower-
frequency counterpart.

##  4 The connection between the parameters of the fitted spectra and the
physical characteristics of their sources

From Eqs. (6), (4) and (1), it follows that the parameters $\kappa$, $\zeta$,
$\rho$ and $\sigma_{0}$ in Eq. (5) are related to the characteristics of the
source of the observed radiation via the quantities $\omega$, $\kappa_{0}$,
$\boldsymbol{\cal P}_{2}^{(2)}$ and $\boldsymbol{\cal Q}_{2}^{(2)}$ that
appear in the expression for the SED ${\cal S}_{\nu}$. In this section we use
the results of the analysis presented in Ardavan (2021) to express these
quantities in terms of the inclination angle of the central neutron star,
$\alpha$, the magnitude of the star’s magnetic field at its magnetic pole,
$B_{0}=10^{12}{\hat{B}}_{0}$ Gauss, the radius of the star, $r_{s0}=10^{6}d$
cm, the rotation frequency of the star, $\omega=10^{2}{\hat{P}}^{-1}$ rad/s,
and the spherical polar coordinates, $R_{P}=D\,\,{\rm kpc}=3.085\times
10^{21}D$ cm, $\varphi_{P}$ and $\theta_{P}$, of the observation point $P$ in
a frame whose centre and $z$-axis coincide with the centre and spin axis of
the star.

When $\sigma_{21}\ll 1$ and so the quantities $\tau_{2{\rm min}}$ and
$\tau_{2{\rm max}}$ that appear in Eqs. (177), (138)–(146) and (136) of
Ardavan (2021) are approximately equal, these equations imply that

$\kappa_{0}=4.15\times
10^{18}w_{1}^{2}w_{3}^{2}{\hat{B}}_{0}^{2}d^{4}D^{-2}{\hat{P}}^{-1}\quad{\rm
Jy},$ (15)

and

$\left|\boldsymbol{\cal P}_{2}^{(2)}\right|=\frac{b}{3}\left(\frac{2}{\pi
a^{3}}\right)^{1/2}\left(\frac{2\sigma_{21}}{\partial^{2}f_{2C}/\partial\tau^{2}|_{\tau=\tau_{2{\rm
min}}}}\right)^{1/2}\left|{\tilde{\bf P}}_{2}\right|,$ (16)

where

$w_{1}=\left|1-2\alpha/\pi\right|,\qquad w_{3}=1+0.2\sin^{2}\alpha,$ (17)
$\displaystyle f_{2C}$ $\displaystyle=$
$\displaystyle({\hat{r}}_{P}^{2}{\hat{r}}_{sC}^{2}\sin^{2}\theta-1)^{1/2}-{\hat{R}}_{P}-\arccos\left({\hat{r}}_{P}^{-1}{\hat{r}}_{sC}^{-1}\csc\theta\right)$
(18)
$\displaystyle-\arccos(\cot\alpha\cot\theta)+{\hat{r}}_{sC}+\varphi_{P}-{\hat{r}}_{s0}-2\pi,$
$\displaystyle a$ $\displaystyle=$
$\displaystyle({\hat{r}}_{P}^{2}{\hat{r}}_{sC}^{2}\sin^{2}\theta-1)\left[\frac{({\hat{r}}_{P}^{2}{\hat{r}}_{sC}^{2}\sin^{2}\theta-1)^{1/2}+{\hat{r}}_{sC}}{{\hat{r}}_{sC}({\hat{r}}_{P}^{2}-1)^{1/2}({\hat{R}}_{P}^{2}\sin^{2}\theta-1)^{1/2}}\right.$
(19)
$\displaystyle\left.-\frac{1}{({\hat{r}}_{P}^{2}{\hat{r}}_{sC}^{2}\sin^{2}\theta-1)^{1/2}}\right],$
$b=\frac{{\hat{r}}_{P}^{2}{\hat{r}}_{sC}^{2}\sin^{2}\theta-1}{{\hat{r}}_{sC}({\hat{r}}_{P}^{2}-1)^{1/2}({\hat{R}}_{P}^{2}\sin^{2}\theta-1)^{1/2}},$
(20)
${\hat{r}}_{sC}=\frac{({\hat{r}}_{P}^{2}-1)^{1/2}({\hat{R}}_{P}^{2}\sin^{2}\theta-1)^{1/2}-{\hat{z}}_{P}\cos\theta}{{\hat{r}}_{P}^{2}\sin^{2}\theta-1},$
(21) $\sigma_{21}=[{\textstyle\frac{3}{4}}(f_{2C}|_{\tau=\tau_{2{\rm
max}}}-f_{2C}|_{\tau=\tau_{2{\rm min}}})]^{1/3},$ (22)
$\theta=\arccos(\sin\alpha\cos\tau),$ (23) $\displaystyle\left|{\tilde{\bf
P}}_{2}\right|$ $\displaystyle=$
$\displaystyle\cos\alpha[{\hat{r}}_{sC}^{2}(1+\cos^{2}\theta_{P}-2\cos\theta\cos\theta_{P})+\tan^{2}\alpha$
$\displaystyle-\cot^{2}\theta+2{\hat{r}}_{sC}\sin\theta\cos\theta_{P}(\tan^{2}\alpha-\cot^{2}\theta)^{1/2}]^{1/2},\quad\,\,$

and the variables $\tau_{2{\rm min}}$ and $\tau_{2{\rm max}}$ stand for the
minimum and maximum of the function $f_{2C}$ (see also Eqs. (7), (9), (174),
(175), (93)–(95), (97) and (88) of Ardavan 2021). The caret on $R_{P}$ and
$r_{s0}$ (and $r_{P}=R_{P}\sin\theta_{P}$, $z_{P}=R_{P}\cos\theta_{P}$) is
used here to designate a variable that is rendered dimensionless by being
measured in units of the light-cylinder radius $c/\omega$. (Note the following
two corrections: the vector ${\bf P}_{2}$ in Eq. (145) and the numerical
coefficient $1.54\times 10^{19}$ in Eq. (177) of Ardavan 2021 have been
corrected to read ${\tilde{\bf P}}_{2}$ and $4.15\times 10^{18}$,
respectively.)

The expression for $|{\tilde{\bf P}}_{2}|$ in Eq. (LABEL:E24) is derived from
that for ${\bf P}_{2}$ given by Eqs. (98), (80), (78) and (62) of Ardavan
(2021). In this derivation, we have set the observation point on the cusp
locus of the bifurcation surface where ${\hat{r}}_{s}={\hat{r}}_{sC}$, have
approximated $(p_{1},p_{2},p_{3})$ by its far-field value
$2^{1/3}{\hat{R}}_{P}^{-1}({\hat{R}}_{P}^{-1},-1,1)$ and have let ${\bf
P}_{2}={\hat{R}}_{P}^{-1}{\tilde{\bf P}}_{2}$. The factor ${\hat{R}}_{P}^{-2}$
that would have otherwise appeared in the resulting expression for
$|\boldsymbol{\cal P}^{(2)}_{2}|^{2}$ is thus incorporated in the factor
$\kappa_{0}$ in Eq. (15).

For certain values of $\theta_{P}$, denoted by $\theta_{P2S}$, the function
$f_{2C}(\tau)$ has an inflection point (see Ardavan 2021, Sect. 4.4, and
Ardavan 2023, Sect. 2). For any given inclination angle $\alpha$, the position
$\tau_{2S}$ of this inflection point and the colatitude $\theta_{P2S}$ of the
observation points for which $f_{2C}(\tau)$ has an inflection point follow
from the solutions to the simultaneous equations $\partial
f_{2C}/\partial\tau=0$ and $\partial^{2}f_{2C}/\partial\tau^{2}=0$. (Explicit
expressions for the derivatives that appear in these equations can be found in
Appendix A of Ardavan 2021.) For values of $\theta_{P}$ close to
$\theta_{P2S}$, the separation between the maximum $\tau=\tau_{2{\rm max}}$
and minimum $\tau=\tau_{2{\rm min}}$ of $f_{2C}$ and hence the value of
$\sigma_{21}$ are vanishingly small. Since the high-frequency radiation that
arises from the current sheet is detectable only in the vicinity of the
latitudinal direction $\theta_{P}=\theta_{P2S}$ (see Fig. 4), here we have
evaluated the variables that appear in Eqs. (16)–(LABEL:E24) for
$\sigma_{21}\ll 1$. The ratio appearing inside the second pair of parentheses
in Eq. (16), which becomes indeterminate as $\sigma_{21}$ tends to zero,
approaches a finite value in this limit.

The dimensionless frequency $k=2\pi\nu/\omega$ of a $1$ GeV photon has the
value $8\times 10^{21}$ in the case of the Crab pulsar. In contrast, as
indicated by the small values of the fit parameter $\sigma_{0}$ encountered in
Sect. 3, the value of the variable $\chi=2\sigma_{0}^{3}k/3$ that appears in
the description of the SED ${\cal S}_{\nu}$ is of the order of $10^{-3}$ for a
$1$ GeV gamma-ray photon. Because any measurement of the spectrum described by
Eq. (5) is carried out by counting the number of detected gamma-ray photons in
a bin covering an energy interval of a few MeV or GeV, the widths of the
equivalent frequency bins involved in such measurements are of the order of
$\Delta k=3\sigma_{0}^{-3}\Delta\chi/2$ with
$10^{-4}\lesssim\Delta\chi\lesssim 10^{-2}$.

Once the the expression for the SED ${\cal S}_{\nu}$ is accordingly multiplied
by $3\sigma_{0}^{-3}\Delta\chi/2$, Eqs. (6), (4), (15) and (16) jointly yield

$\kappa=2.42\times
10^{-6}{\hat{B}}_{0}^{2}d^{4}D^{-2}{\hat{P}}^{-2}\Delta\chi\,{\hat{\kappa}}_{\rm
th}\qquad{\rm erg}\,\,{\rm s}^{-1}{\rm cm}^{-2},$ (25)

in which

${\hat{\kappa}}_{\rm
th}=w_{1}^{2}w_{3}^{2}\frac{b^{2}}{a^{3}}\left(\frac{2\sigma_{21}}{\partial^{2}f_{2C}/\partial\tau^{2}|_{\tau=\tau_{2{\rm
min}}}}\right)\left|{\tilde{\bf P}}_{2}\right|^{2},$ (26)

and Jy has been expressed in terms of erg s-1 cm-2 Hz-1. Equation (25) can be
written as ${\hat{\kappa}}_{\rm th}={\hat{\kappa}}_{\rm obs}$, where

${\hat{\kappa}}_{\rm obs}=4.13\times
10^{5}({\hat{B}}_{0}d^{2})^{-2}(\Delta\chi)^{-1}D^{2}{\hat{P}}^{2}\kappa.$
(27)

While ${\hat{\kappa}}_{\rm obs}$ only contains the observed parameters of the
pulsar and its emission, the value of ${\hat{\kappa}}_{\rm th}$ is determined
by the physical characteristics of the magnetospheric current sheet that acts
as the source of the observed emission. When, as in Fig. 20 of Ardavan (2021),
the value of the colatitude $\theta_{P}$ of the observation point is
sufficiently close to that of the critical angle $\theta_{P2S}$ for
$|\tau_{2{\rm max}}-\tau_{2{\rm min}}|$ and hence $\sigma_{21}$ to be much
smaller than unity, the right-hand side of Eq. (26) is a function of the
inclination angle $\alpha$ and the observer’s distance ${\hat{R}}_{P}$ only.

Equation (25) thus enables us to connect the the parameters of the fitted
spectra to the physical characteristics of their sources.

Figure 4: The dependence on the star’s inclination angle $\alpha$ of the
critical colatitude $\theta_{P2S}$ along which the high-frequency radiation
beam generated by the superluminally moving current sheet propagates.

Figure 5: Dependence of the ${\hat{B}}_{0}d^{2}$ on $\alpha$ for the Crab
pulsar, where ${\hat{B}}_{0}$ is the magnitude of the star’s magnetic field at
its magnetic pole in units of $10^{12}$ Gauss, $d$ is the radius of the star
in units of $10^{6}$ cm and $\alpha$ is the angle between the magnetic and
rotation axes of the star in degrees.

###  4.1 Characteristics of the source of the Crab pulsar’s emission

Once the value of the fit parameter $\kappa$ given in Eq. (10), the period
($0.033$ s) and the distance ($2$ kpc) of the Crab pulsar are inserted in Eq.
(27), the resulting value of ${\hat{\kappa}}_{\rm obs}$ can be equated to
${\hat{\kappa}}_{\rm th}$ to obtain

${\hat{B}}_{0}d^{2}=1.93\times 10^{-1}{\hat{\kappa}}_{\rm th}^{-1/2}$ (28)

for $\Delta\chi=10^{-3}$. The dependence of ${\hat{B}}_{0}d^{2}$ on $\alpha$
implied by this relationship, in conjunction with Eq. (26), is plotted in Fig.
5.

One can infer from the separation between the two peaks in the gamma-ray light
curve of the Crab pulsar ($\sim 0.4$ cycles; Abdo et al., 2010) that the
inclination angle $\alpha$ of this pulsar has a value lying in the range
$60^{\circ}$ to $70^{\circ}$ (Ardavan, 2021, Sect. 5.4). This in turn implies
that the value of ${\hat{B}}_{0}d^{2}$ lies between $1.05$ and $2.61$ in this
case (see Fig. 5), i.e. that the magnitude of the star’s magnetic field at its
magnetic pole lies in the range $1.05\times 10^{12}$ to $2.61\times 10^{12}$
Gauss if we assume that the star’s radius has the value $10^{6}$ cm. These are
of the same order of magnitude as the value of $B_{0}$ obtained from the
conventional formula for magnetic dipole radiation (i.e. $3.79\times 10^{12}$
Gauss; Manchester et al. 2005).

The corresponding value of the privileged latitudinal direction $\theta_{P2S}$
in which the high-frequency radiation from the Crab pulsar can be observed
lies between $90^{\circ}$ and $117^{\circ}$ according to Fig. 4. The range of
values of the inclination angle $\alpha$ determines also the scale factor
$j_{0}$ of the electric current density in the magnetosphere of the pulsar:
from Eq. (22) of Ardavan (2021) it follows that

$j_{0}=1.29\times 10^{8}{\hat{B}}_{0}d^{2}{\hat{P}}^{-3}w_{1}w_{3},\qquad{\rm
statamp\,cm}^{-2}$ (29)

a scale factor whose value lies between $3.55\times 10^{8}$ and $6.01\times
10^{8}$ statamp cm-2 in the case of the Crab pulsar.

The above results supersede those reported in Ardavan (2023, Sect. 4).

###  4.2 Characteristics of the source of the Vela pulsar’s emission

Properties of the central neutron star of the Vela pulsar and its
magnetosphere are more accurately reflected in the parameters of its spectrum
over the photon energies $0.1$ to $10$ GeV, i.e. the interval that encompasses
the majority of the data points plotted in Fig. 2. Insertion of the value of
the fit parameter $\kappa$ given in Eq. (11), the period ($0.089$ s) and the
distance ($287$ pc) of this pulsar in Eq. (27), together with
${\hat{\kappa}}_{\rm obs}={\hat{\kappa}}_{\rm th}$, results in

${\hat{B}}_{0}d^{2}=1.50\times 10^{-1}{\hat{\kappa}}_{\rm th}^{-1/2}$ (30)

for $\Delta\chi=10^{-3}$. Thus the dependence of ${\hat{B}}_{0}d^{2}$ on
$\alpha$ is given by a shifted version of the curve in Fig. 5: one that is
vertically lowered by 0.109 (cf. Eqs. (28) and (30)).

If we assume that the radius of the central neutron star of this pulsar has
the value $10^{6}$ cm, i.e. $d=1$, then the value $3.38\times 10^{12}$ Gauss
of its magnetic field $B_{0}$ (inferred from the conventional formula for
magnetic dipole radiation; Manchester et al. 2005) corresponds, according to
the modified version of Fig. 5, to the inclination angle $\alpha=70^{\circ}$.
This would in turn imply, in conjunction with Eq. (29) and Fig. 4, that the
scale factor $j_{0}$ of the electric current density has the value $4.02\times
10^{7}$ statamp cm-2 and the latitudinal direction along which the pulsar is
observed is given by $\theta_{P2S}=117^{\circ}$.

That the main peaks of the gamma-ray light curve of the Vela pulsar are
separated by $\sim 0.43$ cycles (The H.E.S.S. Collaboration et al., 2023) is
consistent with an inclination angle close to $70^{\circ}$ (see Ardavan, 2021,
Sect. 5.4).

###  4.3 Characteristics of the source of the Geminga pulsar’s emission

The value of $\kappa$ given in Eq. (13), which applies to a wider range of
photon energies than that given in Eq. (14), together with the period ($0.237$
s) and the distance ($190$ pc) of the Geminga pulsar yield

${\hat{B}}_{0}d^{2}=2.24\times 10^{-1}{\hat{\kappa}}_{\rm th}^{-1/2}$ (31)

for $\Delta\chi=10^{-3}$ (see Eq. (27) and recall that ${\hat{\kappa}}_{\rm
obs}={\hat{\kappa}}_{\rm th}$). The dependence of ${\hat{B}}_{0}d^{2}$ on
$\alpha$ is therefore given by a set of points whose only difference with
those plotted in Fig. 5 is that their vertical coordinates are shifted upward
by $0.065$ (cf. Eqs. (28) and (31)).

Given that the two peaks of the gamma-ray light curve of this pulsar are
separated by $\sim 0.5$ cycle (The MAGIC Collaboration et al., 2020), the
inclination angle of its central neutron star must be about $60^{\circ}$
(Ardavan, 2021, Sect. 5.4). This, together with the modified version of Fig.
5, implies the value ${\hat{B}}_{0}\simeq 1.22$, i.e. $B_{0}=1.22\times
10^{12}$ Gauss, for the magnetic field of the star at its pole (when $d=1$), a
value that is not too different from that predicted by the conventional
expression for magnetic dipole radiation, i.e. $1.63\times 10^{12}$ Gauss
(Manchester et al., 2005). Equation (29) and Fig. 4 respectively yield
$j_{0}=1.12\times 10^{6}$ statamp cm-2 and $\theta_{P2S}=90^{\circ}$ for the
scale factor of the electric current density and the colatitude of the
observation point in this case.

##  5 Concluding remarks

An up to date discussion of the emission mechanisms that have so far been
proposed for explaining the data on the observed gamma-ray spectra of pulsars
can be found in The H.E.S.S. Collaboration et al. (2023). Even in the most
favourable circumstances, spectra of the emissions from relativistic charged
particles that are accelerated by synchro-curvature, inverse Compton or
magnetic reconnection processes fit the observed data only over disjoint
limited ranges of photon energies. In particular, the constraints set by the
$20$ TeV component of the radiation from the Vela pulsar “provide
unprecedented challenges to the state-of-the-art models of HE [high energy]
and VHE [very high energy] emission from pulsars” (The H.E.S.S. Collaboration
et al., 2023).

In contrast, the SED of the tightly focused caustics that are generated by the
superluminally moving current sheet in the magnetosphere of a non-aligned
neutron star fits the observed spectra of the Crab, Vela and Geminga pulsars
over the entire range of photon energies so far detected by Fermi-LAT, MAGIC
and H.E.S.S. from them: over $10^{2}$ MeV to $20$ TeV. Not only are these data
fitted by a single spectral distribution function but they are also explained
by a single emission mechanism, an emission mechanism that has already
accounted for a number of other salient features of the radiation received
from pulsars: its brightness temperature, polarization, radio spectrum, and
profile with microstructure and with a phase lag between the radio and gamma-
ray peaks (Ardavan, 2021, 2022).

Two final remarks concerning the analysis that has led to the present SED are
in order:

(i)

It is often presumed that the plasma equations used in the numerical
simulations of the magnetospheric structure of an oblique rotator should, at
the same time, predict any radiation that the resulting structure would be
capable of emitting (Spitkovsky, 2006; Kalapotharakos et al., 2012). This
presumption stems from ignoring the role of boundary conditions in the
solution of Maxwell’s equations. The boundary conditions with which the
structure of the pulsar magnetosphere is computed are radically different from
the boundary conditions with which the retarded solution of these equations
(i.e. the solution describing the radiation from the charges and currents in
the magnetosphere) is derived (see the last paragraph in Sect. 6 of Ardavan
2021).

(ii)

Thickness of the current sheet, which sets an upper limit on the frequency of
the present radiation, is dictated by microphysical processes that are not
well understood: the standard Harris solution of the Vlasov-Maxwell equations
(Harris, 1962) that is commonly used in analysing a current sheet is not
applicable in the present case because the current sheet in the magnetosphere
of a non-aligned neutron star moves faster than light and so has no rest
frame. Even in stationary or subluminally moving cases, there is no consensus
on whether equilibrium current sheets in realistic geometries have finite or
zero thickness (Klimchuk et al., 2023). The fact that the SED described by Eq.
(5) yields such good fits to the gamma-ray spectra of the pulsars analysed
here lends support to treating the magnetospheric current sheet as volume-
distributed but ultra-thin (see Ardavan, 2021, Sect. 4.7).

###### Acknowledgements.

I thank S. Campana for his helpful comments on an earlier version of this
paper.

## References

  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 708, 1254
  * Abdo et al. (2013) Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17
  * Aleksić et al. (2011) Aleksić, J., Alvarez, E. A., Antonelli, L. A., et al. 2011, ApJ, 742, 43
  * Ansoldi et al. (2016) Ansoldi, S., Antonelli, L. A., Antoranz, P., et al. 2016, A&A, 585, A133
  * Ardavan (2021) Ardavan, H. 2021, MNRAS, 507, 4530
  * Ardavan (2022) Ardavan, H. 2022, arXiv:2206.02729
  * Ardavan (2023) Ardavan, H. 2023, A&A, 672, A154
  * Bolotovskii & Bykov (1990) Bolotovskii, B. M. & Bykov, V. P. 1990, Sov. Phys-Usp., 33, 477
  * Harris (1962) Harris, E. G. 1962, Nuovo Cim., 23, 115
  * Kalapotharakos et al. (2012) Kalapotharakos, C., Contopoulos, I., & Kazanas, D. 2012, MNRAS, 420, 2793
  * Klimchuk et al. (2023) Klimchuk, J. A., Leake, J. E., Daldorff, L. K. S., & Johnston, C. D. 2023, Front. Phys., 11:1198194
  * Manchester et al. (2005) Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993
  * Olver et al. (2010) Olver, F. W. J., Lozier, D. W., Boisvert, R. F., & Clark, C. W. 2010, NIST Handbook of Mathematical Functions (Cambridge: Cambridge Univ. Press)
  * Philippov & Kramer (2022) Philippov, A. & Kramer, M. 2022, ARA&A, 60, 495
  * Spitkovsky (2006) Spitkovsky, A. 2006, ApJ, 648, L51
  * Stamnes (1986) Stamnes, J. J. 1986, Waves in Focal Regions (Boston: Hilgar)
  * The H.E.S.S. Collaboration et al. (2023) The H.E.S.S. Collaboration et al. 2023, Nature Astron., https://doi.org/10.1038/s41550-023-02052-3
  * The MAGIC Collaboration et al. (2020) The MAGIC Collaboration et al. 2020, A&A, 643, L14
  * Zanin (2017) Zanin, R. 2017, Eur. Phys. J. Web Conf., 136, 03003